File: /nfs/home/0/users/jenkins/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     !                                                                      C
7     !  Author: M. Syamlal                                 Date: 19-DEC-96  C
8     !  Reviewer:                                          Date:            C
9     !                                                                      C
10     !  Revision Number: 1                                                  C
11     !  Purpose: To incorporate Cartesian grid modifications                C
12     !  Author: Jeff Dietiker                              Date: 01-Jul-09  C
13     !                                                                      C
14     !                                                                      C
15     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
16     
17           SUBROUTINE CALC_TAU_U_S(TAU_U_S)
18     
19     !-----------------------------------------------
20     ! Modules
21     !-----------------------------------------------
22           USE param1, only: zero, half, one
23     
24     ! number of solids phases
25           USE physprop, only: smax, mmax
26     
27     ! x,y,z-components of solids velocity
28           USE fldvar, only: u_s, v_s, w_s
29     ! solids phase particle bulk density and material density
30           USE fldvar, only: rop_s, ro_s, ep_s
31     
32     ! viscous solids transport coefficients
33           USE visc_s, only: mu_s, lambda_s
34     ! trace of rate of strain tensor
35           USE visc_s, only: trd_s
36     
37     ! dilute threshold
38           USE toleranc, only: dil_ep_s
39     
40     ! kinetic theories
41           USE run, only: kt_type_enum, ghd_2007
42     ! runtime flag for treating the system as if shearing
43           USE run, only: shear
44     
45     ! shear velocity
46           Use vshear, only: VSH
47     
48     ! primarily needed for function.inc
49           USE compar
50           USE geometry
51           USE indices
52     
53     ! for sendrecv calls
54           USE sendrecv
55           USE fun_avg
56     
57     ! for cutcell:
58     ! wall velocity for slip conditions
59           USE bc, only: bc_hw_s, bc_uw_s, bc_vw_s, bc_ww_s
60           USE bc, only: bc_type
61           USE quadric
62           USE cutcell
63     
64           USE functions
65     
66           IMPLICIT NONE
67     !-----------------------------------------------
68     ! Dummy arguments
69     !-----------------------------------------------
70     ! TAU_U_s
71           DOUBLE PRECISION, INTENT(OUT) :: TAU_U_s(DIMENSION_3, DIMENSION_M)
72     !-----------------------------------------------
73     ! Local variables
74     !-----------------------------------------------
75     ! Indices
76           INTEGER :: I, J, JM, K, KM, IJK, IJKE, IPJK, IP, IMJK, IJKN,&
77                      IJKNE, IJKS, IJKSE, IPJMK, IJMK, IJKT, IJKTE,&
78                      IJKB, IJKBE, IJKM, IPJKM
79     ! Phase index
80           INTEGER :: M, L
81     ! Average volume fraction
82           DOUBLE PRECISION :: EPSA, EPStmp
83     ! Average EP_s*viscosity
84           DOUBLE PRECISION :: EPMU_ste, EPMU_sbe, EPMUSA
85     ! Average dW/Xdz
86           DOUBLE PRECISION :: dWoXdz
87     ! Source terms (Surface)
88           DOUBLE PRECISION :: Sbv, Ssx, Ssy, Ssz
89     ! Source terms (Volumetric)
90           DOUBLE PRECISION :: Vtzb
91     
92     ! for cartesian grid implementation:
93           DOUBLE PRECISION :: DEL_H,Nx,Ny,Nz
94           LOGICAL :: V_NODE_AT_NE,V_NODE_AT_NW,V_NODE_AT_SE,V_NODE_AT_SW
95           LOGICAL :: W_NODE_AT_TE,W_NODE_AT_TW,W_NODE_AT_BE,W_NODE_AT_BW
96           DOUBLE PRECISION :: dvdx_at_N,dvdx_at_S
97           DOUBLE PRECISION :: dwdx_at_T,dwdx_at_B
98           DOUBLE PRECISION :: Xi,Yi,Zi,Vi,Wi,Sx,Sy,Sz
99           DOUBLE PRECISION :: MU_S_CUT,SSY_CUT,SSZ_CUT
100           DOUBLE PRECISION :: UW_s,VW_s,WW_s
101           INTEGER :: BCV
102           CHARACTER(LEN=9) :: BCT
103     
104     !-----------------------------------------------
105     
106           DO M = 1, MMAX
107             IF(KT_TYPE_ENUM == GHD_2007 .AND. M /= MMAX) CYCLE
108     
109     ! update to true velocity
110             IF (SHEAR) THEN
111     !!$omp  parallel do private(IJK)
112               DO IJK = IJKSTART3, IJKEND3
113                 IF (FLUID_AT(IJK)) THEN
114                   V_S(IJK,m)=V_S(IJK,m)+VSH(IJK)
115                 ENDIF
116               ENDDO
117             ENDIF
118     
119     !!$omp  parallel do private( IJK, I, IJKE, EPSA, EPStmp, IP, J, JM, K, KM,  &
120     !!$omp&  IPJK,IMJK,IJKN,IJKNE,IJKS,IJKSE,IPJMK,IJMK,IJKT,IJKTE,  &
121     !!$omp&  IJKB,IJKBE,IJKM,IPJKM, &
122     !!$omp&  SBV,  SSX,SSY,   SSZ, EPMU_STE,EPMU_SBE, &
123     !!$omp&  EPMUSA,DWOXDZ,VTZB ) &
124     !!$omp&  schedule(static)
125     
126             DO IJK = IJKSTART3, IJKEND3
127     
128     ! Skip walls where some values are undefined.
129               IF(WALL_AT(IJK)) cycle
130     
131               I = I_OF(IJK)
132               IJKE = EAST_OF(IJK)
133     
134               IF (KT_TYPE_ENUM == GHD_2007) THEN
135                 EPStmp = ZERO
136                 DO L = 1, SMAX
137                   EPStmp = EPStmp + AVG_X(EP_S(IJK,L),EP_S(IJKE,L),I)
138                 ENDDO
139                 EPSA = EPStmp
140               ELSE
141                 EPSA = AVG_X(EP_S(IJK,M),EP_S(IJKE,M),I)
142               ENDIF
143     
144               IF ( .NOT.SIP_AT_E(IJK) .AND. EPSA>DIL_EP_S) THEN
145                 IP = IP1(I)
146                 J = J_OF(IJK)
147                 JM = JM1(J)
148                 K = K_OF(IJK)
149                 KM = KM1(K)
150                 IPJK = IP_OF(IJK)
151                 IMJK = IM_OF(IJK)
152                 IJKN = NORTH_OF(IJK)
153                 IJKNE = EAST_OF(IJKN)
154                 IJKS = SOUTH_OF(IJK)
155                 IJKSE = EAST_OF(IJKS)
156                 IPJMK = JM_OF(IPJK)
157                 IJMK = JM_OF(IJK)
158                 IJKT = TOP_OF(IJK)
159                 IJKTE = EAST_OF(IJKT)
160                 IJKB = BOTTOM_OF(IJK)
161                 IJKBE = EAST_OF(IJKB)
162                 IJKM = KM_OF(IJK)
163                 IPJKM = IP_OF(IJKM)
164     
165                 IF((.NOT.CARTESIAN_GRID).OR.(CG_SAFE_MODE(3)==1)) THEN
166     ! NON CARTESIAN GRID CASE
167     ! ---------------------------------------------------------------->>>
168     ! Surface forces
169     
170     ! bulk viscosity term
171                   SBV = (LAMBDA_S(IJKE,M)*TRD_S(IJKE,M)-&
172                          LAMBDA_S(IJK,M)*TRD_S(IJK,M))*AYZ(IJK)
173     
174     ! shear stress terms
175                   SSX = MU_S(IJKE,M)*(U_S(IPJK,M)-U_S(IJK,M))*ODX(IP)*AYZ_U(IJK)&
176                       - MU_S(IJK,M)*(U_S(IJK,M)-U_S(IMJK,M))*ODX(I)*AYZ_U(IMJK)
177                   SSY = AVG_X_H(AVG_Y_H(MU_S(IJK,M),MU_S(IJKN,M),J),&
178                                 AVG_Y_H(MU_S(IJKE,M),MU_S(IJKNE,M),J),I)*&
179                         (V_S(IPJK,M)-V_S(IJK,M))*ODX_E(I)*AXZ_U(IJK) - &
180                         AVG_X_H(AVG_Y_H(MU_S(IJKS,M),MU_S(IJK,M),JM),&
181                                 AVG_Y_H(MU_S(IJKSE,M),MU_S(IJKE,M),JM),I)*&
182                         (V_S(IPJMK,M)-V_S(IJMK,M))*ODX_E(I)*AXZ_U(IJMK)
183     !?              IF(DO_K) THEN
184                     EPMU_STE = AVG_X_H(AVG_Z_H(MU_S(IJK,M),MU_S(IJKT,M),K),&
185                                        AVG_Z_H(MU_S(IJKE,M),MU_S(IJKTE,M),K),I)
186                     EPMU_SBE = AVG_X_H(AVG_Z_H(MU_S(IJKB,M),MU_S(IJK,M),KM),&
187                                        AVG_Z_H(MU_S(IJKBE,M),MU_S(IJKE,M),KM),I)
188                     SSZ = EPMU_STE*(W_S(IPJK,M)-W_S(IJK,M))*ODX_E(I)*AXY_U(IJK) - &
189                           EPMU_SBE*(W_S(IPJKM,M)-W_S(IJKM,M))*ODX_E(I)*AXY_U(IJKM)
190     !?              ELSE
191     !?                SSZ = ZERO
192     !?              ENDIF
193     ! ----------------------------------------------------------------<<<
194     
195                 ELSE
196     
197     ! CARTESIAN GRID CASE
198     ! ---------------------------------------------------------------->>>
199     ! Surface forces
200     
201     ! bulk viscosity term
202                   SBV = (LAMBDA_S(IJKE,M)*TRD_S(IJKE,M)) * AYZ_U(IJK) - &
203                         (LAMBDA_S(IJK,M) *TRD_S(IJK,M) ) * AYZ_U(IMJK)
204     
205     ! shear stress terms
206                   IF(.NOT.CUT_U_CELL_AT(IJK))   THEN
207                     SSX = MU_S(IJKE,M)*(U_S(IPJK,M)-U_S(IJK,M))*ONEoDX_E_U(IJK)*AYZ_U(IJK) - &
208                           MU_S(IJK,M)*(U_S(IJK,M)-U_S(IMJK,M))*ONEoDX_E_U(IMJK)*AYZ_U(IMJK)
209                     SSY = AVG_X_H(AVG_Y_H(MU_S(IJK,M),MU_S(IJKN,M),J),&
210                                   AVG_Y_H(MU_S(IJKE,M),MU_S(IJKNE,M),J),I)*&
211                           (V_S(IPJK,M)-V_S(IJK,M))*ONEoDX_E_V(IJK)*AXZ_U(IJK) - &
212                           AVG_X_H(AVG_Y_H(MU_S(IJKS,M),MU_S(IJK,M),JM),&
213                                   AVG_Y_H(MU_S(IJKSE,M),MU_S(IJKE,M),JM),I)*&
214                           (V_S(IPJMK,M)-V_S(IJMK,M))*ONEoDX_E_V(IJMK)*AXZ_U(IJMK)
215                     IF(DO_K) THEN
216                       EPMU_STE = AVG_X_H(AVG_Z_H(MU_S(IJK,M),MU_S(IJKT,M),K),&
217                                          AVG_Z_H(MU_S(IJKE,M),MU_S(IJKTE,M),K),I)
218                       EPMU_SBE = AVG_X_H(AVG_Z_H(MU_S(IJKB,M),MU_S(IJK,M),KM),&
219                                          AVG_Z_H(MU_S(IJKBE,M),MU_S(IJKE,M),KM),I)
220                       SSZ = EPMU_STE*(W_S(IPJK,M)-W_S(IJK,M))*ONEoDX_E_W(IJK)*AXY_U(IJK) - &
221                             EPMU_SBE*(W_S(IPJKM,M)-W_S(IJKM,M))*ONEoDX_E_W(IJKM)*AXY_U(IJKM)
222                     ELSE
223                       SSZ = ZERO
224                     ENDIF
225     
226     
227                   ELSE   ! CUT_CELL
228                     BCV = BC_U_ID(IJK)
229                     IF(BCV > 0 ) THEN
230                       BCT = BC_TYPE(BCV)
231                     ELSE
232                       BCT = 'NONE'
233                     ENDIF
234     
235                     SELECT CASE (BCT)
236                       CASE ('CG_NSW')
237                          CUT_TAU_US = .TRUE.
238                          NOC_US     = .TRUE.
239                          UW_s = ZERO
240                          VW_s = ZERO
241                          WW_s = ZERO
242                       CASE ('CG_FSW')
243                          CUT_TAU_US = .FALSE.
244                          NOC_US     = .FALSE.
245                          UW_s = ZERO
246                          VW_s = ZERO
247                          WW_s = ZERO
248                       CASE('CG_PSW')
249                          IF(BC_HW_S(BC_U_ID(IJK),M)==UNDEFINED) THEN   ! same as NSW
250                             CUT_TAU_US = .TRUE.
251                             NOC_US     = .TRUE.
252                             UW_s = BC_UW_S(BCV,M)
253                             VW_s = BC_VW_S(BCV,M)
254                             WW_s = BC_WW_S(BCV,M)
255                          ELSEIF(BC_HW_S(BC_U_ID(IJK),M)==ZERO) THEN   ! same as FSW
256                             CUT_TAU_US = .FALSE.
257                             NOC_US     = .FALSE.
258                             UW_s = ZERO
259                             VW_s = ZERO
260                             WW_s = ZERO
261                          ELSE                              ! partial slip
262                             CUT_TAU_US = .FALSE.
263                             NOC_US     = .FALSE.
264                          ENDIF
265                       CASE ('NONE')
266                          TAU_U_S(IJK,M) = ZERO
267                          CYCLE
268                     END SELECT
269     
270                     IF(CUT_TAU_US) THEN
271                       MU_S_CUT =  (VOL(IJK)*MU_S(IJK,M) + &
272                          VOL(IPJK)*MU_S(IJKE,M))/(VOL(IJK) + VOL(IPJK))
273                     ELSE
274                       MU_S_CUT = ZERO
275                     ENDIF
276     
277     ! SSX:
278                     CALL GET_DEL_H(IJK,'U_MOMENTUM',X_U(IJK),Y_U(IJK),Z_U(IJK),Del_H,Nx,Ny,Nz)
279                     SSX = MU_S(IJKE,M)*(U_S(IPJK,M)-U_S(IJK,M))*ONEoDX_E_U(IJK)*AYZ_U(IJK) - &
280                           MU_S(IJK,M)*(U_S(IJK,M)-U_S(IMJK,M))*ONEoDX_E_U(IMJK)*AYZ_U(IMJK) - &
281                           MU_S_CUT * (U_S(IJK,M) - UW_s) / DEL_H * (Nx**2) * Area_U_CUT(IJK)
282     
283     ! SSY:
284                     V_NODE_AT_NE = ((.NOT.BLOCKED_V_CELL_AT(IPJK)).AND.(.NOT.WALL_V_AT(IPJK)))
285                     V_NODE_AT_NW = ((.NOT.BLOCKED_V_CELL_AT(IJK)).AND.(.NOT.WALL_V_AT(IJK)))
286                     V_NODE_AT_SE = ((.NOT.BLOCKED_V_CELL_AT(IPJMK)).AND.(.NOT.WALL_V_AT(IPJMK)))
287                     V_NODE_AT_SW = ((.NOT.BLOCKED_V_CELL_AT(IJMK)).AND.(.NOT.WALL_V_AT(IJMK)))
288     
289                     IF(V_NODE_AT_NE.AND.V_NODE_AT_NW) THEN
290                        Vi = HALF * (V_S(IPJK,M) + V_S(IJK,M))
291                        Xi = HALF * (X_V(IPJK) + X_V(IJK))
292                        Yi = HALF * (Y_V(IPJK) + Y_V(IJK))
293                        Zi = HALF * (Z_V(IPJK) + Z_V(IJK))
294                        Sx = X_V(IPJK) - X_V(IJK)
295                        Sy = Y_V(IPJK) - Y_V(IJK)
296                        Sz = Z_V(IPJK) - Z_V(IJK)
297                        CALL GET_DEL_H(IJK,'U_MOMENTUM',Xi,Yi,Zi,Del_H,Nx,Ny,Nz)
298                        dvdx_at_N =  (V_S(IPJK,M) - V_S(IJK,M)) * ONEoDX_E_V(IJK)
299                        IF(NOC_US) dvdx_at_N = dvdx_at_N - ((Vi-VW_s) *&
300                                      ONEoDX_E_V(IJK) /DEL_H*(Sy*Ny+Sz*Nz))
301                     ELSE
302                        dvdx_at_N =  ZERO
303                     ENDIF
304     
305                     IF(V_NODE_AT_SE.AND.V_NODE_AT_SW) THEN
306                        Vi = HALF * (V_S(IPJMK,M) + V_S(IJMK,M))
307                        Xi = HALF * (X_V(IPJMK) + X_V(IJMK))
308                        Yi = HALF * (Y_V(IPJMK) + Y_V(IJMK))
309                        Zi = HALF * (Z_V(IPJMK) + Z_V(IJMK))
310                        Sx = X_V(IPJMK) - X_V(IJMK)
311                        Sy = Y_V(IPJMK) - Y_V(IJMK)
312                        Sz = Z_V(IPJMK) - Z_V(IJMK)
313                        CALL GET_DEL_H(IJK,'U_MOMENTUM',Xi,Yi,Zi,Del_H,Nx,Ny,Nz)
314                        dvdx_at_S =  (V_S(IPJMK,M)-V_S(IJMK,M))*ONEoDX_E_V(IJMK)
315                        IF(NOC_US) dvdx_at_S = dvdx_at_S - ((Vi-VW_s) * &
316                                      ONEoDX_E_V(IJMK)/DEL_H*(Sy*Ny+Sz*Nz))
317                     ELSE
318                        dvdx_at_S =  ZERO
319                     ENDIF
320     
321                     IF(V_NODE_AT_NW) THEN
322                        CALL GET_DEL_H(IJK,'U_MOMENTUM',X_V(IJK),Y_V(IJK),&
323                                       Z_V(IJK),Del_H,Nx,Ny,Nz)
324                        SSY_CUT = - MU_S_CUT * (V_S(IJK,M) - VW_s) / DEL_H * &
325                                  (Nx*Ny) * Area_U_CUT(IJK)
326                     ELSE
327                        SSY_CUT =  ZERO
328                     ENDIF
329     
330                     SSY = AVG_X_H(AVG_Y_H(MU_S(IJK,M),MU_S(IJKN,M),J),&
331                                   AVG_Y_H(MU_S(IJKE,M),MU_S(IJKNE,M),J),I)*&
332                           dvdx_at_N*AXZ_U(IJK) - &
333                           AVG_X_H(AVG_Y_H(MU_S(IJKS,M),MU_S(IJK,M),JM),&
334                                   AVG_Y_H(MU_S(IJKSE,M),MU_S(IJKE,M),JM),I)*&
335                           dvdx_at_S*AXZ_U(IJMK) + SSY_CUT
336     
337     ! SSZ:
338                     IF(DO_K) THEN
339                        W_NODE_AT_TE = ((.NOT.BLOCKED_W_CELL_AT(IPJK)).AND.(.NOT.WALL_W_AT(IPJK)))
340                        W_NODE_AT_TW = ((.NOT.BLOCKED_W_CELL_AT(IJK)).AND.(.NOT.WALL_W_AT(IJK)))
341                        W_NODE_AT_BE = ((.NOT.BLOCKED_W_CELL_AT(IPJKM)).AND.(.NOT.WALL_W_AT(IPJKM)))
342                        W_NODE_AT_BW = ((.NOT.BLOCKED_W_CELL_AT(IJKM)).AND.(.NOT.WALL_W_AT(IJKM)))
343     
344                        IF(W_NODE_AT_TE.AND.W_NODE_AT_TW) THEN
345                           Wi = HALF * (W_S(IPJK,M) + W_S(IJK,M))
346                           Xi = HALF * (X_W(IPJK) + X_W(IJK))
347                           Yi = HALF * (Y_W(IPJK) + Y_W(IJK))
348                           Zi = HALF * (Z_W(IPJK) + Z_W(IJK))
349                           Sx = X_W(IPJK) - X_W(IJK)
350                           Sy = Y_W(IPJK) - Y_W(IJK)
351                           Sz = Z_W(IPJK) - Z_W(IJK)
352                           CALL GET_DEL_H(IJK,'U_MOMENTUM',Xi,Yi,Zi,Del_H,Nx,Ny,Nz)
353                           dwdx_at_T =  (W_S(IPJK,M) - W_S(IJK,M)) * ONEoDX_E_W(IJK)
354                           IF(NOC_US) dwdx_at_T = dwdx_at_T - ((Wi-WW_s) * &
355                                         ONEoDX_E_W(IJK)/DEL_H*(Sy*Ny+Sz*Nz))
356                        ELSE
357                           dwdx_at_T =  ZERO
358                        ENDIF
359     
360                        IF(W_NODE_AT_BE.AND.W_NODE_AT_BW) THEN
361                           Wi = HALF * (W_S(IPJKM,M) + W_S(IJKM,M))
362                           Xi = HALF * (X_W(IPJKM) + X_W(IJKM))
363                           Yi = HALF * (Y_W(IPJKM) + Y_W(IJKM))
364                           Zi = HALF * (Z_W(IPJKM) + Z_W(IJKM))
365                           Sx = X_W(IPJKM) - X_W(IJKM)
366                           Sy = Y_W(IPJKM) - Y_W(IJKM)
367                           Sz = Z_W(IPJKM) - Z_W(IJKM)
368                           CALL GET_DEL_H(IJK,'U_MOMENTUM',Xi,Yi,Zi,Del_H,Nx,Ny,Nz)
369                           dwdx_at_B =  (W_S(IPJKM,M) - W_S(IJKM,M)) * ONEoDX_E_W(IJKM)
370                           IF(NOC_US) dwdx_at_B = dwdx_at_B  - ((Wi-WW_s) * &
371                                         ONEoDX_E_W(IJKM)/DEL_H*(Sy*Ny+Sz*Nz))
372                        ELSE
373                           dwdx_at_B =  ZERO
374                        ENDIF
375     
376                        IF(W_NODE_AT_TW) THEN
377                           CALL GET_DEL_H(IJK,'U_MOMENTUM',X_W(IJK),Y_W(IJK),&
378                                          Z_W(IJK),Del_H,Nx,Ny,Nz)
379                           SSZ_CUT = - MU_S_CUT * (W_S(IJK,M) - WW_s) / &
380                                     DEL_H * (Nx*Nz) * Area_U_CUT(IJK)
381                        ELSE
382                           SSZ_CUT =  ZERO
383                        ENDIF
384     
385                        EPMU_STE = AVG_X_H(AVG_Z_H(MU_S(IJK,M),MU_S(IJKT,M),K),&
386                                           AVG_Z_H(MU_S(IJKE,M),MU_S(IJKTE,M),K),I)
387                        EPMU_SBE = AVG_X_H(AVG_Z_H(MU_S(IJKB,M),MU_S(IJK,M),KM),&
388                                           AVG_Z_H(MU_S(IJKBE,M),MU_S(IJKE,M),KM),I)
389                        SSZ =   EPMU_STE*dwdx_at_T*AXY_U(IJK)  &
390                              - EPMU_SBE*dwdx_at_B*AXY_U(IJKM) &
391                              + SSZ_CUT
392                     ELSE
393                        SSZ = ZERO
394                     ENDIF  ! end if do_k
395     
396                   ENDIF  ! end if/else cut_cell
397     ! ----------------------------------------------------------------<<<
398                 ENDIF   ! end if/else cartesian grid
399     
400     
401     ! Special terms for cylindrical coordinates
402                 IF (CYLINDRICAL) THEN
403     ! modify Ssz: integral of (1/x) (d/dz) (mu*(-w/x))
404                    SSZ = SSZ - &
405                       (EPMU_STE*(HALF*(W_S(IPJK,M)+W_S(IJK,M))*OX_E(I))*AXY_U(IJK)-&
406                        EPMU_SBE*(HALF*(W_S(IPJKM,M)+W_S(IJKM,M))*OX_E(I))*AXY_U(IJKM))
407     
408     ! -(2mu/x)*(1/x)*dw/dz part of Tau_zz/X
409                     EPMUSA = AVG_X(MU_S(IJK,M),MU_S(IJKE,M),I)
410                     DWOXDZ = HALF*((W_S(IJK,M)-W_S(IJKM,M))*OX(I)*ODZ(K)+(W_S(&
411                          IPJK,M)-W_S(IPJKM,M))*OX(IP)*ODZ(K))
412                     VTZB = -2.d0*EPMUSA*OX_E(I)*DWOXDZ
413                 ELSE
414                     VTZB = ZERO
415                 ENDIF
416     
417     ! Add the terms
418                 TAU_U_S(IJK,M) = SBV + SSX + SSY + SSZ + VTZB*VOL_U(IJK)
419               ELSE
420                 TAU_U_S(IJK,M) = ZERO
421               ENDIF   ! end if/else .not.sip_at_e .and. epsa>dil_ep_s
422             ENDDO   ! end do ijk
423     
424             IF (SHEAR) THEN
425     !!$omp  parallel do private(IJK)
426               DO IJK = IJKSTART3, IJKEND3
427                 IF (FLUID_AT(IJK)) THEN
428                   V_S(IJK,m)=V_S(IJK,m)-VSH(IJK)
429                 ENDIF
430               ENDDO
431             ENDIF
432     
433           ENDDO
434     
435           call send_recv(tau_u_s,2)
436           RETURN
437           END SUBROUTINE CALC_TAU_U_S
438