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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Subroutine: CALC_Tau_V_g                                            C
4     !  Purpose: Cross terms in the gradient of stress in V_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 v-momentum equation that cannot be     C
11     !  cast in the form: mu.grad(v). These components are stored in the    C
12     !  passed variable, which is then used as a source of the v-momentum   C
13     !  equation.                                                           C
14     !                                                                      C
15     !  The following v component viscous stress tensor terms are           C
16     !  calculated here:                                                    C
17     !  > part of d/dy (tau_yy) xdxdydz =>                                  C
18     !            d/dy (lambda.trcD) xdxdydz =>                             C
19     !    delta (lambda.trcD)Ap|N-S                                         C
20     !  > part of 1/x d/dx(x.tau_xy) xdxdydz =>                             C
21     !            1/x d/dx (x.mu.du/dy) xdxdydz =>                          C
22     !    delta (x.mu.du/dy)Ayz |E-W                                        C
23     !  > part of d/dy (tau_xy) xdxdydz =>                                  C
24     !           d/dy (mu.dv/dy) xdxdydz =>                                 C
25     !    delta (mu.dv/dx)Axz |N-S                                          C
26     !  > part of 1/x d/dz (tau_xz) xdxdydz =>                              C
27     !            1/x d/dz (mu.dw/dy) xdxdydz =>                            C
28     !    delta (mu.dw/dx)Axy |T-B                                          C
29     !                                                                      C
30     !  To reconstitute the full v-momentum gas phase viscous stress        C
31     !  tensor would require including the the 'diffusional' components     C
32     !  (i.e., those of the form mu.grad(v)                                 C
33     !                                                                      C
34     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
35           SUBROUTINE CALC_TAU_V_G(lTAU_V_G, lctau_v_g)
36     
37     ! Modules
38     !---------------------------------------------------------------------//
39           USE param
40           USE param1
41           USE constant
42           USE physprop
43           USE fldvar
44           USE visc_g
45           USE run
46           USE toleranc
47           USE geometry
48           USE indices
49           USE is
50           USE sendrecv
51           USE compar
52           USE functions
53           USE cutcell
54           IMPLICIT NONE
55     
56     ! Dummy arguments
57     !---------------------------------------------------------------------//
58     ! TAU_V_g
59           DOUBLE PRECISION, INTENT(OUT) :: lTAU_V_g(DIMENSION_3)
60     ! cTAU_V_g
61           DOUBLE PRECISION, INTENT(OUT) :: lcTAU_V_g(DIMENSION_3)
62     
63     ! Local variables
64     !---------------------------------------------------------------------//
65     ! Indices
66           INTEGER :: I, J, K, IM, JP, KM
67           INTEGER :: IJK, IJKN, IJKE, IJKW, IJKT, IJKB
68           INTEGER :: IJKTN, IJKBN, IJKNE, IJKNW
69           INTEGER :: IJPK, IJMK, IMJK, IJKM
70           INTEGER :: IMJPK, IJPKM
71     ! Average volume fraction
72           DOUBLE PRECISION :: EPGA
73     ! Source terms (Surface)
74           DOUBLE PRECISION :: Sbv, Ssx, Ssy, Ssz
75     !---------------------------------------------------------------------//
76     
77           IF((.NOT.CARTESIAN_GRID).OR.(CG_SAFE_MODE(4)==1)) THEN
78     
79     !$omp  parallel do default(none) &
80     !$omp  private(I, J, K, IJK, IM, JP, KM,                               &
81     !$omp          IJKE, IJKW, IJKN, IJKT, IJKB,                           &
82     !$omp          IJKTN, IJKBN, IJKNW, IJKNE,                             &
83     !$omp          IJPK, IMJK, IJMK, IJKM, IJPKM, IMJPK,                   &
84     !$omp          EPGA, SBV, SSX, SSY, SSZ)                               &
85     !$omp  shared(ijkstart3, ijkend3, i_of, j_of, k_of, im1, jp1, km1,     &
86     !$omp         do_k, ltau_v_g, lctau_v_g,                               &
87     !$omp         ep_g, lambda_gt, trd_g, epmu_gt, mu_g, u_g, v_g, w_g,      &
88     !$omp         ayz_v, axz_v, axz, axy_v,                                &
89     !$omp         ody_n, ody)
90              DO IJK = IJKSTART3, IJKEND3
91                 J = J_OF(IJK)
92                 IJKN = NORTH_OF(IJK)
93                 EPGA = AVG_Y(EP_G(IJK),EP_G(IJKN),J)
94                 IF ( .NOT.IP_AT_N(IJK) .AND. EPGA>DIL_EP_S) THEN
95                    I = I_OF(IJK)
96                    K = K_OF(IJK)
97                    JP = JP1(J)
98                    IM = IM1(I)
99                    KM = KM1(K)
100     
101                    IJPK = JP_OF(IJK)
102                    IJMK = JM_OF(IJK)
103                    IMJK = IM_OF(IJK)
104                    IJKM = KM_OF(IJK)
105                    IJPKM = JP_OF(IJKM)
106                    IMJPK = IM_OF(IJPK)
107     
108                    IJKE = EAST_OF(IJK)
109                    IJKNE = EAST_OF(IJKN)
110                    IJKW = WEST_OF(IJK)
111                    IJKNW = NORTH_OF(IJKW)
112                    IJKT = TOP_OF(IJK)
113                    IJKTN = NORTH_OF(IJKT)
114                    IJKB = BOTTOM_OF(IJK)
115                    IJKBN = NORTH_OF(IJKB)
116     
117     ! Surface forces at i, j+1/2, k
118     ! bulk viscosity term
119     ! part of d/dy (tau_yy) xdxdydz =>
120     !         d/dy (lambda.trcD) xdxdydz =>
121     ! delta (lambda.trcD)Ap|N-S  : at (i, j+1 - j-1, k)
122                    SBV = (LAMBDA_GT(IJKN)*TRD_G(IJKN)-&
123                           LAMBDA_GT(IJK)*TRD_G(IJK))*AXZ(IJK)
124     
125     ! shear stress terms at i, j+1/2, k
126     
127     ! part of 1/x d/dx(x.tau_xy) xdxdydz =>
128     !         1/x d/dx (x.mu.du/dy) xdxdydz =>
129     ! delta (x.mu.du/dy)Ayz |E-W : at (i+1/2 - i-1/2, j+1/2, k)
130                    SSX = AVG_Y_H(AVG_X_H(EPMU_GT(IJK),EPMU_GT(IJKE),I),&
131                                  AVG_X_H(EPMU_GT(IJKN),EPMU_GT(IJKNE),I),J)*&
132                          (U_G(IJPK)-U_G(IJK))*ODY_N(J)*AYZ_V(IJK) - &
133                          AVG_Y_H(AVG_X_H(EPMU_GT(IJKW),EPMU_GT(IJK),IM),&
134                                  AVG_X_H(EPMU_GT(IJKNW),EPMU_GT(IJKN),IM),J)*&
135                          (U_G(IMJPK)-U_G(IMJK))*ODY_N(J)*AYZ_V(IMJK)
136     
137     ! part of d/dy (tau_xy) xdxdydz =>
138     !         d/dy (mu.dv/dy) xdxdydz =>
139     ! delta (mu.dv/dx)Axz |N-S : at (i, j+1 - j-1, k)
140                    SSY = EPMU_GT(IJKN)*(V_G(IJPK)-V_G(IJK))*ODY(JP)*&
141                             AXZ_V(IJK) - &
142                          EPMU_GT(IJK)*(V_G(IJK)-V_G(IJMK))*ODY(J)*&
143                             AXZ_V(IJMK)
144     
145     ! part of 1/x d/dz (tau_xz) xdxdydz =>
146     !         1/x d/dz (mu.dw/dy) xdxdydz =>
147     ! delta (mu.dw/dx)Axy |T-B : at (i, j+1/2, k+1/2 - k-1/2)
148                    SSZ = AVG_Y_H(AVG_Z_H(EPMU_GT(IJK),EPMU_GT(IJKT),K),&
149                                  AVG_Z_H(EPMU_GT(IJKN),EPMU_GT(IJKTN),K),J)*&
150                          (W_G(IJPK)-W_G(IJK))*ODY_N(J)*AXY_V(IJK) - &
151                          AVG_Y_H(AVG_Z_H(EPMU_GT(IJKB),EPMU_GT(IJK),KM),&
152                                  AVG_Z_H(EPMU_GT(IJKBN),EPMU_GT(IJKN),KM),J)*&
153                          (W_G(IJPKM)-W_G(IJKM))*ODY_N(J)*AXY_V(IJKM)
154     
155     ! Add the terms
156                    lTAU_V_G(IJK) = SBV + SSX + SSY + SSZ
157     
158     ! Also calculate and store the full gas phase viscous stress tensor
159                    CALL GET_FULL_TAU_V_G(IJK, ltau_v_g, lctau_v_g)
160     
161                 ELSE
162                    lTAU_V_G(IJK) = ZERO
163                    lctau_v_g(IJK) = ZERO
164                 ENDIF   ! end if (.NOT. IP_AT_N(IJK) .AND. EPGA>DIL_EP_S)
165              ENDDO   ! end do ijk
166     !$omp end parallel do
167     
168           ELSE
169     ! cartesian grid case
170              CALL CALC_CG_TAU_V_G(lTAU_V_G, lctau_v_g)
171           ENDIF
172     
173           call send_recv(ltau_v_g,2)
174           call send_recv(lctau_v_g,2)
175           RETURN
176     
177         CONTAINS
178     
179           INCLUDE 'fun_avg.inc'
180     
181         END SUBROUTINE CALC_TAU_V_G
182     
183     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
184     !                                                                      C
185     !  Subroutine: CALC_CG_Tau_v_g                                         C
186     !  Purpose: Cross terms in the gradient of stress in V_g momentum      C
187     !  based on cartesian grid cut cell.                                   C
188     !                                                                      C
189     !  Author: Jeff Dietiker                              Date: 01-Jul-09  C
190     !                                                                      C
191     !                                                                      C
192     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
193           SUBROUTINE CALC_CG_TAU_V_G(lTAU_V_G, lctau_v_g)
194     
195     ! Modules
196     !---------------------------------------------------------------------//
197           USE param
198           USE param1
199           USE constant
200           USE physprop
201           USE fldvar
202           USE visc_g
203           USE run
204           USE toleranc
205           USE geometry
206           USE indices
207           USE is
208           USE sendrecv
209           USE compar
210           USE fun_avg
211           USE functions
212     
213           USE cutcell
214     
215           USE bc
216           USE quadric
217           USE cutcell
218     
219           IMPLICIT NONE
220     
221     ! Dummy arguments
222     !---------------------------------------------------------------------//
223     ! TAU_V_g
224           DOUBLE PRECISION, INTENT(OUT) :: lTAU_V_g(DIMENSION_3)
225     ! CTAU_V_g
226           DOUBLE PRECISION, INTENT(OUT) :: lcTAU_V_g(DIMENSION_3)
227     
228     ! Local variables
229     !---------------------------------------------------------------------//
230     ! Indices
231           INTEGER :: I, J, K, IM, JP, KM
232           INTEGER :: IJK, IJKN, IJKE, IJKW, IJKT, IJKB
233           INTEGER :: IJKTN, IJKBN, IJKNE, IJKNW
234           INTEGER :: IJPK, IJMK, IMJK, IJKM
235           INTEGER :: IMJPK, IJPKM
236     ! Average volume fraction
237           DOUBLE PRECISION EPGA
238     ! Source terms (Surface)
239           DOUBLE PRECISION Sbv, Ssx, Ssy, Ssz
240     
241     ! Cartesian grid variables
242           DOUBLE PRECISION :: DEL_H, Nx, Ny, Nz
243           LOGICAL :: U_NODE_AT_NE, U_NODE_AT_NW, U_NODE_AT_SE, U_NODE_AT_SW
244           LOGICAL :: W_NODE_AT_TN, W_NODE_AT_TS, W_NODE_AT_BN, W_NODE_AT_BS
245           DOUBLE PRECISION :: dudy_at_E, dudy_at_W
246           DOUBLE PRECISION :: dwdy_at_T, dwdy_at_B
247           DOUBLE PRECISION :: Xi, Yi, Zi, Ui, Wi, Sx, Sy, Sz
248           DOUBLE PRECISION :: MU_GT_CUT, SSX_CUT, SSZ_CUT
249           DOUBLE PRECISION :: UW_g, VW_g, WW_g
250           INTEGER :: BCV
251           INTEGER :: BCT
252     !---------------------------------------------------------------------//
253     
254     !$omp  parallel do default(none)                                       &
255     !$omp  private(I, J, K, IM, JP, KM, IJK,                               &
256     !$omp          IJKE, IJKW, IJKN, IJKT, IJKB,                           &
257     !$omp          IJKTN, IJKBN, IJKNW, IJKNE,                             &
258     !$omp          IJPK, IMJK, IJMK, IJKM, IJPKM, IMJPK,                   &
259     !$omp          EPGA, SBV, SSX, SSY, SSZ,                               &
260     !$omp          BCV, BCT, BC_TYPE_ENUM, NOC_VG, uw_g, vw_g, ww_g,       &
261     !$omp          del_h, nx, ny, nz, xi, yi, zi, sx, sy, sz, wi, ui,      &
262     !$omp          u_node_at_sw, u_node_at_se, u_node_at_nw, u_node_at_ne, &
263     !$omp          w_node_at_bs, w_node_at_bn, w_node_at_ts, w_node_at_tn, &
264     !$omp          dudy_at_W, dudy_at_E, dwdy_at_B, dwdy_at_T,             &
265     !$omp          cut_tau_vg, mu_gt_cut, ssz_cut, ssx_cut)                &
266     !$omp  shared(ijkstart3, ijkend3, i_of, j_of, k_of, do_k,              &
267     !$omp         im1, jm1, jp1, km1, ltau_v_g, lctau_v_g,                 &
268     !$omp         epmu_gt, ep_g, lambda_gt, trd_g, v_g, w_g, u_g,          &
269     !$omp         ayz_v, axz_v, axz, axy_v, vol,                           &
270     !$omp         bc_type, bc_v_id, bc_hw_g, bc_uw_g, bc_vw_g, bc_ww_g,    &
271     !$omp         x_u, y_u, z_u, x_v, y_v, z_v, x_w, y_w, z_w,             &
272     !$omp         oneody_n_v, oneody_n_w, oneody_n_u,                      &
273     !$omp         cut_v_cell_at, wall_u_at, wall_w_at, area_v_cut,         &
274     !$omp         blocked_u_cell_at, blocked_w_cell_at)
275     
276           DO IJK = IJKSTART3, IJKEND3
277              J = J_OF(IJK)
278              IJKN = NORTH_OF(IJK)
279              EPGA = AVG_Y(EP_G(IJK),EP_G(IJKN),J)
280              IF ( .NOT.IP_AT_N(IJK) .AND. EPGA>DIL_EP_S) THEN
281                 I = I_OF(IJK)
282                 K = K_OF(IJK)
283                 JP = JP1(J)
284                 IM = IM1(I)
285                 KM = KM1(K)
286     
287                 IJPK = JP_OF(IJK)
288                 IJMK = JM_OF(IJK)
289                 IMJK = IM_OF(IJK)
290                 IJKM = KM_OF(IJK)
291                 IMJPK = IM_OF(IJPK)
292                 IJPKM = JP_OF(IJKM)
293     
294                 IJKE = EAST_OF(IJK)
295                 IJKNE = EAST_OF(IJKN)
296                 IJKW = WEST_OF(IJK)
297                 IJKNW = NORTH_OF(IJKW)
298                 IJKT = TOP_OF(IJK)
299                 IJKTN = NORTH_OF(IJKT)
300                 IJKB = BOTTOM_OF(IJK)
301                 IJKBN = NORTH_OF(IJKB)
302     
303     
304     ! bulk viscosity term
305                 SBV =  (LAMBDA_GT(IJKN)*TRD_G(IJKN)) * AXZ_V(IJK) - &
306                        (LAMBDA_GT(IJK) *TRD_G(IJK) ) * AXZ_V(IJMK)
307     
308     ! shear stress terms
309                 IF(.NOT.CUT_V_CELL_AT(IJK)) THEN
310                    SSX = AVG_Y_H(AVG_X_H(EPMU_GT(IJK),EPMU_GT(IJKE),I),&
311                                  AVG_X_H(EPMU_GT(IJKN),EPMU_GT(IJKNE),I),J)*&
312                          (U_G(IJPK)-U_G(IJK))*ONEoDY_N_U(IJK)*AYZ_V(IJK) - &
313                          AVG_Y_H(AVG_X_H(EPMU_GT(IJKW),EPMU_GT(IJK),IM),&
314                                  AVG_X_H(EPMU_GT(IJKNW),EPMU_GT(IJKN),IM),J)*&
315                          (U_G(IMJPK)-U_G(IMJK))*ONEoDY_N_U(IMJK)*AYZ_V(IMJK)
316     
317                    SSY = EPMU_GT(IJKN)*(V_G(IJPK)-V_G(IJK))*&
318                             ONEoDY_N_V(IJK)*AXZ_V(IJK) - &
319                          EPMU_GT(IJK)*(V_G(IJK)-V_G(IJMK))*&
320                             ONEoDY_N_V(IJMK)*AXZ_V(IJMK)
321     
322                    IF(DO_K) THEN
323                       SSZ = AVG_Y_H(AVG_Z_H(EPMU_GT(IJK),EPMU_GT(IJKT),K),&
324                                     AVG_Z_H(EPMU_GT(IJKN),EPMU_GT(IJKTN),K),J)*&
325                             (W_G(IJPK)-W_G(IJK))*ONEoDY_N_W(IJK)*AXY_V(IJK) - &
326                             AVG_Y_H(AVG_Z_H(EPMU_GT(IJKB),EPMU_GT(IJK),KM),&
327                                     AVG_Z_H(EPMU_GT(IJKBN),EPMU_GT(IJKN),KM),J)*&
328                             (W_G(IJPKM)-W_G(IJKM))*ONEoDY_N_W(IJKM)*AXY_V(IJKM)
329                    ELSE
330                       SSZ = ZERO
331                    ENDIF
332     
333     ! cut cell modifications
334     !---------------------------------------------------------------------//
335                 ELSE
336     
337                    BCV = BC_V_ID(IJK)
338                    IF(BCV > 0 ) THEN
339                       BCT = BC_TYPE_ENUM(BCV)
340                    ELSE
341                       BCT = NONE
342                    ENDIF
343     
344                    SELECT CASE (BCT)
345                       CASE (CG_NSW,CG_MI)
346                          CUT_TAU_VG = .TRUE.
347                          NOC_VG     = .TRUE.
348                          UW_g = ZERO
349                          VW_g = ZERO
350                          WW_g = ZERO
351                       CASE (CG_FSW)
352                          CUT_TAU_VG = .FALSE.
353                          NOC_VG     = .FALSE.
354                          UW_g = ZERO
355                          VW_g = ZERO
356                          WW_g = ZERO
357                       CASE(CG_PSW)
358                          IF(BC_HW_G(BC_V_ID(IJK))==UNDEFINED) THEN   ! same as NSW
359                             CUT_TAU_VG = .TRUE.
360                             NOC_VG     = .TRUE.
361                             UW_g = BC_UW_G(BCV)
362                             VW_g = BC_VW_G(BCV)
363                             WW_g = BC_WW_G(BCV)
364                          ELSEIF(BC_HW_G(BC_V_ID(IJK))==ZERO) THEN   ! same as FSW
365                             CUT_TAU_VG = .FALSE.
366                             NOC_VG     = .FALSE.
367                             UW_g = ZERO
368                             VW_g = ZERO
369                             WW_g = ZERO
370                          ELSE                              ! partial slip
371                             CUT_TAU_VG = .FALSE.
372                             NOC_VG     = .FALSE.
373                          ENDIF
374                       CASE (NONE)
375                          lTAU_V_G(IJK) = ZERO
376                          CYCLE
377                    END SELECT
378     
379     
380                    IF(CUT_TAU_VG) THEN
381                       MU_GT_CUT = (VOL(IJK)*EPMU_GT(IJK) + &
382                                    VOL(IJPK)*EPMU_GT(IJKN))/&
383                                   (VOL(IJK) + VOL(IJPK))
384                    ELSE
385                       MU_GT_CUT = ZERO
386                    ENDIF
387     
388     ! SSX:
389                    U_NODE_AT_NE = ((.NOT.BLOCKED_U_CELL_AT(IJPK)).AND.&
390                                    (.NOT.WALL_U_AT(IJPK)))
391                    U_NODE_AT_SE = ((.NOT.BLOCKED_U_CELL_AT(IJK)).AND.&
392                                    (.NOT.WALL_U_AT(IJK)))
393                    U_NODE_AT_NW = ((.NOT.BLOCKED_U_CELL_AT(IMJPK)).AND.&
394                                    (.NOT.WALL_U_AT(IMJPK)))
395                    U_NODE_AT_SW = ((.NOT.BLOCKED_U_CELL_AT(IMJK)).AND.&
396                                    (.NOT.WALL_U_AT(IMJK)))
397     
398                    IF(U_NODE_AT_NE.AND.U_NODE_AT_SE) THEN
399                       Ui = HALF * (U_G(IJPK) + U_G(IJK))
400                       Xi = HALF * (X_U(IJPK) + X_U(IJK))
401                       Yi = HALF * (Y_U(IJPK) + Y_U(IJK))
402                       Zi = HALF * (Z_U(IJPK) + Z_U(IJK))
403                       Sx = X_U(IJPK) - X_U(IJK)
404                       Sy = Y_U(IJPK) - Y_U(IJK)
405                       Sz = Z_U(IJPK) - Z_U(IJK)
406                       CALL GET_DEL_H(IJK, 'V_MOMENTUM', Xi, Yi, Zi, &
407                          Del_H, Nx, Ny, Nz)
408     
409                       dudy_at_E =  (U_G(IJPK) - U_G(IJK)) * ONEoDY_N_U(IJK)
410                       IF(NOC_VG) dudy_at_E = dudy_at_E - ((Ui-UW_g) * &
411                          ONEoDY_N_U(IJK)/DEL_H * (Sx*Nx+Sz*Nz) )
412                    ELSE
413                       dudy_at_E =  ZERO
414                    ENDIF
415     
416                    IF(U_NODE_AT_NW.AND.U_NODE_AT_SW) THEN
417                       Ui = HALF * (U_G(IMJPK) + U_G(IMJK))
418                       Xi = HALF * (X_U(IMJPK) + X_U(IMJK))
419                       Yi = HALF * (Y_U(IMJPK) + Y_U(IMJK))
420                       Zi = HALF * (Z_U(IMJPK) + Z_U(IMJK))
421                       Sx = X_U(IMJPK) - X_U(IMJK)
422                       Sy = Y_U(IMJPK) - Y_U(IMJK)
423                       Sz = Z_U(IMJPK) - Z_U(IMJK)
424                       CALL GET_DEL_H(IJK, 'V_MOMENTUM', Xi, Yi, Zi, &
425                          Del_H, Nx, Ny, Nz)
426     
427                       dudy_at_W =  (U_G(IMJPK) - U_G(IMJK)) * ONEoDY_N_U(IMJK)
428                       IF(NOC_VG) dudy_at_W = dudy_at_W - ((Ui-UW_g) * &
429                          ONEoDY_N_U(IMJK)/DEL_H * (Sx*Nx+Sz*Nz) )
430                    ELSE
431                       dudy_at_W =  ZERO
432                    ENDIF
433     
434                    IF(U_NODE_AT_SE) THEN
435                       CALL GET_DEL_H(IJK,'V_MOMENTUM', X_U(IJK), Y_U(IJK), &
436                          Z_U(IJK), Del_H, Nx, Ny, Nz)
437                       SSX_CUT = - MU_GT_CUT * (U_G(IJK) - UW_g) / DEL_H * &
438                          (Ny*Nx) * Area_V_CUT(IJK)
439                    ELSE
440                       SSX_CUT =  ZERO
441                    ENDIF
442     
443                    SSX = AVG_Y_H(AVG_X_H(EPMU_GT(IJK),EPMU_GT(IJKE),I),&
444                                  AVG_X_H(EPMU_GT(IJKN),EPMU_GT(IJKNE),I),J)*&
445                          dudy_at_E*AYZ_V(IJK) - &
446                          AVG_Y_H(AVG_X_H(EPMU_GT(IJKW),EPMU_GT(IJK),IM),&
447                                  AVG_X_H(EPMU_GT(IJKNW),EPMU_GT(IJKN),IM),J)*&
448                          dudy_at_W*AYZ_V(IMJK) + SSX_CUT
449     
450     ! SSY:
451                    CALL GET_DEL_H(IJK, 'V_MOMENTUM', X_V(IJK), Y_V(IJK), &
452                       Z_V(IJK), Del_H, Nx, Ny, Nz)
453     
454                    SSY = EPMU_GT(IJKN)*(V_G(IJPK)-V_G(IJK))*&
455                             ONEoDY_N_V(IJK)*AXZ_V(IJK) - &
456                          EPMU_GT(IJK)*(V_G(IJK)-V_G(IJMK))*&
457                             ONEoDY_N_V(IJMK)*AXZ_V(IJMK) - &
458                          MU_GT_CUT * (V_g(IJK) - VW_g) / DEL_H * &
459                             (Ny**2) * Area_V_CUT(IJK)
460     
461     ! SSZ:
462                    IF(DO_K) THEN
463                       W_NODE_AT_TN = ((.NOT.BLOCKED_W_CELL_AT(IJPK)).AND.&
464                                       (.NOT.WALL_W_AT(IJPK)))
465                       W_NODE_AT_TS = ((.NOT.BLOCKED_W_CELL_AT(IJK)).AND.&
466                                       (.NOT.WALL_W_AT(IJK)))
467                       W_NODE_AT_BN = ((.NOT.BLOCKED_W_CELL_AT(IJPKM)).AND.&
468                                       (.NOT.WALL_W_AT(IJPKM)))
469                       W_NODE_AT_BS = ((.NOT.BLOCKED_W_CELL_AT(IJKM)).AND.&
470                                       (.NOT.WALL_W_AT(IJKM)))
471     
472                       IF(W_NODE_AT_TN.AND.W_NODE_AT_TS) THEN
473                             Wi = HALF * (W_G(IJPK) + W_G(IJK))
474                             Xi = HALF * (X_W(IJPK) + X_W(IJK))
475                             Yi = HALF * (Y_W(IJPK) + Y_W(IJK))
476                             Zi = HALF * (Z_W(IJPK) + Z_W(IJK))
477                             Sx = X_W(IJPK) - X_W(IJK)
478                             Sy = Y_W(IJPK) - Y_W(IJK)
479                             Sz = Z_W(IJPK) - Z_W(IJK)
480                             CALL GET_DEL_H(IJK, 'V_MOMENTUM', Xi, Yi, Zi, &
481                                Del_H, Nx, Ny, Nz)
482     
483                             dwdy_at_T = (W_G(IJPK)-W_G(IJK)) * &
484                                ONEoDY_N_W(IJK)
485                             IF(NOC_VG) dwdy_at_T = dwdy_at_T - ((Wi-WW_g) * &
486                                ONEoDY_N_W(IJK)/DEL_H*(Sx*Nx+Sz*Nz))
487                       ELSE
488                          dwdy_at_T =  ZERO
489                       ENDIF
490     
491                       IF(W_NODE_AT_BN.AND.W_NODE_AT_BS) THEN
492                          Wi = HALF * (W_G(IJPKM) + W_G(IJKM))
493                          Xi = HALF * (X_W(IJPKM) + X_W(IJKM))
494                          Yi = HALF * (Y_W(IJPKM) + Y_W(IJKM))
495                          Zi = HALF * (Z_W(IJPKM) + Z_W(IJKM))
496                          Sx = X_W(IJPKM) - X_W(IJKM)
497                          Sy = Y_W(IJPKM) - Y_W(IJKM)
498                          Sz = Z_W(IJPKM) - Z_W(IJKM)
499                          CALL GET_DEL_H(IJK, 'V_MOMENTUM', Xi, Yi, Zi,&
500                             Del_H, Nx, Ny, Nz)
501     
502                          dwdy_at_B =  (W_G(IJPKM) - W_G(IJKM)) * &
503                             ONEoDY_N_W(IJKM)
504                          IF(NOC_VG) dwdy_at_B = dwdy_at_B - ((Wi-WW_g) * &
505                             ONEoDY_N_W(IJKM)/DEL_H*(Sx*Nx+Sz*Nz))
506                       ELSE
507                          dwdy_at_B =  ZERO
508                       ENDIF
509     
510                       IF(W_NODE_AT_TS) THEN
511                          CALL GET_DEL_H(IJK, 'V_MOMENTUM', X_W(IJK), &
512                             Y_W(IJK), Z_W(IJK), Del_H, Nx, Ny, Nz)
513                          SSZ_CUT = - MU_GT_CUT * (W_G(IJK) - WW_g) / &
514                             DEL_H * (Ny*Nz) * Area_V_CUT(IJK)
515                       ELSE
516                          SSZ_CUT =  ZERO
517                       ENDIF
518     
519                       SSZ = AVG_Y_H(AVG_Z_H(EPMU_GT(IJK),EPMU_GT(IJKT),K),&
520                                     AVG_Z_H(EPMU_GT(IJKN),EPMU_GT(IJKTN),K),J)*&
521                                dwdy_at_T*AXY_V(IJK) - &
522                             AVG_Y_H(AVG_Z_H(EPMU_GT(IJKB),EPMU_GT(IJK),KM),&
523                                     AVG_Z_H(EPMU_GT(IJKBN),EPMU_GT(IJKN),KM),J)*&
524                                dwdy_at_B*AXY_V(IJKM) + SSZ_CUT
525                    ELSE
526                       SSZ = ZERO
527                    ENDIF  ! DO_K
528     
529                 ENDIF   ! end if/else cut_v_cell_at
530     
531     ! Add the terms
532                 lTAU_V_G(IJK) = SBV + SSX + SSY + SSZ
533     
534     ! Also calculate and store the full gas phase viscous stress tensor
535                 CALL GET_FULL_TAU_V_G(IJK, ltau_v_g, lctau_v_g)
536              ELSE
537                 lTAU_V_G(IJK) = ZERO
538                 lctau_v_g(IJK) = ZERO
539              ENDIF   ! end if/else (.NOT. IP_AT_N(IJK) .AND. EPGA>DIL_EP_S)
540           ENDDO   ! end do ijk
541     !$omp end parallel do
542     
543           RETURN
544           END SUBROUTINE CALC_CG_TAU_V_G
545     
546     
547     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
548     !                                                                      C
549     !  Subroutine: GET_FULL_Tau_V_g                                        C
550     !  Purpose: Calculate the divergence of complete gas phase stress      C
551     !  tensor including all terms in v-direction                                            C
552     !                                                                      C
553     !                                                                      C
554     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
555           SUBROUTINE GET_FULL_TAU_V_G(IJK, ltau_v_g, lctau_v_g)
556     
557     ! Modules
558     !---------------------------------------------------------------------//
559           USE fldvar, only: v_g
560     
561           USE functions, only: flow_at_n
562           USE functions, only: im_of, jm_of, km_of
563           USE functions, only: ip_of, jp_of, kp_of
564     
565           USE geometry, only: do_k
566     
567           USE param
568           USE param1, only: zero
569           USE visc_g, only: df_gv
570           IMPLICIT NONE
571     
572     ! Dummy arguments
573     !---------------------------------------------------------------------//
574     ! ijk index
575           INTEGER, INTENT(IN) :: IJK
576     ! TAU_V_g
577           DOUBLE PRECISION, INTENT(IN) :: lTAU_V_g(DIMENSION_3)
578     ! cTAU_V_g
579           DOUBLE PRECISION, INTENT(INOUT) :: lcTAU_V_g(DIMENSION_3)
580     
581     ! Local Variables
582     !---------------------------------------------------------------------//
583     ! indices
584           INTEGER :: IPJK, IJPK, IJKP, IMJK, IJMK, IJKM
585     ! source terms
586           DOUBLE PRECISION :: SSX, SSY, SSZ
587     !---------------------------------------------------------------------//
588     
589           IPJK = IP_OF(IJK)
590           IJPK = JP_OF(IJK)
591           IJKP = KP_OF(IJK)
592           IMJK = IM_OF(IJK)
593           IJMK = JM_OF(IJK)
594           IJKM = KM_OF(IJK)
595     
596     ! convection terms: see conv_dif_v_g
597           SSX = ZERO
598           SSY = ZERO
599           SSZ = ZERO
600           IF (FLOW_AT_N(IJK)) THEN
601     ! part of 1/x d/dx (x.tau_xx) xdxdydz =>
602     !         1/x d/dx (x.mu.dv/dx) xdxdydz =>
603     ! delta (mu.dv/dx.Ayz) |E-W : at (i+1/2 - i-1/2), j+1/2, k
604              SSX = DF_GV(IJK,east)*(V_G(IPJK) - V_G(IJK)) - &
605                    DF_GV(IJK,west)*(V_G(IJK) - V_G(IJKM))
606     
607     ! part of d/dy (tau_xy) xdxdydz =>
608     !         d/dy (mu.dv/dy) xdxdydz =>
609     ! delta (mu.dv/dy.Axz) |N-S : at (i, j+1 - j-1, k)
610              SSY = DF_GV(IJK,north)*(V_G(IJPK)-V_G(IJK)) - &
611                    DF_GV(IJK,south)*(V_G(IJK)-V_G(IJMK))
612     
613              IF (DO_K) THEN
614     ! part of 1/x d/dz (tau_xz) xdxdydz =>
615     !         1/x d/dz (mu/x.dv/dz) xdxdydz =>
616     ! delta (mu/x.dv/dz.Axy) |T-B : at (i, j+1/2, k+1/2 - k-1/2)
617                 SSZ = DF_GV(IJK,top)*(V_G(IJKP)-V_G(IJK)) - &
618                       DF_GV(IJK,bottom)*(V_G(IJK)-V_G(IJKM))
619              ENDIF
620           ENDIF   ! end if flow_at_n
621     
622     ! Add the terms
623           lctau_v_g(IJK) = (lTAU_v_G(IJK) + SSX + SSY + SSZ)
624     
625           RETURN
626           END SUBROUTINE GET_FULL_TAU_V_G
627     
628