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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Subroutine: CALC_Tau_U_s                                            C
4     !  Purpose: Cross terms in the gradient of stress in U_s momentum      C
5     !                                                                      C
6     !  Comments: This routine calculates the components of the solids      C
7     !  phase viscous stress tensor of the u-momentum equation that cannot  C
8     !  be cast in the form: mu.grad(u). These components are stored in the C
9     !  passed variable, which is then used as a source of the u-momentum   C
10     !  equation.                                                           C
11     !  For greater details see calc_tau_u_g.                               C
12     !                                                                      C
13     !                                                                      C
14     !  Author: M. Syamlal                                 Date: 19-DEC-96  C
15     !  Reviewer:                                          Date:            C
16     !                                                                      C
17     !                                                                      C
18     !                                                                      C
19     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
20           SUBROUTINE CALC_TAU_U_S(lTAU_U_S)
21     
22     ! Modules
23     !---------------------------------------------------------------------//
24           USE param, only: dimension_3, dimension_m
25           USE param1, only: zero
26           USE physprop, only: smax, mmax
27           USE fldvar, only: v_s
28           USE fldvar, only: ep_s
29           USE toleranc, only: dil_ep_s
30           USE run, only: kt_type_enum, ghd_2007
31           USE run, only: shear
32           USE vshear, only: VSH
33           USE cutcell, only: cartesian_grid, cg_safe_mode
34     
35           USE compar
36           USE geometry
37           USE indices
38           USE functions
39           USE sendrecv
40           USE fun_avg
41           IMPLICIT NONE
42     
43     ! Dummy arguments
44     !---------------------------------------------------------------------//
45     ! TAU_U_s
46           DOUBLE PRECISION, INTENT(OUT) :: lTAU_U_s(DIMENSION_3, DIMENSION_M)
47     
48     ! Local variables
49     !---------------------------------------------------------------------//
50     ! Indices
51           INTEGER :: IJK, I, IJKE
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     ! Cylindrical source terms
59           DOUBLE PRECISION :: Vtzb
60     !---------------------------------------------------------------------//
61     
62     
63           DO M = 1, MMAX
64              IF(KT_TYPE_ENUM == GHD_2007 .AND. M /= MMAX) CYCLE
65     
66              IF (SHEAR) THEN
67     !!$omp  parallel do private(IJK)
68                 DO IJK = IJKSTART3, IJKEND3
69                    IF (FLUID_AT(IJK)) THEN
70                       V_S(IJK,m)=V_S(IJK,m)+VSH(IJK)
71                    ENDIF
72                 ENDDO
73              ENDIF
74     
75     !!$omp  parallel do &
76     !!$omp  private(I, IJK, IJKE,                                         &
77     !!$omp&         EPSA, EPStmp, SSX, SSY, SSZ, SBV, VTZB)               &
78     !!$omp&  schedule(static)
79              DO IJK = IJKSTART3, IJKEND3
80     
81     ! Skip walls where some values are undefined.
82                 IF(WALL_AT(IJK)) cycle
83                 I = I_OF(IJK)
84                 IJKE = EAST_OF(IJK)
85     
86                 IF (KT_TYPE_ENUM == GHD_2007) THEN
87                    EPStmp = ZERO
88                    DO L = 1, SMAX
89                      EPStmp = EPStmp + AVG_X(EP_S(IJK,L),EP_S(IJKE,L),I)
90                    ENDDO
91                    EPSA = EPStmp
92                 ELSE
93                    EPSA = AVG_X(EP_S(IJK,M),EP_S(IJKE,M),I)
94                 ENDIF
95     
96                 IF ( .NOT.SIP_AT_E(IJK) .AND. EPSA>DIL_EP_S) THEN
97     
98                    IF ((.NOT.CARTESIAN_GRID) .OR. CG_SAFE_MODE(3) ==1) THEN
99                       CALL CALC_REG_TAU_U_S(IJK, M, SSX, SSY, SSZ, SBV, VTZB)
100                    ELSE
101                       VTZB = ZERO   ! no cylindrical terms in CG
102                       CALL CALC_CG_TAU_U_S(IJK, M, SSX, SSY, SSZ, SBV)
103                    ENDIF
104     
105     ! Add the terms
106                    lTAU_U_S(IJK,M) = SBV + SSX + SSY + SSZ + VTZB*VOL_U(IJK)
107                 ELSE
108                    lTAU_U_S(IJK,M) = ZERO
109                 ENDIF   ! end if/else .not.sip_at_e .and. epsa>dil_ep_s
110              ENDDO   ! end do ijk
111     !!$omp end parallel do
112     
113              IF (SHEAR) THEN
114     !!$omp  parallel do private(IJK)
115                 DO IJK = IJKSTART3, IJKEND3
116                    IF (FLUID_AT(IJK)) THEN
117                       V_S(IJK,m)=V_S(IJK,m)-VSH(IJK)
118                    ENDIF
119                 ENDDO
120              ENDIF
121     
122           ENDDO   ! end do m=1,mmax
123     
124           call send_recv(ltau_u_s,2)
125           RETURN
126           END SUBROUTINE CALC_TAU_U_S
127     
128     
129     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
130     !                                                                      C
131     !  Subroutine: calc_reg_tau_u_s                                        C
132     !  Purpose: Cross terms in the gradient of stress in U_s momentum      C
133     !  based on NON cartesian grid cut cell.                               C
134     !                                                                      C
135     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
136           SUBROUTINE CALC_REG_TAU_U_S(IJK, M, SSX, SSY, SSZ, SBV, VTZB)
137     
138     ! Modules
139     !---------------------------------------------------------------------//
140           USE param1, only: zero, half
141           USE fldvar, only: u_s, v_s, w_s
142           USE visc_s, only: epmu_s, eplambda_s
143           USE visc_s, only: trd_s
144     
145           USE compar
146           USE geometry
147           USE indices
148           USE functions
149           USE fun_avg
150     
151           IMPLICIT NONE
152     
153     ! Dummy arguments
154     !---------------------------------------------------------------------//
155     ! current ijk index
156           INTEGER, INTENT(IN) :: IJK
157     ! solids phase index
158           INTEGER, INTENT(IN) :: M
159     ! Source terms (surface)
160           DOUBLE PRECISION, INTENT(OUT) :: SSX
161           DOUBLE PRECISION, INTENT(OUT) :: SSY
162           DOUBLE PRECISION, INTENT(OUT) :: SSZ
163           DOUBLE PRECISION, INTENT(OUT) :: SBV
164     ! Cylindrical source terms
165           DOUBLE PRECISION, INTENT(OUT) :: VTZB
166     
167     ! Local variables
168     !---------------------------------------------------------------------//
169     ! Indices
170           INTEGER :: I, J, K, JM, KM, IP
171           INTEGER :: IJKE, IJKN, IJKS, IJKT, IJKB
172           INTEGER :: IJKNE, IJKSE, IJKTE, IJKBE
173           INTEGER :: IPJK, IMJK, IJKM, IJMK
174           INTEGER :: IPJMK, IPJKM
175     ! Average viscosity
176           DOUBLE PRECISION :: MU_ste, MU_sbe, MUSA
177     ! Average dW/Xdz
178           DOUBLE PRECISION :: dWoXdz
179     !---------------------------------------------------------------------//
180     
181           I = I_OF(IJK)
182           IP = IP1(I)
183           J = J_OF(IJK)
184           JM = JM1(J)
185           K = K_OF(IJK)
186           KM = KM1(K)
187     
188           IPJK = IP_OF(IJK)
189           IMJK = IM_OF(IJK)
190           IJMK = JM_OF(IJK)
191           IJKM = KM_OF(IJK)
192           IPJKM = IP_OF(IJKM)
193           IPJMK = JM_OF(IPJK)
194     
195           IJKE = EAST_OF(IJK)
196           IJKN = NORTH_OF(IJK)
197           IJKNE = EAST_OF(IJKN)
198           IJKS = SOUTH_OF(IJK)
199           IJKSE = EAST_OF(IJKS)
200           IJKT = TOP_OF(IJK)
201           IJKTE = EAST_OF(IJKT)
202           IJKB = BOTTOM_OF(IJK)
203           IJKBE = EAST_OF(IJKB)
204     
205     ! Surface forces at i+1/2, j, k
206     ! bulk viscosity term
207     ! combines part of 1/x d/dx (x.tau_xx) xdxdydz and -tau_zz/x xdxdydz =>
208     ! combines 1/x d/dx (x.lambda.trcD) xdxdydz - (lambda/x.trcD) xdxdydz =>
209     !              d/dx (lambda.trcD) xdxdydz
210     ! delta (lambda.trcD)Ap |E-W : at (i+1 - i-1), j, k
211           SBV = (EPLAMBDA_S(IJKE,M)*TRD_S(IJKE,M)-&
212                  EPLAMBDA_S(IJK,M)*TRD_S(IJK,M))*AYZ(IJK)
213     
214     ! shear stress terms at i+1/2, j, k
215     ! part of 1/x d/dx(x.tau_xx) xdxdydz =>
216     !         1/x d/dx (x.mu.du/dx) xdxdydz =>
217     ! delta (mu du/dx)Ayz |E-W : at (i+1 - i-1), j, k
218           SSX = EPMU_S(IJKE,M)*(U_S(IPJK,M)-U_S(IJK,M))*&
219                    ODX(IP)*AYZ_U(IJK) - &
220                 EPMU_S(IJK,M)*(U_S(IJK,M)-U_S(IMJK,M))*&
221                    ODX(I)*AYZ_U(IMJK)
222     
223     ! part of d/dy (tau_xy) xdxdydz =>
224     !         d/dy (mu.dv/dx) xdxdydz =>
225     ! delta (mu.dv/dx)Axz |N-S : at i+1/2, (j+1/2 - j-1/2), k
226           SSY = AVG_X_H(AVG_Y_H(EPMU_S(IJK,M),EPMU_S(IJKN,M),J),&
227                         AVG_Y_H(EPMU_S(IJKE,M),EPMU_S(IJKNE,M),J),I)*&
228                    (V_S(IPJK,M)-V_S(IJK,M))*ODX_E(I)*AXZ_U(IJK) - &
229                 AVG_X_H(AVG_Y_H(EPMU_S(IJKS,M),EPMU_S(IJK,M),JM),&
230                         AVG_Y_H(EPMU_S(IJKSE,M),EPMU_S(IJKE,M),JM),I)*&
231                    (V_S(IPJMK,M)-V_S(IJMK,M))*ODX_E(I)*AXZ_U(IJMK)
232     
233     ! part of 1/x d/dz (tau_xz) xdxdydz =>
234     !         1/x d/dz (mu.dw/dx) xdxdydz =>
235     ! delta (mu.dw/dx)Axy |T-B : at i+1/2, j, (k+1/2 - k-1/2)
236           MU_STE = AVG_X_H(AVG_Z_H(EPMU_S(IJK,M),EPMU_S(IJKT,M),K),&
237                            AVG_Z_H(EPMU_S(IJKE,M),EPMU_S(IJKTE,M),K),I)
238           MU_SBE = AVG_X_H(AVG_Z_H(EPMU_S(IJKB,M),EPMU_S(IJK,M),KM),&
239                            AVG_Z_H(EPMU_S(IJKBE,M),EPMU_S(IJKE,M),KM),I)
240           SSZ = MU_STE*(W_S(IPJK,M)-W_S(IJK,M))*ODX_E(I)*&
241                             AXY_U(IJK) - &
242                 MU_SBE*(W_S(IPJKM,M)-W_S(IJKM,M))*ODX_E(I)*&
243                             AXY_U(IJKM)
244     
245     
246     ! Special terms for cylindrical coordinates
247           IF (CYLINDRICAL) THEN
248     ! part of 1/x d/dz (tau_xz) xdxdydz =>
249     !         1/x d/dz (mu.(-w/x)) xdxdydz =>
250     ! delta (mu(-w/x))Axy |T-B : at i+1/2, j, (k+1/2- k-1/2)
251              SSZ = SSZ - (&
252                    MU_STE*(HALF*(W_S(IPJK,M)+W_S(IJK,M))*&
253                       OX_E(I))*AXY_U(IJK)-&
254                    MU_SBE*(HALF*(W_S(IPJKM,M)+W_S(IJKM,M))*&
255                       OX_E(I))*AXY_U(IJKM))
256     
257     ! part of -tau_zz/x xdxdydz =>
258     !         -(2.mu/x).(1/x).dw/dz xdxdydz
259     ! delta (2.mu/x.1/x.dw/dz)Vp |p : at i+1/2, j, k
260              MUSA = AVG_X(EPMU_S(IJK,M),EPMU_S(IJKE,M),I)
261              DWOXDZ = HALF*&
262                       ((W_S(IJK,M)-W_S(IJKM,M))*OX(I)*ODZ(K)+&
263                        (W_S(IPJK,M)-W_S(IPJKM,M))*OX(IP)*ODZ(K))
264              VTZB = -2.d0*MUSA*OX_E(I)*DWOXDZ
265           ELSE
266              VTZB = ZERO
267           ENDIF
268     
269           RETURN
270           END SUBROUTINE CALC_REG_TAU_U_S
271     
272     
273     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
274     !                                                                      C
275     !  Subroutine: calc_cg_tau_u_s                                         C
276     !  Purpose: Cross terms in the gradient of stress in U_s momentum      C
277     !  based on cartesian grid cut cell.                                   C
278     !                                                                      C
279     !  Author: Jeff Dietiker                              Date: 01-Jul-09  C
280     !                                                                      C
281     !                                                                      C
282     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
283           SUBROUTINE CALC_CG_TAU_U_S(IJK, M, SSX, SSY, SSZ, SBV)
284     
285     ! Modules
286     !---------------------------------------------------------------------//
287           USE param1, only: zero, half
288           USE fldvar, only: u_s, v_s, w_s
289           USE visc_s, only: epmu_s, eplambda_s
290           USE visc_s, only: trd_s
291     
292           USE compar
293           USE geometry
294           USE indices
295           USE functions
296           USE fun_avg
297     
298     ! for cutcell:
299     ! wall velocity for slip conditions
300           USE bc, only: bc_hw_s, bc_uw_s, bc_vw_s, bc_ww_s
301           USE bc, only: bc_type
302           USE quadric
303           USE cutcell
304     
305     ! Dummy arguments
306     !---------------------------------------------------------------------//
307     ! current ijk index
308           INTEGER, INTENT(IN) :: IJK
309     ! solids phase index
310           INTEGER, INTENT(IN) :: M
311     ! Source terms (surface)
312           DOUBLE PRECISION, INTENT(OUT) :: SSX
313           DOUBLE PRECISION, INTENT(OUT) :: SSY
314           DOUBLE PRECISION, INTENT(OUT) :: SSZ
315           DOUBLE PRECISION, INTENT(OUT) :: SBV
316     
317     ! Local variables
318     !---------------------------------------------------------------------//
319     ! Indices
320           INTEGER :: I, J, K, JM, KM, IP
321           INTEGER :: IJKE, IJKN, IJKS, IJKT, IJKB
322           INTEGER :: IJKNE, IJKSE, IJKTE, IJKBE
323           INTEGER :: IPJK, IMJK, IJKM, IJMK
324           INTEGER :: IPJMK, IPJKM
325     ! Average viscosity
326           DOUBLE PRECISION :: MU_ste, MU_sbe
327     
328     ! for cartesian grid implementation:
329           DOUBLE PRECISION :: DEL_H,Nx,Ny,Nz
330           LOGICAL :: V_NODE_AT_NE, V_NODE_AT_NW, V_NODE_AT_SE, V_NODE_AT_SW
331           LOGICAL :: W_NODE_AT_TE, W_NODE_AT_TW, W_NODE_AT_BE, W_NODE_AT_BW
332           DOUBLE PRECISION :: dvdx_at_N, dvdx_at_S
333           DOUBLE PRECISION :: dwdx_at_T, dwdx_at_B
334           DOUBLE PRECISION :: Xi, Yi, Zi, Vi, Wi, Sx, Sy, Sz
335           DOUBLE PRECISION :: MU_S_CUT, SSY_CUT, SSZ_CUT
336           DOUBLE PRECISION :: UW_s, VW_s, WW_s
337           INTEGER :: BCV
338           CHARACTER(LEN=9) :: BCT
339     !---------------------------------------------------------------------//
340     
341           I = I_OF(IJK)
342           IP = IP1(I)
343           J = J_OF(IJK)
344           JM = JM1(J)
345           K = K_OF(IJK)
346           KM = KM1(K)
347     
348           IPJK = IP_OF(IJK)
349           IMJK = IM_OF(IJK)
350           IJMK = JM_OF(IJK)
351           IJKM = KM_OF(IJK)
352           IPJKM = IP_OF(IJKM)
353           IPJMK = JM_OF(IPJK)
354     
355           IJKE = EAST_OF(IJK)
356           IJKN = NORTH_OF(IJK)
357           IJKNE = EAST_OF(IJKN)
358           IJKS = SOUTH_OF(IJK)
359           IJKSE = EAST_OF(IJKS)
360           IJKT = TOP_OF(IJK)
361           IJKTE = EAST_OF(IJKT)
362           IJKB = BOTTOM_OF(IJK)
363           IJKBE = EAST_OF(IJKB)
364     
365     
366     ! bulk viscosity term
367           SBV = (EPLAMBDA_S(IJKE,M)*TRD_S(IJKE,M)) * AYZ_U(IJK) - &
368                 (EPLAMBDA_S(IJK,M) *TRD_S(IJK,M) ) * AYZ_U(IMJK)
369     
370     ! shear stress terms
371           IF(.NOT.CUT_U_CELL_AT(IJK))   THEN
372              SSX = EPMU_S(IJKE,M)*(U_S(IPJK,M)-U_S(IJK,M))*&
373                       ONEoDX_E_U(IJK)*AYZ_U(IJK) - &
374                    EPMU_S(IJK,M)*(U_S(IJK,M)-U_S(IMJK,M))*&
375                       ONEoDX_E_U(IMJK)*AYZ_U(IMJK)
376              SSY = AVG_X_H(AVG_Y_H(EPMU_S(IJK,M),EPMU_S(IJKN,M),J),&
377                            AVG_Y_H(EPMU_S(IJKE,M),EPMU_S(IJKNE,M),J),I)*&
378                       (V_S(IPJK,M)-V_S(IJK,M))*&
379                       ONEoDX_E_V(IJK)*AXZ_U(IJK) - &
380                    AVG_X_H(AVG_Y_H(EPMU_S(IJKS,M),EPMU_S(IJK,M),JM),&
381                            AVG_Y_H(EPMU_S(IJKSE,M),EPMU_S(IJKE,M),JM),I)*&
382                       (V_S(IPJMK,M)-V_S(IJMK,M))*&
383                       ONEoDX_E_V(IJMK)*AXZ_U(IJMK)
384              IF(DO_K) THEN
385                 MU_STE = AVG_X_H(AVG_Z_H(EPMU_S(IJK,M),EPMU_S(IJKT,M),K),&
386                                  AVG_Z_H(EPMU_S(IJKE,M),EPMU_S(IJKTE,M),K),I)
387                 MU_SBE = AVG_X_H(AVG_Z_H(EPMU_S(IJKB,M),EPMU_S(IJK,M),KM),&
388                                  AVG_Z_H(EPMU_S(IJKBE,M),EPMU_S(IJKE,M),KM),I)
389                 SSZ = MU_STE*(W_S(IPJK,M)-W_S(IJK,M))*&
390                          ONEoDX_E_W(IJK)*AXY_U(IJK) - &
391                       MU_SBE*(W_S(IPJKM,M)-W_S(IJKM,M))*&
392                          ONEoDX_E_W(IJKM)*AXY_U(IJKM)
393              ELSE
394                 SSZ = ZERO
395              ENDIF
396     
397     ! cut cell modifications
398     !---------------------------------------------------------------------//
399           ELSE
400     
401              BCV = BC_U_ID(IJK)
402              IF(BCV > 0 ) THEN
403                 BCT = BC_TYPE(BCV)
404              ELSE
405                 BCT = 'NONE'
406              ENDIF
407     
408              SELECT CASE (BCT)
409                 CASE ('CG_NSW')
410                    CUT_TAU_US = .TRUE.
411                    NOC_US     = .TRUE.
412                    UW_s = ZERO
413                    VW_s = ZERO
414                    WW_s = ZERO
415                 CASE ('CG_FSW')
416                    CUT_TAU_US = .FALSE.
417                    NOC_US     = .FALSE.
418                    UW_s = ZERO
419                    VW_s = ZERO
420                    WW_s = ZERO
421                 CASE('CG_PSW')
422                    IF(BC_HW_S(BC_U_ID(IJK),M)==UNDEFINED) THEN   ! same as NSW
423                       CUT_TAU_US = .TRUE.
424                       NOC_US     = .TRUE.
425                       UW_s = BC_UW_S(BCV,M)
426                       VW_s = BC_VW_S(BCV,M)
427                       WW_s = BC_WW_S(BCV,M)
428                    ELSEIF(BC_HW_S(BC_U_ID(IJK),M)==ZERO) THEN   ! same as FSW
429                       CUT_TAU_US = .FALSE.
430                       NOC_US     = .FALSE.
431                       UW_s = ZERO
432                       VW_s = ZERO
433                       WW_s = ZERO
434                    ELSE                              ! partial slip
435                       CUT_TAU_US = .FALSE.
436                       NOC_US     = .FALSE.
437                    ENDIF
438                 CASE ('NONE')
439     ! equivalent of setting tau_u_s to zero at this ijk cell
440                    SSX = ZERO
441                    SSY = ZERO
442                    SSZ = ZERO
443                    SBV = ZERO
444                    RETURN
445              END SELECT
446     
447              IF(CUT_TAU_US) THEN
448                 MU_S_CUT = (VOL(IJK)*EPMU_S(IJK,M) + &
449                             VOL(IPJK)*EPMU_S(IJKE,M))/&
450                            (VOL(IJK) + VOL(IPJK))
451              ELSE
452                 MU_S_CUT = ZERO
453              ENDIF
454     
455     
456     ! SSX:
457              CALL GET_DEL_H(IJK, 'U_MOMENTUM', X_U(IJK), &
458                 Y_U(IJK), Z_U(IJK), Del_H, Nx, Ny, Nz)
459              SSX = EPMU_S(IJKE,M)*(U_S(IPJK,M)-U_S(IJK,M))*&
460                       ONEoDX_E_U(IJK)*AYZ_U(IJK) - &
461                    EPMU_S(IJK,M)*(U_S(IJK,M)-U_S(IMJK,M))*&
462                       ONEoDX_E_U(IMJK)*AYZ_U(IMJK) - &
463                    MU_S_CUT * (U_S(IJK,M) - UW_s) / DEL_H * &
464                       (Nx**2) * Area_U_CUT(IJK)
465     
466     ! SSY:
467              V_NODE_AT_NE = ((.NOT.BLOCKED_V_CELL_AT(IPJK)).AND.&
468                              (.NOT.WALL_V_AT(IPJK)))
469              V_NODE_AT_NW = ((.NOT.BLOCKED_V_CELL_AT(IJK)).AND.&
470                              (.NOT.WALL_V_AT(IJK)))
471              V_NODE_AT_SE = ((.NOT.BLOCKED_V_CELL_AT(IPJMK)).AND.&
472                              (.NOT.WALL_V_AT(IPJMK)))
473              V_NODE_AT_SW = ((.NOT.BLOCKED_V_CELL_AT(IJMK)).AND.&
474                              (.NOT.WALL_V_AT(IJMK)))
475     
476              IF(V_NODE_AT_NE.AND.V_NODE_AT_NW) THEN
477                 Vi = HALF * (V_S(IPJK,M) + V_S(IJK,M))
478                 Xi = HALF * (X_V(IPJK) + X_V(IJK))
479                 Yi = HALF * (Y_V(IPJK) + Y_V(IJK))
480                 Zi = HALF * (Z_V(IPJK) + Z_V(IJK))
481                 Sx = X_V(IPJK) - X_V(IJK)
482                 Sy = Y_V(IPJK) - Y_V(IJK)
483                 Sz = Z_V(IPJK) - Z_V(IJK)
484                 CALL GET_DEL_H(IJK, 'U_MOMENTUM', Xi, Yi, Zi, &
485                    Del_H, Nx, Ny, Nz)
486                 dvdx_at_N = (V_S(IPJK,M) - V_S(IJK,M)) * &
487                              ONEoDX_E_V(IJK)
488                 IF(NOC_US) dvdx_at_N = dvdx_at_N - ((Vi-VW_s) *&
489                               ONEoDX_E_V(IJK) /DEL_H*(Sy*Ny+Sz*Nz))
490              ELSE
491                 dvdx_at_N =  ZERO
492              ENDIF
493     
494              IF(V_NODE_AT_SE.AND.V_NODE_AT_SW) THEN
495                 Vi = HALF * (V_S(IPJMK,M) + V_S(IJMK,M))
496                 Xi = HALF * (X_V(IPJMK) + X_V(IJMK))
497                 Yi = HALF * (Y_V(IPJMK) + Y_V(IJMK))
498                 Zi = HALF * (Z_V(IPJMK) + Z_V(IJMK))
499                 Sx = X_V(IPJMK) - X_V(IJMK)
500                 Sy = Y_V(IPJMK) - Y_V(IJMK)
501                 Sz = Z_V(IPJMK) - Z_V(IJMK)
502                 CALL GET_DEL_H(IJK, 'U_MOMENTUM', Xi, Yi, Zi, &
503                    Del_H, Nx, Ny, Nz)
504                 dvdx_at_S = (V_S(IPJMK,M)-V_S(IJMK,M))*ONEoDX_E_V(IJMK)
505                 IF(NOC_US) dvdx_at_S = dvdx_at_S - ((Vi-VW_s) * &
506                               ONEoDX_E_V(IJMK)/DEL_H*(Sy*Ny+Sz*Nz))
507              ELSE
508                 dvdx_at_S =  ZERO
509              ENDIF
510     
511              IF(V_NODE_AT_NW) THEN
512                 CALL GET_DEL_H(IJK, 'U_MOMENTUM', X_V(IJK), &
513                    Y_V(IJK), Z_V(IJK), Del_H, Nx, Ny, Nz)
514                 SSY_CUT = -MU_S_CUT * (V_S(IJK,M) - VW_s)/DEL_H*&
515                            (Nx*Ny) * Area_U_CUT(IJK)
516              ELSE
517                 SSY_CUT =  ZERO
518              ENDIF
519     
520              SSY = AVG_X_H(AVG_Y_H(EPMU_S(IJK,M),EPMU_S(IJKN,M),J),&
521                            AVG_Y_H(EPMU_S(IJKE,M),EPMU_S(IJKNE,M),J),I)*&
522                       dvdx_at_N*AXZ_U(IJK) - &
523                    AVG_X_H(AVG_Y_H(EPMU_S(IJKS,M),EPMU_S(IJK,M),JM),&
524                            AVG_Y_H(EPMU_S(IJKSE,M),EPMU_S(IJKE,M),JM),I)*&
525                       dvdx_at_S*AXZ_U(IJMK) + SSY_CUT
526     
527     ! SSZ:
528              IF(DO_K) THEN
529                 W_NODE_AT_TE = ((.NOT.BLOCKED_W_CELL_AT(IPJK)).AND.&
530                                 (.NOT.WALL_W_AT(IPJK)))
531                 W_NODE_AT_TW = ((.NOT.BLOCKED_W_CELL_AT(IJK)).AND.&
532                                 (.NOT.WALL_W_AT(IJK)))
533                 W_NODE_AT_BE = ((.NOT.BLOCKED_W_CELL_AT(IPJKM)).AND.&
534                                 (.NOT.WALL_W_AT(IPJKM)))
535                 W_NODE_AT_BW = ((.NOT.BLOCKED_W_CELL_AT(IJKM)).AND.&
536                                 (.NOT.WALL_W_AT(IJKM)))
537     
538                 IF(W_NODE_AT_TE.AND.W_NODE_AT_TW) THEN
539                    Wi = HALF * (W_S(IPJK,M) + W_S(IJK,M))
540                    Xi = HALF * (X_W(IPJK) + X_W(IJK))
541                    Yi = HALF * (Y_W(IPJK) + Y_W(IJK))
542                    Zi = HALF * (Z_W(IPJK) + Z_W(IJK))
543                    Sx = X_W(IPJK) - X_W(IJK)
544                    Sy = Y_W(IPJK) - Y_W(IJK)
545                    Sz = Z_W(IPJK) - Z_W(IJK)
546                    CALL GET_DEL_H(IJK, 'U_MOMENTUM', Xi, Yi, Zi, &
547                       Del_H, Nx, Ny, Nz)
548                    dwdx_at_T = (W_S(IPJK,M) - W_S(IJK,M)) * &
549                                ONEoDX_E_W(IJK)
550                    IF(NOC_US) dwdx_at_T = dwdx_at_T - ((Wi-WW_s)*&
551                                  ONEoDX_E_W(IJK)/DEL_H*(Sy*Ny+Sz*Nz))
552                 ELSE
553                    dwdx_at_T =  ZERO
554                 ENDIF
555     
556                 IF(W_NODE_AT_BE.AND.W_NODE_AT_BW) THEN
557                    Wi = HALF * (W_S(IPJKM,M) + W_S(IJKM,M))
558                    Xi = HALF * (X_W(IPJKM) + X_W(IJKM))
559                    Yi = HALF * (Y_W(IPJKM) + Y_W(IJKM))
560                    Zi = HALF * (Z_W(IPJKM) + Z_W(IJKM))
561                    Sx = X_W(IPJKM) - X_W(IJKM)
562                    Sy = Y_W(IPJKM) - Y_W(IJKM)
563                    Sz = Z_W(IPJKM) - Z_W(IJKM)
564                    CALL GET_DEL_H(IJK, 'U_MOMENTUM', Xi, Yi, &
565                       Zi, Del_H, Nx, Ny, Nz)
566                    dwdx_at_B = (W_S(IPJKM,M) - W_S(IJKM,M)) * &
567                                ONEoDX_E_W(IJKM)
568                    IF(NOC_US) dwdx_at_B = dwdx_at_B - ((Wi-WW_s)*&
569                                  ONEoDX_E_W(IJKM)/DEL_H*(Sy*Ny+Sz*Nz))
570                 ELSE
571                    dwdx_at_B =  ZERO
572                 ENDIF
573     
574                 IF(W_NODE_AT_TW) THEN
575                    CALL GET_DEL_H(IJK, 'U_MOMENTUM', X_W(IJK), &
576                       Y_W(IJK), Z_W(IJK), Del_H, Nx, Ny, Nz)
577                    SSZ_CUT = -MU_S_CUT * (W_S(IJK,M) - WW_s) / &
578                              DEL_H * (Nx*Nz) * Area_U_CUT(IJK)
579                 ELSE
580                    SSZ_CUT =  ZERO
581                 ENDIF
582     
583                 MU_STE = AVG_X_H(AVG_Z_H(EPMU_S(IJK,M),EPMU_S(IJKT,M),K),&
584                                  AVG_Z_H(EPMU_S(IJKE,M),EPMU_S(IJKTE,M),K),I)
585                 MU_SBE = AVG_X_H(AVG_Z_H(EPMU_S(IJKB,M),EPMU_S(IJK,M),KM),&
586                                  AVG_Z_H(EPMU_S(IJKBE,M),EPMU_S(IJKE,M),KM),I)
587                 SSZ = MU_STE*dwdx_at_T*AXY_U(IJK) -  &
588                       MU_SBE*dwdx_at_B*AXY_U(IJKM) + SSZ_CUT
589              ELSE
590                 SSZ = ZERO
591              ENDIF  ! end if do_k
592     
593           ENDIF  ! end if/else cut_cell
594     
595           RETURN
596           END SUBROUTINE CALC_CG_TAU_U_S
597