File: /nfs/home/0/users/jenkins/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     !                                                                      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_W_S(TAU_W_S)
18     
19     !-----------------------------------------------
20     ! Modules
21     !-----------------------------------------------
22           USE param1, only: zero, half, one
23           USE fun_avg
24     
25     ! number of solids phases
26           USE physprop, only: smax, mmax
27     
28     ! x,y,z-components of solids velocity
29           USE fldvar, only: u_s, v_s, w_s
30     ! solids phase particle bulk density and material density
31           USE fldvar, only: rop_s, ro_s, ep_s
32     
33     ! viscous solids transport coefficients
34           USE visc_s, only: mu_s, lambda_s
35     ! trace of rate of strain tensor
36           USE visc_s, only: trd_s
37     
38     ! dilute threshold
39           USE toleranc, only: dil_ep_s
40     
41     ! kinetic theories
42           USE run, only: kt_type_enum, ghd_2007
43     ! runtime flag for treating the system as if shearing
44           USE run, only: shear
45     
46     ! primarily needed for function.inc
47           USE compar
48           USE geometry
49           USE indices
50     
51     ! for sendrecv calls
52           USE sendrecv
53     
54     ! for cutcell:
55     ! wall velocity for slip conditions
56           USE bc, only: bc_hw_s, bc_uw_s, bc_vw_s, bc_ww_s
57           USE bc, only: bc_type
58           USE quadric
59           USE cutcell
60     
61           USE functions
62     
63           IMPLICIT NONE
64     !-----------------------------------------------
65     ! Dummy arguments
66     !-----------------------------------------------
67     ! TAU_W_s
68           DOUBLE PRECISION, INTENT(OUT) :: TAU_W_s(DIMENSION_3, DIMENSION_M)
69     !-----------------------------------------------
70     ! Local Variables
71     !-----------------------------------------------
72     ! Indices
73           INTEGER :: IJK, J, I, IM, IJKP, IMJK, IJKN, IJKNT, IJKS,&
74                      IJKST, IJMKP, IJMK, IJKE, IJKTE, IJKW, IJKTW,&
75                      IMJKP, K, IJKT, JM, KP, IJKM
76     ! Phase index
77           INTEGER :: M, L
78     ! Average volume fraction
79           DOUBLE PRECISION :: EPSA, EPStmp
80     ! Average velocity gradients
81           DOUBLE PRECISION :: duodz
82     ! Source terms (Surface)
83           DOUBLE PRECISION :: Sbv, Ssx, Ssy, Ssz
84     ! Source terms (Volumetric)
85           DOUBLE PRECISION :: Vxz
86     
87     ! for cartesian grid implementation:
88           DOUBLE PRECISION :: DEL_H,Nx,Ny,Nz
89           LOGICAL :: U_NODE_AT_ET,U_NODE_AT_EB,U_NODE_AT_WT,U_NODE_AT_WB
90           LOGICAL :: V_NODE_AT_NT,V_NODE_AT_NB,V_NODE_AT_ST,V_NODE_AT_SB
91     
92           DOUBLE PRECISION :: dudz_at_E,dudz_at_W
93           DOUBLE PRECISION :: dvdz_at_N,dvdz_at_S
94           DOUBLE PRECISION :: MU_S_CUT,SSX_CUT,SSY_CUT
95           DOUBLE PRECISION :: Xi,Yi,Zi,Ui,Vi,Sx,Sy,Sz
96           DOUBLE PRECISION :: UW_s,VW_s,WW_s
97           INTEGER :: BCV
98           CHARACTER(LEN=9) :: BCT
99     !-----------------------------------------------
100     
101           DO M = 1, MMAX
102             IF(KT_TYPE_ENUM == GHD_2007 .AND. M /= MMAX) CYCLE
103     
104     
105     !!$omp  parallel do private( IJK, I, IJKE, EPSA, EPStmp  J,  K,   &
106     !!$omp&  JM,IJKP,IJKNT,IJKS,IJKST,IJMKP,IJKTE,IJKTW,IMJKP,KP, &
107     !!$omp&  DUODZ,VXZ, &
108     !!$omp&  IM,IJKW, IPJK,&
109     !!$omp&  IMJK,IJKN,IJMK,IJKT,  &
110     !!$omp&  IJKM, &
111     !!$omp&  SBV,  SSX,SSY,   SSZ) &
112     !!$omp&  schedule(static)
113     
114             DO IJK = IJKSTART3, IJKEND3
115     
116     ! Skip walls where some values are undefined.
117               IF(WALL_AT(IJK)) cycle
118     
119               K = K_OF(IJK)
120               IJKT = TOP_OF(IJK)
121     
122               IF (KT_TYPE_ENUM == GHD_2007) THEN
123                  EPStmp = ZERO
124                  DO L = 1, SMAX
125                     EPStmp = EPStmp + AVG_Z(EP_S(IJK,L),EP_S(IJKT,L),K)
126                  ENDDO
127                  EPSA = EPStmp
128               ELSE
129                  EPSA = AVG_Z(EP_S(IJK,M),EP_S(IJKT,M),K)
130               ENDIF
131     
132               IF ( .NOT.SIP_AT_T(IJK) .AND. EPSA>DIL_EP_S) THEN
133                 J = J_OF(IJK)
134                 I = I_OF(IJK)
135                 IM = IM1(I)
136                 JM = JM1(J)
137                 IJKP = KP_OF(IJK)
138                 IMJK = IM_OF(IJK)
139                 IJKN = NORTH_OF(IJK)
140                 IJKNT = TOP_OF(IJKN)
141                 IJKS = SOUTH_OF(IJK)
142                 IJKST = TOP_OF(IJKS)
143                 IJMKP = JM_OF(IJKP)
144                 IJMK = JM_OF(IJK)
145                 IJKE = EAST_OF(IJK)
146                 IJKTE = EAST_OF(IJKT)
147                 IJKW = WEST_OF(IJK)
148                 IJKTW = WEST_OF(IJKT)
149                 IMJKP = KP_OF(IMJK)
150                 KP = KP1(K)
151                 IJKM = KM_OF(IJK)
152     
153                 IF((.NOT.CARTESIAN_GRID).OR.(CG_SAFE_MODE(5)==1)) THEN
154     ! NON CARTESIAN GRID CASE
155     ! ---------------------------------------------------------------->>>
156     ! Surface forces
157     
158     ! bulk viscosity term
159                   SBV = (LAMBDA_S(IJKT,M)*TRD_S(IJKT,M)-&
160                          LAMBDA_S(IJK,M)*TRD_S(IJK,M))*AXY(IJK)
161     
162     ! shear stress terms
163                   SSX = AVG_Z_H(AVG_X_H(MU_S(IJK,M),MU_S(IJKE,M),I),&
164                                 AVG_X_H(MU_S(IJKT,M),MU_S(IJKTE,M),I),K)*&
165                         (U_S(IJKP,M)-U_S(IJK,M))*OX_E(I)*ODZ_T(K)*AYZ_W(IJK) - &
166                         AVG_Z_H(AVG_X_H(MU_S(IJKW,M),MU_S(IJK,M),IM),&
167                                 AVG_X_H(MU_S(IJKTW,M),MU_S(IJKT,M),IM),K)*&
168                         (U_S(IMJKP,M)-U_S(IMJK,M))*ODZ_T(K)*DY(J)*(HALF*(DZ(K)+DZ(KP)))
169                         !same as oX_E(IM)*AYZ_W(IMJK), but avoids singularity
170                   SSY = AVG_Z_H(AVG_Y_H(MU_S(IJK,M),MU_S(IJKN,M),J),&
171                                 AVG_Y_H(MU_S(IJKT,M),MU_S(IJKNT,M),J),K)*&
172                         (V_S(IJKP,M)-V_S(IJK,M))*OX(I)*ODZ_T(K)*AXZ_W(IJK) - &
173                         AVG_Z_H(AVG_Y_H(MU_S(IJKS,M),MU_S(IJK,M),JM),&
174                                 AVG_Y_H(MU_S(IJKST,M),MU_S(IJKT,M),JM),K)*&
175                         (V_S(IJMKP,M)-V_S(IJMK,M))*OX(I)*ODZ_T(K)*AXZ_W(IJMK)
176                   SSZ = MU_S(IJKT,M)*(W_S(IJKP,M)-W_S(IJK,M))*&
177                         OX(I)*ODZ(KP)*AXY_W(IJK) - &
178                         MU_S(IJK,M)*(W_S(IJK,M)-W_S(IJKM,M))*&
179                         OX(I)*ODZ(K)* AXY_W(IJKM)
180     ! ----------------------------------------------------------------<<<
181     
182                 ELSE
183     
184     ! CARTESIAN GRID CASE
185     ! ---------------------------------------------------------------->>>
186     ! Surface forces
187     
188     ! bulk viscosity term
189                   SBV =  (LAMBDA_S(IJKT,M)*TRD_S(IJKT,M)) * AXY_W(IJK) - &
190                          (LAMBDA_S(IJK,M) *TRD_S(IJK,M) ) * AXY_W(IJKM)
191     
192     ! shear stress terms
193                   IF(.NOT.CUT_W_CELL_AT(IJK)) THEN
194                     SSX = AVG_Z_H(AVG_X_H(MU_S(IJK,M),MU_S(IJKE,M),I),&
195                                   AVG_X_H(MU_S(IJKT,M),MU_S(IJKTE,M),I),K)*&
196                           (U_S(IJKP,M)-U_S(IJK,M))*ONEoDZ_T_U(IJK)*AYZ_W(IJK) - &
197                           AVG_Z_H(AVG_X_H(MU_S(IJKW,M),MU_S(IJK,M),IM),&
198                                   AVG_X_H(MU_S(IJKTW,M),MU_S(IJKT,M),IM),K)*&
199                           (U_S(IMJKP,M)-U_S(IMJK,M))*ONEoDZ_T_U(IMJK)*AYZ_W(IMJK)
200                     SSY = AVG_Z_H(AVG_Y_H(MU_S(IJK,M),MU_S(IJKN,M),J),&
201                                   AVG_Y_H(MU_S(IJKT,M),MU_S(IJKNT,M),J),K)*&
202                           (V_S(IJKP,M)-V_S(IJK,M))*ONEoDZ_T_V(IJK)*AXZ_W(IJK) - &
203                           AVG_Z_H(AVG_Y_H(MU_S(IJKS,M),MU_S(IJK,M),JM),&
204                                   AVG_Y_H(MU_S(IJKST,M),MU_S(IJKT,M),JM),K)*&
205                           (V_S(IJMKP,M)-V_S(IJMK,M))*ONEoDZ_T_V(IJMK)*AXZ_W(IJMK)
206                     SSZ = MU_S(IJKT,M)*(W_S(IJKP,M)-W_S(IJK,M))*&
207                           ONEoDZ_T_W(IJK)*AXY_W(IJK) - &
208                           MU_S(IJK,M)*(W_S(IJK,M)-W_S(IJKM,M))*&
209                           ONEoDZ_T_W(IJKM)*AXY_W(IJKM)
210     
211                   ELSE   ! CUT_CELL
212                     BCV = BC_W_ID(IJK)
213                     IF(BCV > 0 ) THEN
214                        BCT = BC_TYPE(BCV)
215                     ELSE
216                        BCT = 'NONE'
217                     ENDIF
218     
219                     SELECT CASE (BCT)
220                       CASE ('CG_NSW')
221                          CUT_TAU_WS = .TRUE.
222                          NOC_WS     = .TRUE.
223                          UW_s = ZERO
224                          VW_s = ZERO
225                          WW_s = ZERO
226                       CASE ('CG_FSW')
227                          CUT_TAU_WS = .FALSE.
228                          NOC_WS     = .FALSE.
229                          UW_s = ZERO
230                          VW_s = ZERO
231                          WW_s = ZERO
232                       CASE('CG_PSW')
233                          IF(BC_HW_S(BC_W_ID(IJK),M)==UNDEFINED) THEN   ! same as NSW
234                             CUT_TAU_WS = .TRUE.
235                             NOC_WS     = .TRUE.
236                             UW_s = BC_UW_S(BCV,M)
237                             VW_s = BC_VW_S(BCV,M)
238                             WW_s = BC_WW_S(BCV,M)
239                          ELSEIF(BC_HW_S(BC_W_ID(IJK),M)==ZERO) THEN   ! same as FSW
240                             CUT_TAU_WS = .FALSE.
241                             NOC_WS     = .FALSE.
242                             UW_s = ZERO
243                             VW_s = ZERO
244                             WW_s = ZERO
245                          ELSE                              ! partial slip
246                             CUT_TAU_WS = .FALSE.
247                             NOC_WS     = .FALSE.
248                          ENDIF
249                       CASE ('NONE')
250                          TAU_W_S(IJK,M) = ZERO
251                          CYCLE
252                     END SELECT
253     
254                     IF(CUT_TAU_WS) THEN
255                        MU_S_CUT = (VOL(IJK)*MU_S(IJK,M) + &
256                           VOL(IJKP)*MU_S(IJKT,M))/(VOL(IJK) + VOL(IJKP))
257                     ELSE
258                        MU_S_CUT = ZERO
259                     ENDIF
260     
261     ! SSX:
262                     U_NODE_AT_ET = ((.NOT.BLOCKED_U_CELL_AT(IJKP)).AND.(.NOT.WALL_U_AT(IJKP)))
263                     U_NODE_AT_EB = ((.NOT.BLOCKED_U_CELL_AT(IJK)).AND.(.NOT.WALL_U_AT(IJK)))
264                     U_NODE_AT_WT = ((.NOT.BLOCKED_U_CELL_AT(IMJKP)).AND.(.NOT.WALL_U_AT(IMJKP)))
265                     U_NODE_AT_WB = ((.NOT.BLOCKED_U_CELL_AT(IMJK)).AND.(.NOT.WALL_U_AT(IMJK)))
266     
267                     IF(U_NODE_AT_ET.AND.U_NODE_AT_EB) THEN
268                        Ui = HALF * (U_S(IJKP,M) + U_S(IJK,M))
269                        Xi = HALF * (X_U(IJKP) + X_U(IJK))
270                        Yi = HALF * (Y_U(IJKP) + Y_U(IJK))
271                        Zi = HALF * (Z_U(IJKP) + Z_U(IJK))
272                        Sx = X_U(IJKP) - X_U(IJK)
273                        Sy = Y_U(IJKP) - Y_U(IJK)
274                        Sz = Z_U(IJKP) - Z_U(IJK)
275                        CALL GET_DEL_H(IJK,'W_MOMENTUM',Xi,Yi,Zi,Del_H,Nx,Ny,Nz)
276                        dudz_at_E =  (U_S(IJKP,M) - U_S(IJK,M)) * ONEoDZ_T_U(IJK)
277                        IF(NOC_WS) dudz_at_E = dudz_at_E - ((Ui-UW_s) * &
278                                      ONEoDZ_T_U(IJK)/DEL_H*(Sx*Nx+Sy*Ny))
279                     ELSE
280                        dudz_at_E =  ZERO
281                     ENDIF
282     
283                     IF(U_NODE_AT_WT.AND.U_NODE_AT_WB) THEN
284                        Ui = HALF * (U_S(IMJKP,M) + U_S(IMJK,M))
285                        Xi = HALF * (X_U(IMJKP) + X_U(IMJK))
286                        Yi = HALF * (Y_U(IMJKP) + Y_U(IMJK))
287                        Zi = HALF * (Z_U(IMJKP) + Z_U(IMJK))
288                        Sx = X_U(IMJKP) - X_U(IMJK)
289                        Sy = Y_U(IMJKP) - Y_U(IMJK)
290                        Sz = Z_U(IMJKP) - Z_U(IMJK)
291                        CALL GET_DEL_H(IJK,'W_MOMENTUM',Xi,Yi,Zi,Del_H,Nx,Ny,Nz)
292                        dudz_at_W =  (U_S(IMJKP,M) - U_S(IMJK,M)) * ONEoDZ_T_U(IMJK)
293                        IF(NOC_WS) dudz_at_W = dudz_at_W - ((Ui-UW_s) * &
294                                      ONEoDZ_T_U(IMJK)/DEL_H*(Sx*Nx+Sy*Ny))
295                     ELSE
296                        dudz_at_W =  ZERO
297                     ENDIF
298     
299                     IF(U_NODE_AT_EB) THEN
300                        CALL GET_DEL_H(IJK,'W_MOMENTUM',X_U(IJK),Y_U(IJK),&
301                                       Z_U(IJK),Del_H,Nx,Ny,Nz)
302                        SSX_CUT = - MU_S_CUT * (U_S(IJK,M) - UW_s) / DEL_H * &
303                                  (Nz*Nx) * Area_W_CUT(IJK)
304                     ELSE
305                        SSX_CUT =  ZERO
306                     ENDIF
307     
308                     SSX = AVG_Z_H(AVG_X_H(MU_S(IJK,M),MU_S(IJKE,M),I),&
309                                   AVG_X_H(MU_S(IJKT,M),MU_S(IJKTE,M),I),K)*&
310                           dudz_at_E*AYZ_W(IJK) - &
311                           AVG_Z_H(AVG_X_H(MU_S(IJKW,M),MU_S(IJK,M),IM),&
312                                   AVG_X_H(MU_S(IJKTW,M),MU_S(IJKT,M),IM),K)*&
313                           dudz_at_W*AYZ_W(IMJK) + SSX_CUT
314     
315     ! SSY:
316                     V_NODE_AT_NT = ((.NOT.BLOCKED_V_CELL_AT(IJKP)).AND.(.NOT.WALL_V_AT(IJKP)))
317                     V_NODE_AT_NB = ((.NOT.BLOCKED_V_CELL_AT(IJK)).AND.(.NOT.WALL_V_AT(IJK)))
318                     V_NODE_AT_ST = ((.NOT.BLOCKED_V_CELL_AT(IJMKP)).AND.(.NOT.WALL_V_AT(IJMKP)))
319                     V_NODE_AT_SB = ((.NOT.BLOCKED_V_CELL_AT(IJMK)).AND.(.NOT.WALL_V_AT(IJMK)))
320     
321                     IF(V_NODE_AT_NT.AND.V_NODE_AT_NB) THEN
322                        Vi = HALF * (V_S(IJKP,M) + V_S(IJK,M))
323                        Xi = HALF * (X_V(IJKP) + X_V(IJK))
324                        Yi = HALF * (Y_V(IJKP) + Y_V(IJK))
325                        Zi = HALF * (Z_V(IJKP) + Z_V(IJK))
326                        Sx = X_V(IJKP) - X_V(IJK)
327                        Sy = Y_V(IJKP) - Y_V(IJK)
328                        Sz = Z_V(IJKP) - Z_V(IJK)
329                        CALL GET_DEL_H(IJK,'W_MOMENTUM',Xi,Yi,Zi,Del_H,Nx,Ny,Nz)
330                        dvdz_at_N =  (V_S(IJKP,M) - V_S(IJK,M)) * ONEoDZ_T_V(IJK)
331                        IF(NOC_WS) dvdz_at_N = dvdz_at_N - ((Vi-VW_s) * &
332                                      ONEoDZ_T_V(IJK)/DEL_H*(Sx*Nx+Sy*Ny))
333                     ELSE
334                        dvdz_at_N =  ZERO
335                     ENDIF
336     
337                     IF(V_NODE_AT_ST.AND.V_NODE_AT_SB) THEN
338                        Vi = HALF * (V_S(IJMKP,M) + V_S(IJMK,M))
339                        Xi = HALF * (X_V(IJMKP) + X_V(IJMK))
340                        Yi = HALF * (Y_V(IJMKP) + Y_V(IJMK))
341                        Zi = HALF * (Z_V(IJMKP) + Z_V(IJMK))
342                        Sx = X_V(IJMKP) - X_V(IJMK)
343                        Sy = Y_V(IJMKP) - Y_V(IJMK)
344                        Sz = Z_V(IJMKP) - Z_V(IJMK)
345                        CALL GET_DEL_H(IJK,'W_MOMENTUM',Xi,Yi,Zi,Del_H,Nx,Ny,Nz)
346                        dvdz_at_S =  (V_S(IJMKP,M) - V_S(IJMK,M)) * ONEoDZ_T_V(IJMK)
347                        IF(NOC_WS) dvdz_at_S = dvdz_at_S - ((Vi-VW_s) * &
348                                      ONEoDZ_T_V(IJMK)/DEL_H*(Sx*Nx+Sy*Ny))
349                     ELSE
350                        dvdz_at_S =  ZERO
351                     ENDIF
352     
353                     IF(V_NODE_AT_NB) THEN
354                        CALL GET_DEL_H(IJK,'W_MOMENTUM',X_V(IJK),Y_V(IJK),&
355                                       Z_V(IJK),Del_H,Nx,Ny,Nz)
356                        SSY_CUT = - MU_S_CUT * (V_S(IJK,M) - VW_s) / DEL_H * &
357                                  (Nz*Ny) * Area_W_CUT(IJK)
358                     ELSE
359                        SSY_CUT =  ZERO
360                     ENDIF
361     
362                     SSY = AVG_Z_H(AVG_Y_H(MU_S(IJK,M),MU_S(IJKN,M),J),&
363                                   AVG_Y_H(MU_S(IJKT,M),MU_S(IJKNT,M),J),K)*&
364                           dvdz_at_N*AXZ_W(IJK) - &
365                           AVG_Z_H(AVG_Y_H(MU_S(IJKS,M),MU_S(IJK,M),JM),&
366                                   AVG_Y_H(MU_S(IJKST,M),MU_S(IJKT,M),JM),K)*&
367                           dvdz_at_S*AXZ_W(IJMK) + SSY_CUT
368     
369     ! SSZ:
370                     CALL GET_DEL_H(IJK,'W_MOMENTUM',X_W(IJK),Y_W(IJK),&
371                                    Z_W(IJK),Del_H,Nx,Ny,Nz)
372                     SSZ = MU_S(IJKT,M)*(W_S(IJKP,M)-W_S(IJK,M))*&
373                           ONEoDZ_T_W(IJK)*AXY_W(IJK) - &
374                           MU_S(IJK,M)*(W_S(IJK,M)-W_S(IJKM,M))*&
375                           OX(I)*ONEoDZ_T_W(IJKM)*AXY_W(IJKM) &
376                           - MU_S_CUT * (W_S(IJK,M) - WW_s) / DEL_H * &
377                           (Nz**2) * Area_W_CUT(IJK)
378     
379                   ENDIF  ! end if/else cut_cell
380     ! ----------------------------------------------------------------<<<
381                 ENDIF   ! end if/else cartesian grid
382     
383     
384     ! Special terms for cylindrical coordinates
385                 IF (CYLINDRICAL) THEN
386     ! modify Szz to include integral of (1/x)*(d/dz)(2*u/x)
387                    SSZ = SSZ + &
388                       MU_S(IJKT,M)*(U_S(IJKP,M)+U_S(IMJKP,M))*OX(I)*AXY_W(IJK) - &
389                       MU_S(IJK,M)*(U_S(IJK,M)+U_S(IMJK,M))*OX(I)*AXY_W(IJKM)
390     
391     ! (mu_s/x)* part of tau_xz/X
392                    IF (OX_E(IM) /= UNDEFINED) THEN
393                       DUODZ = (U_S(IMJKP,M)-U_S(IMJK,M))*OX_E(IM)*ODZ_T(K)
394                    ELSE
395                       DUODZ = ZERO
396                    ENDIF
397                    VXZ = AVG_Z(MU_S(IJK,M),MU_S(IJKT,M),K)*OX(I)*HALF*&
398                       ((U_S(IJKP,M)-U_S(IJK,M))*OX_E(I)*ODZ_T(K)+DUODZ)
399                 ELSE
400                    VXZ = ZERO
401                 ENDIF
402     
403     ! Add the terms
404                   TAU_W_S(IJK,M) = SBV + SSX + SSY + SSZ + VXZ*VOL_W(IJK)
405               ELSE
406                 TAU_W_S(IJK,M) = ZERO
407               ENDIF   ! end if/else .not.sip_at_t .and. epsa>dil_ep_s
408     
409             ENDDO   ! end do ijk
410           ENDDO
411     
412           call send_recv(tau_w_s,2)
413           RETURN
414           END SUBROUTINE CALC_TAU_W_S
415