File: N:\mfix\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     
288           USE compar
289           USE fldvar, only: u_s, v_s, w_s
290           USE fun_avg
291           USE functions
292           USE geometry
293           USE indices
294           USE param1, only: zero, half, undefined
295           USE visc_s, only: epmu_s, eplambda_s
296           USE visc_s, only: trd_s
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, cg_nsw, cg_fsw, cg_psw, cg_mi, none, bc_type_enum
301           USE quadric
302           USE cutcell
303     
304     ! Dummy arguments
305     !---------------------------------------------------------------------//
306     ! current ijk index
307           INTEGER, INTENT(IN) :: IJK
308     ! solids phase index
309           INTEGER, INTENT(IN) :: M
310     ! Source terms (surface)
311           DOUBLE PRECISION, INTENT(OUT) :: SSX
312           DOUBLE PRECISION, INTENT(OUT) :: SSY
313           DOUBLE PRECISION, INTENT(OUT) :: SSZ
314           DOUBLE PRECISION, INTENT(OUT) :: SBV
315     
316     ! Local variables
317     !---------------------------------------------------------------------//
318     ! Indices
319           INTEGER :: I, J, K, JM, KM, IP
320           INTEGER :: IJKE, IJKN, IJKS, IJKT, IJKB
321           INTEGER :: IJKNE, IJKSE, IJKTE, IJKBE
322           INTEGER :: IPJK, IMJK, IJKM, IJMK
323           INTEGER :: IPJMK, IPJKM
324     ! Average viscosity
325           DOUBLE PRECISION :: MU_ste, MU_sbe
326     
327     ! for cartesian grid implementation:
328           DOUBLE PRECISION :: DEL_H,Nx,Ny,Nz
329           LOGICAL :: V_NODE_AT_NE, V_NODE_AT_NW, V_NODE_AT_SE, V_NODE_AT_SW
330           LOGICAL :: W_NODE_AT_TE, W_NODE_AT_TW, W_NODE_AT_BE, W_NODE_AT_BW
331           DOUBLE PRECISION :: dvdx_at_N, dvdx_at_S
332           DOUBLE PRECISION :: dwdx_at_T, dwdx_at_B
333           DOUBLE PRECISION :: Xi, Yi, Zi, Vi, Wi, Sx, Sy, Sz
334           DOUBLE PRECISION :: MU_S_CUT, SSY_CUT, SSZ_CUT
335           DOUBLE PRECISION :: UW_s, VW_s, WW_s
336           INTEGER :: BCV
337           INTEGER :: BCT
338     !---------------------------------------------------------------------//
339     
340           I = I_OF(IJK)
341           IP = IP1(I)
342           J = J_OF(IJK)
343           JM = JM1(J)
344           K = K_OF(IJK)
345           KM = KM1(K)
346     
347           IPJK = IP_OF(IJK)
348           IMJK = IM_OF(IJK)
349           IJMK = JM_OF(IJK)
350           IJKM = KM_OF(IJK)
351           IPJKM = IP_OF(IJKM)
352           IPJMK = JM_OF(IPJK)
353     
354           IJKE = EAST_OF(IJK)
355           IJKN = NORTH_OF(IJK)
356           IJKNE = EAST_OF(IJKN)
357           IJKS = SOUTH_OF(IJK)
358           IJKSE = EAST_OF(IJKS)
359           IJKT = TOP_OF(IJK)
360           IJKTE = EAST_OF(IJKT)
361           IJKB = BOTTOM_OF(IJK)
362           IJKBE = EAST_OF(IJKB)
363     
364     
365     ! bulk viscosity term
366           SBV = (EPLAMBDA_S(IJKE,M)*TRD_S(IJKE,M)) * AYZ_U(IJK) - &
367                 (EPLAMBDA_S(IJK,M) *TRD_S(IJK,M) ) * AYZ_U(IMJK)
368     
369     ! shear stress terms
370           IF(.NOT.CUT_U_CELL_AT(IJK))   THEN
371              SSX = EPMU_S(IJKE,M)*(U_S(IPJK,M)-U_S(IJK,M))*&
372                       ONEoDX_E_U(IJK)*AYZ_U(IJK) - &
373                    EPMU_S(IJK,M)*(U_S(IJK,M)-U_S(IMJK,M))*&
374                       ONEoDX_E_U(IMJK)*AYZ_U(IMJK)
375              SSY = AVG_X_H(AVG_Y_H(EPMU_S(IJK,M),EPMU_S(IJKN,M),J),&
376                            AVG_Y_H(EPMU_S(IJKE,M),EPMU_S(IJKNE,M),J),I)*&
377                       (V_S(IPJK,M)-V_S(IJK,M))*&
378                       ONEoDX_E_V(IJK)*AXZ_U(IJK) - &
379                    AVG_X_H(AVG_Y_H(EPMU_S(IJKS,M),EPMU_S(IJK,M),JM),&
380                            AVG_Y_H(EPMU_S(IJKSE,M),EPMU_S(IJKE,M),JM),I)*&
381                       (V_S(IPJMK,M)-V_S(IJMK,M))*&
382                       ONEoDX_E_V(IJMK)*AXZ_U(IJMK)
383              IF(DO_K) THEN
384                 MU_STE = AVG_X_H(AVG_Z_H(EPMU_S(IJK,M),EPMU_S(IJKT,M),K),&
385                                  AVG_Z_H(EPMU_S(IJKE,M),EPMU_S(IJKTE,M),K),I)
386                 MU_SBE = AVG_X_H(AVG_Z_H(EPMU_S(IJKB,M),EPMU_S(IJK,M),KM),&
387                                  AVG_Z_H(EPMU_S(IJKBE,M),EPMU_S(IJKE,M),KM),I)
388                 SSZ = MU_STE*(W_S(IPJK,M)-W_S(IJK,M))*&
389                          ONEoDX_E_W(IJK)*AXY_U(IJK) - &
390                       MU_SBE*(W_S(IPJKM,M)-W_S(IJKM,M))*&
391                          ONEoDX_E_W(IJKM)*AXY_U(IJKM)
392              ELSE
393                 SSZ = ZERO
394              ENDIF
395     
396     ! cut cell modifications
397     !---------------------------------------------------------------------//
398           ELSE
399     
400              BCV = BC_U_ID(IJK)
401              IF(BCV > 0 ) THEN
402                 BCT = BC_TYPE_ENUM(BCV)
403              ELSE
404                 BCT = NONE
405              ENDIF
406     
407              SELECT CASE (BCT)
408                 CASE (CG_NSW,CG_MI)
409                    CUT_TAU_US = .TRUE.
410                    NOC_US     = .TRUE.
411                    UW_s = ZERO
412                    VW_s = ZERO
413                    WW_s = ZERO
414                 CASE (CG_FSW)
415                    CUT_TAU_US = .FALSE.
416                    NOC_US     = .FALSE.
417                    UW_s = ZERO
418                    VW_s = ZERO
419                    WW_s = ZERO
420                 CASE(CG_PSW)
421                    IF(BC_HW_S(BC_U_ID(IJK),M)==UNDEFINED) THEN   ! same as NSW
422                       CUT_TAU_US = .TRUE.
423                       NOC_US     = .TRUE.
424                       UW_s = BC_UW_S(BCV,M)
425                       VW_s = BC_VW_S(BCV,M)
426                       WW_s = BC_WW_S(BCV,M)
427                    ELSEIF(BC_HW_S(BC_U_ID(IJK),M)==ZERO) THEN   ! same as FSW
428                       CUT_TAU_US = .FALSE.
429                       NOC_US     = .FALSE.
430                       UW_s = ZERO
431                       VW_s = ZERO
432                       WW_s = ZERO
433                    ELSE                              ! partial slip
434                       CUT_TAU_US = .FALSE.
435                       NOC_US     = .FALSE.
436                    ENDIF
437                 CASE (NONE)
438     ! equivalent of setting tau_u_s to zero at this ijk cell
439                    SSX = ZERO
440                    SSY = ZERO
441                    SSZ = ZERO
442                    SBV = ZERO
443                    RETURN
444              END SELECT
445     
446              IF(CUT_TAU_US) THEN
447                 MU_S_CUT = (VOL(IJK)*EPMU_S(IJK,M) + &
448                             VOL(IPJK)*EPMU_S(IJKE,M))/&
449                            (VOL(IJK) + VOL(IPJK))
450              ELSE
451                 MU_S_CUT = ZERO
452              ENDIF
453     
454     
455     ! SSX:
456              CALL GET_DEL_H(IJK, 'U_MOMENTUM', X_U(IJK), &
457                 Y_U(IJK), Z_U(IJK), Del_H, Nx, Ny, Nz)
458              SSX = EPMU_S(IJKE,M)*(U_S(IPJK,M)-U_S(IJK,M))*&
459                       ONEoDX_E_U(IJK)*AYZ_U(IJK) - &
460                    EPMU_S(IJK,M)*(U_S(IJK,M)-U_S(IMJK,M))*&
461                       ONEoDX_E_U(IMJK)*AYZ_U(IMJK) - &
462                    MU_S_CUT * (U_S(IJK,M) - UW_s) / DEL_H * &
463                       (Nx**2) * Area_U_CUT(IJK)
464     
465     ! SSY:
466              V_NODE_AT_NE = ((.NOT.BLOCKED_V_CELL_AT(IPJK)).AND.&
467                              (.NOT.WALL_V_AT(IPJK)))
468              V_NODE_AT_NW = ((.NOT.BLOCKED_V_CELL_AT(IJK)).AND.&
469                              (.NOT.WALL_V_AT(IJK)))
470              V_NODE_AT_SE = ((.NOT.BLOCKED_V_CELL_AT(IPJMK)).AND.&
471                              (.NOT.WALL_V_AT(IPJMK)))
472              V_NODE_AT_SW = ((.NOT.BLOCKED_V_CELL_AT(IJMK)).AND.&
473                              (.NOT.WALL_V_AT(IJMK)))
474     
475              IF(V_NODE_AT_NE.AND.V_NODE_AT_NW) THEN
476                 Vi = HALF * (V_S(IPJK,M) + V_S(IJK,M))
477                 Xi = HALF * (X_V(IPJK) + X_V(IJK))
478                 Yi = HALF * (Y_V(IPJK) + Y_V(IJK))
479                 Zi = HALF * (Z_V(IPJK) + Z_V(IJK))
480                 Sx = X_V(IPJK) - X_V(IJK)
481                 Sy = Y_V(IPJK) - Y_V(IJK)
482                 Sz = Z_V(IPJK) - Z_V(IJK)
483                 CALL GET_DEL_H(IJK, 'U_MOMENTUM', Xi, Yi, Zi, &
484                    Del_H, Nx, Ny, Nz)
485                 dvdx_at_N = (V_S(IPJK,M) - V_S(IJK,M)) * &
486                              ONEoDX_E_V(IJK)
487                 IF(NOC_US) dvdx_at_N = dvdx_at_N - ((Vi-VW_s) *&
488                               ONEoDX_E_V(IJK) /DEL_H*(Sy*Ny+Sz*Nz))
489              ELSE
490                 dvdx_at_N =  ZERO
491              ENDIF
492     
493              IF(V_NODE_AT_SE.AND.V_NODE_AT_SW) THEN
494                 Vi = HALF * (V_S(IPJMK,M) + V_S(IJMK,M))
495                 Xi = HALF * (X_V(IPJMK) + X_V(IJMK))
496                 Yi = HALF * (Y_V(IPJMK) + Y_V(IJMK))
497                 Zi = HALF * (Z_V(IPJMK) + Z_V(IJMK))
498                 Sx = X_V(IPJMK) - X_V(IJMK)
499                 Sy = Y_V(IPJMK) - Y_V(IJMK)
500                 Sz = Z_V(IPJMK) - Z_V(IJMK)
501                 CALL GET_DEL_H(IJK, 'U_MOMENTUM', Xi, Yi, Zi, &
502                    Del_H, Nx, Ny, Nz)
503                 dvdx_at_S = (V_S(IPJMK,M)-V_S(IJMK,M))*ONEoDX_E_V(IJMK)
504                 IF(NOC_US) dvdx_at_S = dvdx_at_S - ((Vi-VW_s) * &
505                               ONEoDX_E_V(IJMK)/DEL_H*(Sy*Ny+Sz*Nz))
506              ELSE
507                 dvdx_at_S =  ZERO
508              ENDIF
509     
510              IF(V_NODE_AT_NW) THEN
511                 CALL GET_DEL_H(IJK, 'U_MOMENTUM', X_V(IJK), &
512                    Y_V(IJK), Z_V(IJK), Del_H, Nx, Ny, Nz)
513                 SSY_CUT = -MU_S_CUT * (V_S(IJK,M) - VW_s)/DEL_H*&
514                            (Nx*Ny) * Area_U_CUT(IJK)
515              ELSE
516                 SSY_CUT =  ZERO
517              ENDIF
518     
519              SSY = AVG_X_H(AVG_Y_H(EPMU_S(IJK,M),EPMU_S(IJKN,M),J),&
520                            AVG_Y_H(EPMU_S(IJKE,M),EPMU_S(IJKNE,M),J),I)*&
521                       dvdx_at_N*AXZ_U(IJK) - &
522                    AVG_X_H(AVG_Y_H(EPMU_S(IJKS,M),EPMU_S(IJK,M),JM),&
523                            AVG_Y_H(EPMU_S(IJKSE,M),EPMU_S(IJKE,M),JM),I)*&
524                       dvdx_at_S*AXZ_U(IJMK) + SSY_CUT
525     
526     ! SSZ:
527              IF(DO_K) THEN
528                 W_NODE_AT_TE = ((.NOT.BLOCKED_W_CELL_AT(IPJK)).AND.&
529                                 (.NOT.WALL_W_AT(IPJK)))
530                 W_NODE_AT_TW = ((.NOT.BLOCKED_W_CELL_AT(IJK)).AND.&
531                                 (.NOT.WALL_W_AT(IJK)))
532                 W_NODE_AT_BE = ((.NOT.BLOCKED_W_CELL_AT(IPJKM)).AND.&
533                                 (.NOT.WALL_W_AT(IPJKM)))
534                 W_NODE_AT_BW = ((.NOT.BLOCKED_W_CELL_AT(IJKM)).AND.&
535                                 (.NOT.WALL_W_AT(IJKM)))
536     
537                 IF(W_NODE_AT_TE.AND.W_NODE_AT_TW) THEN
538                    Wi = HALF * (W_S(IPJK,M) + W_S(IJK,M))
539                    Xi = HALF * (X_W(IPJK) + X_W(IJK))
540                    Yi = HALF * (Y_W(IPJK) + Y_W(IJK))
541                    Zi = HALF * (Z_W(IPJK) + Z_W(IJK))
542                    Sx = X_W(IPJK) - X_W(IJK)
543                    Sy = Y_W(IPJK) - Y_W(IJK)
544                    Sz = Z_W(IPJK) - Z_W(IJK)
545                    CALL GET_DEL_H(IJK, 'U_MOMENTUM', Xi, Yi, Zi, &
546                       Del_H, Nx, Ny, Nz)
547                    dwdx_at_T = (W_S(IPJK,M) - W_S(IJK,M)) * &
548                                ONEoDX_E_W(IJK)
549                    IF(NOC_US) dwdx_at_T = dwdx_at_T - ((Wi-WW_s)*&
550                                  ONEoDX_E_W(IJK)/DEL_H*(Sy*Ny+Sz*Nz))
551                 ELSE
552                    dwdx_at_T =  ZERO
553                 ENDIF
554     
555                 IF(W_NODE_AT_BE.AND.W_NODE_AT_BW) THEN
556                    Wi = HALF * (W_S(IPJKM,M) + W_S(IJKM,M))
557                    Xi = HALF * (X_W(IPJKM) + X_W(IJKM))
558                    Yi = HALF * (Y_W(IPJKM) + Y_W(IJKM))
559                    Zi = HALF * (Z_W(IPJKM) + Z_W(IJKM))
560                    Sx = X_W(IPJKM) - X_W(IJKM)
561                    Sy = Y_W(IPJKM) - Y_W(IJKM)
562                    Sz = Z_W(IPJKM) - Z_W(IJKM)
563                    CALL GET_DEL_H(IJK, 'U_MOMENTUM', Xi, Yi, &
564                       Zi, Del_H, Nx, Ny, Nz)
565                    dwdx_at_B = (W_S(IPJKM,M) - W_S(IJKM,M)) * &
566                                ONEoDX_E_W(IJKM)
567                    IF(NOC_US) dwdx_at_B = dwdx_at_B - ((Wi-WW_s)*&
568                                  ONEoDX_E_W(IJKM)/DEL_H*(Sy*Ny+Sz*Nz))
569                 ELSE
570                    dwdx_at_B =  ZERO
571                 ENDIF
572     
573                 IF(W_NODE_AT_TW) THEN
574                    CALL GET_DEL_H(IJK, 'U_MOMENTUM', X_W(IJK), &
575                       Y_W(IJK), Z_W(IJK), Del_H, Nx, Ny, Nz)
576                    SSZ_CUT = -MU_S_CUT * (W_S(IJK,M) - WW_s) / &
577                              DEL_H * (Nx*Nz) * Area_U_CUT(IJK)
578                 ELSE
579                    SSZ_CUT =  ZERO
580                 ENDIF
581     
582                 MU_STE = AVG_X_H(AVG_Z_H(EPMU_S(IJK,M),EPMU_S(IJKT,M),K),&
583                                  AVG_Z_H(EPMU_S(IJKE,M),EPMU_S(IJKTE,M),K),I)
584                 MU_SBE = AVG_X_H(AVG_Z_H(EPMU_S(IJKB,M),EPMU_S(IJK,M),KM),&
585                                  AVG_Z_H(EPMU_S(IJKBE,M),EPMU_S(IJKE,M),KM),I)
586                 SSZ = MU_STE*dwdx_at_T*AXY_U(IJK) -  &
587                       MU_SBE*dwdx_at_B*AXY_U(IJKM) + SSZ_CUT
588              ELSE
589                 SSZ = ZERO
590              ENDIF  ! end if do_k
591     
592           ENDIF  ! end if/else cut_cell
593     
594           RETURN
595           END SUBROUTINE CALC_CG_TAU_U_S
596