File: RELATIVE:/../../../mfix.git/model/tau_w_g.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Subroutine: CALC_Tau_W_g                                            C
4     !  Purpose: Cross terms in the gradient of stress in W_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 w-momentum equation that cannot be     C
11     !  cast in the form: mu.grad(w). These components are stored in the    C
12     !  passed variable, which is then used as a source of the w-momentum   C
13     !  equation.                                                           C
14     !                                                                      C
15     !  The following w component viscous stress tensor terms are           C
16     !  calculated here:                                                    C
17     !  > part of 1/x d/dz (tau_zz) xdxdydz =>                              C
18     !            1/x d/dz (lambda.trcD) xdxdydz=>                          C
19     !    delta (lambda.trcD)Ap |T-B                                        C
20     !  > part of 1/x^2 d/dx (x^2 tau_xz) xdxdydz => or equivalently        C
21     !    part of (tau_xz/x + 1/x d/dx (x tau_xz) ) xdxdydz =>              C
22     !            1/x d/dx(mu.du/dz) xdxdydz =>                             C
23     !    delta (mu/x du/dz)Ayz |E-W                                        C
24     !  > part of 1/x^2 d/dx (x^2 tau_xz) xdxdydz => or equivalently        C
25     !    part of (tau_xz/x + 1/x d/dx (x tau_xz) ) xdxdydz =>              C
26     !            1/x d/dx(mu.du/dz) xdxdydz =>                             C
27     !    delta (mu/x du/dz)Ayz |E-W                                        C
28     !  > part of 1/x d/dz (tau_zz) xdxdydz =>                              C
29     !            1/x d/dz (mu/x dw/dz) xdxdydz =>                          C
30     !    delta (mu/x dw/dz)Axy |T-B                                        C
31     !  CYLINDRICAL TERMS:                                                  C
32     !  > part of 1/x d/dz (tau_zz) xdxdydz =>                              C
33     !            1/x d/dz (2.mu/x u) xdxdydz =>                            C
34     !    delta (2.mu/x u)Axy |T-B                                          C
35     !  > part of 1/x^2 d/dx (x^2 tau_xz) xdxdydz => or equivalently        C
36     !    part of (tau_xz/x + 1/x d/dx (x tau_xz)) xdxdydz =>               C
37     !            1/x (mu/x du/dz) xdxdydz =>                               C
38     !    delta (1/x mu/x du/dz)Vp                                          C
39     !                                                                      C
40     !  The following w-component CYLINDRICAL viscous stress tensor terms   C
41     !  are not calculated here. They are handled in source_w_g:            C
42     !  > part of 1/x^2 d/dx (x^2 tau_xz) xdxdydz => or equivalently        C
43     !    part of (tau_xz/x + 1/x d/dx (x tau_xz)) xdxdydz =>               C
44     !            1/x d/dx (x.mu.(-w/x)) xdxdydz =>                         C
45     !    delta (mu/x.(-w))Ayz |E-W                                         C
46     !  > part of 1/x^2 d/dx (x^2 tau_xz) xdxdydz => or equivalently        C
47     !    part of (tau_xz/x + 1/x d/dx (x tau_xz)) xdxdydz =>               C
48     !            mu/x dw/dx xdxdydz =>                                     C
49     !    delta (mu/x.(dw/dx))Vp |p                                         C
50     !  > part of 1/x^2 d/dx (x^2 tau_xz) xdxdydz => or equivalently        C
51     !    part of (tau_xz/x + 1/x d/dx (x tau_xz)) xdxdydz =>               C
52     !            1/x d/dx (x.mu.(-w/x)) xdxdydz =>                         C
53     !    delta (mu/x.(-w/x))Vp |p                                          C
54     !                                                                      C
55     !  To reconstitute the complete w-momentum gas phase viscous stress    C
56     !  tensor requires including the terms calculated in source_w_g        C
57     !  and the 'diffusional' components (i.e., those of the form           C
58     !  mu.grad(w)                                                          C
59     !                                                                      C
60     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
61           SUBROUTINE CALC_TAU_W_G(lTAU_W_G, lctau_w_g)
62     
63     ! Modules
64     !---------------------------------------------------------------------//
65           USE param
66           USE param1
67           USE constant
68           USE physprop
69           USE fldvar
70           USE visc_g
71           USE run
72           USE toleranc
73           USE geometry
74           USE indices
75           USE is
76           USE sendrecv
77           USE compar
78           USE functions
79           USE cutcell
80           IMPLICIT NONE
81     
82     ! Dummy arguments
83     !---------------------------------------------------------------------//
84     ! TAU_W_g
85           DOUBLE PRECISION, INTENT(OUT) :: lTAU_w_g(DIMENSION_3)
86     ! cTAU_W_g
87           DOUBLE PRECISION, INTENT(OUT) :: lcTAU_w_g(DIMENSION_3)
88     
89     ! Local variables
90     !---------------------------------------------------------------------//
91     ! Indices
92           INTEGER :: I, J, K, IM, JM, KP
93           INTEGER :: IJK, IJKE, IJKW, IJKN, IJKS, IJKT
94           INTEGER :: IJKNT, IJKST, IJKTE, IJKTW
95           INTEGER :: IJKP, IMJK, IJMK, IJKM
96           INTEGER :: IMJKP, IJMKP
97     ! Average volume fraction
98           DOUBLE PRECISION :: EPGA
99     ! Average gradients
100           DOUBLE PRECISION :: duodz
101     ! Source terms (Surface)
102           DOUBLE PRECISION :: Sbv, Ssx, Ssy, Ssz
103     ! Source terms (Volumetric)
104           DOUBLE PRECISION :: Vxz
105     !---------------------------------------------------------------------//
106     
107           IF((.NOT.CARTESIAN_GRID).OR.(CG_SAFE_MODE(5)==1)) THEN
108     
109     !$omp  parallel do default(none) &
110     !$omp  private(I, J, K, IM, JM, KP, IJK,                               &
111     !$omp          IJKE, IJKW, IJKN, IJKS, IJKT,                           &
112     !$omp          IJKTE, IJKNT, IJKTW, IJKST,                             &
113     !$omp          IMJK, IJMK, IJKM, IJKP, IJMKP, IMJKP,                   &
114     !$omp          EPGA, SBV, SSX, SSY, SSZ, vxz, duodz)                   &
115     !$omp  shared(ijkstart3, ijkend3, i_of, j_of, k_of, im1, jm1, kp1,     &
116     !$omp         do_k, cylindrical, ltau_w_g, lctau_w_g,                  &
117     !$omp         ep_g, epmu_gt, lambda_gt, trd_g, v_g, w_g, u_g,          &
118     !$omp         axy,  axy_w, ayz_w, axz_w, vol_w,                        &
119     !$omp         dy, dz, ox, ox_e, odz, odz_t, eplambda_gt)
120              DO IJK = IJKSTART3, IJKEND3
121                 K = K_OF(IJK)
122                 IJKT = TOP_OF(IJK)
123                 EPGA = AVG_Z(EP_G(IJK),EP_G(IJKT),K)
124                 IF ( .NOT.IP_AT_T(IJK) .AND. EPGA>DIL_EP_S) THEN
125                    J = J_OF(IJK)
126                    I = I_OF(IJK)
127                    IM = IM1(I)
128                    JM = JM1(J)
129                    KP = KP1(K)
130     
131                    IJKP = KP_OF(IJK)
132                    IMJK = IM_OF(IJK)
133                    IJMK = JM_OF(IJK)
134                    IJKM = KM_OF(IJK)
135                    IMJKP = KP_OF(IMJK)
136                    IJMKP = JM_OF(IJKP)
137     
138                    IJKN = NORTH_OF(IJK)
139                    IJKS = SOUTH_OF(IJK)
140                    IJKE = EAST_OF(IJK)
141                    IJKW = WEST_OF(IJK)
142                    IJKNT = TOP_OF(IJKN)
143                    IJKST = TOP_OF(IJKS)
144                    IJKTE = EAST_OF(IJKT)
145                    IJKTW = WEST_OF(IJKT)
146     
147     
148     ! Surface forces
149     ! Bulk viscosity term
150     ! part of 1/x d/dz (tau_zz) xdxdydz =>
151     !         1/x d/dz (lambda.trcD) xdxdydz=>
152     ! delta (lambda.trcD)Ap |T-B : at (i, j, k+1 - k-1)
153                    SBV = (EPLAMBDA_GT(IJKT)*TRD_G(IJKT)-&
154                           EPLAMBDA_GT(IJK)*TRD_G(IJK))*AXY(IJK)
155     
156     ! shear stress terms
157     ! part of 1/x^2 d/dx (x^2 tau_xz) xdxdydz => or equivalently
158     ! part of (tau_xz/x + 1/x d/dx (x tau_xz) ) xdxdydz =>
159     !         1/x d/dx(mu.du/dz) xdxdydz =>
160     ! delta (mu/x du/dz)Ayz |E-W : at (i+1/2-i-1/2, j, k+1/2)
161                    SSX = AVG_Z_H(AVG_X_H(EPMU_GT(IJK),EPMU_GT(IJKE),I),&
162                                  AVG_X_H(EPMU_GT(IJKT),EPMU_GT(IJKTE),I),K)*&
163                             (U_G(IJKP)-U_G(IJK))*OX_E(I)*ODZ_T(K)*&
164                             AYZ_W(IJK) - &
165                          AVG_Z_H(AVG_X_H(EPMU_GT(IJKW),EPMU_GT(IJK),IM),&
166                                  AVG_X_H(EPMU_GT(IJKTW),EPMU_GT(IJKT),IM),K)*&
167                             (U_G(IMJKP)-U_G(IMJK))*ODZ_T(K)*DY(J)*&
168                             (HALF*(DZ(K)+DZ(KP)))
169     ! DY(J)*HALF(DZ(k)+DZ(kp)) = oX_E(IM)*AYZ_W(IMJK), but avoids singularity
170     
171     ! part of d/dy (tau_zy) xdxdydz =>
172     !         d/dy (mu/x dv/dz) xdxdydz =>
173     ! delta (mu/x dv/dz)Axz |N-S : at (i, j+1/2 - j-1/2, k+1/2)
174                    SSY = AVG_Z_H(AVG_Y_H(EPMU_GT(IJK),EPMU_GT(IJKN),J),&
175                                  AVG_Y_H(EPMU_GT(IJKT),EPMU_GT(IJKNT),J),K)*&
176                             (V_G(IJKP)-V_G(IJK))*OX(I)*ODZ_T(K)*AXZ_W(IJK) -&
177                          AVG_Z_H(AVG_Y_H(EPMU_GT(IJKS),EPMU_GT(IJK),JM),&
178                                  AVG_Y_H(EPMU_GT(IJKST),EPMU_GT(IJKT),JM),K)*&
179                             (V_G(IJMKP)-V_G(IJMK))*OX(I)*ODZ_T(K)*AXZ_W(IJMK)
180     
181     ! part of 1/x d/dz (tau_zz) xdxdydz =>
182     !         1/x d/dz (mu/x dw/dz) xdxdydz =>
183     ! delta (mu/x dw/dz)Axy |T-B : at (i, j, k+1 - k-1)
184                    SSZ = EPMU_GT(IJKT)*(W_G(IJKP)-W_G(IJK))*OX(I)*ODZ(KP)*&
185                             AXY_W(IJK) - &
186                          EPMU_GT(IJK)*(W_G(IJK)-W_G(IJKM))*OX(I)*ODZ(K)*&
187                             AXY_W(IJKM)
188     
189     
190     ! Special terms for cylindrical coordinates
191                    IF (CYLINDRICAL) THEN
192     
193     ! part of 1/x d/dz (tau_zz) xdxdydz =>
194     !         1/x d/dz (2.mu/x u) xdxdydz =>
195     ! delta (2.mu/x u)Axy |T-B : at (i, j, k+1 - k-1)
196                       SSZ = SSZ + EPMU_GT(IJKT)*(U_G(IJKP)+U_G(IMJKP))*&
197                                      OX(I)*AXY_W(IJK) - &
198                                   EPMU_GT(IJK)*(U_G(IJK)+U_G(IMJK))*&
199                                      OX(I)*AXY_W(IJKM)
200     
201     ! part of 1/x^2 d/dx (x^2 tau_xz) xdxdydz => or equivalently
202     ! part of (tau_xz/x + 1/x d/dx (x tau_xz)) xdxdydz =>
203     !         1/x (mu/x du/dz) xdxdydz =>
204     ! delta (1/x mu/x du/dz)Vp : at (i, j, k+1/2)
205                       IF (OX_E(IM) /= UNDEFINED) THEN
206                          DUODZ = (U_G(IMJKP)-U_G(IMJK))*OX_E(IM)*ODZ_T(K)
207                       ELSE
208                          DUODZ = ZERO
209                       ENDIF
210                       VXZ = AVG_Z(EPMU_GT(IJK),EPMU_GT(IJKT),K)*OX(I)*HALF*&
211                             ( (U_G(IJKP)-U_G(IJK))*OX_E(I)*ODZ_T(K) + &
212                               DUODZ )
213                    ELSE
214                       VXZ = ZERO
215                    ENDIF
216     
217     ! Add the terms
218                    lTAU_W_G(IJK) =  SBV + SSX + SSY + SSZ + VXZ*VOL_W(IJK)
219     
220     ! Also calculate and store the full gas phase viscous stress tensor
221                    CALL GET_FULL_TAU_W_G(IJK, ltau_w_g, lctau_w_g)
222     
223                 ELSE
224                    lTAU_W_G(IJK) = ZERO
225                    lcTAU_W_G(IJK) = ZERO
226                 ENDIF   ! end if (.NOT. IP_AT_T(IJK) .AND. EPGA>DIL_EP_S)
227              ENDDO   ! end do ijk
228     !$omp end parallel do
229     
230           ELSE
231     ! if cartesian grid
232              CALL CALC_CG_TAU_W_G(lTAU_W_G, lctau_w_g)
233           ENDIF
234     
235           call send_recv(ltau_w_g,2)
236           call send_recv(lctau_w_g,2)
237     
238           RETURN
239     
240         CONTAINS
241     
242           INCLUDE 'fun_avg.inc'
243     
244         END SUBROUTINE CALC_TAU_W_G
245     
246     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
247     !                                                                      C
248     !  Subroutine: CALC_CG_Tau_W_g                                         C
249     !  Purpose: Cross terms in the gradient of stress in w_g momentum      C
250     !  based on cartesian grid cut cell.                                   C
251     !                                                                      C
252     !  Author: Jeff Dietiker                              Date: 01-Jul-09  C
253     !                                                                      C
254     !                                                                      C
255     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
256           SUBROUTINE CALC_CG_TAU_w_G(lTAu_w_G, lctau_w_g)
257     
258     ! Modules
259     !---------------------------------------------------------------------//
260           USE param
261           USE param1
262           USE constant
263           USE physprop
264           USE fldvar
265           USE visc_g
266           USE run
267           USE toleranc
268           USE geometry
269           USE indices
270           USE is
271           USE sendrecv
272           USE compar
273           USE fun_avg
274           USE functions
275     
276           USE bc
277           USE quadric
278           USE cutcell
279           IMPLICIT NONE
280     
281     ! Dummy arguments
282     !---------------------------------------------------------------------//
283     ! TAU_W_g
284           DOUBLE PRECISION, INTENT(OUT) :: lTAU_w_g(DIMENSION_3)
285     ! cTAU_W_g
286           DOUBLE PRECISION, INTENT(OUT) :: lcTAU_w_g(DIMENSION_3)
287     
288     ! Local variables
289     !---------------------------------------------------------------------//
290     ! Indices
291           INTEGER :: I, J, K, IM, JM, KP
292           INTEGER :: IJK, IJKE, IJKW, IJKN, IJKS, IJKT
293           INTEGER :: IJKNT, IJKST, IJKTE, IJKTW
294           INTEGER :: IJKP, IMJK, IJMK, IJKM
295           INTEGER :: IMJKP, IJMKP
296     ! Average volume fraction
297           DOUBLE PRECISION :: EPGA
298     ! Source terms (Surface)
299           DOUBLE PRECISION :: Sbv, Ssx, Ssy, Ssz
300     ! cartesian grid
301           DOUBLE PRECISION :: DEL_H, Nx, Ny, Nz
302           LOGICAL :: U_NODE_AT_ET, U_NODE_AT_EB, U_NODE_AT_WT, U_NODE_AT_WB
303           LOGICAL :: V_NODE_AT_NT, V_NODE_AT_NB, V_NODE_AT_ST, V_NODE_AT_SB
304           DOUBLE PRECISION :: dudz_at_E, dudz_at_W
305           DOUBLE PRECISION :: dvdz_at_N, dvdz_at_S
306           DOUBLE PRECISION :: Xi, Yi, Zi, Ui, Vi, Sx, Sy, Sz
307           DOUBLE PRECISION :: MU_GT_CUT, SSX_CUT, SSY_CUT
308           DOUBLE PRECISION :: UW_g, VW_g, WW_g
309           INTEGER :: BCV
310           CHARACTER(LEN=9) :: BCT
311     
312     !---------------------------------------------------------------------//
313     
314     !$omp  parallel do default(none) &
315     !$omp  private(I, J, K, IM, JM, KP, IJK,                               &
316     !$omp          IJKE, IJKW, IJKN, IJKS, IJKT,                           &
317     !$omp          IJKTE, IJKNT, IJKTW, IJKST,                             &
318     !$omp          IMJK, IJMK, IJKM, IJKP, IJMKP, IMJKP,                   &
319     !$omp          EPGA, SBV, SSX, SSY, SSZ,                               &
320     !$omp          BCV, BCT, noc_wg, uw_g, vw_g, ww_g,                     &
321     !$omp          del_h, nx, ny, nz, xi, yi, zi, sx, sy, sz, vi, ui,      &
322     !$omp          cut_tau_wg, mu_gt_cut, ssy_cut, ssx_cut,                &
323     !$omp          dudz_at_e, dudz_at_w, dvdz_at_S, dvdz_at_N,             &
324     !$omp          u_node_at_wb, u_node_at_wt, u_node_at_eb, u_node_at_et, &
325     !$omp          v_node_at_sb, v_node_at_st, v_node_at_nb, v_node_at_nt) &
326     !$omp  shared(ijkstart3, ijkend3, i_of, j_of, k_of, im1, jm1, kp1,     &
327     !$omp         do_k, cylindrical, ltau_w_g, lctau_w_g,                  &
328     !$omp         ep_g, mu_g, epmu_gt, lambda_gt, trd_g, v_g, w_g, u_g,    &
329     !$omp         axy,  axy_w, ayz_w, axz_w, vol, ox,                      &
330     !$omp         bc_type, bc_w_id, bc_hw_g, bc_uw_g, bc_vw_g, bc_ww_g,    &
331     !$omp         oneodz_t_u, oneodz_t_v, oneodz_t_w,                      &
332     !$omp         x_u, y_u, z_u, x_v, y_v, z_v, x_w, y_w, z_w,             &
333     !$omp         wall_u_at, wall_v_at, area_w_cut, cut_w_cell_at,         &
334     !$omp         blocked_u_cell_at, blocked_v_cell_at, eplambda_gt)
335           DO IJK = IJKSTART3, IJKEND3
336              K = K_OF(IJK)
337              IJKT = TOP_OF(IJK)
338              EPGA = AVG_Z(EP_G(IJK),EP_G(IJKT),K)
339              IF ( .NOT.IP_AT_T(IJK) .AND. EPGA>DIL_EP_S) THEN
340                 J = J_OF(IJK)
341                 I = I_OF(IJK)
342                 IM = IM1(I)
343                 JM = JM1(J)
344                 KP = KP1(K)
345     
346                 IJKP = KP_OF(IJK)
347                 IMJK = IM_OF(IJK)
348                 IJMK = JM_OF(IJK)
349                 IJKM = KM_OF(IJK)
350                 IMJKP = KP_OF(IMJK)
351                 IJMKP = JM_OF(IJKP)
352     
353                 IJKN = NORTH_OF(IJK)
354                 IJKS = SOUTH_OF(IJK)
355                 IJKE = EAST_OF(IJK)
356                 IJKW = WEST_OF(IJK)
357                 IJKNT = TOP_OF(IJKN)
358                 IJKST = TOP_OF(IJKS)
359                 IJKTE = EAST_OF(IJKT)
360                 IJKTW = WEST_OF(IJKT)
361     
362     ! bulk viscosity term
363                 SBV =  (EPLAMBDA_GT(IJKT)*TRD_G(IJKT)) * AXY_W(IJK) &
364                       -(EPLAMBDA_GT(IJK) *TRD_G(IJK) ) * AXY_W(IJKM)
365     
366     ! shear stress terms
367                 IF(.NOT.CUT_W_CELL_AT(IJK))   THEN
368                    SSX = AVG_Z_H(AVG_X_H(EPMU_GT(IJK),EPMU_GT(IJKE),I),&
369                                  AVG_X_H(EPMU_GT(IJKT),EPMU_GT(IJKTE),I),K)*&
370                             (U_G(IJKP)-U_G(IJK))*ONEoDZ_T_U(IJK)*AYZ_W(IJK) -&
371                          AVG_Z_H(AVG_X_H(EPMU_GT(IJKW),EPMU_GT(IJK),IM),&
372                                  AVG_X_H(EPMU_GT(IJKTW),EPMU_GT(IJKT),IM),K)*&
373                             (U_G(IMJKP)-U_G(IMJK))*ONEoDZ_T_U(IMJK)*AYZ_W(IMJK)
374     
375                    SSY = AVG_Z_H(AVG_Y_H(EPMU_GT(IJK),EPMU_GT(IJKN),J),&
376                                  AVG_Y_H(EPMU_GT(IJKT),EPMU_GT(IJKNT),J),K)*&
377                             (V_G(IJKP)-V_G(IJK))*ONEoDZ_T_V(IJK)*AXZ_W(IJK) - &
378                          AVG_Z_H(AVG_Y_H(EPMU_GT(IJKS),EPMU_GT(IJK),JM),&
379                                  AVG_Y_H(EPMU_GT(IJKST),EPMU_GT(IJKT),JM),K)*&
380                             (V_G(IJMKP)-V_G(IJMK))*ONEoDZ_T_V(IJMK)*AXZ_W(IJMK)
381     
382                    SSZ = EPMU_GT(IJKT)*(W_G(IJKP)-W_G(IJK))*&
383                             ONEoDZ_T_W(IJK)*AXY_W(IJK) - &
384                          EPMU_GT(IJK)*(W_G(IJK)-W_G(IJKM))*&
385                             ONEoDZ_T_W(IJKM)*AXY_W(IJKM)
386     
387     ! cut cell modifications
388     !---------------------------------------------------------------------//
389                 ELSE
390                    BCV = BC_W_ID(IJK)
391     
392                    IF(BCV > 0 ) THEN
393                       BCT = BC_TYPE(BCV)
394                    ELSE
395                       BCT = 'NONE'
396                    ENDIF
397     
398                    SELECT CASE (BCT)
399                       CASE ('CG_NSW')
400                          CUT_TAU_WG = .TRUE.
401                          NOC_WG     = .TRUE.
402                          UW_g = ZERO
403                          VW_g = ZERO
404                          WW_g = ZERO
405                       CASE ('CG_FSW')
406                          CUT_TAU_WG = .FALSE.
407                          NOC_WG     = .FALSE.
408                          UW_g = ZERO
409                          VW_g = ZERO
410                          WW_g = ZERO
411                       CASE('CG_PSW')
412                          IF(BC_HW_G(BC_W_ID(IJK))==UNDEFINED) THEN   ! same as NSW
413                             CUT_TAU_WG = .TRUE.
414                             NOC_WG     = .TRUE.
415                             UW_g = BC_UW_G(BCV)
416                             VW_g = BC_VW_G(BCV)
417                             WW_g = BC_WW_G(BCV)
418                          ELSEIF(BC_HW_G(BC_W_ID(IJK))==ZERO) THEN   ! same as FSW
419                             CUT_TAU_WG = .FALSE.
420                             NOC_WG     = .FALSE.
421                             UW_g = ZERO
422                             VW_g = ZERO
423                             WW_g = ZERO
424                          ELSE                              ! partial slip
425                             CUT_TAU_WG = .FALSE.
426                             NOC_WG     = .FALSE.
427                          ENDIF
428                       CASE ('NONE')
429                          lTAU_W_G(IJK) = ZERO
430                          CYCLE
431                    END SELECT
432     
433                    IF(CUT_TAU_WG) THEN
434                       MU_GT_CUT = (VOL(IJK)*EPMU_GT(IJK) + &
435                          VOL(IJKP)*EPMU_GT(IJKT))/(VOL(IJK) + VOL(IJKP))
436                    ELSE
437                       MU_GT_CUT = ZERO
438                    ENDIF
439     
440     ! SSX:
441                    U_NODE_AT_ET = ((.NOT.BLOCKED_U_CELL_AT(IJKP)).AND.&
442                                    (.NOT.WALL_U_AT(IJKP)))
443                    U_NODE_AT_EB = ((.NOT.BLOCKED_U_CELL_AT(IJK)).AND.&
444                                    (.NOT.WALL_U_AT(IJK)))
445                    U_NODE_AT_WT = ((.NOT.BLOCKED_U_CELL_AT(IMJKP)).AND.&
446                                    (.NOT.WALL_U_AT(IMJKP)))
447                    U_NODE_AT_WB = ((.NOT.BLOCKED_U_CELL_AT(IMJK)).AND.&
448                                    (.NOT.WALL_U_AT(IMJK)))
449     
450                    IF(U_NODE_AT_ET.AND.U_NODE_AT_EB) THEN
451                       Ui = HALF * (U_G(IJKP) + U_G(IJK))
452                       Xi = HALF * (X_U(IJKP) + X_U(IJK))
453                       Yi = HALF * (Y_U(IJKP) + Y_U(IJK))
454                       Zi = HALF * (Z_U(IJKP) + Z_U(IJK))
455                       Sx = X_U(IJKP) - X_U(IJK)
456                       Sy = Y_U(IJKP) - Y_U(IJK)
457                       Sz = Z_U(IJKP) - Z_U(IJK)
458                       CALL GET_DEL_H(IJK, 'W_MOMENTUM', Xi, Yi, Zi, &
459                          Del_H, Nx, Ny, Nz)
460     
461                       dudz_at_E =  (U_G(IJKP) - U_G(IJK)) * &
462                          ONEoDZ_T_U(IJK)
463                       IF(NOC_WG) dudz_at_E = dudz_at_E - ((Ui-UW_g) * &
464                          ONEoDZ_T_U(IJK)/DEL_H*(Sx*Nx+Sy*Ny))
465                    ELSE
466                       dudz_at_E =  ZERO
467                    ENDIF
468     
469                    IF(U_NODE_AT_WT.AND.U_NODE_AT_WB) THEN
470                       Ui = HALF * (U_G(IMJKP) + U_G(IMJK))
471                       Xi = HALF * (X_U(IMJKP) + X_U(IMJK))
472                       Yi = HALF * (Y_U(IMJKP) + Y_U(IMJK))
473                       Zi = HALF * (Z_U(IMJKP) + Z_U(IMJK))
474                       Sx = X_U(IMJKP) - X_U(IMJK)
475                       Sy = Y_U(IMJKP) - Y_U(IMJK)
476                       Sz = Z_U(IMJKP) - Z_U(IMJK)
477                       CALL GET_DEL_H(IJK, 'W_MOMENTUM', Xi, Yi, Zi, &
478                          Del_H, Nx, Ny, Nz)
479     
480                       dudz_at_W =  (U_G(IMJKP) - U_G(IMJK)) * &
481                          ONEoDZ_T_U(IMJK)
482                       IF(NOC_WG) dudz_at_W = dudz_at_W - ((Ui-UW_g) * &
483                          ONEoDZ_T_U(IMJK)/DEL_H*(Sx*Nx+Sy*Ny))
484                    ELSE
485                       dudz_at_W =  ZERO
486                    ENDIF
487     
488                    IF(U_NODE_AT_EB) THEN
489                       CALL GET_DEL_H(IJK, 'W_MOMENTUM', X_U(IJK), &
490                          Y_U(IJK), Z_U(IJK), Del_H, Nx, Ny, Nz)
491                       SSX_CUT = - MU_GT_CUT * (U_G(IJK) - UW_g) / &
492                          DEL_H * (Nz*Nx) * Area_W_CUT(IJK)
493                    ELSE
494                       SSX_CUT =  ZERO
495                    ENDIF
496                    SSX = AVG_Z_H(AVG_X_H(EPMU_GT(IJK),EPMU_GT(IJKE),I),&
497                                  AVG_X_H(EPMU_GT(IJKT),EPMU_GT(IJKTE),I),K)*&
498                             dudz_at_E*AYZ_W(IJK) - &
499                          AVG_Z_H(AVG_X_H(EPMU_GT(IJKW),EPMU_GT(IJK),IM),&
500                                  AVG_X_H(EPMU_GT(IJKTW),EPMU_GT(IJKT),IM),K)*&
501                             dudz_at_W*AYZ_W(IMJK) + SSX_CUT
502     
503     ! SSY:
504                    V_NODE_AT_NT = ((.NOT.BLOCKED_V_CELL_AT(IJKP)).AND.&
505                                    (.NOT.WALL_V_AT(IJKP)))
506                    V_NODE_AT_NB = ((.NOT.BLOCKED_V_CELL_AT(IJK)).AND.&
507                                    (.NOT.WALL_V_AT(IJK)))
508                    V_NODE_AT_ST = ((.NOT.BLOCKED_V_CELL_AT(IJMKP)).AND.&
509                                    (.NOT.WALL_V_AT(IJMKP)))
510                    V_NODE_AT_SB = ((.NOT.BLOCKED_V_CELL_AT(IJMK)).AND.&
511                                    (.NOT.WALL_V_AT(IJMK)))
512     
513                    IF(V_NODE_AT_NT.AND.V_NODE_AT_NB) THEN
514                       Vi = HALF * (V_G(IJKP) + V_G(IJK))
515                       Xi = HALF * (X_V(IJKP) + X_V(IJK))
516                       Yi = HALF * (Y_V(IJKP) + Y_V(IJK))
517                       Zi = HALF * (Z_V(IJKP) + Z_V(IJK))
518                       Sx = X_V(IJKP) - X_V(IJK)
519                       Sy = Y_V(IJKP) - Y_V(IJK)
520                       Sz = Z_V(IJKP) - Z_V(IJK)
521                       CALL GET_DEL_H(IJK, 'W_MOMENTUM', Xi, Yi, Zi, &
522                          Del_H, Nx, Ny, Nz)
523     
524                       dvdz_at_N =  (V_G(IJKP) - V_G(IJK)) * &
525                          ONEoDZ_T_V(IJK)
526                       IF(NOC_WG) dvdz_at_N = dvdz_at_N - ((Vi-VW_g) * &
527                          ONEoDZ_T_V(IJK)/DEL_H*(Sx*Nx+Sy*Ny))
528                    ELSE
529                          dvdz_at_N =  ZERO
530                    ENDIF
531     
532                    IF(V_NODE_AT_ST.AND.V_NODE_AT_SB) THEN
533                       Vi = HALF * (V_G(IJMKP) + V_G(IJMK))
534                       Xi = HALF * (X_V(IJMKP) + X_V(IJMK))
535                       Yi = HALF * (Y_V(IJMKP) + Y_V(IJMK))
536                       Zi = HALF * (Z_V(IJMKP) + Z_V(IJMK))
537                       Sx = X_V(IJMKP) - X_V(IJMK)
538                       Sy = Y_V(IJMKP) - Y_V(IJMK)
539                       Sz = Z_V(IJMKP) - Z_V(IJMK)
540                       CALL GET_DEL_H(IJK, 'W_MOMENTUM', Xi, Yi, Zi, &
541                          Del_H, Nx, Ny, Nz)
542     
543                       dvdz_at_S =  (V_G(IJMKP) - V_G(IJMK)) * &
544                          ONEoDZ_T_V(IJMK)
545                       IF(NOC_WG) dvdz_at_S = dvdz_at_S - ((Vi-VW_g) * &
546                          ONEoDZ_T_V(IJMK)/DEL_H*(Sx*Nx+Sy*Ny))
547                    ELSE
548                       dvdz_at_S =  ZERO
549                    ENDIF
550     
551                    IF(V_NODE_AT_NB) THEN
552                       CALL GET_DEL_H(IJK, 'W_MOMENTUM', X_V(IJK), &
553                          Y_V(IJK), Z_V(IJK), Del_H, Nx, Ny, Nz)
554                       SSY_CUT = - MU_GT_CUT * (V_G(IJK) - VW_g) / &
555                          DEL_H * (Nz*Ny) * Area_W_CUT(IJK)
556                    ELSE
557                       SSY_CUT =  ZERO
558                    ENDIF
559     
560                    SSY = AVG_Z_H(AVG_Y_H(EPMU_GT(IJK),EPMU_GT(IJKN),J),&
561                                  AVG_Y_H(EPMU_GT(IJKT),EPMU_GT(IJKNT),J),K)*&
562                             dvdz_at_N*AXZ_W(IJK) - &
563                          AVG_Z_H(AVG_Y_H(EPMU_GT(IJKS),EPMU_GT(IJK),JM),&
564                                  AVG_Y_H(EPMU_GT(IJKST),EPMU_GT(IJKT),JM),K)*&
565                             dvdz_at_S*AXZ_W(IJMK) + SSY_CUT
566     
567     ! SSZ:
568                    CALL GET_DEL_H(IJK, 'W_MOMENTUM', X_W(IJK), &
569                       Y_W(IJK), Z_W(IJK), Del_H, Nx, Ny, Nz)
570     
571                    SSZ = EPMU_GT(IJKT)*(W_G(IJKP)-W_G(IJK))*&
572                             ONEoDZ_T_W(IJK)*AXY_W(IJK) - &
573                          EPMU_GT(IJK)*(W_G(IJK)-W_G(IJKM))*&
574                             OX(I)*ONEoDZ_T_W(IJKM)*AXY_W(IJKM) - &
575                          MU_GT_CUT * (W_g(IJK) - WW_g) / DEL_H * &
576                             (Nz**2) * Area_W_CUT(IJK)
577     
578                 ENDIF  ! end if/else cut_w_cell_at
579     
580     ! Add the terms
581                 lTAU_W_G(IJK) =  SBV + SSX + SSY + SSZ
582     
583     ! Also calculate and store the full gas phase viscous stress tensor
584                 CALL GET_FULL_TAU_W_G(IJK, ltau_w_g, lctau_w_g)
585     
586              ELSE
587                 lTAU_W_G(IJK) = ZERO
588                 lcTAU_W_G(IJK) = ZERO
589              ENDIF   ! end if (.NOT. IP_AT_T(IJK) .AND. EPGA>DIL_EP_S)
590           ENDDO   ! end do ijk
591     !$omp end parallel do
592     
593           RETURN
594           END SUBROUTINE CALC_CG_TAU_W_G
595     
596     
597     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
598     !                                                                      C
599     !  Subroutine: GET_FULL_Tau_W_g                                        C
600     !  Purpose: Calculate the divergence of the complete gas phase stress  C
601     !  tensor including all terms in the w-direction.                      C
602     !                                                                      C
603     !                                                                      C
604     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
605           SUBROUTINE GET_FULL_TAU_W_G(IJK, ltau_w_g, lctau_w_g)
606     
607     ! Modules
608     !---------------------------------------------------------------------//
609           USE fldvar, only: w_g
610     
611           USE functions, only: east_of, west_of, top_of, flow_at_t
612           USE functions, only: im_of, jm_of, km_of
613           USE functions, only: ip_of, jp_of, kp_of
614     
615           USE fun_avg, only: avg_x_h, avg_z_h, avg_z
616     
617           USE geometry, only: cylindrical, ox_e, ox, odx_e
618           USE geometry, only: dy, dz, vol_w, ayz_w
619           USE indices, only: i_of, j_of, k_of, im1, kp1
620     
621           USE matrix, only: e, w, n, s, t, b
622           use param, only: dimension_3
623           USE param1, only: zero, half
624           USE visc_g, only: mu_gt, df_gw
625           IMPLICIT NONE
626     
627     ! Dummy arguments
628     !---------------------------------------------------------------------//
629     ! ijk index
630           INTEGER, INTENT(IN) :: IJK
631     ! TAU_W_g
632           DOUBLE PRECISION, INTENT(IN) :: lTAU_W_g(DIMENSION_3)
633     ! cTAU_W_g
634           DOUBLE PRECISION, INTENT(INOUT) :: lcTAU_W_g(DIMENSION_3)
635     
636     ! Local Variables
637     !---------------------------------------------------------------------//
638     ! indices
639           INTEGER :: I, J, K, IM
640           INTEGER :: IJKT, IJKE, IJKW, IJKTE, IJKTW
641           INTEGER :: IPJK, IJPK, IJKP, IMJK, IJMK, IJKM
642     ! average viscosity
643           DOUBLE PRECISION :: MUOX
644     ! source terms
645           DOUBLE PRECISION :: SSX, SSY, SSZ
646     ! cylindrical terms
647           DOUBLE PRECISION :: vxza
648           DOUBLE PRECISION :: cte, ctw, vx14
649           DOUBLE PRECISION :: cpe, cpw, vx15
650     !---------------------------------------------------------------------//
651     
652           IPJK = IP_OF(IJK)
653           IJPK = JP_OF(IJK)
654           IJKP = KP_OF(IJK)
655           IMJK = IM_OF(IJK)
656           IJMK = JM_OF(IJK)
657           IJKM = KM_OF(IJK)
658     
659     ! additional terms (found in source_W_g)
660           VXZA = ZERO
661           VX14 = ZERO
662           VX15 = ZERO
663           IF (CYLINDRICAL) THEN
664              I = I_OF(IJK)
665              J = J_OF(IJK)
666              K = K_OF(IJK)
667              IM = IM1(I)
668              IJKT = TOP_OF(IJK)
669              IJKE = EAST_OF(IJK)
670              IJKW = WEST_OF(IJK)
671              IJKTE = TOP_OF(IJKE)
672              IJKTW = TOP_OF(IJKW)
673     
674     ! part of 1/x^2 d/dx (x^2 tau_xz) xdxdydz => or equivalently
675     ! part of (tau_xz/x + 1/x d/dx (x tau_xz)) xdxdydz =>
676     !         1/x d/dx (x.mu.(-w/x)) xdxdydz =>
677     ! delta (mu/x.(-w))Ayz |E-W : at (i+1/2 - i-1/2, j, k+1/2)
678              CTE = HALF*AVG_Z_H(AVG_X_H(MU_GT(IJK),MU_GT(IJKE),I),&
679                                 AVG_X_H(MU_GT(IJKT),MU_GT(IJKTE),I),K)*&
680                    OX_E(I)*AYZ_W(IJK)
681              CTW = HALF*AVG_Z_H(AVG_X_H(MU_GT(IJKW),MU_GT(IJK),IM),&
682                           AVG_X_H(MU_GT(IJKTW),MU_GT(IJKT),IM),K)*&
683                    DY(J)*(HALF*(DZ(K)+DZ(KP1(K))))
684     ! DY(J)*HALF(DZ(k)+DZ(kp)) = oX_E(IM)*AYZ_W(IMJK), but avoids singularity
685              VX14 = -CTE*(W_G(IPJK)+W_G(IJK))+CTW*(W_G(IJK)+W_G(IMJK))
686     
687     ! part of 1/x^2 d/dx (x^2 tau_xz) xdxdydz => or equivalently
688     ! part of (tau_xz/x + 1/x d/dx (x tau_xz)) xdxdydz =>
689     !         mu/x dw/dx xdxdydz =>
690     ! delta (mu/x.(dw/dx))Vp |p : at (i, j, k+1/2)
691              MUOX = AVG_Z(MU_GT(IJK),MU_GT(IJKT),K)*OX(I)
692              CPE = HALF*MUOX*ODX_E(I)*VOL_W(IJK)
693              CPW = HALF*MUOX*ODX_E(IM)*VOL_W(IJK)
694              VX15 = CPE*(W_G(IPJK)-W_G(IJK))+CPW*(W_G(IJK)-W_G(IMJK))
695     
696     ! part of 1/x^2 d/dx (x^2 tau_xz) xdxdydz => or equivalently
697     ! part of (tau_xz/x + 1/x d/dx (x tau_xz)) xdxdydz =>
698     !         1/x d/dx (x.mu.(-w/x)) xdxdydz =>
699     ! delta (mu/x.(-w/x))Vp |p : at (i, j, k+1/2)
700              VXZA = -MUOX*OX(I)*W_G(IJK)*VOL_W(IJK)
701           ENDIF
702     
703     
704     ! convection terms
705           SSX = ZERO
706           SSY = ZERO
707           SSZ = ZERO
708           IF (FLOW_AT_T(IJK)) THEN
709     ! part of 1/x^2 d/dx (x^2 tau_xz) xdxdydz => or equivalently
710     ! part of (tau_xz/x + 1/x d/dx (x tau_xz)) xdxdydz =>
711     !         1/x d/dx (x.mu.dw/dx) xdxdydz =>
712     ! delta (mu.dw/dx.Ayz) |E-W : at (i+1/2 - i-1/2), j, k+1/2
713              SSX = DF_GW(IJK,E)*(W_G(IPJK) - W_G(IJK)) - &
714                    DF_GW(IJK,W)*(W_G(IJK) - W_G(IJKM))
715     
716     ! part of d/dy (tau_zy) xdxdydz =>
717     !         d/dy (mu.dw/dy) xdxdydz =>
718     ! delta (mu.dw/dy.Axz) |N-S : at (i, j+1/2 - j-1/2, k+1/2)
719              SSY = DF_GW(IJK,N)*(W_G(IJPK)-W_G(IJK)) - &
720                    DF_GW(IJK,S)*(W_G(IJK)-W_G(IJMK))
721     
722     ! part of 1/x d/dz (tau_zz) xdxdydz =>
723     !         1/x d/dz (mu/x.dw/dz) xdxdydz =>
724     ! delta (mu/x.dw/dz.Axy) |T-B : at (i, j, k+1 - k-1)
725              SSZ = DF_GW(IJK,T)*(W_G(IJKP)-W_G(IJK)) - &
726                    DF_GW(IJK,B)*(W_G(IJK)-W_G(IJKM))
727     
728           ENDIF
729     
730     ! Add the terms
731           lctau_W_g(IJK) = (lTAU_W_G(IJK) + SSX + SSY + SSZ + VXZA + &
732                             VX14 + VX15)
733     
734           RETURN
735           END SUBROUTINE GET_FULL_TAU_W_G
736     
737