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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Module name: CALC_Tau_V_g(TAU_V_g, IER)                             C
4     !  Purpose: Cross terms in the gradient of stress in V_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_V_G(TAU_V_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     
62           IMPLICIT NONE
63     !-----------------------------------------------
64     !   G l o b a l   P a r a m e t e r s
65     !-----------------------------------------------
66     !-----------------------------------------------
67     !   D u m m y   A r g u m e n t s
68     !-----------------------------------------------
69     !
70     !                      TAU_V_g
71           DOUBLE PRECISION TAU_V_g(DIMENSION_3)
72     !
73     !                      Indices
74           INTEGER          I, J, K, IJK, IJKN, JP, IM,  KM, IJPK, IJMK,&
75                            IJKE, IJKNE, IJKW, IJKNW, IMJPK, IMJK, IJKT,&
76                            IJKTN, IJKB, IJKBN, IJKM, IJPKM
77     !
78     !                      Average volume fraction
79           DOUBLE PRECISION EPGA
80     !
81     !                      Source terms (Surface)
82           DOUBLE PRECISION Sbv, Ssx, Ssy, Ssz
83     !
84     !=======================================================================
85     ! JFD: START MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
86     !=======================================================================
87           DOUBLE PRECISION :: DEL_H,Nx,Ny,Nz
88           LOGICAL :: U_NODE_AT_NE,U_NODE_AT_NW,U_NODE_AT_SE,U_NODE_AT_SW
89           LOGICAL :: W_NODE_AT_TN,W_NODE_AT_TS,W_NODE_AT_BN,W_NODE_AT_BS
90           DOUBLE PRECISION :: U_SUM,X_SUM,Y_SUM,Z_SUM
91           DOUBLE PRECISION :: dudy_at_E,dudy_at_W
92           DOUBLE PRECISION :: dwdy_at_T,dwdy_at_B
93           DOUBLE PRECISION :: Xi,Yi,Zi,Ui,Wi,Sx,Sy,Sz
94           DOUBLE PRECISION :: MU_GT_CUT,SSX_CUT,SSZ_CUT
95           DOUBLE PRECISION :: UW_g,VW_g,WW_g
96           INTEGER :: N_SUM
97           INTEGER :: BCV
98           CHARACTER(LEN=9) :: BCT
99     !=======================================================================
100     ! JFD: END MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
101     !=======================================================================
102     
103     !$omp  parallel do default(none) &
104     !$omp              private( IJK, I, IM, IJPK, IJKW, IJKTN, IJPKM, IJKE, IJKNW, IMJPK, IJKBN, EPGA, J, JP,K, KM,  &
105     !$omp              IMJK,IJKN,IJKNE,IJMK,IJKT,  &
106     !$omp              IJKB,IJKM, &
107     !$omp              SBV,  SSX,SSY, SSZ, &
108     !$omp              BCV,NOC_UG, uw_g, vw_g, ww_g, cut_tau_vg, bct,mu_gt_cut,del_h,nx,ny,nz,xi,yi,zi,sx,sy,sz,wi,ssz_cut,ui,noc_vg,u_node_at_sw,u_node_at_se,u_node_at_nw,u_node_at_ne,u_sum,x_sum,y_sum,z_sum,n_sum,dudy_at_W,dudy_at_E,dwdy_at_B,dwdy_at_T,w_node_at_bs,w_node_at_bn,w_node_at_ts,w_node_at_tn,ssx_cut) &
109     !$omp  shared(ijkstart3, ijkend3, do_k, bc_type, cut_u_cell_at, bc_u_id, tau_v_g, bc_uw_g, bc_vw_g, bc_ww_g, x_u, y_u, z_u, vol, x_v, y_v,z_v, wall_v_at, area_u_cut,oneodx_e_u,oneodx_e_v,oneodx_e_w,x_w,y_w,z_w,bc_hw_g,blocked_v_cell_at,blocked_w_cell_at,cylindrical,ox_e,mu_g,v_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,blocked_u_cell_at,oneody_n_v,ayz_v,axz_v,cut_v_cell_at,oneody_n_w,bc_v_id,wall_u_at,area_v_cut,axz,ody_n,ody,axy_v,oneody_n_u,jp1,im1)
110           DO IJK = IJKSTART3, IJKEND3
111              J = J_OF(IJK)
112              IJKN = NORTH_OF(IJK)
113              EPGA = AVG_Y(EP_G(IJK),EP_G(IJKN),J)
114              IF ( .NOT.IP_AT_N(IJK) .AND. EPGA>DIL_EP_S) THEN
115                 JP = JP1(J)
116                 I = I_OF(IJK)
117                 IM = IM1(I)
118                 K = K_OF(IJK)
119                 KM = KM1(K)
120                 IJPK = JP_OF(IJK)
121                 IJMK = JM_OF(IJK)
122                 IJKE = EAST_OF(IJK)
123                 IJKNE = EAST_OF(IJKN)
124                 IJKW = WEST_OF(IJK)
125                 IJKNW = NORTH_OF(IJKW)
126                 IMJPK = IM_OF(IJPK)
127                 IMJK = IM_OF(IJK)
128                 IJKT = TOP_OF(IJK)
129                 IJKTN = NORTH_OF(IJKT)
130                 IJKB = BOTTOM_OF(IJK)
131                 IJKBN = NORTH_OF(IJKB)
132                 IJKM = KM_OF(IJK)
133                 IJPKM = JP_OF(IJKM)
134     !=======================================================================
135     ! JFD: START MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
136     !=======================================================================
137                 IF((.NOT.CARTESIAN_GRID).OR.(CG_SAFE_MODE(4)==1)) THEN
138     !
139     !       Surface forces
140     !
141     !         bulk viscosity term
142                    SBV = (LAMBDA_GT(IJKN)*TRD_G(IJKN)-LAMBDA_GT(IJK)*TRD_G(IJK))*AXZ(&
143                       IJK)
144     !
145     !         shear stress terms
146                    SSX = AVG_Y_H(AVG_X_H(MU_GT(IJK),MU_GT(IJKE),I),AVG_X_H(MU_GT(IJKN)&
147                       ,MU_GT(IJKNE),I),J)*(U_G(IJPK)-U_G(IJK))*ODY_N(J)*AYZ_V(IJK) - &
148                       AVG_Y_H(AVG_X_H(MU_GT(IJKW),MU_GT(IJK),IM),AVG_X_H(MU_GT(IJKNW),&
149                       MU_GT(IJKN),IM),J)*(U_G(IMJPK)-U_G(IMJK))*ODY_N(J)*AYZ_V(IMJK)
150                    SSY = MU_GT(IJKN)*(V_G(IJPK)-V_G(IJK))*ODY(JP)*AXZ_V(IJK) - MU_GT(&
151                       IJK)*(V_G(IJK)-V_G(IJMK))*ODY(J)*AXZ_V(IJMK)
152                    SSZ = AVG_Y_H(AVG_Z_H(MU_GT(IJK),MU_GT(IJKT),K),AVG_Z_H(MU_GT(IJKN)&
153                       ,MU_GT(IJKTN),K),J)*(W_G(IJPK)-W_G(IJK))*ODY_N(J)*AXY_V(IJK) - &
154                       AVG_Y_H(AVG_Z_H(MU_GT(IJKB),MU_GT(IJK),KM),AVG_Z_H(MU_GT(IJKBN),&
155                       MU_GT(IJKN),KM),J)*(W_G(IJPKM)-W_G(IJKM))*ODY_N(J)*AXY_V(IJKM)
156     
157                 ELSE  ! CARTESIAN GRID CASE
158     !
159     !       Surface forces
160     !
161     !         bulk viscosity term
162                    SBV =  (LAMBDA_GT(IJKN)*TRD_G(IJKN)) * AXZ_V(IJK) &
163                          -(LAMBDA_GT(IJK) *TRD_G(IJK) ) * AXZ_V(IJMK)
164     !
165     !         shear stress terms
166     
167                    IF(.NOT.CUT_V_CELL_AT(IJK))   THEN
168     
169                       SSX = AVG_Y_H(AVG_X_H(MU_GT(IJK),MU_GT(IJKE),I),AVG_X_H(MU_GT(IJKN)&
170                          ,MU_GT(IJKNE),I),J)*(U_G(IJPK)-U_G(IJK))*ONEoDY_N_U(IJK)*AYZ_V(IJK) - &
171                          AVG_Y_H(AVG_X_H(MU_GT(IJKW),MU_GT(IJK),IM),AVG_X_H(MU_GT(IJKNW),&
172                          MU_GT(IJKN),IM),J)*(U_G(IMJPK)-U_G(IMJK))*ONEoDY_N_U(IMJK)*AYZ_V(IMJK)
173     
174                       SSY = MU_GT(IJKN)*(V_G(IJPK)-V_G(IJK))*ONEoDY_N_V(IJK)*AXZ_V(IJK) - MU_GT(&
175                          IJK)*(V_G(IJK)-V_G(IJMK))*ONEoDY_N_V(IJMK)*AXZ_V(IJMK)
176     
177     
178     
179                       IF(DO_K) THEN
180     
181                          SSZ = AVG_Y_H(AVG_Z_H(MU_GT(IJK),MU_GT(IJKT),K),AVG_Z_H(MU_GT(IJKN)&
182                             ,MU_GT(IJKTN),K),J)*(W_G(IJPK)-W_G(IJK))*ONEoDY_N_W(IJK)*AXY_V(IJK) - &
183                             AVG_Y_H(AVG_Z_H(MU_GT(IJKB),MU_GT(IJK),KM),AVG_Z_H(MU_GT(IJKBN),&
184                             MU_GT(IJKN),KM),J)*(W_G(IJPKM)-W_G(IJKM))*ONEoDY_N_W(IJKM)*AXY_V(IJKM)
185     
186                       ELSE
187                          SSZ = ZERO
188                       ENDIF
189     
190                    ELSE  ! CUT CELL
191     
192                       BCV = BC_V_ID(IJK)
193     
194                       IF(BCV > 0 ) THEN
195                          BCT = BC_TYPE(BCV)
196                       ELSE
197                          BCT = 'NONE'
198                       ENDIF
199     
200                       SELECT CASE (BCT)
201                          CASE ('CG_NSW')
202                             CUT_TAU_VG = .TRUE.
203                             NOC_VG     = .TRUE.
204                             UW_g = ZERO
205                             VW_g = ZERO
206                             WW_g = ZERO
207                          CASE ('CG_FSW')
208                             CUT_TAU_VG = .FALSE.
209                             NOC_VG     = .FALSE.
210                             UW_g = ZERO
211                             VW_g = ZERO
212                             WW_g = ZERO
213                          CASE('CG_PSW')
214                             IF(BC_HW_G(BC_V_ID(IJK))==UNDEFINED) THEN   ! same as NSW
215                                CUT_TAU_VG = .TRUE.
216                                NOC_VG     = .TRUE.
217                                UW_g = BC_UW_G(BCV)
218                                VW_g = BC_VW_G(BCV)
219                                WW_g = BC_WW_G(BCV)
220                             ELSEIF(BC_HW_G(BC_V_ID(IJK))==ZERO) THEN   ! same as FSW
221                                CUT_TAU_VG = .FALSE.
222                                NOC_VG     = .FALSE.
223                                UW_g = ZERO
224                                VW_g = ZERO
225                                WW_g = ZERO
226                             ELSE                              ! partial slip
227                                CUT_TAU_VG = .FALSE.
228                                NOC_VG     = .FALSE.
229                             ENDIF
230                          CASE ('NONE')
231                             TAU_V_G(IJK) = ZERO
232                             CYCLE
233                       END SELECT
234     
235     
236                       IF(CUT_TAU_VG) THEN
237                          MU_GT_CUT = (VOL(IJK)*MU_GT(IJK) + VOL(IJPK)*MU_GT(IJKN))/(VOL(IJK) + VOL(IJPK))
238                       ELSE
239                          MU_GT_CUT = ZERO
240                       ENDIF
241     
242     !           SSX:
243     
244                       U_NODE_AT_NE = ((.NOT.BLOCKED_U_CELL_AT(IJPK)).AND.(.NOT.WALL_U_AT(IJPK)))
245                       U_NODE_AT_SE = ((.NOT.BLOCKED_U_CELL_AT(IJK)).AND.(.NOT.WALL_U_AT(IJK)))
246                       U_NODE_AT_NW = ((.NOT.BLOCKED_U_CELL_AT(IMJPK)).AND.(.NOT.WALL_U_AT(IMJPK)))
247                       U_NODE_AT_SW = ((.NOT.BLOCKED_U_CELL_AT(IMJK)).AND.(.NOT.WALL_U_AT(IMJK)))
248     
249     
250                       IF(U_NODE_AT_NE.AND.U_NODE_AT_SE) THEN
251     
252                          Ui = HALF * (U_G(IJPK) + U_G(IJK))
253                          Xi = HALF * (X_U(IJPK) + X_U(IJK))
254                          Yi = HALF * (Y_U(IJPK) + Y_U(IJK))
255                          Zi = HALF * (Z_U(IJPK) + Z_U(IJK))
256                          Sx = X_U(IJPK) - X_U(IJK)
257                          Sy = Y_U(IJPK) - Y_U(IJK)
258                          Sz = Z_U(IJPK) - Z_U(IJK)
259     
260     
261                          CALL GET_DEL_H(IJK,'V_MOMENTUM',Xi,Yi,Zi,Del_H,Nx,Ny,Nz)
262     
263                          dudy_at_E =  (U_G(IJPK) - U_G(IJK)) * ONEoDY_N_U(IJK)
264     
265                          IF(NOC_VG) dudy_at_E = dudy_at_E - ((Ui-UW_g) * ONEoDY_N_U(IJK)/DEL_H*(Sx*Nx+Sz*Nz))
266     
267                       ELSE
268                          dudy_at_E =  ZERO
269                       ENDIF
270     
271                       IF(U_NODE_AT_NW.AND.U_NODE_AT_SW) THEN
272     
273                          Ui = HALF * (U_G(IMJPK) + U_G(IMJK))
274                          Xi = HALF * (X_U(IMJPK) + X_U(IMJK))
275                          Yi = HALF * (Y_U(IMJPK) + Y_U(IMJK))
276                          Zi = HALF * (Z_U(IMJPK) + Z_U(IMJK))
277                          Sx = X_U(IMJPK) - X_U(IMJK)
278                          Sy = Y_U(IMJPK) - Y_U(IMJK)
279                          Sz = Z_U(IMJPK) - Z_U(IMJK)
280     
281                          CALL GET_DEL_H(IJK,'V_MOMENTUM',Xi,Yi,Zi,Del_H,Nx,Ny,Nz)
282     
283                          dudy_at_W =  (U_G(IMJPK) - U_G(IMJK)) * ONEoDY_N_U(IMJK)
284     
285                          IF(NOC_VG) dudy_at_W = dudy_at_W - ((Ui-UW_g) * ONEoDY_N_U(IMJK)/DEL_H*(Sx*Nx+Sz*Nz))
286     
287                       ELSE
288                          dudy_at_W =  ZERO
289                       ENDIF
290     
291                       U_SUM = ZERO
292                       X_SUM = ZERO
293                       Y_SUM = ZERO
294                       Z_SUM = ZERO
295                       N_SUM = 0
296     
297                      IF(U_NODE_AT_SE) THEN
298                         CALL GET_DEL_H(IJK,'V_MOMENTUM',X_U(IJK),Y_U(IJK),Z_U(IJK),Del_H,Nx,Ny,Nz)
299                         SSX_CUT = - MU_GT_CUT * (U_G(IJK) - UW_g) / DEL_H * (Ny*Nx) * Area_V_CUT(IJK)
300                      ELSE
301                         SSX_CUT =  ZERO
302                      ENDIF
303     
304                       SSX = AVG_Y_H(AVG_X_H(MU_GT(IJK),MU_GT(IJKE),I),AVG_X_H(MU_GT(IJKN)&
305                          ,MU_GT(IJKNE),I),J)*dudy_at_E*AYZ_V(IJK) - &
306                          AVG_Y_H(AVG_X_H(MU_GT(IJKW),MU_GT(IJK),IM),AVG_X_H(MU_GT(IJKNW),&
307                          MU_GT(IJKN),IM),J)*dudy_at_W*AYZ_V(IMJK) &
308                         + SSX_CUT
309     
310     !           SSY:
311                       CALL GET_DEL_H(IJK,'V_MOMENTUM',X_V(IJK),Y_V(IJK),Z_V(IJK),Del_H,Nx,Ny,Nz)
312     
313                       SSY = MU_GT(IJKN)*(V_G(IJPK)-V_G(IJK))*ONEoDY_N_V(IJK)*AXZ_V(IJK) - MU_GT(&
314                             IJK)*(V_G(IJK)-V_G(IJMK))*ONEoDY_N_V(IJMK)*AXZ_V(IJMK) &
315                           - MU_GT_CUT * (V_g(IJK) - VW_g) / DEL_H * (Ny**2) * Area_V_CUT(IJK)
316     
317     !           SSZ:
318                       IF(DO_K) THEN
319     
320                          W_NODE_AT_TN = ((.NOT.BLOCKED_W_CELL_AT(IJPK)).AND.(.NOT.WALL_W_AT(IJPK)))
321                          W_NODE_AT_TS = ((.NOT.BLOCKED_W_CELL_AT(IJK)).AND.(.NOT.WALL_W_AT(IJK)))
322                          W_NODE_AT_BN = ((.NOT.BLOCKED_W_CELL_AT(IJPKM)).AND.(.NOT.WALL_W_AT(IJPKM)))
323                          W_NODE_AT_BS = ((.NOT.BLOCKED_W_CELL_AT(IJKM)).AND.(.NOT.WALL_W_AT(IJKM)))
324     
325                          IF(W_NODE_AT_TN.AND.W_NODE_AT_TS) THEN
326     
327                             Wi = HALF * (W_G(IJPK) + W_G(IJK))
328                             Xi = HALF * (X_W(IJPK) + X_W(IJK))
329                             Yi = HALF * (Y_W(IJPK) + Y_W(IJK))
330                             Zi = HALF * (Z_W(IJPK) + Z_W(IJK))
331                             Sx = X_W(IJPK) - X_W(IJK)
332                             Sy = Y_W(IJPK) - Y_W(IJK)
333                             Sz = Z_W(IJPK) - Z_W(IJK)
334     
335                             CALL GET_DEL_H(IJK,'V_MOMENTUM',Xi,Yi,Zi,Del_H,Nx,Ny,Nz)
336     
337                             dwdy_at_T =  (W_G(IJPK) - W_G(IJK)) * ONEoDY_N_W(IJK)
338     
339                             IF(NOC_VG) dwdy_at_T = dwdy_at_T - ((Wi-WW_g) * ONEoDY_N_W(IJK)/DEL_H*(Sx*Nx+Sz*Nz))
340     
341                          ELSE
342                             dwdy_at_T =  ZERO
343                          ENDIF
344     
345                          IF(W_NODE_AT_BN.AND.W_NODE_AT_BS) THEN
346     
347                             Wi = HALF * (W_G(IJPKM) + W_G(IJKM))
348                             Xi = HALF * (X_W(IJPKM) + X_W(IJKM))
349                             Yi = HALF * (Y_W(IJPKM) + Y_W(IJKM))
350                             Zi = HALF * (Z_W(IJPKM) + Z_W(IJKM))
351                             Sx = X_W(IJPKM) - X_W(IJKM)
352                             Sy = Y_W(IJPKM) - Y_W(IJKM)
353                             Sz = Z_W(IJPKM) - Z_W(IJKM)
354     
355                             CALL GET_DEL_H(IJK,'V_MOMENTUM',Xi,Yi,Zi,Del_H,Nx,Ny,Nz)
356     
357                             dwdy_at_B =  (W_G(IJPKM) - W_G(IJKM)) * ONEoDY_N_W(IJKM)
358     
359                             IF(NOC_VG) dwdy_at_B = dwdy_at_B - ((Wi-WW_g) * ONEoDY_N_W(IJKM)/DEL_H*(Sx*Nx+Sz*Nz))
360     
361                          ELSE
362                             dwdy_at_B =  ZERO
363                          ENDIF
364     
365                          IF(W_NODE_AT_TS) THEN
366                             CALL GET_DEL_H(IJK,'V_MOMENTUM',X_W(IJK),Y_W(IJK),Z_W(IJK),Del_H,Nx,Ny,Nz)
367                             SSZ_CUT = - MU_GT_CUT * (W_G(IJK) - WW_g) / DEL_H * (Ny*Nz) * Area_V_CUT(IJK)
368                          ELSE
369                             SSZ_CUT =  ZERO
370                          ENDIF
371     
372                          SSZ = AVG_Y_H(AVG_Z_H(MU_GT(IJK),MU_GT(IJKT),K),AVG_Z_H(MU_GT(IJKN)&
373                             ,MU_GT(IJKTN),K),J)*dwdy_at_T*AXY_V(IJK) - &
374                             AVG_Y_H(AVG_Z_H(MU_GT(IJKB),MU_GT(IJK),KM),AVG_Z_H(MU_GT(IJKBN),&
375                             MU_GT(IJKN),KM),J)*dwdy_at_B*AXY_V(IJKM) &
376                            + SSZ_CUT
377     
378                       ELSE
379     
380                          SSZ = ZERO
381     
382                       ENDIF  ! DO_K
383     
384     
385                    END IF  ! CUT CELL
386     
387     
388     ! Original terms
389     !            SSX = AVG_Y_H(AVG_X_H(MU_GT(IJK),MU_GT(IJKE),I),AVG_X_H(MU_GT(IJKN)&
390     !              ,MU_GT(IJKNE),I),J)*(U_G(IJPK)-U_G(IJK))*ODY_N(J)*AYZ_V(IJK) - &
391     !               AVG_Y_H(AVG_X_H(MU_GT(IJKW),MU_GT(IJK),IM),AVG_X_H(MU_GT(IJKNW),&
392     !               MU_GT(IJKN),IM),J)*(U_G(IMJPK)-U_G(IMJK))*ODY_N(J)*AYZ_V(IMJK)
393     
394     !            SSY = MU_GT(IJKN)*(V_G(IJPK)-V_G(IJK))*ODY(JP)*AXZ_V(IJK) - MU_GT(&
395     !               IJK)*(V_G(IJK)-V_G(IJMK))*ODY(J)*AXZ_V(IJMK)
396     
397     !            SSZ = AVG_Y_H(AVG_Z_H(MU_GT(IJK),MU_GT(IJKT),K),AVG_Z_H(MU_GT(IJKN)&
398     !              ,MU_GT(IJKTN),K),J)*(W_G(IJPK)-W_G(IJK))*ODY_N(J)*AXY_V(IJK) - &
399     !               AVG_Y_H(AVG_Z_H(MU_GT(IJKB),MU_GT(IJK),KM),AVG_Z_H(MU_GT(IJKBN),&
400     !               MU_GT(IJKN),KM),J)*(W_G(IJPKM)-W_G(IJKM))*ODY_N(J)*AXY_V(IJKM)
401     
402                 ENDIF  ! CARTESIAN GRID
403     !=======================================================================
404     ! JFD: END MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
405     !=======================================================================
406     !
407     !
408     !         Add the terms
409                 TAU_V_G(IJK) = SBV + SSX + SSY + SSZ
410              ELSE
411                 TAU_V_G(IJK) = ZERO
412              ENDIF
413           END DO
414           call send_recv(tau_v_g,2)
415           RETURN
416           END SUBROUTINE CALC_TAU_V_G
417     
418     !// Comments on the modifications for DMP version implementation
419     !// 001 Include header file and common declarations for parallelization
420     !// 350 Changed do loop limits: 1,ijkmax2-> ijkstart3, ijkend3
421     !// 400 Added sendrecv module and send_recv calls for COMMunication
422