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

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