File: /nfs/home/0/users/jenkins/mfix.git/model/tau_w_g.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Module name: CALC_Tau_W_g(TAU_W_g, IER)                             C
4     !  Purpose: Cross terms in the gradient of stress in W_g momentum      c
5     !                                                                      C
6     !                                                                      C
7     !  Author: M. Syamlal                                 Date: 18-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     !  Literature/Document References:                                     C
15     !                                                                      C
16     !  Variables referenced:                                               C
17     !  Variables modified:                                                 C
18     !                                                                      C
19     !  Local variables:                                                    C
20     !                                                                      C
21     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
22     !
23           SUBROUTINE CALC_TAU_W_G(TAU_W_G)
24     !...Translated by Pacific-Sierra Research VAST-90 2.06G5  12:17:31  12/09/98
25     !...Switches: -xf
26     !
27     !  Include param.inc file to specify parameter values
28     !
29     !-----------------------------------------------
30     !   M o d u l e s
31     !-----------------------------------------------
32           USE param
33           USE param1
34           USE parallel
35           USE matrix
36           USE scales
37           USE constant
38           USE physprop
39           USE fldvar
40           USE visc_g
41           USE rxns
42           USE run
43           USE toleranc
44           USE geometry
45           USE indices
46           USE is
47           USE sendrecv
48           USE compar
49           USE fun_avg
50           USE functions
51     
52     !=======================================================================
53     ! JFD: START MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
54     !=======================================================================
55           USE bc
56           USE quadric
57           USE cutcell
58     !=======================================================================
59     ! JFD: END MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
60     !=======================================================================
61           IMPLICIT NONE
62     !-----------------------------------------------
63     !   G l o b a l   P a r a m e t e r s
64     !-----------------------------------------------
65     !-----------------------------------------------
66     !   D u m m y   A r g u m e n t s
67     !-----------------------------------------------
68     !
69     !
70     !                      Error index
71           INTEGER          IER
72     !
73     !                      TAU_W_g
74           DOUBLE PRECISION TAU_W_g(DIMENSION_3)
75     !
76     !                      Indices
77           INTEGER          IJK, J, I, IM, IJKP, IMJK, IJKN, IJKNT, IJKS,&
78                            IJKST, IJMKP, IJMK, IJKE, IJKTE, IJKW, IJKTW,&
79                            IMJKP, K, IJKT, JM, KP, IJKM
80     !
81     !                      Average volume fraction
82           DOUBLE PRECISION EPGA
83     !
84     !                      Average gradients
85           DOUBLE PRECISION duodz
86     !
87     !                      Source terms (Surface)
88           DOUBLE PRECISION Sbv, Ssx, Ssy, Ssz
89     !
90     !                      Source terms (Volumetric)
91           DOUBLE PRECISION Vxz
92     !
93     !=======================================================================
94     ! JFD: START MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
95     !=======================================================================
96           DOUBLE PRECISION :: DEL_H,Nx,Ny,Nz
97           LOGICAL :: U_NODE_AT_ET,U_NODE_AT_EB,U_NODE_AT_WT,U_NODE_AT_WB
98           LOGICAL :: V_NODE_AT_NT,V_NODE_AT_NB,V_NODE_AT_ST,V_NODE_AT_SB
99     
100           DOUBLE PRECISION :: dudz_at_E,dudz_at_W
101           DOUBLE PRECISION :: dvdz_at_N,dvdz_at_S
102           DOUBLE PRECISION :: Xi,Yi,Zi,Ui,Vi,Sx,Sy,Sz
103           DOUBLE PRECISION :: MU_GT_CUT,SSX_CUT,SSY_CUT
104           DOUBLE PRECISION :: UW_g,VW_g,WW_g
105           INTEGER :: BCV
106           CHARACTER(LEN=9) :: BCT
107     !=======================================================================
108     ! JFD: END MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
109     !=======================================================================
110     
111     !$omp  parallel do default(none) &
112     !$omp              private( IJK, I, IJKE, EPGA, J, JM, K,  &
113     !$omp              IMJK,IJKN,IJKS,IJMK,IJKT,IJKTE,IJKNT,  &
114     !$omp              IJKM, IJMKP, IJKTW, IMJKP, IJKP, IM, IJKW, IJKST, KP, &
115     !$omp              SBV,  SSX,SSY, SSZ, &
116     !$omp              BCV,NOC_UG, uw_g, vw_g, ww_g, cut_tau_ug, bct,mu_gt_cut,del_h,nx,ny,nz,xi, vi,yi,zi,sx,sy,sz,ssy_cut,noc_wg,cut_tau_wg,ui,dudz_at_e,dudz_at_w,ssx_cut,u_node_at_wb,u_node_at_wt,u_node_at_eb,u_node_at_et,vxz,duodz,dvdz_at_S,dvdz_at_N,v_node_at_sb,v_node_at_st,v_node_at_nb,v_node_at_nt) &
117     !$omp  shared(ijkstart3, ijkend3, do_k, bc_type, cut_u_cell_at, bc_u_id,bc_uw_g, bc_vw_g, bc_ww_g, x_u, y_u, z_u, vol, x_v, y_v,z_v,area_u_cut,oneodx_e_u,oneodx_e_w,x_w,y_w,z_w,bc_hw_g,blocked_w_cell_at,cylindrical,ox_e,mu_g,w_g,odx_e,axy_u,ox,vol_u,odx,odz,axz_u,wall_w_at,jm1,k_of,cartesian_grid,lambda_gt,trd_g,ayz,u_g,ayz_u,i_of,ip1,j_of,km1,cg_safe_mode,mu_gt,ep_g,tau_w_g,bc_w_id,oneodz_t_u,oneodz_t_v,wall_u_at,blocked_u_cell_at,oneodz_t_w,cut_w_cell_at,area_w_cut,kp1,ayz_w,axz_w,dy,wall_v_at,dz,axy_w,vol_w,blocked_v_cell_at,im1,axy,v_g,odz_t)
118           DO IJK = IJKSTART3, IJKEND3
119              K = K_OF(IJK)
120              IJKT = TOP_OF(IJK)
121              EPGA = AVG_Z(EP_G(IJK),EP_G(IJKT),K)
122              IF ( .NOT.IP_AT_T(IJK) .AND. EPGA>DIL_EP_S) THEN
123                 J = J_OF(IJK)
124                 I = I_OF(IJK)
125                 IM = IM1(I)
126                 JM = JM1(J)
127                 IJKP = KP_OF(IJK)
128                 IMJK = IM_OF(IJK)
129                 IJKN = NORTH_OF(IJK)
130                 IJKNT = TOP_OF(IJKN)
131                 IJKS = SOUTH_OF(IJK)
132                 IJKST = TOP_OF(IJKS)
133                 IJMKP = JM_OF(IJKP)
134                 IJMK = JM_OF(IJK)
135                 IJKE = EAST_OF(IJK)
136                 IJKTE = EAST_OF(IJKT)
137                 IJKW = WEST_OF(IJK)
138                 IJKTW = WEST_OF(IJKT)
139                 IMJKP = KP_OF(IMJK)
140                 KP = KP1(K)
141                 IJKM = KM_OF(IJK)
142     
143     
144     !=======================================================================
145     ! JFD: START MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
146     !=======================================================================
147                 IF((.NOT.CARTESIAN_GRID).OR.(CG_SAFE_MODE(5)==1)) THEN
148     !
149     !       Surface forces
150     !
151     !         bulk viscosity term
152                    SBV = (LAMBDA_GT(IJKT)*TRD_G(IJKT)-LAMBDA_GT(IJK)*TRD_G(IJK))*AXY(&
153                       IJK)
154     !
155     !         shear stress terms
156                    SSX = AVG_Z_H(AVG_X_H(MU_GT(IJK),MU_GT(IJKE),I),AVG_X_H(MU_GT(IJKT)&
157                       ,MU_GT(IJKTE),I),K)*(U_G(IJKP)-U_G(IJK))*OX_E(I)*ODZ_T(K)*AYZ_W(&
158                       IJK) - AVG_Z_H(AVG_X_H(MU_GT(IJKW),MU_GT(IJK),IM),AVG_X_H(MU_GT(&
159                       IJKTW),MU_GT(IJKT),IM),K)*(U_G(IMJKP)-U_G(IMJK))*ODZ_T(K)*DY(J)*&
160                       (HALF*(DZ(K)+DZ(KP)))
161                       !same as oX_E(IM)*AYZ_W(IMJK), but avoids singularity
162                    SSY = AVG_Z_H(AVG_Y_H(MU_GT(IJK),MU_GT(IJKN),J),AVG_Y_H(MU_GT(IJKT)&
163                       ,MU_GT(IJKNT),J),K)*(V_G(IJKP)-V_G(IJK))*OX(I)*ODZ_T(K)*AXZ_W(&
164                       IJK) - AVG_Z_H(AVG_Y_H(MU_GT(IJKS),MU_GT(IJK),JM),AVG_Y_H(MU_GT(&
165                       IJKST),MU_GT(IJKT),JM),K)*(V_G(IJMKP)-V_G(IJMK))*OX(I)*ODZ_T(K)*&
166                       AXZ_W(IJMK)
167     !
168                    SSZ = MU_GT(IJKT)*(W_G(IJKP)-W_G(IJK))*OX(I)*ODZ(KP)*AXY_W(IJK) - &
169                       MU_GT(IJK)*(W_G(IJK)-W_G(IJKM))*OX(I)*ODZ(K)*AXY_W(IJKM)
170     
171                 ELSE  !CARTESIAN GRID CASE
172     
173     !       Surface forces
174     !
175     !         bulk viscosity term
176                    SBV =  (LAMBDA_GT(IJKT)*TRD_G(IJKT)) * AXY_W(IJK) &
177                          -(LAMBDA_GT(IJK) *TRD_G(IJK) ) * AXY_W(IJKM)
178     !
179     !         shear stress terms
180     
181                    IF(.NOT.CUT_W_CELL_AT(IJK))   THEN
182     
183                       SSX = AVG_Z_H(AVG_X_H(MU_GT(IJK),MU_GT(IJKE),I),AVG_X_H(MU_GT(IJKT)&
184                          ,MU_GT(IJKTE),I),K)*(U_G(IJKP)-U_G(IJK))*ONEoDZ_T_U(IJK)*AYZ_W(&
185                          IJK) - AVG_Z_H(AVG_X_H(MU_GT(IJKW),MU_GT(IJK),IM),AVG_X_H(MU_GT(&
186                          IJKTW),MU_GT(IJKT),IM),K)*(U_G(IMJKP)-U_G(IMJK))*ONEoDZ_T_U(IMJK)*AYZ_W(IMJK)
187     
188                       SSY = AVG_Z_H(AVG_Y_H(MU_GT(IJK),MU_GT(IJKN),J),AVG_Y_H(MU_GT(IJKT)&
189                          ,MU_GT(IJKNT),J),K)*(V_G(IJKP)-V_G(IJK))*ONEoDZ_T_V(IJK)*AXZ_W(&
190                          IJK) - AVG_Z_H(AVG_Y_H(MU_GT(IJKS),MU_GT(IJK),JM),AVG_Y_H(MU_GT(&
191                          IJKST),MU_GT(IJKT),JM),K)*(V_G(IJMKP)-V_G(IJMK))*ONEoDZ_T_V(IJMK)*&
192                          AXZ_W(IJMK)
193     !
194                       SSZ = MU_GT(IJKT)*(W_G(IJKP)-W_G(IJK))*ONEoDZ_T_W(IJK)*AXY_W(IJK) - &
195                          MU_GT(IJK)*(W_G(IJK)-W_G(IJKM))*ONEoDZ_T_W(IJKM)*AXY_W(IJKM)
196     
197                    ELSE  ! CUT CELL
198     
199                       BCV = BC_W_ID(IJK)
200     
201                       IF(BCV > 0 ) THEN
202                          BCT = BC_TYPE(BCV)
203                       ELSE
204                          BCT = 'NONE'
205                       ENDIF
206     
207                       SELECT CASE (BCT)
208                          CASE ('CG_NSW')
209                             CUT_TAU_WG = .TRUE.
210                             NOC_WG     = .TRUE.
211                             UW_g = ZERO
212                             VW_g = ZERO
213                             WW_g = ZERO
214                          CASE ('CG_FSW')
215                             CUT_TAU_WG = .FALSE.
216                             NOC_WG     = .FALSE.
217                             UW_g = ZERO
218                             VW_g = ZERO
219                             WW_g = ZERO
220                          CASE('CG_PSW')
221                             IF(BC_HW_G(BC_W_ID(IJK))==UNDEFINED) THEN   ! same as NSW
222                                CUT_TAU_WG = .TRUE.
223                                NOC_WG     = .TRUE.
224                                UW_g = BC_UW_G(BCV)
225                                VW_g = BC_VW_G(BCV)
226                                WW_g = BC_WW_G(BCV)
227                             ELSEIF(BC_HW_G(BC_W_ID(IJK))==ZERO) THEN   ! same as FSW
228                                CUT_TAU_WG = .FALSE.
229                                NOC_WG     = .FALSE.
230                                UW_g = ZERO
231                                VW_g = ZERO
232                                WW_g = ZERO
233                             ELSE                              ! partial slip
234                                CUT_TAU_WG = .FALSE.
235                                NOC_WG     = .FALSE.
236                             ENDIF
237                          CASE ('NONE')
238                             TAU_W_G(IJK) = ZERO
239                             CYCLE
240                       END SELECT
241     
242                       IF(CUT_TAU_WG) THEN
243                          MU_GT_CUT = (VOL(IJK)*MU_GT(IJK) + VOL(IJKP)*MU_GT(IJKT))/(VOL(IJK) + VOL(IJKP))
244                       ELSE
245                          MU_GT_CUT = ZERO
246                       ENDIF
247     
248     !           SSX:
249                       U_NODE_AT_ET = ((.NOT.BLOCKED_U_CELL_AT(IJKP)).AND.(.NOT.WALL_U_AT(IJKP)))
250                       U_NODE_AT_EB = ((.NOT.BLOCKED_U_CELL_AT(IJK)).AND.(.NOT.WALL_U_AT(IJK)))
251                       U_NODE_AT_WT = ((.NOT.BLOCKED_U_CELL_AT(IMJKP)).AND.(.NOT.WALL_U_AT(IMJKP)))
252                       U_NODE_AT_WB = ((.NOT.BLOCKED_U_CELL_AT(IMJK)).AND.(.NOT.WALL_U_AT(IMJK)))
253     
254                       IF(U_NODE_AT_ET.AND.U_NODE_AT_EB) THEN
255     
256                          Ui = HALF * (U_G(IJKP) + U_G(IJK))
257                          Xi = HALF * (X_U(IJKP) + X_U(IJK))
258                          Yi = HALF * (Y_U(IJKP) + Y_U(IJK))
259                          Zi = HALF * (Z_U(IJKP) + Z_U(IJK))
260                          Sx = X_U(IJKP) - X_U(IJK)
261                          Sy = Y_U(IJKP) - Y_U(IJK)
262                          Sz = Z_U(IJKP) - Z_U(IJK)
263     
264                          CALL GET_DEL_H(IJK,'W_MOMENTUM',Xi,Yi,Zi,Del_H,Nx,Ny,Nz)
265     
266                          dudz_at_E =  (U_G(IJKP) - U_G(IJK)) * ONEoDZ_T_U(IJK)
267     
268                          IF(NOC_WG) dudz_at_E = dudz_at_E - ((Ui-UW_g) * ONEoDZ_T_U(IJK)/DEL_H*(Sx*Nx+Sy*Ny))
269     
270                       ELSE
271                          dudz_at_E =  ZERO
272                       ENDIF
273     
274                       IF(U_NODE_AT_WT.AND.U_NODE_AT_WB) THEN
275     
276                          Ui = HALF * (U_G(IMJKP) + U_G(IMJK))
277                          Xi = HALF * (X_U(IMJKP) + X_U(IMJK))
278                          Yi = HALF * (Y_U(IMJKP) + Y_U(IMJK))
279                          Zi = HALF * (Z_U(IMJKP) + Z_U(IMJK))
280                          Sx = X_U(IMJKP) - X_U(IMJK)
281                          Sy = Y_U(IMJKP) - Y_U(IMJK)
282                          Sz = Z_U(IMJKP) - Z_U(IMJK)
283     
284                          CALL GET_DEL_H(IJK,'W_MOMENTUM',Xi,Yi,Zi,Del_H,Nx,Ny,Nz)
285     
286                          dudz_at_W =  (U_G(IMJKP) - U_G(IMJK)) * ONEoDZ_T_U(IMJK)
287     
288                          IF(NOC_WG) dudz_at_W = dudz_at_W - ((Ui-UW_g) * ONEoDZ_T_U(IMJK)/DEL_H*(Sx*Nx+Sy*Ny))
289     
290                       ELSE
291                          dudz_at_W =  ZERO
292                       ENDIF
293     
294                       IF(U_NODE_AT_EB) THEN
295                          CALL GET_DEL_H(IJK,'W_MOMENTUM',X_U(IJK),Y_U(IJK),Z_U(IJK),Del_H,Nx,Ny,Nz)
296                          SSX_CUT = - MU_GT_CUT * (U_G(IJK) - UW_g) / DEL_H * (Nz*Nx) * Area_W_CUT(IJK)
297                       ELSE
298                          SSX_CUT =  ZERO
299                       ENDIF
300     
301                       SSX = AVG_Z_H(AVG_X_H(MU_GT(IJK),MU_GT(IJKE),I),AVG_X_H(MU_GT(IJKT)&
302                          ,MU_GT(IJKTE),I),K)*dudz_at_E*AYZ_W(&
303                          IJK) - AVG_Z_H(AVG_X_H(MU_GT(IJKW),MU_GT(IJK),IM),AVG_X_H(MU_GT(&
304                          IJKTW),MU_GT(IJKT),IM),K)*dudz_at_W*AYZ_W(IMJK) &
305                         + SSX_CUT
306     
307     !           SSY:
308     
309                       V_NODE_AT_NT = ((.NOT.BLOCKED_V_CELL_AT(IJKP)).AND.(.NOT.WALL_V_AT(IJKP)))
310                       V_NODE_AT_NB = ((.NOT.BLOCKED_V_CELL_AT(IJK)).AND.(.NOT.WALL_V_AT(IJK)))
311                       V_NODE_AT_ST = ((.NOT.BLOCKED_V_CELL_AT(IJMKP)).AND.(.NOT.WALL_V_AT(IJMKP)))
312                       V_NODE_AT_SB = ((.NOT.BLOCKED_V_CELL_AT(IJMK)).AND.(.NOT.WALL_V_AT(IJMK)))
313     
314                       IF(V_NODE_AT_NT.AND.V_NODE_AT_NB) THEN
315     
316                          Vi = HALF * (V_G(IJKP) + V_G(IJK))
317                          Xi = HALF * (X_V(IJKP) + X_V(IJK))
318                          Yi = HALF * (Y_V(IJKP) + Y_V(IJK))
319                          Zi = HALF * (Z_V(IJKP) + Z_V(IJK))
320                          Sx = X_V(IJKP) - X_V(IJK)
321                          Sy = Y_V(IJKP) - Y_V(IJK)
322                          Sz = Z_V(IJKP) - Z_V(IJK)
323     
324     
325                          CALL GET_DEL_H(IJK,'W_MOMENTUM',Xi,Yi,Zi,Del_H,Nx,Ny,Nz)
326     
327                          dvdz_at_N =  (V_G(IJKP) - V_G(IJK)) * ONEoDZ_T_V(IJK)
328     
329                          IF(NOC_WG) dvdz_at_N = dvdz_at_N - ((Vi-VW_g) * ONEoDZ_T_V(IJK)/DEL_H*(Sx*Nx+Sy*Ny))
330     
331                       ELSE
332                             dvdz_at_N =  ZERO
333                       ENDIF
334     
335                       IF(V_NODE_AT_ST.AND.V_NODE_AT_SB) THEN
336     
337                          Vi = HALF * (V_G(IJMKP) + V_G(IJMK))
338                          Xi = HALF * (X_V(IJMKP) + X_V(IJMK))
339                          Yi = HALF * (Y_V(IJMKP) + Y_V(IJMK))
340                          Zi = HALF * (Z_V(IJMKP) + Z_V(IJMK))
341                          Sx = X_V(IJMKP) - X_V(IJMK)
342                          Sy = Y_V(IJMKP) - Y_V(IJMK)
343                          Sz = Z_V(IJMKP) - Z_V(IJMK)
344     
345                          CALL GET_DEL_H(IJK,'W_MOMENTUM',Xi,Yi,Zi,Del_H,Nx,Ny,Nz)
346     
347                          dvdz_at_S =  (V_G(IJMKP) - V_G(IJMK)) * ONEoDZ_T_V(IJMK)
348     
349                          IF(NOC_WG) dvdz_at_S = dvdz_at_S - ((Vi-VW_g) * ONEoDZ_T_V(IJMK)/DEL_H*(Sx*Nx+Sy*Ny))
350     
351                       ELSE
352                          dvdz_at_S =  ZERO
353                       ENDIF
354     
355                       IF(V_NODE_AT_NB) THEN
356                          CALL GET_DEL_H(IJK,'W_MOMENTUM',X_V(IJK),Y_V(IJK),Z_V(IJK),Del_H,Nx,Ny,Nz)
357                          SSY_CUT = - MU_GT_CUT * (V_G(IJK) - VW_g) / DEL_H * (Nz*Ny) * Area_W_CUT(IJK)
358                       ELSE
359                          SSY_CUT =  ZERO
360                       ENDIF
361     
362                       SSY = AVG_Z_H(AVG_Y_H(MU_GT(IJK),MU_GT(IJKN),J),AVG_Y_H(MU_GT(IJKT)&
363                          ,MU_GT(IJKNT),J),K)*dvdz_at_N*AXZ_W(&
364                          IJK) - AVG_Z_H(AVG_Y_H(MU_GT(IJKS),MU_GT(IJK),JM),AVG_Y_H(MU_GT(&
365                          IJKST),MU_GT(IJKT),JM),K)*dvdz_at_S*AXZ_W(IJMK) &
366                         + SSY_CUT
367     
368     !           SSZ:
369                       CALL GET_DEL_H(IJK,'W_MOMENTUM',X_W(IJK),Y_W(IJK),Z_W(IJK),Del_H,Nx,Ny,Nz)
370     
371                       SSZ = MU_GT(IJKT)*(W_G(IJKP)-W_G(IJK))*ONEoDZ_T_W(IJK)*AXY_W(IJK) - &
372                             MU_GT(IJK)*(W_G(IJK)-W_G(IJKM))*OX(I)*ONEoDZ_T_W(IJKM)*AXY_W(IJKM) &
373                           - MU_GT_CUT * (W_g(IJK) - WW_g) / DEL_H * (Nz**2) * Area_W_CUT(IJK)
374     
375     
376     
377                    END IF  ! CUT CELL
378     
379     ! Original terms
380     
381     !            SSX = AVG_Z_H(AVG_X_H(MU_GT(IJK),MU_GT(IJKE),I),AVG_X_H(MU_GT(IJKT)&
382     !              ,MU_GT(IJKTE),I),K)*(U_G(IJKP)-U_G(IJK))*OX_E(I)*ODZ_T(K)*AYZ_W(&
383     !               IJK) - AVG_Z_H(AVG_X_H(MU_GT(IJKW),MU_GT(IJK),IM),AVG_X_H(MU_GT(&
384     !               IJKTW),MU_GT(IJKT),IM),K)*(U_G(IMJKP)-U_G(IMJK))*ODZ_T(K)*DY(J)*&
385     !               (HALF*(DZ(K)+DZ(KP)))
386                  !same as oX_E(IM)*AYZ_W(IMJK), but avoids singularity
387     !            SSY = AVG_Z_H(AVG_Y_H(MU_GT(IJK),MU_GT(IJKN),J),AVG_Y_H(MU_GT(IJKT)&
388     !               ,MU_GT(IJKNT),J),K)*(V_G(IJKP)-V_G(IJK))*OX(I)*ODZ_T(K)*AXZ_W(&
389     !               IJK) - AVG_Z_H(AVG_Y_H(MU_GT(IJKS),MU_GT(IJK),JM),AVG_Y_H(MU_GT(&
390     !               IJKST),MU_GT(IJKT),JM),K)*(V_G(IJMKP)-V_G(IJMK))*OX(I)*ODZ_T(K)*&
391     !               AXZ_W(IJMK)
392     !
393     !            SSZ = MU_GT(IJKT)*(W_G(IJKP)-W_G(IJK))*OX(I)*ODZ(KP)*AXY_W(IJK) - &
394     !               MU_GT(IJK)*(W_G(IJK)-W_G(IJKM))*OX(I)*ODZ(K)*AXY_W(IJKM)
395     !
396                 ENDIF  ! CARTESIAN GRID
397     !=======================================================================
398     ! JFD: END MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
399     !=======================================================================
400     !
401     !
402     !         Special terms for cylindrical coordinates
403                 IF (CYLINDRICAL) THEN
404     !
405     !              modify Szz to include integral of (1/x)*(d/dz)(2*u/x)
406                    SSZ = SSZ + MU_GT(IJKT)*(U_G(IJKP)+U_G(IMJKP))*OX(I)*AXY_W(IJK)&
407                        - MU_GT(IJK)*(U_G(IJK)+U_G(IMJK))*OX(I)*AXY_W(IJKM)
408     !
409     !              (mu_s/x)* part of (du/dz*x) tau_xz/X
410                    IF (OX_E(IM) /= UNDEFINED) THEN
411                       DUODZ = (U_G(IMJKP)-U_G(IMJK))*OX_E(IM)*ODZ_T(K)
412                    ELSE
413                       DUODZ = ZERO
414                    ENDIF
415                    VXZ = AVG_Z(MU_GT(IJK),MU_GT(IJKT),K)*OX(I)*HALF*((U_G(IJKP)-U_G&
416                       (IJK))*OX_E(I)*ODZ_T(K)+DUODZ)
417     !
418                 ELSE
419                    VXZ = ZERO
420                 ENDIF
421     !
422     !         Add the terms
423                 TAU_W_G(IJK) =  SBV + SSX + SSY + SSZ + VXZ*VOL_W(IJK)
424              ELSE
425                 TAU_W_G(IJK) = ZERO
426              ENDIF
427           END DO
428           call send_recv(tau_w_g,2)
429           RETURN
430           END SUBROUTINE CALC_TAU_W_G
431     
432     !// Comments on the modifications for DMP version implementation
433     !// 001 Include header file and common declarations for parallelization
434     !// 350 Changed do loop limits: 1,ijkmax2-> ijkstart3, ijkend3
435     !// 400 Added sendrecv module and send_recv calls for COMMunication
436