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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Subroutine: CALC_Tau_U_g                                            C
4     !  Purpose: Cross terms in the gradient of stress in U_g momentum      C
5     !                                                                      C
6     !  Author: M. Syamlal                                 Date: 18-DEC-96  C
7     !                                                                      C
8     !                                                                      C
9     !  Comments: This routine calculates the components of the gas phase   C
10     !  viscous stress tensor of the u-momentum equation that cannot be     C
11     !  cast in the form: mu.grad(u). These components are stored in the    C
12     !  passed variable, which is then used as a source of the u-momentum   C
13     !  equation.                                                           C
14     !                                                                      C
15     !  The following u component viscous stress tensor terms are           C
16     !  calculated here:                                                    C
17     !  > Combined part of 1/x d/dx (x.tau_xx) xdxdydz and                  C
18     !                              -tau_zz/x xdxdydz =>                    C
19     !    1/x d/dx (x.lambda.trcD) xdxdydz - (lambda/x.trcD) xdxdydz =>     C
20     !        d/dx (lambda.trcD) xdxdydz =>                                 C
21     !    delta (lambda.trcD)Ap|E-W                                         C
22     !  > part of 1/x d/dx(x.tau_xx) xdxdydz =>                             C
23     !            1/x d/dx (x.mu.du/dx) xdxdydz =>                          C
24     !    delta (mu du/dx)Ayz |E-W                                          C
25     !  > part of d/dy (tau_xy) xdxdydz =>                                  C
26     !           d/dy (mu.dv/dx) xdxdydz =>                                 C
27     !    delta (mu.dv/dx)Axz |N-S                                          C
28     !  > part of 1/x d/dz (tau_xz) xdxdydz =>                              C
29     !            1/x d/dz (mu.dw/dx) xdxdydz =>                            C
30     !    delta (mu.dw/dx)Axy |T-B                                          C
31     !  CYLINDRICAL TERMS:                                                  C
32     !  > part of 1/x d/dz (tau_xz) xdxdydz =>                              C
33     !            1/x d/dz (mu.(-w/x)) xdxdydz =>                           C
34     !    delta (mu(-w/x))Axy |T-B                                          C
35     !  > part of -tau_zz/x xdxdydz =>                                      C
36     !            -(2.mu/x).(1/x).dw/dz xdxdydz                             C
37     !    delta (2.mu/x.1/x.dw/dz)Vp |E-W                                   C
38     !                                                                      C
39     !  The following u-component CYLINDRICAL viscous stress tensor terms   C
40     !  are not calculated here. They are handled in source_u_g:            C
41     !  > part of the term -tau_zz/x xdxdydz =>                             C
42     !                     -2.mu/x. u/x xdxdydz                             C
43     !    delta(-2.mu.u/xp^2)Vp |p                                          C
44     !                                                                      C
45     !  To reconstitute the complete u-momentum gas phase viscous stress    C
46     !  tensor would require including the term calculated in source_u_g    C
47     !  and the 'diffusional' components (i.e., those of the form           C
48     !  mu.grad(u)                                                          C
49     !                                                                      C
50     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
51           SUBROUTINE CALC_TAU_U_G(lTAU_U_G, lctau_u_g)
52     
53     ! Modules
54     !---------------------------------------------------------------------//
55           USE param, only: dimension_3
56           USE param1, only: zero, half
57           USE constant
58           USE physprop
59           USE fldvar
60           USE visc_g
61           USE run
62           USE toleranc, only: dil_ep_s
63           USE geometry
64           USE indices
65           USE is
66           USE compar
67           USE sendrecv
68           USE functions
69           USE cutcell
70           IMPLICIT NONE
71     
72     ! Dummy arguments
73     !---------------------------------------------------------------------//
74     ! TAU_U_g
75           DOUBLE PRECISION, INTENT(OUT) :: lTAU_U_g(DIMENSION_3)
76     ! cTAU_U_g
77           DOUBLE PRECISION, INTENT(OUT) :: lcTAU_U_g(DIMENSION_3)
78     
79     ! Local variables
80     !---------------------------------------------------------------------//
81     ! Indices
82           INTEGER :: I, J, K, IP, JM, KM
83           INTEGER :: IJK, IJKE, IJKN, IJKS, IJKT, IJKB
84           INTEGER :: IJKNE, IJKSE, IJKTE, IJKBE
85           INTEGER :: IPJK, IMJK, IJMK, IJKM
86           INTEGER :: IPJMK, IPJKM
87     ! Average volume fraction
88           DOUBLE PRECISION :: EPGA
89     ! Average viscosity
90           DOUBLE PRECISION :: MU_gte, MU_gbe, MUGA
91     ! Average dW/Xdz
92           DOUBLE PRECISION :: dWoXdz
93     ! Source terms (Surface)
94           DOUBLE PRECISION :: Sbv, Ssx, Ssy, Ssz
95     ! Source terms (Volumetric)
96           DOUBLE PRECISION :: Vtzb
97     !---------------------------------------------------------------------//
98     
99           IF((.NOT.CARTESIAN_GRID).OR.(CG_SAFE_MODE(3)==1)) THEN
100     
101     !$omp  parallel do default(none) &
102     !$omp  private(I, J, K, IP, JM, KM, IJK,                               &
103     !$omp          IJKE, IJKN, IJKS, IJKT, IJKB,                           &
104     !$omp          IJKNE, IJKSE, IJKTE, IJKBE,                             &
105     !$omp          IPJK, IMJK, IJMK, IJKM, IPJMK, IPJKM,                   &
106     !$omp          EPGA, MUGA, MU_GTE, MU_GBE, SBV, SSX, SSY, SSZ,         &
107     !$omp          DWOXDZ, VTZB)                                           &
108     !$omp  shared(ijkstart3, ijkend3, i_of, j_of, k_of, ip1, jm1, km1,     &
109     !$omp         do_k, cylindrical, ltau_u_g, lctau_u_g,                  &
110     !$omp         ep_g, epmu_gt, lambda_gt, trd_g, u_g, v_g, w_g,          &
111     !$omp         axy_u, axz_u, ayz, ayz_u, vol_u,                         &
112     !$omp         ox_e, odx_e, ox, odx, odz, eplambda_gt)
113              DO IJK = IJKSTART3, IJKEND3
114                 I = I_OF(IJK)
115                 IJKE = EAST_OF(IJK)
116                 EPGA = AVG_X(EP_G(IJK),EP_G(IJKE),I)
117                 IF (.NOT.IP_AT_E(IJK) .AND. EPGA>DIL_EP_S) THEN
118                    J = J_OF(IJK)
119                    K = K_OF(IJK)
120                    IP = IP1(I)
121                    JM = JM1(J)
122                    KM = KM1(K)
123     
124                    IPJK = IP_OF(IJK)
125                    IMJK = IM_OF(IJK)
126                    IJMK = JM_OF(IJK)
127                    IJKM = KM_OF(IJK)
128                    IPJMK = JM_OF(IPJK)
129                    IPJKM = IP_OF(IJKM)
130     
131                    IJKN = NORTH_OF(IJK)
132                    IJKNE = EAST_OF(IJKN)
133                    IJKS = SOUTH_OF(IJK)
134                    IJKSE = EAST_OF(IJKS)
135                    IJKT = TOP_OF(IJK)
136                    IJKTE = EAST_OF(IJKT)
137                    IJKB = BOTTOM_OF(IJK)
138                    IJKBE = EAST_OF(IJKB)
139     
140     
141     ! Surface forces at i+1/2, j, k
142     ! bulk viscosity term
143     ! combines part of 1/x d/dx (x.tau_xx) xdxdydz and -tau_zz/x xdxdydz =>
144     ! combines 1/x d/dx (x.lambda.trcD) xdxdydz - (lambda/x.trcD) xdxdydz =>
145     !              d/dx (lambda.trcD) xdxdydz
146     ! delta (lambda.trcD)Ap |E-W : at (i+1 - i-1), j, k
147                    SBV = (EPLAMBDA_GT(IJKE)*TRD_G(IJKE)-&
148                           EPLAMBDA_GT(IJK)*TRD_G(IJK))*AYZ(IJK)
149     
150     ! shear stress terms at i+1/2, j, k
151     ! part of 1/x d/dx(x.tau_xx) xdxdydz =>
152     !         1/x d/dx (x.mu.du/dx) xdxdydz =>
153     ! delta (mu du/dx)Ayz |E-W : at (i+1 - i-1), j, k
154                    SSX = EPMU_GT(IJKE)*(U_G(IPJK)-U_G(IJK))*ODX(IP)*&
155                             AYZ_U(IJK) - &
156                          EPMU_GT(IJK)*(U_G(IJK)-U_G(IMJK))*ODX(I)*&
157                             AYZ_U(IMJK)
158     
159     ! part of d/dy (tau_xy) xdxdydz =>
160     !         d/dy (mu.dv/dx) xdxdydz =>
161     ! delta (mu.dv/dx)Axz |N-S : at i+1/2, (j+1/2 - j-1/2), k
162                    SSY = AVG_X_H(AVG_Y_H(EPMU_GT(IJK),EPMU_GT(IJKN),J),&
163                                  AVG_Y_H(EPMU_GT(IJKE),EPMU_GT(IJKNE),J),I)*&
164                             (V_G(IPJK)-V_G(IJK))*ODX_E(I)*AXZ_U(IJK) - &
165                          AVG_X_H(AVG_Y_H(EPMU_GT(IJKS),EPMU_GT(IJK),JM),&
166                                  AVG_Y_H(EPMU_GT(IJKSE),EPMU_GT(IJKE),JM),I)*&
167                             (V_G(IPJMK)-V_G(IJMK))*ODX_E(I)*AXZ_U(IJMK)
168     
169     ! part of 1/x d/dz (tau_xz) xdxdydz =>
170     !         1/x d/dz (mu.dw/dx) xdxdydz =>
171     ! delta (mu.dw/dx)Axy |T-B : at i+1/2, j, (k+1/2 - k-1/2)
172                    MU_GTE = AVG_X_H(AVG_Z_H(EPMU_GT(IJK),EPMU_GT(IJKT),K),&
173                                     AVG_Z_H(EPMU_GT(IJKE),EPMU_GT(IJKTE),K),I)
174                    MU_GBE = AVG_X_H(AVG_Z_H(EPMU_GT(IJKB),EPMU_GT(IJK),KM),&
175                                     AVG_Z_H(EPMU_GT(IJKBE),EPMU_GT(IJKE),KM),I)
176                    SSZ = MU_GTE*(W_G(IPJK)-W_G(IJK))*ODX_E(I)*&
177                             AXY_U(IJK) - &
178                          MU_GBE*(W_G(IPJKM)-W_G(IJKM))*ODX_E(I)*&
179                             AXY_U(IJKM)
180     
181     
182     ! Special terms for cylindrical coordinates
183                    IF (CYLINDRICAL) THEN
184     ! part of 1/x d/dz (tau_xz) xdxdydz =>
185     !         1/x d/dz (mu.(-w/x)) xdxdydz =>
186     ! delta (mu(-w/x))Axy |T-B : at i+1/2, j, (k+1/2- k-1/2)
187                       SSZ = SSZ - ( &
188                             MU_GTE*(HALF*(W_G(IPJK)+W_G(IJK))*OX_E(I))*&
189                                AXY_U(IJK)-&
190                             MU_GBE*(HALF*(W_G(IPJKM)+W_G(IJKM))*OX_E(I))*&
191                                AXY_U(IJKM) )
192     
193     ! part of -tau_zz/x xdxdydz =>
194     !         -(2.mu/x).(1/x).dw/dz xdxdydz
195     ! delta (2.mu/x.1/x.dw/dz)Vp |p : at i+1/2, j, k
196                       MUGA = AVG_X(EPMU_GT(IJK),EPMU_GT(IJKE),I)
197                       DWOXDZ = HALF*((W_G(IJK)-W_G(IJKM))*OX(I)*ODZ(K)+&
198                                      (W_G(IPJK)-W_G(IPJKM))*OX(IP)*ODZ(K))
199                       VTZB = -2.d0*MUGA*OX_E(I)*DWOXDZ
200                    ELSE
201                       VTZB = ZERO
202                    ENDIF
203     
204     ! Add the terms
205                    lTAU_U_G(IJK) = SBV + SSX + SSY + SSZ + VTZB*VOL_U(IJK)
206     
207     ! Also calculate and store the full gas phase viscous stress tensor
208                    CALL GET_FULL_TAU_U_G(IJK, ltau_u_g, lctau_u_g)
209     
210                 ELSE
211                    lTAU_U_G(IJK) = ZERO
212                    lctau_u_g(IJK) = zero
213                 ENDIF   ! end if (.NOT.IP_AT_E(IJK) .AND. EPGA>DIL_EP_S)
214              ENDDO   ! end do ijk
215     !$omp end parallel do
216     
217           ELSE
218     ! if cartesian grid
219              CALL CALC_CG_TAU_U_G(lTAU_U_G, lctau_u_g)
220           ENDIF
221     
222           call send_recv(ltau_u_g,2)
223           call send_recv(lctau_u_g,2)
224     
225           RETURN
226     
227         CONTAINS
228     
229           INCLUDE 'fun_avg.inc'
230     
231         END SUBROUTINE CALC_TAU_U_G
232     
233     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
234     !                                                                      C
235     !  Subroutine: CALC_CG_Tau_U_g                                         C
236     !  Purpose: Cross terms in the gradient of stress in U_g momentum      C
237     !  based on cartesian grid cut cell.                                   C
238     !                                                                      C
239     !  Author: Jeff Dietiker                              Date: 01-Jul-09  C
240     !                                                                      C
241     !                                                                      C
242     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
243           SUBROUTINE CALC_CG_TAU_U_G(lTAu_U_G, lctau_u_g)
244     
245     ! Modules
246     !---------------------------------------------------------------------//
247     
248           USE bc
249           USE compar
250           USE constant
251           USE cutcell
252           USE fldvar
253           USE fun_avg
254           USE functions
255           USE geometry
256           USE indices
257           USE is
258           USE param, only: dimension_3
259           USE param1, only: zero, half, undefined
260           USE physprop
261           USE quadric
262           USE run
263           USE rxns
264           USE toleranc, only: dil_ep_s
265           USE visc_g
266     
267           IMPLICIT NONE
268     ! Dummy arguments
269     !---------------------------------------------------------------------//
270     ! TAU_U_g
271           DOUBLE PRECISION, INTENT(OUT) :: lTAU_U_g(DIMENSION_3)
272     ! cTAU_U_g
273           DOUBLE PRECISION, INTENT(OUT) :: lcTAU_U_g(DIMENSION_3)
274     
275     ! Local variables
276     !---------------------------------------------------------------------//
277     ! Indices
278           INTEGER :: I, J, K, IP, JM, KM
279           INTEGER :: IJK, IJKE, IJKN, IJKS, IJKT, IJKB
280           INTEGER :: IJKNE, IJKSE, IJKTE, IJKBE
281           INTEGER :: IPJK, IMJK, IJMK, IJKM
282           INTEGER :: IPJMK, IPJKM
283     ! Average volume fraction
284           DOUBLE PRECISION :: EPGA
285     ! Average viscosity
286           DOUBLE PRECISION :: MU_gte, MU_gbe
287     ! Source terms (Surface)
288           DOUBLE PRECISION :: Sbv, Ssx, Ssy, Ssz
289     ! Cartesian grid variables
290           DOUBLE PRECISION :: DEL_H, Nx, Ny, Nz
291           LOGICAL :: V_NODE_AT_NE, V_NODE_AT_NW, V_NODE_AT_SE, V_NODE_AT_SW
292           LOGICAL :: W_NODE_AT_TE, W_NODE_AT_TW, W_NODE_AT_BE, W_NODE_AT_BW
293           DOUBLE PRECISION :: dvdx_at_N, dvdx_at_S
294           DOUBLE PRECISION :: dwdx_at_T, dwdx_at_B
295           DOUBLE PRECISION :: Xi, Yi, Zi, Vi, Wi, Sx, Sy, Sz
296           DOUBLE PRECISION :: MU_GT_CUT, SSY_CUT, SSZ_CUT
297           DOUBLE PRECISION :: UW_g, VW_g, WW_g
298           INTEGER :: BCV
299           INTEGER :: BCT
300     !---------------------------------------------------------------------//
301     
302     !$omp  parallel do default(none) &
303     !$omp  private(I, J, K, JM, KM, IP, IJK,                               &
304     !$omp          IJKE, IJKN, IJKS, IJKT, IJKB,                           &
305     !$omp          IJKNE, IJKSE, IJKTE, IJKBE,                             &
306     !$omp          IPJK, IMJK, IJMK, IJKM, IPJMK, IPJKM,                   &
307     !$omp          EPGA, MU_GTE, MU_GBE, SBV, SSX, SSY, SSZ,               &
308     !$omp          BCV, BCT, BC_TYPE_ENUM, NOC_UG, uw_g, vw_g, ww_g,       &
309     !$omp          del_h, nx, ny, nz, xi, yi, zi, sx, sy, sz, vi, wi,      &
310     !$omp          v_node_at_ne, v_node_at_sw, v_node_at_se, v_node_at_nw, &
311     !$omp          w_node_at_bw, w_node_at_be, w_node_at_tw, w_node_at_te, &
312     !$omp          dvdx_at_S, dvdx_at_N, dwdx_at_B, dwdx_at_T,             &
313     !$omp          mu_gt_cut, cut_tau_ug, ssy_cut, ssz_cut)                &
314     !$omp  shared(ijkstart3, ijkend3, i_of, j_of, k_of,  ip1, jm1, km1,    &
315     !$omp         do_k, ltau_u_g, lctau_u_g,                               &
316     !$omp         ep_g, epmu_gt, lambda_gt, trd_g, u_g, v_g, w_g,          &
317     !$omp         axy_u, axz_u, ayz, ayz_u, vol,                           &
318     !$omp         bc_type, bc_u_id, bc_uw_g, bc_vw_g, bc_ww_g, bc_hw_g,    &
319     !$omp         oneodx_e_u, oneodx_e_v, oneodx_e_w,                      &
320     !$omp         x_u, y_u, z_u, x_v, y_v, z_v, x_w, y_w, z_w,             &
321     !$omp         wall_v_at, wall_w_at, area_u_cut, cut_u_cell_at,         &
322     !$omp         blocked_v_cell_at, blocked_w_cell_at, eplambda_gt)
323     
324           DO IJK = IJKSTART3, IJKEND3
325              I = I_OF(IJK)
326              IJKE = EAST_OF(IJK)
327              EPGA = AVG_X(EP_G(IJK),EP_G(IJKE),I)
328              IF ( .NOT.IP_AT_E(IJK) .AND. EPGA>DIL_EP_S) THEN
329                 J = J_OF(IJK)
330                 K = K_OF(IJK)
331                 IP = IP1(I)
332                 JM = JM1(J)
333                 KM = KM1(K)
334     
335                 IPJK = IP_OF(IJK)
336                 IMJK = IM_OF(IJK)
337                 IJMK = JM_OF(IJK)
338                 IJKM = KM_OF(IJK)
339                 IPJMK = JM_OF(IPJK)
340                 IPJKM = IP_OF(IJKM)
341     
342                 IJKN = NORTH_OF(IJK)
343                 IJKNE = EAST_OF(IJKN)
344                 IJKS = SOUTH_OF(IJK)
345                 IJKSE = EAST_OF(IJKS)
346                 IJKT = TOP_OF(IJK)
347                 IJKTE = EAST_OF(IJKT)
348                 IJKB = BOTTOM_OF(IJK)
349                 IJKBE = EAST_OF(IJKB)
350     
351     
352     ! bulk viscosity term
353                 SBV = (EPLAMBDA_GT(IJKE)*TRD_G(IJKE)) * AYZ_U(IJK) - &
354                       (EPLAMBDA_GT(IJK) *TRD_G(IJK) ) * AYZ_U(IMJK)
355     
356     ! shear stress terms
357                 IF(.NOT.CUT_U_CELL_AT(IJK)) THEN
358                    SSX = EPMU_GT(IJKE)*(U_G(IPJK)-U_G(IJK))*&
359                             ONEoDX_E_U(IJK)*AYZ_U(IJK) - &
360                          EPMU_GT(IJK)*(U_G(IJK)-U_G(IMJK))*&
361                             ONEoDX_E_U(IMJK)*AYZ_U(IMJK)
362     
363                    SSY = AVG_X_H(AVG_Y_H(EPMU_GT(IJK),EPMU_GT(IJKN),J),&
364                                  AVG_Y_H(EPMU_GT(IJKE),EPMU_GT(IJKNE),J),I)*&
365                          (V_G(IPJK)-V_G(IJK))*ONEoDX_E_V(IJK)*AXZ_U(IJK) - &
366                          AVG_X_H(AVG_Y_H(EPMU_GT(IJKS),EPMU_GT(IJK),JM),&
367                                  AVG_Y_H(EPMU_GT(IJKSE),EPMU_GT(IJKE),JM),I)*&
368                          (V_G(IPJMK)-V_G(IJMK))*ONEoDX_E_V(IJMK)*AXZ_U(IJMK)
369     
370                    IF(DO_K) THEN
371                       MU_GTE = AVG_X_H(AVG_Z_H(EPMU_GT(IJK),EPMU_GT(IJKT),K),&
372                                          AVG_Z_H(EPMU_GT(IJKE),EPMU_GT(IJKTE),K),I)
373                       MU_GBE = AVG_X_H(AVG_Z_H(EPMU_GT(IJKB),EPMU_GT(IJK),KM),&
374                                        AVG_Z_H(EPMU_GT(IJKBE),EPMU_GT(IJKE),KM),I)
375                       SSZ = MU_GTE*(W_G(IPJK)-W_G(IJK))*&
376                                ONEoDX_E_W(IJK)*AXY_U(IJK) - &
377                             MU_GBE*(W_G(IPJKM)-W_G(IJKM))*&
378                                ONEoDX_E_W(IJKM)*AXY_U(IJKM)
379                    ELSE
380                       SSZ = ZERO
381                    ENDIF
382     
383     ! cut cell modifications
384     !---------------------------------------------------------------------//
385                 ELSE
386     
387                    BCV = BC_U_ID(IJK)
388                    IF(BCV > 0 ) THEN
389                       BCT = BC_TYPE_ENUM(BCV)
390                    ELSE
391                       BCT = NONE
392                    ENDIF
393     
394                    SELECT CASE (BCT)
395                       CASE (CG_NSW,CG_MI)
396                          CUT_TAU_UG = .TRUE.
397                          NOC_UG     = .TRUE.
398                          UW_g = ZERO
399                          VW_g = ZERO
400                          WW_g = ZERO
401                       CASE (CG_FSW)
402                          CUT_TAU_UG = .FALSE.
403                          NOC_UG     = .FALSE.
404                          UW_g = ZERO
405                          VW_g = ZERO
406                          WW_g = ZERO
407                       CASE(CG_PSW)
408                          IF(BC_HW_G(BC_U_ID(IJK))==UNDEFINED) THEN   ! same as NSW
409                             CUT_TAU_UG = .TRUE.
410                             NOC_UG     = .TRUE.
411                             UW_g = BC_UW_G(BCV)
412                             VW_g = BC_VW_G(BCV)
413                             WW_g = BC_WW_G(BCV)
414                          ELSEIF(BC_HW_G(BC_U_ID(IJK))==ZERO) THEN   ! same as FSW
415                             CUT_TAU_UG = .FALSE.
416                             NOC_UG     = .FALSE.
417                          ELSE                              ! partial slip
418                             CUT_TAU_UG = .FALSE.
419                             NOC_UG     = .FALSE.
420                             UW_g = ZERO
421                             VW_g = ZERO
422                             WW_g = ZERO
423                          ENDIF
424                       CASE (NONE)
425                          lTAU_U_G(IJK) = ZERO
426                          CYCLE
427                    END SELECT
428     
429                    IF(CUT_TAU_UG) THEN
430                       MU_GT_CUT = (VOL(IJK)*EPMU_GT(IJK) + &
431                                    VOL(IPJK)*EPMU_GT(IJKE))/&
432                                   (VOL(IJK) + VOL(IPJK))
433                    ELSE
434                       MU_GT_CUT = ZERO
435                    ENDIF
436     
437     ! SSX:
438                    CALL GET_DEL_H(IJK,'U_MOMENTUM', X_U(IJK), &
439                       Y_U(IJK), Z_U(IJK), Del_H, Nx, Ny, Nz)
440     
441                    SSX = EPMU_GT(IJKE)*(U_G(IPJK)-U_G(IJK))*&
442                             ONEoDX_E_U(IJK)*AYZ_U(IJK) - &
443                          EPMU_GT(IJK)*(U_G(IJK)-U_G(IMJK))*&
444                             ONEoDX_E_U(IMJK)*AYZ_U(IMJK) - &
445                          MU_GT_CUT * (U_g(IJK) - UW_g) / DEL_H * &
446                             (Nx**2) * Area_U_CUT(IJK)
447     
448     ! SSY:
449                    V_NODE_AT_NE = ((.NOT.BLOCKED_V_CELL_AT(IPJK)).AND.&
450                                    (.NOT.WALL_V_AT(IPJK)))
451                    V_NODE_AT_NW = ((.NOT.BLOCKED_V_CELL_AT(IJK)).AND.&
452                                    (.NOT.WALL_V_AT(IJK)))
453                    V_NODE_AT_SE = ((.NOT.BLOCKED_V_CELL_AT(IPJMK)).AND.&
454                                    (.NOT.WALL_V_AT(IPJMK)))
455                    V_NODE_AT_SW = ((.NOT.BLOCKED_V_CELL_AT(IJMK)).AND.&
456                                    (.NOT.WALL_V_AT(IJMK)))
457     
458                    IF(V_NODE_AT_NE.AND.V_NODE_AT_NW) THEN
459                       Vi = HALF * (V_G(IPJK) + V_G(IJK))
460                       Xi = HALF * (X_V(IPJK) + X_V(IJK))
461                       Yi = HALF * (Y_V(IPJK) + Y_V(IJK))
462                       Zi = HALF * (Z_V(IPJK) + Z_V(IJK))
463                       Sx = X_V(IPJK) - X_V(IJK)
464                       Sy = Y_V(IPJK) - Y_V(IJK)
465                       Sz = Z_V(IPJK) - Z_V(IJK)
466                       CALL GET_DEL_H(IJK,'U_MOMENTUM',Xi, Yi, Zi, &
467                          Del_H, Nx, Ny, Nz)
468     
469                       dvdx_at_N = (V_G(IPJK)-V_G(IJK)) * ONEoDX_E_V(IJK)
470                       IF(NOC_UG) dvdx_at_N = dvdx_at_N - ( (Vi-VW_g)*&
471                          ONEoDX_E_V(IJK)/DEL_H * (Sy*Ny+Sz*Nz) )
472                    ELSE
473                       dvdx_at_N =  ZERO
474                    ENDIF
475     
476                    IF(V_NODE_AT_SE.AND.V_NODE_AT_SW) THEN
477                       Vi = HALF * (V_G(IPJMK) + V_G(IJMK))
478                       Xi = HALF * (X_V(IPJMK) + X_V(IJMK))
479                       Yi = HALF * (Y_V(IPJMK) + Y_V(IJMK))
480                       Zi = HALF * (Z_V(IPJMK) + Z_V(IJMK))
481                       Sx = X_V(IPJMK) - X_V(IJMK)
482                       Sy = Y_V(IPJMK) - Y_V(IJMK)
483                       Sz = Z_V(IPJMK) - Z_V(IJMK)
484                       CALL GET_DEL_H(IJK,'U_MOMENTUM', Xi, Yi, Zi, &
485                          Del_H, Nx, Ny, Nz)
486     
487                       dvdx_at_S = (V_G(IPJMK)-V_G(IJMK)) * ONEoDX_E_V(IJMK)
488                       IF(NOC_UG) dvdx_at_S = dvdx_at_S - ( (Vi-VW_g)*&
489                          ONEoDX_E_V(IJMK)/DEL_H * (Sy*Ny+Sz*Nz) )
490                    ELSE
491                       dvdx_at_S =  ZERO
492                    ENDIF
493     
494                    IF(V_NODE_AT_NW) THEN
495                       CALL GET_DEL_H(IJK,'U_MOMENTUM', X_V(IJK), &
496                          Y_V(IJK), Z_V(IJK), Del_H, Nx, Ny, Nz)
497                       SSY_CUT = -MU_GT_CUT * (V_G(IJK)-VW_g)/DEL_H * &
498                          (Nx*Ny) * Area_U_CUT(IJK)
499                    ELSE
500                       SSY_CUT =  ZERO
501                    ENDIF
502     
503                    SSY = AVG_X_H(AVG_Y_H(EPMU_GT(IJK),EPMU_GT(IJKN),J),&
504                                  AVG_Y_H(EPMU_GT(IJKE),EPMU_GT(IJKNE),J),I)*&
505                             dvdx_at_N*AXZ_U(IJK) - &
506                          AVG_X_H(AVG_Y_H(EPMU_GT(IJKS),EPMU_GT(IJK),JM),&
507                                  AVG_Y_H(EPMU_GT(IJKSE),EPMU_GT(IJKE),JM),I)*&
508                             dvdx_at_S*AXZ_U(IJMK) + SSY_CUT
509     
510     ! SSZ:
511                    IF(DO_K) THEN
512                       W_NODE_AT_TE = ((.NOT.BLOCKED_W_CELL_AT(IPJK)).AND.&
513                                       (.NOT.WALL_W_AT(IPJK)))
514                       W_NODE_AT_TW = ((.NOT.BLOCKED_W_CELL_AT(IJK)).AND.&
515                                       (.NOT.WALL_W_AT(IJK)))
516                       W_NODE_AT_BE = ((.NOT.BLOCKED_W_CELL_AT(IPJKM)).AND.&
517                                       (.NOT.WALL_W_AT(IPJKM)))
518                       W_NODE_AT_BW = ((.NOT.BLOCKED_W_CELL_AT(IJKM)).AND.&
519                                       (.NOT.WALL_W_AT(IJKM)))
520     
521                       IF(W_NODE_AT_TE.AND.W_NODE_AT_TW) THEN
522                          Wi = HALF * (W_G(IPJK) + W_G(IJK))
523                          Xi = HALF * (X_W(IPJK) + X_W(IJK))
524                          Yi = HALF * (Y_W(IPJK) + Y_W(IJK))
525                          Zi = HALF * (Z_W(IPJK) + Z_W(IJK))
526                          Sx = X_W(IPJK) - X_W(IJK)
527                          Sy = Y_W(IPJK) - Y_W(IJK)
528                          Sz = Z_W(IPJK) - Z_W(IJK)
529                          CALL GET_DEL_H(IJK, 'U_MOMENTUM', Xi, Yi, Zi, &
530                             Del_H, Nx, Ny, Nz)
531     
532                          dwdx_at_T = (W_G(IPJK)-W_G(IJK)) * ONEoDX_E_W(IJK)
533                          IF(NOC_UG) dwdx_at_T = dwdx_at_T - ( (Wi-WW_g) * &
534                             ONEoDX_E_W(IJK)/DEL_H * (Sy*Ny+Sz*Nz) )
535                       ELSE
536                          dwdx_at_T =  ZERO
537                       ENDIF
538     
539                       IF(W_NODE_AT_BE.AND.W_NODE_AT_BW) THEN
540                          Wi = HALF * (W_G(IPJKM) + W_G(IJKM))
541                          Xi = HALF * (X_W(IPJKM) + X_W(IJKM))
542                          Yi = HALF * (Y_W(IPJKM) + Y_W(IJKM))
543                          Zi = HALF * (Z_W(IPJKM) + Z_W(IJKM))
544                          Sx = X_W(IPJKM) - X_W(IJKM)
545                          Sy = Y_W(IPJKM) - Y_W(IJKM)
546                          Sz = Z_W(IPJKM) - Z_W(IJKM)
547     
548                          CALL GET_DEL_H(IJK,'U_MOMENTUM', Xi, Yi, Zi,&
549                             Del_H, Nx, Ny, Nz)
550     
551                          dwdx_at_B = (W_G(IPJKM)-W_G(IJKM)) * ONEoDX_E_W(IJKM)
552                          IF(NOC_UG) dwdx_at_B = dwdx_at_B - ( (Wi-WW_g) *&
553                             ONEoDX_E_W(IJKM)/DEL_H * (Sy*Ny+Sz*Nz))
554                       ELSE
555                          dwdx_at_B =  ZERO
556                       ENDIF
557     
558                       IF(W_NODE_AT_TW) THEN
559                          CALL GET_DEL_H(IJK,'U_MOMENTUM', X_W(IJK), &
560                             Y_W(IJK), Z_W(IJK), Del_H, Nx, Ny, Nz)
561                          SSZ_CUT = -MU_GT_CUT * (W_G(IJK)-WW_g)/DEL_H * &
562                             (Nx*Nz) * Area_U_CUT(IJK)
563                       ELSE
564                          SSZ_CUT =  ZERO
565                       ENDIF
566     
567                       MU_GTE = AVG_X_H(AVG_Z_H(EPMU_GT(IJK),EPMU_GT(IJKT),K),&
568                                        AVG_Z_H(EPMU_GT(IJKE),EPMU_GT(IJKTE),K),I)
569                       MU_GBE = AVG_X_H(AVG_Z_H(EPMU_GT(IJKB),EPMU_GT(IJK),KM),&
570                                        AVG_Z_H(EPMU_GT(IJKBE),EPMU_GT(IJKE),KM),I)
571                       SSZ = MU_GTE*dwdx_at_T*AXY_U(IJK) - &
572                             MU_GBE*dwdx_at_B*AXY_U(IJKM)  + SSZ_CUT
573                    ELSE
574                      SSZ = ZERO
575                    ENDIF ! DO_K
576     
577                 ENDIF  ! end if/else cut_u_cell_at
578     
579     ! Add the terms
580                 lTAU_U_G(IJK) = SBV + SSX + SSY + SSZ
581     
582     ! Also calculate and store the full gas phase viscous stress tensor
583                 CALL GET_FULL_TAU_U_G(IJK, ltau_u_g, lctau_u_g)
584              ELSE
585                 lTAU_U_G(IJK) = ZERO
586                 lctau_u_g(IJK) = zero
587              ENDIF
588           ENDDO
589     !$omp end parallel do
590     
591           RETURN
592           END SUBROUTINE CALC_CG_TAU_U_G
593     
594     
595     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
596     !                                                                      C
597     !  Subroutine: GET_FULL_Tau_U_g                                        C
598     !  Purpose: Calculate the divergence of the complete gas phase stress  C
599     !  tensor including all terms in the u-direction.                      C
600     !                                                                      C
601     !                                                                      C
602     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
603           SUBROUTINE GET_FULL_TAU_U_G(IJK, ltau_u_g, lctau_u_g)
604     
605     ! Modules
606     !---------------------------------------------------------------------//
607           USE fldvar, only: u_g
608     
609           USE functions, only: east_of, flow_at_e
610           USE functions, only: im_of, jm_of, km_of
611           USE functions, only: ip_of, jp_of, kp_of
612     
613           USE fun_avg, only: avg_x
614     
615           USE geometry, only: ox_e, do_k, cylindrical, vol_u
616           USE indices, only: i_of
617     
618           USE param
619           USE param1, only: zero
620           USE visc_g, only: mu_gt, df_gu
621           IMPLICIT NONE
622     
623     ! Dummy arguments
624     !---------------------------------------------------------------------//
625     ! ijk index
626           INTEGER, INTENT(IN) :: IJK
627     ! TAU_U_g
628           DOUBLE PRECISION, INTENT(IN) :: lTAU_U_g(DIMENSION_3)
629     ! cTAU_U_g
630           DOUBLE PRECISION, INTENT(INOUT) :: lcTAU_U_g(DIMENSION_3)
631     
632     ! Local Variables
633     !---------------------------------------------------------------------//
634     ! indices
635           INTEGER :: I, IJKE
636           INTEGER :: IPJK, IJPK, IJKP, IMJK, IJMK, IJKM
637     ! average viscosity
638           DOUBLE PRECISION :: MUGA
639     ! source terms
640           DOUBLE PRECISION :: SSX, SSY, SSZ
641     ! cylindrical terms
642           DOUBLE PRECISION :: vtza
643     !---------------------------------------------------------------------//
644     
645           IPJK = IP_OF(IJK)
646           IJPK = JP_OF(IJK)
647           IJKP = KP_OF(IJK)
648           IMJK = IM_OF(IJK)
649           IJMK = JM_OF(IJK)
650           IJKM = KM_OF(IJK)
651     
652     ! additional terms (found in source_u_g)
653           VTZA = ZERO
654           IF (CYLINDRICAL) THEN
655              I = I_OF(IJK)
656              IJKE = EAST_OF(IJK)
657     ! part of -tau_zz/x xdxdydz =>
658     !         -(2mu/x)*(u/x) xdxdydz =>
659     ! delta(-2.mu.u/x^2)V |p
660              MUGA = AVG_X(MU_GT(IJK),MU_GT(IJKE),I)
661              VTZA = 2.d0*MUGA*OX_E(I)*OX_E(I)*VOL_U(IJK)*U_G(IJK)
662           ENDIF
663     
664     ! convection terms: see conv_dif_u_g
665           SSX = ZERO
666           SSY = ZERO
667           SSZ = ZERO
668           IF (FLOW_AT_E(IJK)) THEN
669     ! part of 1/x d/dx (x.tau_xx) xdxdydz =>
670     !         1/x d/dx (x.mu.du/dx) xdxdydz =>
671     ! delta (mu.du/dx.Ayz) |E-W : at (i+1 - i-1), j, k
672              SSX = DF_GU(IJK,east)*(U_G(IPJK) - U_G(IJK)) - &
673                    DF_GU(IJK,west)*(U_G(IJK) - U_G(IJKM))
674     
675     ! part of d/dy (tau_xy) xdxdydz =>
676     !         d/dy (mu.du/dy) xdxdydz =>
677     ! delta (mu.du/dy.Axz) |N-S : at (i+1/2, j+1/2 - j-1/2, k)
678              SSY = DF_GU(IJK,north)*(U_G(IJPK)-U_G(IJK)) - &
679                    DF_GU(IJK,south)*(U_G(IJK)-U_G(IJMK))
680     
681              IF (DO_K) THEN
682     ! part of 1/x d/dz (tau_xz) xdxdydz =>
683     !         1/x d/dz (mu/x.du/dz) xdxdydz =>
684     ! delta (mu/x.du/dz.Axy) |T-B : at (i+1/2, j, k+1/2 - k-1/2)
685                 SSZ = DF_GU(IJK,top)*(U_G(IJKP)-U_G(IJK)) - &
686                       DF_GU(IJK,bottom)*(U_G(IJK)-U_G(IJKM))
687              ENDIF
688           ENDIF   ! end if flow_at_e
689     
690     ! Add the terms
691           lctau_u_g(IJK) = (lTAU_U_G(IJK) + SSX + SSY + SSZ + VTZA)
692     
693           RETURN
694           END SUBROUTINE GET_FULL_TAU_U_G
695