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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Subroutine: CALC_Tau_V_s                                            C
4     !  Purpose: Cross terms in the gradient of stress in V_s momentum      C
5     !                                                                      C
6     !                                                                      C
7     !  Comments: This routine calculates the components of the solids      C
8     !  phase viscous stress tensor of the v-momentum equation that cannot  C
9     !  be cast in the form: mu.grad(v). These components are stored in the C
10     !  passed variable, which is then used as a source of the v-momentum   C
11     !  equation.                                                           C
12     !  For greater details see calc_tau_v_g.                               C
13     !                                                                      C
14     !  Author: M. Syamlal                                 Date: 19-DEC-96  C
15     !  Reviewer:                                          Date:            C
16     !                                                                      C
17     !                                                                      C
18     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
19           SUBROUTINE CALC_TAU_V_S(lTAU_V_S)
20     
21     ! Modules
22     !---------------------------------------------------------------------//
23           USE param, only: dimension_3, dimension_m
24           USE param1, only: zero
25           USE physprop, only: smax, mmax
26           USE fldvar, only: ep_s
27           USE visc_s, only: epmu_s
28           USE toleranc, only: dil_ep_s
29           USE run, only: kt_type_enum, ghd_2007
30           USE run, only: shear, V_sh
31           USE cutcell, only: cartesian_grid, cg_safe_mode
32     
33           USE compar
34           USE geometry
35           USE indices
36           USE functions
37           USE fun_avg
38           USE sendrecv
39           IMPLICIT NONE
40     
41     ! Dummy arguments
42     !---------------------------------------------------------------------//
43     ! TAU_V_s
44           DOUBLE PRECISION, INTENT(OUT) :: lTAU_V_s(DIMENSION_3, DIMENSION_M)
45     
46     ! Local Variables
47     !---------------------------------------------------------------------//
48     ! Indices
49           INTEGER :: IJK, J, IJKN
50     ! Additional indices needed for shear case
51           INTEGER :: I, IJKE, IJKW, IJKNE, IJKNW
52     ! Phase index
53           INTEGER :: M, L
54     ! Average volume fraction
55           DOUBLE PRECISION :: EPSA, EPStmp
56     ! Source terms (Surface)
57           DOUBLE PRECISION :: Sbv, Ssx, Ssy, Ssz
58     ! Shearing variables
59           DOUBLE PRECISION :: Source_diff, Diffco_e, Diffco_w
60     !---------------------------------------------------------------------//
61     
62           DO M = 1, MMAX
63              IF(KT_TYPE_ENUM == GHD_2007 .AND. M /= MMAX) CYCLE
64     
65     !!$omp  parallel do &
66     !!$omp  private(J, IJK, IJKN, I, IJKE, IJKW, IJKNE, IJKNW,            &
67     !!$omp&         EPSA, EPStmp, SSX, SSY, SSZ, SBV,                     &
68     !!$omp&         Source_diff, Diffco_e, Diffco_w)
69     !!$omp&  schedule(static)
70              DO IJK = IJKSTART3, IJKEND3
71     
72     ! Skip walls where some values are undefined.
73                 IF(WALL_AT(IJK)) cycle
74                 J = J_OF(IJK)
75                 IJKN = NORTH_OF(IJK)
76     
77                 IF (KT_TYPE_ENUM == GHD_2007) THEN
78                    EPStmp = ZERO
79                    DO L = 1, SMAX
80                       EPStmp = EPStmp + AVG_Y(EP_S(IJK,L),EP_S(IJKN,L),J)
81                    ENDDO
82                    EPSA = EPStmp
83                 ELSE
84                    EPSA = AVG_Y(EP_S(IJK,M),EP_S(IJKN,M),J)
85                 ENDIF
86     
87                 IF ( .NOT.SIP_AT_N(IJK) .AND. EPSA>DIL_EP_S) THEN
88     
89                    IF((.NOT.CARTESIAN_GRID).OR.(CG_SAFE_MODE(4)==1)) THEN
90                       CALL CALC_REG_TAU_V_S(IJK, M, SSX, SSY, SSZ, SBV)
91                    ELSE
92                       CALL CALC_CG_TAU_V_S(IJK, M, SSX, SSY, SSZ, SBV)
93                    ENDIF
94     
95     
96     ! Source terms from shear stress
97                    IF (SHEAR) THEN
98                       I = I_OF(IJK)
99                       IJKE = EAST_OF(IJK)
100                       IJKN = NORTH_OF(IJK)
101                       IJKNE = EAST_OF(IJKN)
102                       IJKW = WEST_OF(IJK)
103                       IJKNW = NORTH_OF(IJKW)
104                       Diffco_e = AVG_Y_H((AVG_X_H(EPMU_S(IJK,m),EPMU_S(IJKE,m),I)),&
105                                          (AVG_X_H(EPMU_S(IJKN,m),EPMU_S(IJKNE,m),I)),J)*AYZ_V(IJK)
106                       Diffco_w=AVG_Y_H((AVG_X_H(EPMU_S(IJK,m),EPMU_S(IJKW,m),I)),&
107                                        (AVG_X_H(EPMU_S(IJKN,m),EPMU_S(IJKNW,m),I)),J)*AYZ_V(IJKW)
108                       Source_diff=(2d0*V_sh/XLENGTH)*(Diffco_e-Diffco_w)
109                    ELSE
110                       Source_diff=0d0
111                    ENDIF
112     
113     ! Add the terms
114                    lTAU_V_S(IJK,M) = SBV + SSX + SSY + SSZ + Source_diff
115                 ELSE
116                    lTAU_V_S(IJK,M) = ZERO
117                 ENDIF   ! end if/else .not.sip_at_n .and. epsa>dil_ep_s
118              ENDDO   ! end do ijk
119     !!$omp end parallel do
120     
121           ENDDO   ! end do m=1,mmax
122     
123           call send_recv(ltau_v_s,2)
124           RETURN
125           END SUBROUTINE CALC_TAU_V_S
126     
127     
128     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
129     !                                                                      C
130     !  Subroutine: calc_reg_tau_v_s                                        C
131     !  Purpose: Cross terms in the gradient of stress in V_s momentum      C
132     !  based on NON cartesian grid cut cell.                               C
133     !                                                                      C
134     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
135           SUBROUTINE CALC_REG_TAU_V_S(IJK, M, SSX, SSY, SSZ, SBV)
136     
137     ! Modules
138     !---------------------------------------------------------------------//
139           USE fldvar, only: u_s, v_s, w_s
140           USE visc_s, only: epmu_s, eplambda_s
141           USE visc_s, only: trd_s
142     
143           USE fun_avg
144           USE compar
145           USE geometry
146           USE indices
147           USE functions
148           IMPLICIT NONE
149     
150     ! Dummy arguments
151     !---------------------------------------------------------------------//
152     ! current ijk index
153           INTEGER, INTENT(IN) :: IJK
154     ! solids phase index
155           INTEGER, INTENT(IN) :: M
156     ! Source terms (surface)
157           DOUBLE PRECISION, INTENT(OUT) :: SSX
158           DOUBLE PRECISION, INTENT(OUT) :: SSY
159           DOUBLE PRECISION, INTENT(OUT) :: SSZ
160           DOUBLE PRECISION, INTENT(OUT) :: SBV
161     
162     ! Local variables
163     !---------------------------------------------------------------------//
164     ! Indices
165           INTEGER :: I, J, K, IM, KM, JP
166           INTEGER :: IJKE, IJKW, IJKN, IJKT, IJKB
167           INTEGER :: IJKNE, IJKNW, IJKTN, IJKBN
168           INTEGER :: IJPK, IJMK, IMJK, IJKM
169           INTEGER :: IMJPK, IJPKM
170     !---------------------------------------------------------------------//
171     
172           J = J_OF(IJK)
173           JP = JP1(J)
174           I = I_OF(IJK)
175           IM = IM1(I)
176           K = K_OF(IJK)
177           KM = KM1(K)
178     
179           IJPK = JP_OF(IJK)
180           IJMK = JM_OF(IJK)
181           IJKM = KM_OF(IJK)
182           IMJK = IM_OF(IJK)
183           IMJPK = IM_OF(IJPK)
184           IJPKM = JP_OF(IJKM)
185     
186           IJKE = EAST_OF(IJK)
187           IJKN = NORTH_OF(IJK)
188           IJKNE = EAST_OF(IJKN)
189           IJKW = WEST_OF(IJK)
190           IJKNW = NORTH_OF(IJKW)
191           IJKT = TOP_OF(IJK)
192           IJKTN = NORTH_OF(IJKT)
193           IJKB = BOTTOM_OF(IJK)
194           IJKBN = NORTH_OF(IJKB)
195     
196     ! Surface forces
197     
198     ! bulk viscosity term
199     ! part of d/dy (tau_yy) xdxdydz =>
200     !         d/dy (lambda.trcD) xdxdydz =>
201     ! delta (lambda.trcD)Ap|N-S  : at (i, j+1 - j-1, k)
202           SBV = (EPLAMBDA_S(IJKN,M)*TRD_S(IJKN,M)-&
203                  EPLAMBDA_S(IJK,M)*TRD_S(IJK,M))*AXZ(IJK)
204     
205     ! shear stress terms at i, j+1/2, k
206     ! part of 1/x d/dx(x.tau_xy) xdxdydz =>
207     !         1/x d/dx (x.mu.du/dy) xdxdydz =>
208     ! delta (x.mu.du/dy)Ayz |E-W : at (i+1/2 - i-1/2, j+1/2, k)
209           SSX = AVG_Y_H(AVG_X_H(EPMU_S(IJK,M),EPMU_S(IJKE,M),I),&
210                         AVG_X_H(EPMU_S(IJKN,M),EPMU_S(IJKNE,M),I),J)*&
211                    (U_S(IJPK,M)-U_S(IJK,M))*ODY_N(J)*AYZ_V(IJK) - &
212                 AVG_Y_H(AVG_X_H(EPMU_S(IJKW,M),EPMU_S(IJK,M),IM),&
213                         AVG_X_H(EPMU_S(IJKNW,M),EPMU_S(IJKN,M),IM),J)*&
214                    (U_S(IMJPK,M)-U_S(IMJK,M))*ODY_N(J)*AYZ_V(IMJK)
215     
216     ! part of d/dy (tau_xy) xdxdydz =>
217     !         d/dy (mu.dv/dy) xdxdydz =>
218     ! delta (mu.dv/dx)Axz |N-S : at (i, j+1 - j-1, k)
219           SSY = EPMU_S(IJKN,M)*(V_S(IJPK,M)-V_S(IJK,M))*ODY(JP)*AXZ_V(IJK) - &
220                 EPMU_S(IJK,M)*(V_S(IJK,M)-V_S(IJMK,M))*ODY(J)*AXZ_V(IJMK)
221     
222     ! part of 1/x d/dz (tau_xz) xdxdydz =>
223     !         1/x d/dz (mu.dw/dy) xdxdydz =>
224     ! delta (mu.dw/dx)Axy |T-B : at (i, j+1/2, k+1/2 - k-1/2)
225           SSZ = AVG_Y_H(AVG_Z_H(EPMU_S(IJK,M),EPMU_S(IJKT,M),K),&
226                         AVG_Z_H(EPMU_S(IJKN,M),EPMU_S(IJKTN,M),K),J)*&
227                    (W_S(IJPK,M)-W_S(IJK,M))*ODY_N(J)*AXY_V(IJK) - &
228                 AVG_Y_H(AVG_Z_H(EPMU_S(IJKB,M),EPMU_S(IJK,M),KM),&
229                         AVG_Z_H(EPMU_S(IJKBN,M),EPMU_S(IJKN,M),KM),J)*&
230                    (W_S(IJPKM,M)-W_S(IJKM,M))*ODY_N(J)*AXY_V(IJKM)
231     
232           RETURN
233           END SUBROUTINE CALC_REG_TAU_V_S
234     
235     
236     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
237     !                                                                      C
238     !  Subroutine: calc_cg_tau_v_s                                         C
239     !  Purpose: Cross terms in the gradient of stress in V_s momentum      C
240     !  based on cartesian grid cut cell.                                   C
241     !                                                                      C
242     !  Author: Jeff Dietiker                              Date: 01-Jul-09  C
243     !                                                                      C
244     !                                                                      C
245     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
246           SUBROUTINE CALC_CG_TAU_V_S(IJK, M, SSX, SSY, SSZ, SBV)
247     
248     ! Modules
249     !---------------------------------------------------------------------//
250           USE param1, only: zero, half
251           USE fldvar, only: u_s, v_s, w_s
252           USE visc_s, only: epmu_s, eplambda_s
253           USE visc_s, only: trd_s
254     
255           USE fun_avg
256           USE compar
257           USE geometry
258           USE indices
259           USE functions
260     
261     ! for cutcell:
262     ! wall velocity for slip conditions
263           USE bc, only: bc_hw_s, bc_uw_s, bc_vw_s, bc_ww_s
264           USE bc, only: bc_type
265           USE quadric
266           USE cutcell
267           IMPLICIT NONE
268     
269     ! Dummy arguments
270     !---------------------------------------------------------------------//
271     ! current ijk index
272           INTEGER, INTENT(IN) :: IJK
273     ! solids phase index
274           INTEGER, INTENT(IN) :: M
275     ! Source terms (surface)
276           DOUBLE PRECISION, INTENT(OUT) :: SSX
277           DOUBLE PRECISION, INTENT(OUT) :: SSY
278           DOUBLE PRECISION, INTENT(OUT) :: SSZ
279           DOUBLE PRECISION, INTENT(OUT) :: SBV
280     
281     ! Local variables
282     !---------------------------------------------------------------------//
283     ! Indices
284           INTEGER :: I, J, K, IM, KM, JP
285           INTEGER :: IJKE, IJKW, IJKN, IJKT, IJKB
286           INTEGER :: IJKNE, IJKNW, IJKTN, IJKBN
287           INTEGER :: IJPK, IJMK, IMJK, IJKM
288           INTEGER :: IMJPK, IJPKM
289     
290     ! for cartesian grid implementation:
291           DOUBLE PRECISION :: DEL_H,Nx,Ny,Nz
292           LOGICAL :: U_NODE_AT_NE,U_NODE_AT_NW,U_NODE_AT_SE,U_NODE_AT_SW
293           LOGICAL :: W_NODE_AT_TN,W_NODE_AT_TS,W_NODE_AT_BN,W_NODE_AT_BS
294           DOUBLE PRECISION :: dudy_at_E,dudy_at_W
295           DOUBLE PRECISION :: dwdy_at_T,dwdy_at_B
296           DOUBLE PRECISION :: Xi,Yi,Zi,Ui,Wi,Sx,Sy,Sz
297           DOUBLE PRECISION :: MU_S_CUT,SSX_CUT,SSZ_CUT
298           DOUBLE PRECISION :: UW_s,VW_s,WW_s
299           INTEGER :: BCV
300           CHARACTER(LEN=9) :: BCT
301     !---------------------------------------------------------------------//
302     
303           J = J_OF(IJK)
304           JP = JP1(J)
305           I = I_OF(IJK)
306           IM = IM1(I)
307           K = K_OF(IJK)
308           KM = KM1(K)
309     
310           IJPK = JP_OF(IJK)
311           IJMK = JM_OF(IJK)
312           IJKM = KM_OF(IJK)
313           IMJK = IM_OF(IJK)
314           IMJPK = IM_OF(IJPK)
315           IJPKM = JP_OF(IJKM)
316     
317           IJKE = EAST_OF(IJK)
318           IJKN = NORTH_OF(IJK)
319           IJKNE = EAST_OF(IJKN)
320           IJKW = WEST_OF(IJK)
321           IJKNW = NORTH_OF(IJKW)
322           IJKT = TOP_OF(IJK)
323           IJKTN = NORTH_OF(IJKT)
324           IJKB = BOTTOM_OF(IJK)
325           IJKBN = NORTH_OF(IJKB)
326     
327     ! Surface forces
328     ! bulk viscosity term
329           SBV = (EPLAMBDA_S(IJKN,M)*TRD_S(IJKN,M)) * AXZ_V(IJK) - &
330                 (EPLAMBDA_S(IJK,M) *TRD_S(IJK,M) ) * AXZ_V(IJMK)
331     
332     ! shear stress terms
333           IF(.NOT.CUT_V_CELL_AT(IJK)) THEN
334              SSX = AVG_Y_H(AVG_X_H(EPMU_S(IJK,M),EPMU_S(IJKE,M),I),&
335                            AVG_X_H(EPMU_S(IJKN,M),EPMU_S(IJKNE,M),I),J)*&
336                       (U_S(IJPK,M)-U_S(IJK,M))*ONEoDY_N_U(IJK)*AYZ_V(IJK) - &
337                    AVG_Y_H(AVG_X_H(EPMU_S(IJKW,M),EPMU_S(IJK,M),IM),&
338                            AVG_X_H(EPMU_S(IJKNW,M),EPMU_S(IJKN,M),IM),J)*&
339                       (U_S(IMJPK,M)-U_S(IMJK,M))*ONEoDY_N_U(IMJK)*AYZ_V(IMJK)
340              SSY = EPMU_S(IJKN,M)*(V_S(IJPK,M)-V_S(IJK,M))*&
341                       ONEoDY_N_V(IJK)*AXZ_V(IJK) - &
342                    EPMU_S(IJK,M)*(V_S(IJK,M)-V_S(IJMK,M))*&
343                       ONEoDY_N_V(IJMK)*AXZ_V(IJMK)
344              IF(DO_K) THEN
345                SSZ = AVG_Y_H(AVG_Z_H(EPMU_S(IJK,M),EPMU_S(IJKT,M),K),&
346                              AVG_Z_H(EPMU_S(IJKN,M),EPMU_S(IJKTN,M),K),J)*&
347                         (W_S(IJPK,M)-W_S(IJK,M))*ONEoDY_N_W(IJK)*AXY_V(IJK) - &
348                      AVG_Y_H(AVG_Z_H(EPMU_S(IJKB,M),EPMU_S(IJK,M),KM),&
349                              AVG_Z_H(EPMU_S(IJKBN,M),EPMU_S(IJKN,M),KM),J)*&
350                         (W_S(IJPKM,M)-W_S(IJKM,M))*ONEoDY_N_W(IJKM)*AXY_V(IJKM)
351              ELSE
352                SSZ = ZERO
353              ENDIF
354     
355     ! cut cell modifications
356     !---------------------------------------------------------------------//
357           ELSE
358              BCV = BC_V_ID(IJK)
359              IF(BCV > 0 ) THEN
360                BCT = BC_TYPE(BCV)
361              ELSE
362                BCT = 'NONE'
363              ENDIF
364     
365              SELECT CASE (BCT)
366                 CASE ('CG_NSW')
367                    CUT_TAU_VS = .TRUE.
368                    NOC_VS     = .TRUE.
369                    UW_s = ZERO
370                    VW_s = ZERO
371                    WW_s = ZERO
372                 CASE ('CG_FSW')
373                    CUT_TAU_VS = .FALSE.
374                    NOC_VS     = .FALSE.
375                    UW_s = ZERO
376                    VW_s = ZERO
377                    WW_s = ZERO
378                 CASE('CG_PSW')
379                    IF(BC_HW_S(BC_V_ID(IJK),M)==UNDEFINED) THEN   ! same as NSW
380                       CUT_TAU_VS = .TRUE.
381                       NOC_VS     = .TRUE.
382                       UW_s = BC_UW_S(BCV,M)
383                       VW_s = BC_VW_S(BCV,M)
384                       WW_s = BC_WW_S(BCV,M)
385                    ELSEIF(BC_HW_S(BC_V_ID(IJK),M)==ZERO) THEN   ! same as FSW
386                       CUT_TAU_VS = .FALSE.
387                       NOC_VS     = .FALSE.
388                       UW_s = ZERO
389                       VW_s = ZERO
390                       WW_s = ZERO
391                    ELSE                              ! partial slip
392                       CUT_TAU_VS = .FALSE.
393                       NOC_VS     = .FALSE.
394                    ENDIF
395                 CASE ('NONE')
396     ! equivalent of setting tau_v_s to zero at this ijk cell
397                    SSX = ZERO
398                    SSY = ZERO
399                    SSZ = ZERO
400                    SBV = ZERO
401                    RETURN
402              END SELECT
403     
404              IF(CUT_TAU_VS) THEN
405                 MU_S_CUT = (VOL(IJK)*EPMU_S(IJK,M) + &
406                             VOL(IJPK)*EPMU_S(IJKN,M))/&
407                            (VOL(IJK) + VOL(IJPK))
408              ELSE
409                 MU_S_CUT = ZERO
410              ENDIF
411     
412     ! SSX:
413              U_NODE_AT_NE = ((.NOT.BLOCKED_U_CELL_AT(IJPK)).AND.&
414                              (.NOT.WALL_U_AT(IJPK)))
415              U_NODE_AT_SE = ((.NOT.BLOCKED_U_CELL_AT(IJK)).AND.&
416                              (.NOT.WALL_U_AT(IJK)))
417              U_NODE_AT_NW = ((.NOT.BLOCKED_U_CELL_AT(IMJPK)).AND.&
418                              (.NOT.WALL_U_AT(IMJPK)))
419              U_NODE_AT_SW = ((.NOT.BLOCKED_U_CELL_AT(IMJK)).AND.&
420                              (.NOT.WALL_U_AT(IMJK)))
421     
422              IF(U_NODE_AT_NE.AND.U_NODE_AT_SE) THEN
423                 Ui = HALF * (U_S(IJPK,M) + U_S(IJK,M))
424                 Xi = HALF * (X_U(IJPK) + X_U(IJK))
425                 Yi = HALF * (Y_U(IJPK) + Y_U(IJK))
426                 Zi = HALF * (Z_U(IJPK) + Z_U(IJK))
427                 Sx = X_U(IJPK) - X_U(IJK)
428                 Sy = Y_U(IJPK) - Y_U(IJK)
429                 Sz = Z_U(IJPK) - Z_U(IJK)
430                 CALL GET_DEL_H(IJK, 'V_MOMENTUM', Xi, Yi, Zi, Del_H, &
431                    Nx, Ny, Nz)
432                 dudy_at_E = (U_S(IJPK,M) - U_S(IJK,M)) * ONEoDY_N_U(IJK)
433                 IF(NOC_VS) dudy_at_E = dudy_at_E - ((Ui-UW_s) * &
434                               ONEoDY_N_U(IJK)/DEL_H*(Sx*Nx+Sz*Nz))
435              ELSE
436                 dudy_at_E =  ZERO
437              ENDIF
438     
439              IF(U_NODE_AT_NW.AND.U_NODE_AT_SW) THEN
440                 Ui = HALF * (U_S(IMJPK,M) + U_S(IMJK,M))
441                 Xi = HALF * (X_U(IMJPK) + X_U(IMJK))
442                 Yi = HALF * (Y_U(IMJPK) + Y_U(IMJK))
443                 Zi = HALF * (Z_U(IMJPK) + Z_U(IMJK))
444                 Sx = X_U(IMJPK) - X_U(IMJK)
445                 Sy = Y_U(IMJPK) - Y_U(IMJK)
446                 Sz = Z_U(IMJPK) - Z_U(IMJK)
447                 CALL GET_DEL_H(IJK, 'V_MOMENTUM', Xi, Yi, Zi, Del_H, &
448                    Nx, Ny, Nz)
449                 dudy_at_W =  (U_S(IMJPK,M)-U_S(IMJK,M))*ONEoDY_N_U(IMJK)
450                 IF(NOC_VS) dudy_at_W = dudy_at_W - ((Ui-UW_s) * &
451                               ONEoDY_N_U(IMJK)/DEL_H*(Sx*Nx+Sz*Nz))
452              ELSE
453                 dudy_at_W =  ZERO
454              ENDIF
455     
456              IF(U_NODE_AT_SE) THEN
457                 CALL GET_DEL_H(IJK, 'V_MOMENTUM', X_U(IJK), Y_U(IJK), &
458                    Z_U(IJK), Del_H, Nx, Ny, Nz)
459                 SSX_CUT = -MU_S_CUT * (U_S(IJK,M) - UW_s) / &
460                            DEL_H * (Ny*Nx) * Area_V_CUT(IJK)
461              ELSE
462                 SSX_CUT =  ZERO
463              ENDIF
464     
465              SSX = AVG_Y_H(AVG_X_H(EPMU_S(IJK,M),EPMU_S(IJKE,M),I),&
466                            AVG_X_H(EPMU_S(IJKN,M),EPMU_S(IJKNE,M),I),J)*&
467                       dudy_at_E*AYZ_V(IJK) - &
468                    AVG_Y_H(AVG_X_H(EPMU_S(IJKW,M),EPMU_S(IJK,M),IM),&
469                            AVG_X_H(EPMU_S(IJKNW,M),EPMU_S(IJKN,M),IM),J)*&
470                       dudy_at_W*AYZ_V(IMJK) + SSX_CUT
471     
472     ! SSY:
473              CALL GET_DEL_H(IJK,'V_MOMENTUM', X_V(IJK), Y_V(IJK),&
474                 Z_V(IJK), Del_H, Nx, Ny, Nz)
475              SSY = EPMU_S(IJKN,M)*(V_S(IJPK,M)-V_S(IJK,M))*&
476                       ONEoDY_N_V(IJK)*AXZ_V(IJK) - &
477                    EPMU_S(IJK,M)*(V_S(IJK,M)-V_S(IJMK,M))*&
478                       ONEoDY_N_V(IJMK)*AXZ_V(IJMK) -&
479                    MU_S_CUT * (V_S(IJK,M) - VW_s)/DEL_H * &
480                       (Ny**2) * Area_V_CUT(IJK)
481     
482     ! SSZ:
483              IF(DO_K) THEN
484                 W_NODE_AT_TN = ((.NOT.BLOCKED_W_CELL_AT(IJPK)).AND.&
485                                 (.NOT.WALL_W_AT(IJPK)))
486                 W_NODE_AT_TS = ((.NOT.BLOCKED_W_CELL_AT(IJK)).AND.&
487                                 (.NOT.WALL_W_AT(IJK)))
488                 W_NODE_AT_BN = ((.NOT.BLOCKED_W_CELL_AT(IJPKM)).AND.&
489                                 (.NOT.WALL_W_AT(IJPKM)))
490                 W_NODE_AT_BS = ((.NOT.BLOCKED_W_CELL_AT(IJKM)).AND.&
491                                 (.NOT.WALL_W_AT(IJKM)))
492     
493                 IF(W_NODE_AT_TN.AND.W_NODE_AT_TS) THEN
494                    Wi = HALF * (W_S(IJPK,M) + W_S(IJK,M))
495                    Xi = HALF * (X_W(IJPK) + X_W(IJK))
496                    Yi = HALF * (Y_W(IJPK) + Y_W(IJK))
497                    Zi = HALF * (Z_W(IJPK) + Z_W(IJK))
498                    Sx = X_W(IJPK) - X_W(IJK)
499                    Sy = Y_W(IJPK) - Y_W(IJK)
500                    Sz = Z_W(IJPK) - Z_W(IJK)
501                    CALL GET_DEL_H(IJK, 'V_MOMENTUM', Xi, Yi, Zi, Del_H, &
502                       Nx, Ny, Nz)
503                    dwdy_at_T = (W_S(IJPK,M) - W_S(IJK,M)) * ONEoDY_N_W(IJK)
504                    IF(NOC_VS) dwdy_at_T = dwdy_at_T - ((Wi-WW_s) * &
505                                  ONEoDY_N_W(IJK)/DEL_H*(Sx*Nx+Sz*Nz))
506                 ELSE
507                    dwdy_at_T =  ZERO
508                 ENDIF
509     
510                 IF(W_NODE_AT_BN.AND.W_NODE_AT_BS) THEN
511                    Wi = HALF * (W_S(IJPKM,M) + W_S(IJKM,M))
512                    Xi = HALF * (X_W(IJPKM) + X_W(IJKM))
513                    Yi = HALF * (Y_W(IJPKM) + Y_W(IJKM))
514                    Zi = HALF * (Z_W(IJPKM) + Z_W(IJKM))
515                    Sx = X_W(IJPKM) - X_W(IJKM)
516                    Sy = Y_W(IJPKM) - Y_W(IJKM)
517                    Sz = Z_W(IJPKM) - Z_W(IJKM)
518                    CALL GET_DEL_H(IJK, 'V_MOMENTUM', Xi, Yi, Zi, Del_H, &
519                       Nx, Ny, Nz)
520                    dwdy_at_B = (W_S(IJPKM,M)-W_S(IJKM,M)) * ONEoDY_N_W(IJKM)
521                    IF(NOC_VS) dwdy_at_B = dwdy_at_B - ((Wi-WW_s) * &
522                                  ONEoDY_N_W(IJKM)/DEL_H*(Sx*Nx+Sz*Nz))
523                 ELSE
524                    dwdy_at_B =  ZERO
525                 ENDIF
526     
527                 IF(W_NODE_AT_TS) THEN
528                    CALL GET_DEL_H(IJK, 'V_MOMENTUM', X_W(IJK), Y_W(IJK), &
529                       Z_W(IJK), Del_H, Nx, Ny, Nz)
530                    SSZ_CUT = -MU_S_CUT * (W_S(IJK,M) - WW_s) / &
531                               DEL_H * (Ny*Nz) * Area_V_CUT(IJK)
532                 ELSE
533                    SSZ_CUT =  ZERO
534                 ENDIF
535     
536                 SSZ = AVG_Y_H(AVG_Z_H(EPMU_S(IJK,M),EPMU_S(IJKT,M),K),&
537                               AVG_Z_H(EPMU_S(IJKN,M),EPMU_S(IJKTN,M),K),J)*&
538                          dwdy_at_T*AXY_V(IJK) - &
539                       AVG_Y_H(AVG_Z_H(EPMU_S(IJKB,M),EPMU_S(IJK,M),KM),&
540                               AVG_Z_H(EPMU_S(IJKBN,M),EPMU_S(IJKN,M),KM),J)*&
541                          dwdy_at_B*AXY_V(IJKM) + SSZ_CUT
542              ELSE
543                 SSZ = ZERO
544              ENDIF  ! end if do_k
545     
546           ENDIF  ! end if/else cut_cell
547     
548           RETURN
549           END SUBROUTINE CALC_CG_TAU_V_S
550     
551     
552     
553