File: RELATIVE:/../../../mfix.git/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           USE param, only: dimension_3
248           USE param1, only: zero, half
249           USE constant
250           USE physprop
251           USE fldvar
252           USE visc_g
253           USE rxns
254           USE run
255           USE toleranc, only: dil_ep_s
256           USE geometry
257           USE indices
258           USE is
259           USE compar
260           USE fun_avg
261           USE functions
262     
263           USE bc
264           USE quadric
265           USE cutcell
266           IMPLICIT NONE
267     ! Dummy arguments
268     !---------------------------------------------------------------------//
269     ! TAU_U_g
270           DOUBLE PRECISION, INTENT(OUT) :: lTAU_U_g(DIMENSION_3)
271     ! cTAU_U_g
272           DOUBLE PRECISION, INTENT(OUT) :: lcTAU_U_g(DIMENSION_3)
273     
274     ! Local variables
275     !---------------------------------------------------------------------//
276     ! Indices
277           INTEGER :: I, J, K, IP, JM, KM
278           INTEGER :: IJK, IJKE, IJKN, IJKS, IJKT, IJKB
279           INTEGER :: IJKNE, IJKSE, IJKTE, IJKBE
280           INTEGER :: IPJK, IMJK, IJMK, IJKM
281           INTEGER :: IPJMK, IPJKM
282     ! Average volume fraction
283           DOUBLE PRECISION :: EPGA
284     ! Average viscosity
285           DOUBLE PRECISION :: MU_gte, MU_gbe
286     ! Source terms (Surface)
287           DOUBLE PRECISION :: Sbv, Ssx, Ssy, Ssz
288     ! Cartesian grid variables
289           DOUBLE PRECISION :: DEL_H, Nx, Ny, Nz
290           LOGICAL :: V_NODE_AT_NE, V_NODE_AT_NW, V_NODE_AT_SE, V_NODE_AT_SW
291           LOGICAL :: W_NODE_AT_TE, W_NODE_AT_TW, W_NODE_AT_BE, W_NODE_AT_BW
292           DOUBLE PRECISION :: dvdx_at_N, dvdx_at_S
293           DOUBLE PRECISION :: dwdx_at_T, dwdx_at_B
294           DOUBLE PRECISION :: Xi, Yi, Zi, Vi, Wi, Sx, Sy, Sz
295           DOUBLE PRECISION :: MU_GT_CUT, SSY_CUT, SSZ_CUT
296           DOUBLE PRECISION :: UW_g, VW_g, WW_g
297           INTEGER :: BCV
298           CHARACTER(LEN=9) :: BCT
299     !---------------------------------------------------------------------//
300     
301     !$omp  parallel do default(none) &
302     !$omp  private(I, J, K, JM, KM, IP, IJK,                               &
303     !$omp          IJKE, IJKN, IJKS, IJKT, IJKB,                           &
304     !$omp          IJKNE, IJKSE, IJKTE, IJKBE,                             &
305     !$omp          IPJK, IMJK, IJMK, IJKM, IPJMK, IPJKM,                   &
306     !$omp          EPGA, MU_GTE, MU_GBE, SBV, SSX, SSY, SSZ,               &
307     !$omp          BCV, BCT, NOC_UG, uw_g, vw_g, ww_g,                     &
308     !$omp          del_h, nx, ny, nz, xi, yi, zi, sx, sy, sz, vi, wi,      &
309     !$omp          v_node_at_ne, v_node_at_sw, v_node_at_se, v_node_at_nw, &
310     !$omp          w_node_at_bw, w_node_at_be, w_node_at_tw, w_node_at_te, &
311     !$omp          dvdx_at_S, dvdx_at_N, dwdx_at_B, dwdx_at_T,             &
312     !$omp          mu_gt_cut, cut_tau_ug, ssy_cut, ssz_cut)                &
313     !$omp  shared(ijkstart3, ijkend3, i_of, j_of, k_of,  ip1, jm1, km1,    &
314     !$omp         do_k, ltau_u_g, lctau_u_g,                               &
315     !$omp         ep_g, epmu_gt, lambda_gt, trd_g, u_g, v_g, w_g,          &
316     !$omp         axy_u, axz_u, ayz, ayz_u, vol,                           &
317     !$omp         bc_type, bc_u_id, bc_uw_g, bc_vw_g, bc_ww_g, bc_hw_g,    &
318     !$omp         oneodx_e_u, oneodx_e_v, oneodx_e_w,                      &
319     !$omp         x_u, y_u, z_u, x_v, y_v, z_v, x_w, y_w, z_w,             &
320     !$omp         wall_v_at, wall_w_at, area_u_cut, cut_u_cell_at,         &
321     !$omp         blocked_v_cell_at, blocked_w_cell_at, eplambda_gt)
322     
323           DO IJK = IJKSTART3, IJKEND3
324              I = I_OF(IJK)
325              IJKE = EAST_OF(IJK)
326              EPGA = AVG_X(EP_G(IJK),EP_G(IJKE),I)
327              IF ( .NOT.IP_AT_E(IJK) .AND. EPGA>DIL_EP_S) THEN
328                 J = J_OF(IJK)
329                 K = K_OF(IJK)
330                 IP = IP1(I)
331                 JM = JM1(J)
332                 KM = KM1(K)
333     
334                 IPJK = IP_OF(IJK)
335                 IMJK = IM_OF(IJK)
336                 IJMK = JM_OF(IJK)
337                 IJKM = KM_OF(IJK)
338                 IPJMK = JM_OF(IPJK)
339                 IPJKM = IP_OF(IJKM)
340     
341                 IJKN = NORTH_OF(IJK)
342                 IJKNE = EAST_OF(IJKN)
343                 IJKS = SOUTH_OF(IJK)
344                 IJKSE = EAST_OF(IJKS)
345                 IJKT = TOP_OF(IJK)
346                 IJKTE = EAST_OF(IJKT)
347                 IJKB = BOTTOM_OF(IJK)
348                 IJKBE = EAST_OF(IJKB)
349     
350     
351     ! bulk viscosity term
352                 SBV = (EPLAMBDA_GT(IJKE)*TRD_G(IJKE)) * AYZ_U(IJK) - &
353                       (EPLAMBDA_GT(IJK) *TRD_G(IJK) ) * AYZ_U(IMJK)
354     
355     ! shear stress terms
356                 IF(.NOT.CUT_U_CELL_AT(IJK)) THEN
357                    SSX = EPMU_GT(IJKE)*(U_G(IPJK)-U_G(IJK))*&
358                             ONEoDX_E_U(IJK)*AYZ_U(IJK) - &
359                          EPMU_GT(IJK)*(U_G(IJK)-U_G(IMJK))*&
360                             ONEoDX_E_U(IMJK)*AYZ_U(IMJK)
361     
362                    SSY = AVG_X_H(AVG_Y_H(EPMU_GT(IJK),EPMU_GT(IJKN),J),&
363                                  AVG_Y_H(EPMU_GT(IJKE),EPMU_GT(IJKNE),J),I)*&
364                          (V_G(IPJK)-V_G(IJK))*ONEoDX_E_V(IJK)*AXZ_U(IJK) - &
365                          AVG_X_H(AVG_Y_H(EPMU_GT(IJKS),EPMU_GT(IJK),JM),&
366                                  AVG_Y_H(EPMU_GT(IJKSE),EPMU_GT(IJKE),JM),I)*&
367                          (V_G(IPJMK)-V_G(IJMK))*ONEoDX_E_V(IJMK)*AXZ_U(IJMK)
368     
369                    IF(DO_K) THEN
370                       MU_GTE = AVG_X_H(AVG_Z_H(EPMU_GT(IJK),EPMU_GT(IJKT),K),&
371                                          AVG_Z_H(EPMU_GT(IJKE),EPMU_GT(IJKTE),K),I)
372                       MU_GBE = AVG_X_H(AVG_Z_H(EPMU_GT(IJKB),EPMU_GT(IJK),KM),&
373                                        AVG_Z_H(EPMU_GT(IJKBE),EPMU_GT(IJKE),KM),I)
374                       SSZ = MU_GTE*(W_G(IPJK)-W_G(IJK))*&
375                                ONEoDX_E_W(IJK)*AXY_U(IJK) - &
376                             MU_GBE*(W_G(IPJKM)-W_G(IJKM))*&
377                                ONEoDX_E_W(IJKM)*AXY_U(IJKM)
378                    ELSE
379                       SSZ = ZERO
380                    ENDIF
381     
382     ! cut cell modifications
383     !---------------------------------------------------------------------//
384                 ELSE
385     
386                    BCV = BC_U_ID(IJK)
387                    IF(BCV > 0 ) THEN
388                       BCT = BC_TYPE(BCV)
389                    ELSE
390                       BCT = 'NONE'
391                    ENDIF
392     
393                    SELECT CASE (BCT)
394                       CASE ('CG_NSW')
395                          CUT_TAU_UG = .TRUE.
396                          NOC_UG     = .TRUE.
397                          UW_g = ZERO
398                          VW_g = ZERO
399                          WW_g = ZERO
400                       CASE ('CG_FSW')
401                          CUT_TAU_UG = .FALSE.
402                          NOC_UG     = .FALSE.
403                          UW_g = ZERO
404                          VW_g = ZERO
405                          WW_g = ZERO
406                       CASE('CG_PSW')
407                          IF(BC_HW_G(BC_U_ID(IJK))==UNDEFINED) THEN   ! same as NSW
408                             CUT_TAU_UG = .TRUE.
409                             NOC_UG     = .TRUE.
410                             UW_g = BC_UW_G(BCV)
411                             VW_g = BC_VW_G(BCV)
412                             WW_g = BC_WW_G(BCV)
413                          ELSEIF(BC_HW_G(BC_U_ID(IJK))==ZERO) THEN   ! same as FSW
414                             CUT_TAU_UG = .FALSE.
415                             NOC_UG     = .FALSE.
416                          ELSE                              ! partial slip
417                             CUT_TAU_UG = .FALSE.
418                             NOC_UG     = .FALSE.
419                             UW_g = ZERO
420                             VW_g = ZERO
421                             WW_g = ZERO
422                          ENDIF
423                       CASE ('NONE')
424                          lTAU_U_G(IJK) = ZERO
425                          CYCLE
426                    END SELECT
427     
428                    IF(CUT_TAU_UG) THEN
429                       MU_GT_CUT = (VOL(IJK)*EPMU_GT(IJK) + &
430                                    VOL(IPJK)*EPMU_GT(IJKE))/&
431                                   (VOL(IJK) + VOL(IPJK))
432                    ELSE
433                       MU_GT_CUT = ZERO
434                    ENDIF
435     
436     ! SSX:
437                    CALL GET_DEL_H(IJK,'U_MOMENTUM', X_U(IJK), &
438                       Y_U(IJK), Z_U(IJK), Del_H, Nx, Ny, Nz)
439     
440                    SSX = EPMU_GT(IJKE)*(U_G(IPJK)-U_G(IJK))*&
441                             ONEoDX_E_U(IJK)*AYZ_U(IJK) - &
442                          EPMU_GT(IJK)*(U_G(IJK)-U_G(IMJK))*&
443                             ONEoDX_E_U(IMJK)*AYZ_U(IMJK) - &
444                          MU_GT_CUT * (U_g(IJK) - UW_g) / DEL_H * &
445                             (Nx**2) * Area_U_CUT(IJK)
446     
447     ! SSY:
448                    V_NODE_AT_NE = ((.NOT.BLOCKED_V_CELL_AT(IPJK)).AND.&
449                                    (.NOT.WALL_V_AT(IPJK)))
450                    V_NODE_AT_NW = ((.NOT.BLOCKED_V_CELL_AT(IJK)).AND.&
451                                    (.NOT.WALL_V_AT(IJK)))
452                    V_NODE_AT_SE = ((.NOT.BLOCKED_V_CELL_AT(IPJMK)).AND.&
453                                    (.NOT.WALL_V_AT(IPJMK)))
454                    V_NODE_AT_SW = ((.NOT.BLOCKED_V_CELL_AT(IJMK)).AND.&
455                                    (.NOT.WALL_V_AT(IJMK)))
456     
457                    IF(V_NODE_AT_NE.AND.V_NODE_AT_NW) THEN
458                       Vi = HALF * (V_G(IPJK) + V_G(IJK))
459                       Xi = HALF * (X_V(IPJK) + X_V(IJK))
460                       Yi = HALF * (Y_V(IPJK) + Y_V(IJK))
461                       Zi = HALF * (Z_V(IPJK) + Z_V(IJK))
462                       Sx = X_V(IPJK) - X_V(IJK)
463                       Sy = Y_V(IPJK) - Y_V(IJK)
464                       Sz = Z_V(IPJK) - Z_V(IJK)
465                       CALL GET_DEL_H(IJK,'U_MOMENTUM',Xi, Yi, Zi, &
466                          Del_H, Nx, Ny, Nz)
467     
468                       dvdx_at_N = (V_G(IPJK)-V_G(IJK)) * ONEoDX_E_V(IJK)
469                       IF(NOC_UG) dvdx_at_N = dvdx_at_N - ( (Vi-VW_g)*&
470                          ONEoDX_E_V(IJK)/DEL_H * (Sy*Ny+Sz*Nz) )
471                    ELSE
472                       dvdx_at_N =  ZERO
473                    ENDIF
474     
475                    IF(V_NODE_AT_SE.AND.V_NODE_AT_SW) THEN
476                       Vi = HALF * (V_G(IPJMK) + V_G(IJMK))
477                       Xi = HALF * (X_V(IPJMK) + X_V(IJMK))
478                       Yi = HALF * (Y_V(IPJMK) + Y_V(IJMK))
479                       Zi = HALF * (Z_V(IPJMK) + Z_V(IJMK))
480                       Sx = X_V(IPJMK) - X_V(IJMK)
481                       Sy = Y_V(IPJMK) - Y_V(IJMK)
482                       Sz = Z_V(IPJMK) - Z_V(IJMK)
483                       CALL GET_DEL_H(IJK,'U_MOMENTUM', Xi, Yi, Zi, &
484                          Del_H, Nx, Ny, Nz)
485     
486                       dvdx_at_S = (V_G(IPJMK)-V_G(IJMK)) * ONEoDX_E_V(IJMK)
487                       IF(NOC_UG) dvdx_at_S = dvdx_at_S - ( (Vi-VW_g)*&
488                          ONEoDX_E_V(IJMK)/DEL_H * (Sy*Ny+Sz*Nz) )
489                    ELSE
490                       dvdx_at_S =  ZERO
491                    ENDIF
492     
493                    IF(V_NODE_AT_NW) THEN
494                       CALL GET_DEL_H(IJK,'U_MOMENTUM', X_V(IJK), &
495                          Y_V(IJK), Z_V(IJK), Del_H, Nx, Ny, Nz)
496                       SSY_CUT = -MU_GT_CUT * (V_G(IJK)-VW_g)/DEL_H * &
497                          (Nx*Ny) * Area_U_CUT(IJK)
498                    ELSE
499                       SSY_CUT =  ZERO
500                    ENDIF
501     
502                    SSY = AVG_X_H(AVG_Y_H(EPMU_GT(IJK),EPMU_GT(IJKN),J),&
503                                  AVG_Y_H(EPMU_GT(IJKE),EPMU_GT(IJKNE),J),I)*&
504                             dvdx_at_N*AXZ_U(IJK) - &
505                          AVG_X_H(AVG_Y_H(EPMU_GT(IJKS),EPMU_GT(IJK),JM),&
506                                  AVG_Y_H(EPMU_GT(IJKSE),EPMU_GT(IJKE),JM),I)*&
507                             dvdx_at_S*AXZ_U(IJMK) + SSY_CUT
508     
509     ! SSZ:
510                    IF(DO_K) THEN
511                       W_NODE_AT_TE = ((.NOT.BLOCKED_W_CELL_AT(IPJK)).AND.&
512                                       (.NOT.WALL_W_AT(IPJK)))
513                       W_NODE_AT_TW = ((.NOT.BLOCKED_W_CELL_AT(IJK)).AND.&
514                                       (.NOT.WALL_W_AT(IJK)))
515                       W_NODE_AT_BE = ((.NOT.BLOCKED_W_CELL_AT(IPJKM)).AND.&
516                                       (.NOT.WALL_W_AT(IPJKM)))
517                       W_NODE_AT_BW = ((.NOT.BLOCKED_W_CELL_AT(IJKM)).AND.&
518                                       (.NOT.WALL_W_AT(IJKM)))
519     
520                       IF(W_NODE_AT_TE.AND.W_NODE_AT_TW) THEN
521                          Wi = HALF * (W_G(IPJK) + W_G(IJK))
522                          Xi = HALF * (X_W(IPJK) + X_W(IJK))
523                          Yi = HALF * (Y_W(IPJK) + Y_W(IJK))
524                          Zi = HALF * (Z_W(IPJK) + Z_W(IJK))
525                          Sx = X_W(IPJK) - X_W(IJK)
526                          Sy = Y_W(IPJK) - Y_W(IJK)
527                          Sz = Z_W(IPJK) - Z_W(IJK)
528                          CALL GET_DEL_H(IJK, 'U_MOMENTUM', Xi, Yi, Zi, &
529                             Del_H, Nx, Ny, Nz)
530     
531                          dwdx_at_T = (W_G(IPJK)-W_G(IJK)) * ONEoDX_E_W(IJK)
532                          IF(NOC_UG) dwdx_at_T = dwdx_at_T - ( (Wi-WW_g) * &
533                             ONEoDX_E_W(IJK)/DEL_H * (Sy*Ny+Sz*Nz) )
534                       ELSE
535                          dwdx_at_T =  ZERO
536                       ENDIF
537     
538                       IF(W_NODE_AT_BE.AND.W_NODE_AT_BW) THEN
539                          Wi = HALF * (W_G(IPJKM) + W_G(IJKM))
540                          Xi = HALF * (X_W(IPJKM) + X_W(IJKM))
541                          Yi = HALF * (Y_W(IPJKM) + Y_W(IJKM))
542                          Zi = HALF * (Z_W(IPJKM) + Z_W(IJKM))
543                          Sx = X_W(IPJKM) - X_W(IJKM)
544                          Sy = Y_W(IPJKM) - Y_W(IJKM)
545                          Sz = Z_W(IPJKM) - Z_W(IJKM)
546     
547                          CALL GET_DEL_H(IJK,'U_MOMENTUM', Xi, Yi, Zi,&
548                             Del_H, Nx, Ny, Nz)
549     
550                          dwdx_at_B = (W_G(IPJKM)-W_G(IJKM)) * ONEoDX_E_W(IJKM)
551                          IF(NOC_UG) dwdx_at_B = dwdx_at_B - ( (Wi-WW_g) *&
552                             ONEoDX_E_W(IJKM)/DEL_H * (Sy*Ny+Sz*Nz))
553                       ELSE
554                          dwdx_at_B =  ZERO
555                       ENDIF
556     
557                       IF(W_NODE_AT_TW) THEN
558                          CALL GET_DEL_H(IJK,'U_MOMENTUM', X_W(IJK), &
559                             Y_W(IJK), Z_W(IJK), Del_H, Nx, Ny, Nz)
560                          SSZ_CUT = -MU_GT_CUT * (W_G(IJK)-WW_g)/DEL_H * &
561                             (Nx*Nz) * Area_U_CUT(IJK)
562                       ELSE
563                          SSZ_CUT =  ZERO
564                       ENDIF
565     
566                       MU_GTE = AVG_X_H(AVG_Z_H(EPMU_GT(IJK),EPMU_GT(IJKT),K),&
567                                        AVG_Z_H(EPMU_GT(IJKE),EPMU_GT(IJKTE),K),I)
568                       MU_GBE = AVG_X_H(AVG_Z_H(EPMU_GT(IJKB),EPMU_GT(IJK),KM),&
569                                        AVG_Z_H(EPMU_GT(IJKBE),EPMU_GT(IJKE),KM),I)
570                       SSZ = MU_GTE*dwdx_at_T*AXY_U(IJK) - &
571                             MU_GBE*dwdx_at_B*AXY_U(IJKM)  + SSZ_CUT
572                    ELSE
573                      SSZ = ZERO
574                    ENDIF ! DO_K
575     
576                 ENDIF  ! end if/else cut_u_cell_at
577     
578     ! Add the terms
579                 lTAU_U_G(IJK) = SBV + SSX + SSY + SSZ
580     
581     ! Also calculate and store the full gas phase viscous stress tensor
582                 CALL GET_FULL_TAU_U_G(IJK, ltau_u_g, lctau_u_g)
583              ELSE
584                 lTAU_U_G(IJK) = ZERO
585                 lctau_u_g(IJK) = zero
586              ENDIF
587           ENDDO
588     !$omp end parallel do
589     
590           RETURN
591           END SUBROUTINE CALC_CG_TAU_U_G
592     
593     
594     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
595     !                                                                      C
596     !  Subroutine: GET_FULL_Tau_U_g                                        C
597     !  Purpose: Calculate the divergence of the complete gas phase stress  C
598     !  tensor including all terms in the u-direction.                      C
599     !                                                                      C
600     !                                                                      C
601     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
602           SUBROUTINE GET_FULL_TAU_U_G(IJK, ltau_u_g, lctau_u_g)
603     
604     ! Modules
605     !---------------------------------------------------------------------//
606           USE fldvar, only: u_g
607     
608           USE functions, only: east_of, flow_at_e
609           USE functions, only: im_of, jm_of, km_of
610           USE functions, only: ip_of, jp_of, kp_of
611     
612           USE fun_avg, only: avg_x
613     
614           USE geometry, only: ox_e, do_k, cylindrical, vol_u
615           USE indices, only: i_of
616     
617           USE matrix, only: e, w, n, s, t, b
618           USE param, only: dimension_3
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,E)*(U_G(IPJK) - U_G(IJK)) - &
673                    DF_GU(IJK,W)*(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,N)*(U_G(IJPK)-U_G(IJK)) - &
679                    DF_GU(IJK,S)*(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,T)*(U_G(IJKP)-U_G(IJK)) - &
686                       DF_GU(IJK,B)*(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     
696