File: N:\mfix\model\calc_trd_g.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Subroutine: CALC_trD_g                                              C
4     !  Purpose: Calculate the trace of the gas phase rate of strain        C
5     !  tensor at i, j, k                                                   C
6     !                                                                      C
7     !  Author: M. Syamlal                                 Date: 19-DEC-96  C
8     !                                                                      C
9     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
10           SUBROUTINE CALC_TRD_G(lTRD_G)
11     
12     ! Modules
13     !-----------------------------------------------
14           USE param, only: DIMENSION_3
15           USE param1, only: ZERO
16           USE parallel
17           USE geometry
18           USE fldvar
19           USE indices
20           USE compar
21           USE sendrecv
22           USE functions
23           USE cutcell
24           IMPLICIT NONE
25     
26     ! Dummy arguments
27     !-----------------------------------------------
28     ! Strain rate tensor components for mth solids phase
29           DOUBLE PRECISION, INTENT(OUT) :: ltrD_g(DIMENSION_3)
30     
31     ! Local variables
32     !-----------------------------------------------
33     ! Indices
34           INTEGER :: I, J, K, IJK, IMJK, IJMK, IJKM, IM
35     
36     !-----------------------------------------------
37     
38     
39     !!$omp  parallel do private( IJK, I, J, K, IM, IMJK, IJMK, IJKM ) &
40     !!$omp& schedule(dynamic,chunk_size)
41           DO IJK = ijkstart3, ijkend3
42              IF (.NOT.WALL_AT(IJK)) THEN
43                 I = I_OF(IJK)
44                 J = J_OF(IJK)
45                 K = K_OF(IJK)
46                 IM = IM1(I)
47                 IMJK = IM_OF(IJK)
48                 IJMK = JM_OF(IJK)
49                 IJKM = KM_OF(IJK)
50     
51                 IF (.NOT. CUT_CELL_AT(IJK)) THEN
52     ! at i, j, k:
53     ! du/dx + dv/dy + 1/x dw/dz + u/x =
54     ! 1/x d(xu)/dx + dv/dy + 1/x dw/dz  =
55                    lTRD_G(IJK) = &
56                       (X_E(I)*U_G(IJK)-X_E(IM)*U_G(IMJK))*OX(I)*ODX(I) +&
57                       (V_G(IJK)-V_G(IJMK))*ODY(J) + &
58                       (W_G(IJK)-W_G(IJKM))*(OX(I)*ODZ(K))
59                 ELSE
60                    CALL CALC_CG_TRD_G(IJK,lTRD_G)
61                 ENDIF
62     
63              ELSE
64                 lTRD_G(IJK) = ZERO
65              ENDIF
66           ENDDO
67     
68           RETURN
69           END SUBROUTINE CALC_TRD_G
70     
71     
72     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
73     !                                                                      C
74     !  Subroutine: CALC_CG_trD_g                                           C
75     !  Purpose: Calculate the trace of the gas phase rate of strain        C
76     !  tensor at i, j, k with cartesian grid modifications                 C
77     !                                                                      C
78     !  Author: Jeff Dietiker                              Date: 01-Jul-09  C
79     !                                                                      C
80     !                                                                      C
81     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
82           SUBROUTINE CALC_CG_TRD_G(IJK,lTRD_G)
83     
84     ! Modules
85     !-----------------------------------------------
86           USE bc
87           USE compar
88           USE cutcell
89           USE fldvar
90           USE functions
91           USE geometry
92           USE indices
93           USE parallel
94           USE param, only: DIMENSION_3
95           USE param1, only: ZERO, HALF, UNDEFINED
96           USE quadric
97           USE sendrecv
98           IMPLICIT NONE
99     
100     ! Dummy arguments
101     !-----------------------------------------------
102     ! index
103           INTEGER, INTENT(IN) :: IJK
104     ! Strain rate tensor components for mth solids phase
105           DOUBLE PRECISION, INTENT(INOUT) :: ltrD_g(DIMENSION_3)
106     
107     ! Local variables
108     !-----------------------------------------------
109     ! Indices
110           INTEGER :: I, J, K, IMJK, IJMK, IJKM, IM
111     
112           DOUBLE PRECISION :: DEL_H,Nx,Ny,Nz
113           DOUBLE PRECISION :: dudx,dvdy,dwdz
114           DOUBLE PRECISION :: Xi,Yi,Zi,Ui,Vi,Wi,Sx,Sy,Sz
115           DOUBLE PRECISION :: UW_g,VW_g,WW_g
116           LOGICAL :: U_NODE_AT_E, U_NODE_AT_W
117           LOGICAL :: V_NODE_AT_N, V_NODE_AT_S
118           LOGICAL :: W_NODE_AT_T, W_NODE_AT_B
119           INTEGER :: BCV
120           INTEGER :: BCT
121     !-----------------------------------------------
122     
123     
124           I = I_OF(IJK)
125           J = J_OF(IJK)
126           K = K_OF(IJK)
127           IM = IM1(I)
128           IMJK = IM_OF(IJK)
129           IJMK = JM_OF(IJK)
130           IJKM = KM_OF(IJK)
131     
132           IF(FLOW_AT(IJK)) THEN
133              lTRD_G(IJK) = ZERO
134              RETURN
135           ENDIF
136     
137           BCV = BC_ID(IJK)
138           IF(BCV > 0 ) THEN
139              BCT = BC_TYPE_ENUM(BCV)
140           ELSE
141              BCT = NONE
142           ENDIF
143           SELECT CASE (BCT)
144              CASE (CG_NSW)
145                 NOC_TRDG = .TRUE.
146                 UW_g = ZERO
147                 VW_g = ZERO
148                 WW_g = ZERO
149              CASE (CG_FSW)
150                 NOC_TRDG = .FALSE.
151                 UW_g = ZERO
152                 VW_g = ZERO
153                 WW_g = ZERO
154              CASE(CG_PSW)
155                 IF(BC_HW_G(BCV)==UNDEFINED) THEN   ! same as NSW
156                    NOC_TRDG = .TRUE.
157                    UW_g = BC_UW_G(BCV)
158                    VW_g = BC_VW_G(BCV)
159                    WW_g = BC_WW_G(BCV)
160                 ELSEIF(BC_HW_G(BCV)==ZERO) THEN   ! same as FSW
161                    NOC_TRDG = .FALSE.
162                    UW_g = ZERO
163                    VW_g = ZERO
164                    WW_g = ZERO
165                 ELSE                              ! partial slip
166                    NOC_TRDG = .FALSE.
167                 ENDIF
168              CASE (CG_MI)
169                 lTRD_G(IJK) = ZERO
170                 RETURN
171              CASE (CG_PO)
172                 lTRD_G(IJK) = ZERO
173                 RETURN
174              CASE (NONE)
175                 lTRD_G(IJK) = ZERO
176                 RETURN
177           END SELECT
178     
179     
180     ! du/dx
181     !=======================================================================
182           U_NODE_AT_E = ((.NOT.BLOCKED_U_CELL_AT(IJK)) .AND.&
183                          (.NOT.WALL_U_AT(IJK)))
184           U_NODE_AT_W = ((.NOT.BLOCKED_U_CELL_AT(IMJK)).AND.&
185                          (.NOT.WALL_U_AT(IMJK)))
186           IF(U_NODE_AT_E.AND.U_NODE_AT_W) THEN
187              Ui = HALF * (U_G(IJK) + U_G(IMJK))
188              Xi = HALF * (X_U(IJK) + X_U(IMJK))
189              Yi = HALF * (Y_U(IJK) + Y_U(IMJK))
190              Zi = HALF * (Z_U(IJK) + Z_U(IMJK))
191              Sx = X_U(IJK) - X_U(IMJK)
192              Sy = Y_U(IJK) - Y_U(IMJK)
193              Sz = Z_U(IJK) - Z_U(IMJK)
194              CALL GET_DEL_H(IJK,'SCALAR',Xi,Yi,Zi,Del_H,Nx,Ny,Nz)
195              IF(abs(Sx) > ZERO) THEN
196                 dudx =  (U_G(IJK) - U_G(IMJK))/Sx
197                 IF(NOC_TRDG) dudx = dudx - ((Ui-UW_g)/(Sx*DEL_H)*&
198                    (Sy*Ny+Sz*Nz))
199              ELSE
200                 dudx = ZERO
201              ENDIF
202           ELSE
203              dudx = ZERO
204           ENDIF
205     
206     ! dv/dy
207     !=======================================================================
208           V_NODE_AT_N = ((.NOT.BLOCKED_V_CELL_AT(IJK)) .AND.&
209                          (.NOT.WALL_V_AT(IJK)))
210           V_NODE_AT_S = ((.NOT.BLOCKED_V_CELL_AT(IJMK)).AND.&
211                          (.NOT.WALL_V_AT(IJMK)))
212           IF(V_NODE_AT_N.AND.V_NODE_AT_S) THEN
213              Vi = HALF * (V_G(IJK) + V_G(IJMK))
214              Xi = HALF * (X_V(IJK) + X_V(IJMK))
215              Yi = HALF * (Y_V(IJK) + Y_V(IJMK))
216              Zi = HALF * (Z_V(IJK) + Z_V(IJMK))
217              Sx = X_V(IJK) - X_V(IJMK)
218              Sy = Y_V(IJK) - Y_V(IJMK)
219              Sz = Z_V(IJK) - Z_V(IJMK)
220              CALL GET_DEL_H(IJK,'SCALAR',Xi,Yi,Zi,Del_H,Nx,Ny,Nz)
221              IF(abs(Sy) > ZERO) THEN
222                 dvdy =  (V_G(IJK) - V_G(IJMK))/Sy
223                 IF(NOC_TRDG) dvdy = dvdy - ((Vi-VW_g)/(Sy*DEL_H)*&
224                    (Sx*Nx+Sz*Nz))
225              ELSE
226                 dvdy =  ZERO
227              ENDIF
228           ELSEIF (V_NODE_AT_N.AND.(.NOT.V_NODE_AT_S).AND.NOC_TRDG) THEN
229              CALL GET_DEL_H(IJK,'SCALAR',X_V(IJK),Y_V(IJK),&
230                             Z_V(IJK),DEL_H,Nx,Ny,Nz)
231              dvdy = (V_g(IJK) - VW_g) / DEL_H * Ny
232           ELSEIF ((.NOT.V_NODE_AT_N).AND.V_NODE_AT_S.AND.NOC_TRDG) THEN
233              CALL GET_DEL_H(IJK,'SCALAR',X_V(IJMK),Y_V(IJMK),&
234                             Z_V(IJMK),DEL_H,Nx,Ny,Nz)
235              dvdy = (V_g(IJMK) - VW_g) / DEL_H * Ny
236           ELSE
237              dvdy = ZERO
238           ENDIF
239     
240     ! dw/dz
241     !=======================================================================
242           IF(NO_K) THEN
243              dwdz = ZERO
244           ELSE
245              W_NODE_AT_T = ((.NOT.BLOCKED_W_CELL_AT(IJK)) .AND.&
246                 (.NOT.WALL_W_AT(IJK)))
247              W_NODE_AT_B = ((.NOT.BLOCKED_W_CELL_AT(IJKM)).AND.&
248                 (.NOT.WALL_W_AT(IJKM)))
249              IF(W_NODE_AT_T.AND.W_NODE_AT_B) THEN
250                 Wi = HALF * (W_G(IJK) + W_G(IJKM))
251                 Xi = HALF * (X_W(IJK) + X_W(IJKM))
252                 Yi = HALF * (Y_W(IJK) + Y_W(IJKM))
253                 Zi = HALF * (Z_W(IJK) + Z_W(IJKM))
254                 Sx = X_W(IJK) - X_W(IJKM)
255                 Sy = Y_W(IJK) - Y_W(IJKM)
256                 Sz = Z_W(IJK) - Z_W(IJKM)
257                 CALL GET_DEL_H(IJK,'SCALAR',Xi,Yi,Zi,Del_H,Nx,Ny,Nz)
258                 IF(abs(Sz) > ZERO) THEN
259                    dwdz =  (W_G(IJK) - W_G(IJKM))/Sz
260                    IF(NOC_TRDG) dwdz = dwdz - ((Wi-WW_g)/(Sz*DEL_H)*&
261                       (Sx*Nx+Sy*Ny))
262                 ELSE
263                    dwdz = ZERO
264                 ENDIF
265              ELSEIF (W_NODE_AT_T.AND.(.NOT.W_NODE_AT_B).AND.NOC_TRDG) THEN
266                 CALL GET_DEL_H(IJK,'SCALAR',X_W(IJK),Y_W(IJK),&
267                                Z_W(IJK),DEL_H,Nx,Ny,Nz)
268                 dwdz = (W_g(IJK) - WW_g) / DEL_H * Nz
269              ELSEIF ((.NOT.W_NODE_AT_T).AND.W_NODE_AT_B.AND.NOC_TRDG) THEN
270                 CALL GET_DEL_H(IJK,'SCALAR',X_W(IJKM),Y_W(IJKM),&
271                                Z_W(IJKM),DEL_H,Nx,Ny,Nz)
272                 dwdz = (W_g(IJKM) - WW_g) / DEL_H * Nz
273              ELSE
274                 dwdz = ZERO
275              ENDIF
276           ENDIF   ! end if/else no_k
277     
278           lTRD_G(IJK) = dudx + dvdy + dwdz
279     
280           RETURN
281           END SUBROUTINE CALC_CG_TRD_G
282     
283     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
284     !                                                                      C
285     !  Subroutine: GET_FACE_VEL_GAS                                        C
286     !  Purpose: Evaluate the velocity components at each of the faces of   C
287     !  a scalar cell                                                       C
288     !                                                                      C
289     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
290           SUBROUTINE GET_FACE_VEL_GAS(IJK, uge, ugw, ugn, ugs, ugt, ugb, ugcc, &
291                                            vge, vgw, vgn, vgs, vgt, vgb, &
292                                            wge, wgw, wgn, wgs, wgt, wgb, wgcc)
293     
294     ! Modules
295     !-----------------------------------------------------------------------
296           use fldvar, only: u_g, v_g, w_g
297           use functions, only: im_of, jm_of, km_of
298           use functions, only: ip_of, jp_of, kp_of
299           use fun_avg, only: avg_x, avg_x_e
300           use fun_avg, only: avg_y, avg_y_n
301           use fun_avg, only: avg_z, avg_z_t
302           use geometry, only: cylindrical
303           use indices, only: i_of, j_of, k_of
304           use indices, only: im1, jm1, km1
305           use param1, only: zero
306           IMPLICIT NONE
307     
308     ! Dummy arguments
309     !-----------------------------------------------------------------------
310     ! index
311           INTEGER, INTENT(IN) :: IJK
312     ! U_g at the east (i+1/2, j, k) and west face (i-1/2, j, k)
313           DOUBLE PRECISION, INTENT(OUT) :: UgE, UgW
314     ! U_g at the north (i, j+1/2, k) and south face (i, j-1/2, k) 
315           DOUBLE PRECISION, INTENT(OUT) :: UgN, UgS
316     ! U_g at the top (i, j, k+1/2) and bottom face (i, j, k-1/2)
317           DOUBLE PRECISION, INTENT(OUT) :: UgT, UgB
318     ! U_g at the center of a scalar cell (i, j, k)
319     ! Calculated for Cylindrical coordinates only.
320           DOUBLE PRECISION, INTENT(OUT) :: UgcC
321     
322     ! V_g at the east (i+1/2, j, k) and west face (i-1/2, j, k)
323           DOUBLE PRECISION, INTENT(OUT) :: VgE, VgW
324     ! V_g at the north (i, j+1/2, k) and south face (i, j-1/2, k)
325           DOUBLE PRECISION, INTENT(OUT) :: VgN, VgS
326     ! V_g at the top (i, j, k+1/2) and bottom face (i, j, k-1/2)
327           DOUBLE PRECISION, INTENT(OUT) :: VgT, VgB
328     
329     ! W_g at the east (i+1/2, j, k) and west face (i-1/2, j, k)
330           DOUBLE PRECISION, INTENT(OUT) :: WgE, WgW
331     ! W_g at the north (i, j+1/2, k) and south face (i, j-1/2, k)
332           DOUBLE PRECISION, INTENT(OUT) :: WgN, WgS
333     ! W_g at the top (i, j, k+1/2) and bottom face (i, j, k-1/2)
334           DOUBLE PRECISION, INTENT(OUT) :: WgT, WgB
335     ! W_g at the center of a scalar cell (i, j, k).
336     ! Calculated for Cylindrical coordinates only.
337           DOUBLE PRECISION, INTENT(OUT) :: WgcC
338     
339     ! Local variables
340     !-----------------------------------------------------------------------
341     ! Cell indices
342           INTEGER :: I, J, K, IM, JM, KM
343           INTEGER :: IMJK, IPJK, IJMK, IJPK, IJKM, IJKP
344           INTEGER :: IMJPK, IMJMK, IMJKP, IMJKM, IPJKM, IPJMK
345           INTEGER :: IJMKP, IJMKM, IJPKM
346     !-----------------------------------------------------------------------
347     
348           I = I_OF(IJK)
349           J = J_OF(IJK)
350           K = K_OF(IJK)
351           IM = IM1(I)
352           JM = JM1(J)
353           KM = KM1(K)
354           IMJK = IM_OF(IJK)
355           IPJK = IP_OF(IJK)
356           IJMK = JM_OF(IJK)
357           IJPK = JP_OF(IJK)
358           IJKM = KM_OF(IJK)
359           IJKP = KP_OF(IJK)
360           IMJPK = IM_OF(IJPK)
361           IMJMK = IM_OF(IJMK)
362           IMJKP = IM_OF(IJKP)
363           IMJKM = IM_OF(IJKM)
364           IPJKM = IP_OF(IJKM)
365           IPJMK = IP_OF(IJMK)
366           IJMKP = JM_OF(IJKP)
367           IJMKM = JM_OF(IJKM)
368           IJPKM = JP_OF(IJKM)
369     
370     ! Find fluid velocity values at faces of the cell
371           UgN = AVG_Y(AVG_X_E(U_G(IMJK),U_G(IJK),I), &
372                       AVG_X_E(U_G(IMJPK),U_G(IJPK),I),J)     !i, j+1/2, k
373           UgS = AVG_Y(AVG_X_E(U_G(IMJMK),U_G(IJMK),I),&
374                       AVG_X_E(U_G(IMJK),U_G(IJK),I),JM)      !i, j-1/2, k
375           UgT = AVG_Z(AVG_X_E(U_G(IMJK),U_G(IJK),I),&
376                       AVG_X_E(U_G(IMJKP),U_G(IJKP),I),K)     !i, j, k+1/2
377           UgB = AVG_Z(AVG_X_E(U_G(IMJKM),U_G(IJKM),I),&
378                       AVG_X_E(U_G(IMJK),U_G(IJK),I),KM)      !i, j, k-1/2
379           UgE = U_G(IJK)
380           UgW = U_G(IMJK)
381     
382           VgE = AVG_X(AVG_Y_N(V_G(IJMK),V_G(IJK)),&
383                       AVG_Y_N(V_G(IPJMK),V_G(IPJK)),I)       !i+1/2, j, k
384           VgW = AVG_X(AVG_Y_N(V_G(IMJMK),V_G(IMJK)),&
385                       AVG_Y_N(V_G(IJMK),V_G(IJK)),IM)        !i-1/2, j, k
386           VgT = AVG_Z(AVG_Y_N(V_G(IJMK),V_G(IJK)),&
387                       AVG_Y_N(V_G(IJMKP),V_G(IJKP)),K)       !i, j, k+1/2
388           VgB = AVG_Z(AVG_Y_N(V_G(IJMKM),V_G(IJKM)),&
389                       AVG_Y_N(V_G(IJMK),V_G(IJK)),KM)        !i, j, k-1/2
390           VgN = V_G(IJK)
391           VgS = V_G(IJMK)
392     
393           WgN = AVG_Y(AVG_Z_T(W_G(IJKM),W_G(IJK)),&
394                       AVG_Z_T(W_G(IJPKM),W_G(IJPK)),J)       !i, j+1/2, k
395           WgS = AVG_Y(AVG_Z_T(W_G(IJMKM),W_G(IJMK)),&
396                       AVG_Z_T(W_G(IJKM),W_G(IJK)),JM)        !i, j-1/2, k
397           WgE = AVG_X(AVG_Z_T(W_G(IJKM),W_G(IJK)),&
398                       AVG_Z_T(W_G(IPJKM),W_G(IPJK)),I)       !i+1/2, j, k
399           WgW = AVG_X(AVG_Z_T(W_G(IMJKM),W_G(IMJK)),&
400                       AVG_Z_T(W_G(IJKM),W_G(IJK)),IM)        !i-1/2, j, k
401           WgT = W_G(IJK)
402           WgB = W_G(IJKM)
403     
404           IF (CYLINDRICAL) THEN
405              UgcC = AVG_X_E(U_G(IMJK),U_G(IJK),I)             !i, j, k
406              WgcC = AVG_Z_T(W_G(IJKM),W_G(IJK))               !i, j, k
407           ELSE
408              UgcC = ZERO
409              WgcC = ZERO
410           ENDIF
411     
412           RETURN
413           END SUBROUTINE GET_FACE_VEL_GAS
414     
415     
416     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
417     !                                                                      C
418     !  Subroutine: CALC_DERIV_VEL_GAS                                      C
419     !  Purpose: Calculate the gradient of the gas phase velocity and       C
420     !  the rate of strain tensor at i, j, k                                C
421     !                                                                      C
422     !         |  du/dx    du/dy   du/dz  |                                 C
423     ! DELV =  |  dv/dx    dv/dy   dv/dz  |  =  dVi/dxj                     C
424     !         |  dw/dx    dw/dy   dw/dz  |                                 C
425     !                                                                      C
426     ! 1/2(DELV+DELV^T) = 1/2(dVi/dxj+dVj/dxi)                              C
427     !                                                                      C
428     !                                                                      C
429     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
430           SUBROUTINE CALC_DERIV_VEL_GAS(IJK, lVelGradG, lRateStrainG)
431     
432     ! Modules
433     !-----------------------------------------------------------------------
434           use cutcell, only: cut_cell_at
435           use geometry, only: odx, ody, odz
436           use geometry, only: ox
437           use indices, only: i_of, j_of, k_of
438           use param1, only: half
439     ! Dummy arguments
440     !-----------------------------------------------------------------------
441     ! index
442           INTEGER, INTENT(IN) :: IJK
443     ! gradient of velocity
444           DOUBLE PRECISION, INTENT(OUT) :: lVelGradG(3,3)    ! delV
445     ! rate of strain tensor
446           DOUBLE PRECISION, INTENT(OUT) :: lRateStrainG(3,3) ! D_g
447     
448     ! Local variables
449     !-----------------------------------------------------------------------
450     ! Cell indices
451           INTEGER :: I, J, K
452     ! face values of velocity
453           DOUBLE PRECISION :: uge, ugw, ugn, ugs, ugt, ugb, ugcc
454           DOUBLE PRECISION :: vge, vgw, vgn, vgs, vgt, vgb
455           DOUBLE PRECISION :: wge, wgw, wgn, wgs, wgt, wgb, wgcc
456     !-----------------------------------------------------------------------
457     
458           IF(.NOT.CUT_CELL_AT(IJK)) THEN
459              I = I_OF(IJK)
460              J = J_OF(IJK)
461              K = K_OF(IJK)
462     
463     ! Set the velocity components at the scalar faces
464              CALL GET_FACE_VEL_GAS(IJK, uge, ugw, ugn, ugs, ugt, ugb, ugcc, &
465                               vge, vgw, vgn, vgs, vgt, vgb, &
466                               wge, wgw, wgn, wgs, wgt, wgb, wgcc)
467     
468     ! Gradient of gas velocity  at cell center (i, j, k)
469              lVelGradG(1,1) = (UgE-UgW)*ODX(I)
470              lVelGradG(1,2) = (UgN-UgS)*ODY(J)
471              lVelGradG(1,3) = (UgT-UgB)*(OX(I)*ODZ(K))-WgcC*OX(I)
472     
473              lVelGradG(2,1) = (VgE-VgW)*ODX(I)
474              lVelGradG(2,2) = (VgN-VgS)*ODY(J)
475              lVelGradG(2,3) = (VgT-VgB)*(OX(I)*ODZ(K))
476     
477              lVelGradG(3,1) = (WgE-WgW)*ODX(I)
478              lVelGradG(3,2) = (WgN-WgS)*ODY(J)
479              lVelGradG(3,3) = (WgT-WgB)*(OX(I)*ODZ(K)) + UgcC*OX(I)
480     
481     
482     ! Strain rate tensor, D_G, at cell center (i, j, k)
483     ! or calculated as 0.5*(DelG+DelG^T)
484              lRateStrainG(1,1) = (UgE-UgW)*ODX(I)
485              lRateStrainG(1,2) = HALF*((UgN-UgS)*ODY(J)+&
486                                        (VgE-VgW)*ODX(I))
487              lRateStrainG(1,3) = HALF*((WgE-WgW)*ODX(I)+&
488                                        (UgT-UgB)*(OX(I)*ODZ(K))-&
489                                  WgcC*OX(I))
490     
491              lRateStrainG(2,1) = lRateStrainG(1,2)
492              lRateStrainG(2,2) = (VgN-VgS)*ODY(J)
493              lRateStrainG(2,3) = HALF*((VgT-VgB)*(OX(I)*ODZ(K))+&
494                                        (WgN-WgS)*ODY(J))
495     
496              lRateStrainG(3,1) = lRateStrainG(1,3)
497              lRateStrainG(3,2) = lRateStrainG(2,3)
498              lRateStrainG(3,3) = (WgT-WgB)*(OX(I)*ODZ(K)) +&
499                                  UgcC*OX(I)
500     
501           ELSE   ! cut cell
502              CALL CALC_CG_DERIV_VEL_GAS(IJK,lVelGradG,lRateStrainG)
503           ENDIF
504     
505           RETURN
506           END SUBROUTINE CALC_DERIV_VEL_GAS
507     
508     
509     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
510     !                                                                      C
511     !  Module name: CALC_CG_DERIV_VEL_GAS                                  C
512     !  Purpose: Calculate velocity derivatives in scalar cut-cell          C
513     !           gas phase                                                  C
514     !                                                                      C
515     !  Author: Jeff Dietiker                              Date: 25-JAN-96  C
516     !                                                                      C
517     !                                                                      C
518     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
519           SUBROUTINE CALC_CG_DERIV_VEL_GAS(IJK, lVelGradG, lRateStrainG)
520     
521     ! Modules
522     !-----------------------------------------------
523           USE param
524           USE param1
525           USE parallel
526           USE geometry
527           USE fldvar
528           USE indices
529           USE compar
530           USE sendrecv
531           USE bc
532           USE cutcell
533           USE quadric
534           USE functions
535           IMPLICIT NONE
536     
537     ! Dummy arguments
538     !-----------------------------------------------
539     ! index
540           INTEGER, INTENT(IN) :: IJK
541     ! velocity gradient gas phase
542     !         |  du/dx    du/dy   du/dz  |
543     ! DELV =  |  dv/dx    dv/dy   dv/dz  |  =  dUi/dxj
544     !         |  dw/dx    dw/dy   dw/dz  |
545           DOUBLE PRECISION, INTENT(OUT) :: lVelGradG(3,3)
546     ! rate of strain tensor gas phase
547           DOUBLE PRECISION, INTENT(OUT) :: lRateStrainG(3,3)
548     
549     ! Local variables
550     !-----------------------------------------------
551     ! Indices
552           INTEGER :: I, J, K, IMJK, IJMK, IJKM
553     ! loop counters
554           INTEGER :: P, Q
555     !
556           DOUBLE PRECISION :: DEL_H,Nx,Ny,Nz
557           DOUBLE PRECISION :: dudx,dudy,dudz
558           DOUBLE PRECISION :: dvdx,dvdy,dvdz
559           DOUBLE PRECISION :: dwdx,dwdy,dwdz
560           DOUBLE PRECISION :: Xi,Yi,Zi,Ui,Vi,Wi,Sx,Sy,Sz
561           DOUBLE PRECISION :: UW_g,VW_g,WW_g
562     
563           LOGICAL :: U_NODE_AT_E, U_NODE_AT_W
564           LOGICAL :: V_NODE_AT_N, V_NODE_AT_S
565           LOGICAL :: W_NODE_AT_T, W_NODE_AT_B
566           INTEGER :: BCV
567           INTEGER :: BCT
568     !-----------------------------------------------
569     
570     ! initialize
571           lVelGradG = ZERO
572           lRateStrainG = ZERO
573     
574           I = I_OF(IJK)
575           J = J_OF(IJK)
576           K = K_OF(IJK)
577           IMJK = IM_OF(IJK)
578           IJMK = JM_OF(IJK)
579           IJKM = KM_OF(IJK)
580     
581           IF(FLOW_AT(IJK)) THEN
582              lVelGradG = ZERO
583              RETURN
584           ENDIF
585     
586           BCV = BC_ID(IJK)
587           IF(BCV > 0 ) THEN
588              BCT = BC_TYPE_ENUM(BCV)
589           ELSE
590              BCT = NONE
591           ENDIF
592           SELECT CASE (BCT)
593              CASE (CG_NSW)
594                 NOC_TRDG = .TRUE.
595                 UW_g = ZERO
596                 VW_g = ZERO
597                 WW_g = ZERO
598              CASE (CG_FSW)
599                 NOC_TRDG = .FALSE.
600                 UW_g = ZERO
601                 VW_g = ZERO
602                 WW_g = ZERO
603              CASE(CG_PSW)
604                 IF(BC_HW_G(BCV)==UNDEFINED) THEN   ! same as NSW
605                    NOC_TRDG = .TRUE.
606                    UW_g = BC_UW_G(BCV)
607                    VW_g = BC_VW_G(BCV)
608                    WW_g = BC_WW_G(BCV)
609                 ELSEIF(BC_HW_G(BCV)==ZERO) THEN   ! same as FSW
610                    NOC_TRDG = .FALSE.
611                    UW_g = ZERO
612                    VW_g = ZERO
613                    WW_g = ZERO
614                 ELSE                              ! partial slip
615                    NOC_TRDG = .FALSE.
616                 ENDIF
617              CASE (CG_MI)
618                 lVelGradG = ZERO
619                 RETURN
620              CASE (CG_PO)
621                 lVelGradG = ZERO
622                 RETURN
623              CASE (NONE)
624                 lVelGradG = ZERO
625                 RETURN
626           END SELECT
627     
628     
629     ! du/dx, du/dy, du/dz
630     !=======================================================================
631           U_NODE_AT_E = ((.NOT.BLOCKED_U_CELL_AT(IJK)) .AND.&
632                          (.NOT.WALL_U_AT(IJK)))
633           U_NODE_AT_W = ((.NOT.BLOCKED_U_CELL_AT(IMJK)).AND.&
634                          (.NOT.WALL_U_AT(IMJK)))
635     
636           IF(U_NODE_AT_E.AND.U_NODE_AT_W) THEN
637              Ui = HALF * (U_G(IJK) + U_G(IMJK))
638              Xi = HALF * (X_U(IJK) + X_U(IMJK))
639              Yi = HALF * (Y_U(IJK) + Y_U(IMJK))
640              Zi = HALF * (Z_U(IJK) + Z_U(IMJK))
641              Sx = X_U(IJK) - X_U(IMJK)
642              Sy = Y_U(IJK) - Y_U(IMJK)
643              Sz = Z_U(IJK) - Z_U(IMJK)
644              CALL GET_DEL_H(IJK,'SCALAR',Xi,Yi,Zi,Del_H,Nx,Ny,Nz)
645              IF(abs(Sx) > ZERO) THEN
646                 dudx =  (U_G(IJK) - U_G(IMJK))/Sx
647                 dudy =  ZERO
648                 dudz =  ZERO
649                 IF(NOC_TRDG) THEN
650                    dudx = dudx - ((Ui-UW_g)/(Sx*DEL_H)*(Sy*Ny+Sz*Nz))
651                    dudy = (Ui-UW_g) / DEL_H * Ny
652                    dudz = (Ui-UW_g) / DEL_H * Nz
653                 ENDIF
654              ELSE
655                 dudx = ZERO
656                 dudy = ZERO
657                 dudz = ZERO
658              ENDIF
659           ELSEIF (U_NODE_AT_E.AND.(.NOT.U_NODE_AT_W).AND.NOC_TRDG) THEN
660              CALL GET_DEL_H(IJK,'SCALAR',X_U(IJK),Y_U(IJK),Z_U(IJK),&
661                             DEL_H,Nx,Ny,Nz)
662              dudx = (U_g(IJK) - UW_g) / DEL_H * Nx
663              dudy = (U_g(IJK) - UW_g) / DEL_H * Ny
664              dudz = (U_g(IJK) - UW_g) / DEL_H * Nz
665           ELSEIF ((.NOT.U_NODE_AT_E).AND.U_NODE_AT_W.AND.NOC_TRDG) THEN
666              CALL GET_DEL_H(IJK,'SCALAR',X_U(IMJK),Y_U(IMJK),Z_U(IMJK),&
667                             DEL_H,Nx,Ny,Nz)
668              dudx = (U_g(IMJK) - UW_g) / DEL_H * Nx
669              dudy = (U_g(IMJK) - UW_g) / DEL_H * Ny
670              dudz = (U_g(IMJK) - UW_g) / DEL_H * Nz
671           ELSE
672              dudx = ZERO
673              dudy = ZERO
674              dudz = ZERO
675           ENDIF
676           lVelGradG(1,1) = dudx
677           lVelGradG(1,2) = dudy
678           lVelGradG(1,3) = dudz
679     
680     ! dv/dx, dv/dy, dv/dz
681     !=======================================================================
682           V_NODE_AT_N = ((.NOT.BLOCKED_V_CELL_AT(IJK)) .AND.&
683              (.NOT.WALL_V_AT(IJK)))
684           V_NODE_AT_S = ((.NOT.BLOCKED_V_CELL_AT(IJMK)).AND.&
685              (.NOT.WALL_V_AT(IJMK)))
686     
687           IF(V_NODE_AT_N.AND.V_NODE_AT_S) THEN
688              Vi = HALF * (V_G(IJK) + V_G(IJMK))
689              Xi = HALF * (X_V(IJK) + X_V(IJMK))
690              Yi = HALF * (Y_V(IJK) + Y_V(IJMK))
691              Zi = HALF * (Z_V(IJK) + Z_V(IJMK))
692              Sx = X_V(IJK) - X_V(IJMK)
693              Sy = Y_V(IJK) - Y_V(IJMK)
694              Sz = Z_V(IJK) - Z_V(IJMK)
695              CALL GET_DEL_H(IJK,'SCALAR',Xi,Yi,Zi,Del_H,Nx,Ny,Nz)
696              IF(abs(Sy) > ZERO) THEN
697                 dvdx =  ZERO
698                 dvdy =  (V_G(IJK) - V_G(IJMK))/Sy
699                 dvdz =  ZERO
700                 IF(NOC_TRDG) THEN
701                    dvdx = (Vi-VW_g) / DEL_H * Nx
702                    dvdy = dvdy - ((Vi-VW_g)/(Sy*DEL_H)*(Sx*Nx+Sz*Nz))
703                    dvdz = (Vi-VW_g) / DEL_H * Nz
704                 ENDIF
705              ELSE
706                 dvdx =  ZERO
707                 dvdy =  ZERO
708                 dvdz =  ZERO
709              ENDIF
710           ELSEIF (V_NODE_AT_N.AND.(.NOT.V_NODE_AT_S).AND.NOC_TRDG) THEN
711              CALL GET_DEL_H(IJK,'SCALAR',X_V(IJK),Y_V(IJK),Z_V(IJK),&
712                             DEL_H,Nx,Ny,Nz)
713              dvdx = (V_g(IJK) - VW_g) / DEL_H * Nx
714              dvdy = (V_g(IJK) - VW_g) / DEL_H * Ny
715              dvdz = (V_g(IJK) - VW_g) / DEL_H * Nz
716           ELSEIF ((.NOT.V_NODE_AT_N).AND.V_NODE_AT_S.AND.NOC_TRDG) THEN
717              CALL GET_DEL_H(IJK,'SCALAR',X_V(IJMK),Y_V(IJMK),Z_V(IJMK),&
718                             DEL_H,Nx,Ny,Nz)
719              dvdx = (V_g(IJMK) - VW_g) / DEL_H * Nx
720              dvdy = (V_g(IJMK) - VW_g) / DEL_H * Ny
721              dvdz = (V_g(IJMK) - VW_g) / DEL_H * Nz
722           ELSE
723              dvdx =  ZERO
724              dvdy =  ZERO
725              dvdz =  ZERO
726           ENDIF
727           lVelGradG(2,1) = dvdx
728           lVelGradG(2,2) = dvdy
729           lVelGradG(2,3) = dvdz
730     
731     ! dw/dx, dw/dy, dw/dz
732     !=======================================================================
733           IF(NO_K) THEN
734              dwdx = ZERO
735              dwdy = ZERO
736              dwdz = ZERO
737           ELSE
738              W_NODE_AT_T = ((.NOT.BLOCKED_W_CELL_AT(IJK)) .AND.&
739                 (.NOT.WALL_W_AT(IJK)))
740              W_NODE_AT_B = ((.NOT.BLOCKED_W_CELL_AT(IJKM)).AND.&
741                 (.NOT.WALL_W_AT(IJKM)))
742              IF(W_NODE_AT_T.AND.W_NODE_AT_B) THEN
743                 Wi = HALF * (W_G(IJK) + W_G(IJKM))
744                 Xi = HALF * (X_W(IJK) + X_W(IJKM))
745                 Yi = HALF * (Y_W(IJK) + Y_W(IJKM))
746                 Zi = HALF * (Z_W(IJK) + Z_W(IJKM))
747                 Sx = X_W(IJK) - X_W(IJKM)
748                 Sy = Y_W(IJK) - Y_W(IJKM)
749                 Sz = Z_W(IJK) - Z_W(IJKM)
750                 CALL GET_DEL_H(IJK,'SCALAR',Xi,Yi,Zi,Del_H,Nx,Ny,Nz)
751                 IF(abs(Sz) > ZERO) THEN
752                    dwdx =  ZERO
753                    dwdy =  ZERO
754                    dwdz =  (W_G(IJK) - W_G(IJKM))/Sz
755                    IF(NOC_TRDG) THEN
756                       dwdx = (Wi-WW_g) / DEL_H * Nx
757                       dwdy = (Wi-WW_g) / DEL_H * Ny
758                       dwdz = dwdz - ((Wi-WW_g)/(Sz*DEL_H)*(Sx*Nx+Sy*Ny))
759                    ENDIF
760                 ELSE
761                    dwdx = ZERO
762                    dwdy = ZERO
763                    dwdz = ZERO
764                 ENDIF
765              ELSEIF (W_NODE_AT_T.AND.(.NOT.W_NODE_AT_B).AND.NOC_TRDG) THEN
766                 CALL GET_DEL_H(IJK,'SCALAR',X_W(IJK),Y_W(IJK),Z_W(IJK),&
767                                DEL_H,Nx,Ny,Nz)
768                 dwdx = (W_g(IJK) - WW_g) / DEL_H * Nx
769                 dwdy = (W_g(IJK) - WW_g) / DEL_H * Ny
770                 dwdz = (W_g(IJK) - WW_g) / DEL_H * Nz
771              ELSEIF ((.NOT.W_NODE_AT_T).AND.W_NODE_AT_B.AND.NOC_TRDG) THEN
772                 CALL GET_DEL_H(IJK,'SCALAR',X_W(IJKM),Y_W(IJKM),Z_W(IJKM),&
773                                DEL_H,Nx,Ny,Nz)
774                 dwdx = (W_g(IJKM) - WW_g) / DEL_H * Nx
775                 dwdy = (W_g(IJKM) - WW_g) / DEL_H * Ny
776                 dwdz = (W_g(IJKM) - WW_g) / DEL_H * Nz
777              ELSE
778                 dwdx = ZERO
779                 dwdy = ZERO
780                 dwdz = ZERO
781              ENDIF
782           ENDIF
783           lVelGradG(3,1) = dwdx
784           lVelGradG(3,2) = dwdy
785           lVelGradG(3,3) = dwdz
786     
787           DO P = 1,3
788              DO Q = 1,3
789                 lRateStrainG(P,Q) = HALF * (lVelGradG(P,Q)+lVelGradG(Q,P))
790              ENDDO
791           ENDDO
792     
793     
794           RETURN
795           END SUBROUTINE CALC_CG_DERIV_VEL_GAS
796