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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Module name: CALC_Tau_U_g(TAU_U_g, IER)                             C
4     !  Purpose: Cross terms in the gradient of stress in U_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_U_G(TAU_U_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 compar
48           USE sendrecv
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     !                      TAU_U_g
70           DOUBLE PRECISION TAU_U_g(DIMENSION_3)
71     !
72     !                      Indices
73           INTEGER          I, J, JM, K, KM, IJK, IJKE, IPJK, IP, IMJK, IJKN,&
74                            IJKNE, IJKS, IJKSE, IPJMK, IJMK, IJKT, IJKTE,&
75                            IJKB, IJKBE, IJKM, IPJKM
76     !
77     !                      Average volume fraction
78           DOUBLE PRECISION EPGA
79     !
80     !                      Average viscosity
81           DOUBLE PRECISION EPMU_gte, EPMU_gbe, EPMUGA
82     !
83     !                      Average dW/Xdz
84           DOUBLE PRECISION dWoXdz
85     !
86     !                      Source terms (Surface)
87           DOUBLE PRECISION Sbv, Ssx, Ssy, Ssz
88     !
89     !                      Source terms (Volumetric)
90           DOUBLE PRECISION Vtzb
91     !
92     !=======================================================================
93     ! JFD: START MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
94     !=======================================================================
95           DOUBLE PRECISION :: DEL_H,Nx,Ny,Nz
96           LOGICAL :: V_NODE_AT_NE,V_NODE_AT_NW,V_NODE_AT_SE,V_NODE_AT_SW
97           LOGICAL :: W_NODE_AT_TE,W_NODE_AT_TW,W_NODE_AT_BE,W_NODE_AT_BW
98           DOUBLE PRECISION :: dvdx_at_N,dvdx_at_S
99           DOUBLE PRECISION :: dwdx_at_T,dwdx_at_B
100           DOUBLE PRECISION :: Xi,Yi,Zi,Vi,Wi,Sx,Sy,Sz
101           DOUBLE PRECISION :: MU_GT_CUT,SSY_CUT,SSZ_CUT
102           DOUBLE PRECISION :: UW_g,VW_g,WW_g
103           INTEGER :: BCV
104           CHARACTER(LEN=9) :: BCT
105     !=======================================================================
106     ! JFD: END MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
107     !=======================================================================
108     
109     !$omp  parallel do default(none) &
110     !$omp              private( IJK, I, IJKE, EPGA, IP, J, JM, K, KM,  &
111     !$omp              IPJK,IMJK,IJKN,IJKNE,IJKS,IJKSE,IPJMK,IJMK,IJKT,IJKTE,  &
112     !$omp              IJKB,IJKBE,IJKM,IPJKM, &
113     !$omp              SBV,  SSX,SSY,  EPMU_GTE,EPMU_GBE, SSZ, &
114     !$omp              EPMUGA,DWOXDZ,VTZB,BCV,NOC_UG, uw_g, vw_g, ww_g, cut_tau_ug, bct,mu_gt_cut,del_h,nx,ny,nz,v_node_at_ne, v_node_at_sw,v_node_at_se, v_node_at_nw, w_node_at_bw, w_node_at_be,w_node_at_tw,w_node_at_te,xi, vi,yi,zi,sx,sy,sz,dvdx_at_S,dvdx_at_N,ssy_cut,dwdx_at_B,dwdx_at_T,wi,ssz_cut) &
115     !$omp  shared(ijkstart3, ijkend3, do_k, bc_type, cut_u_cell_at, bc_u_id, tau_u_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)
116           DO IJK = IJKSTART3, IJKEND3
117              I = I_OF(IJK)
118              IJKE = EAST_OF(IJK)
119              EPGA = AVG_X(EP_G(IJK),EP_G(IJKE),I)
120              IF ( .NOT.IP_AT_E(IJK) .AND. EPGA>DIL_EP_S) THEN
121                 IP = IP1(I)
122                 J = J_OF(IJK)
123                 JM = JM1(J)
124                 K = K_OF(IJK)
125                 KM = KM1(K)
126                 IPJK = IP_OF(IJK)
127                 IMJK = IM_OF(IJK)
128                 IJKN = NORTH_OF(IJK)
129                 IJKNE = EAST_OF(IJKN)
130                 IJKS = SOUTH_OF(IJK)
131                 IJKSE = EAST_OF(IJKS)
132                 IPJMK = JM_OF(IPJK)
133                 IJMK = JM_OF(IJK)
134                 IJKT = TOP_OF(IJK)
135                 IJKTE = EAST_OF(IJKT)
136                 IJKB = BOTTOM_OF(IJK)
137                 IJKBE = EAST_OF(IJKB)
138                 IJKM = KM_OF(IJK)
139                 IPJKM = IP_OF(IJKM)
140     !=======================================================================
141     ! JFD: START MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
142     !=======================================================================
143                 IF((.NOT.CARTESIAN_GRID).OR.(CG_SAFE_MODE(3)==1)) THEN
144     !
145     !       Surface forces
146     !
147     !         bulk viscosity term
148                    SBV = (LAMBDA_GT(IJKE)*TRD_G(IJKE)-LAMBDA_GT(IJK)*TRD_G(IJK))*AYZ(&
149                       IJK)
150     !
151     !         shear stress terms
152                    SSX = MU_GT(IJKE)*(U_G(IPJK)-U_G(IJK))*ODX(IP)*AYZ_U(IJK) - MU_GT(&
153                       IJK)*(U_G(IJK)-U_G(IMJK))*ODX(I)*AYZ_U(IMJK)
154                    SSY = AVG_X_H(AVG_Y_H(MU_GT(IJK),MU_GT(IJKN),J),AVG_Y_H(MU_GT(IJKE)&
155                       ,MU_GT(IJKNE),J),I)*(V_G(IPJK)-V_G(IJK))*ODX_E(I)*AXZ_U(IJK) - &
156                       AVG_X_H(AVG_Y_H(MU_GT(IJKS),MU_GT(IJK),JM),AVG_Y_H(MU_GT(IJKSE),&
157                       MU_GT(IJKE),JM),I)*(V_G(IPJMK)-V_G(IJMK))*ODX_E(I)*AXZ_U(IJMK)
158                    EPMU_GTE = AVG_X_H(AVG_Z_H(MU_GT(IJK),MU_GT(IJKT),K),AVG_Z_H(MU_GT(&
159                       IJKE),MU_GT(IJKTE),K),I)
160                    EPMU_GBE = AVG_X_H(AVG_Z_H(MU_GT(IJKB),MU_GT(IJK),KM),AVG_Z_H(MU_GT&
161                       (IJKBE),MU_GT(IJKE),KM),I)
162                    SSZ = EPMU_GTE*(W_G(IPJK)-W_G(IJK))*ODX_E(I)*AXY_U(IJK) - EPMU_GBE*&
163                       (W_G(IPJKM)-W_G(IJKM))*ODX_E(I)*AXY_U(IJKM)
164     
165                 ELSE ! CARTESIAN GRID CASE
166     
167     !         bulk viscosity term
168                    SBV =  (LAMBDA_GT(IJKE)*TRD_G(IJKE)) * AYZ_U(IJK) &
169                          -(LAMBDA_GT(IJK) *TRD_G(IJK) ) * AYZ_U(IMJK)
170     
171     !         shear stress terms
172                    IF(.NOT.CUT_U_CELL_AT(IJK))   THEN
173     
174                       SSX = MU_GT(IJKE)*(U_G(IPJK)-U_G(IJK))*ONEoDX_E_U(IJK)*AYZ_U(IJK) - MU_GT(&
175                       IJK)*(U_G(IJK)-U_G(IMJK))*ONEoDX_E_U(IMJK)*AYZ_U(IMJK)
176     
177                       SSY = AVG_X_H(AVG_Y_H(MU_GT(IJK),MU_GT(IJKN),J),AVG_Y_H(MU_GT(IJKE)&
178                       ,MU_GT(IJKNE),J),I)*(V_G(IPJK)-V_G(IJK))*ONEoDX_E_V(IJK)*AXZ_U(IJK) - &
179                       AVG_X_H(AVG_Y_H(MU_GT(IJKS),MU_GT(IJK),JM),AVG_Y_H(MU_GT(IJKSE),&
180                       MU_GT(IJKE),JM),I)*(V_G(IPJMK)-V_G(IJMK))*ONEoDX_E_V(IJMK)*AXZ_U(IJMK)
181     
182                       IF(DO_K) THEN
183                          EPMU_GTE = AVG_X_H(AVG_Z_H(MU_GT(IJK),MU_GT(IJKT),K),AVG_Z_H(MU_GT(&
184                          IJKE),MU_GT(IJKTE),K),I)
185                          EPMU_GBE = AVG_X_H(AVG_Z_H(MU_GT(IJKB),MU_GT(IJK),KM),AVG_Z_H(MU_GT&
186                          (IJKBE),MU_GT(IJKE),KM),I)
187     
188                          SSZ = EPMU_GTE*(W_G(IPJK)-W_G(IJK))*ONEoDX_E_W(IJK)*AXY_U(IJK) - EPMU_GBE*&
189                          (W_G(IPJKM)-W_G(IJKM))*ONEoDX_E_W(IJKM)*AXY_U(IJKM)
190     
191                       ELSE
192                          SSZ = ZERO
193                       ENDIF
194     
195                    ELSE
196     
197                       BCV = BC_U_ID(IJK)
198     
199                       IF(BCV > 0 ) THEN
200                          BCT = BC_TYPE(BCV)
201                       ELSE
202                          BCT = 'NONE'
203                       ENDIF
204     
205                       SELECT CASE (BCT)
206     
207                          CASE ('CG_NSW')
208                             CUT_TAU_UG = .TRUE.
209                             NOC_UG     = .TRUE.
210                             UW_g = ZERO
211                             VW_g = ZERO
212                             WW_g = ZERO
213                          CASE ('CG_FSW')
214                             CUT_TAU_UG = .FALSE.
215                             NOC_UG     = .FALSE.
216                             UW_g = ZERO
217                             VW_g = ZERO
218                             WW_g = ZERO
219                          CASE('CG_PSW')
220                             IF(BC_HW_G(BC_U_ID(IJK))==UNDEFINED) THEN   ! same as NSW
221                                CUT_TAU_UG = .TRUE.
222                                NOC_UG     = .TRUE.
223                                UW_g = BC_UW_G(BCV)
224                                VW_g = BC_VW_G(BCV)
225                                WW_g = BC_WW_G(BCV)
226                             ELSEIF(BC_HW_G(BC_U_ID(IJK))==ZERO) THEN   ! same as FSW
227                                CUT_TAU_UG = .FALSE.
228                                NOC_UG     = .FALSE.
229                             ELSE                              ! partial slip
230                                CUT_TAU_UG = .FALSE.
231                                NOC_UG     = .FALSE.
232                                UW_g = ZERO
233                                VW_g = ZERO
234                                WW_g = ZERO
235                             ENDIF
236                          CASE ('NONE')
237                             TAU_U_G(IJK) = ZERO
238                             CYCLE
239                       END SELECT
240     
241                       IF(CUT_TAU_UG) THEN
242                          MU_GT_CUT =  (VOL(IJK)*MU_GT(IJK) + VOL(IPJK)*MU_GT(IJKE))/(VOL(IJK) + VOL(IPJK))
243                       ELSE
244                          MU_GT_CUT = ZERO
245                       ENDIF
246     
247     !           SSX:
248     
249                       CALL GET_DEL_H(IJK,'U_MOMENTUM',X_U(IJK),Y_U(IJK),Z_U(IJK),Del_H,Nx,Ny,Nz)
250     
251                       SSX = MU_GT(IJKE)*(U_G(IPJK)-U_G(IJK))*ONEoDX_E_U(IJK)*AYZ_U(IJK) - MU_GT(&
252                             IJK)*(U_G(IJK)-U_G(IMJK))*ONEoDX_E_U(IMJK)*AYZ_U(IMJK) &
253                           - MU_GT_CUT * (U_g(IJK) - UW_g) / DEL_H * (Nx**2) * Area_U_CUT(IJK)
254     
255     !           SSY:
256     
257                       V_NODE_AT_NE = ((.NOT.BLOCKED_V_CELL_AT(IPJK)).AND.(.NOT.WALL_V_AT(IPJK)))
258                       V_NODE_AT_NW = ((.NOT.BLOCKED_V_CELL_AT(IJK)).AND.(.NOT.WALL_V_AT(IJK)))
259                       V_NODE_AT_SE = ((.NOT.BLOCKED_V_CELL_AT(IPJMK)).AND.(.NOT.WALL_V_AT(IPJMK)))
260                       V_NODE_AT_SW = ((.NOT.BLOCKED_V_CELL_AT(IJMK)).AND.(.NOT.WALL_V_AT(IJMK)))
261     
262                       IF(V_NODE_AT_NE.AND.V_NODE_AT_NW) THEN
263     
264                          Vi = HALF * (V_G(IPJK) + V_G(IJK))
265                          Xi = HALF * (X_V(IPJK) + X_V(IJK))
266                          Yi = HALF * (Y_V(IPJK) + Y_V(IJK))
267                          Zi = HALF * (Z_V(IPJK) + Z_V(IJK))
268                          Sx = X_V(IPJK) - X_V(IJK)
269                          Sy = Y_V(IPJK) - Y_V(IJK)
270                          Sz = Z_V(IPJK) - Z_V(IJK)
271     
272                          CALL GET_DEL_H(IJK,'U_MOMENTUM',Xi,Yi,Zi,Del_H,Nx,Ny,Nz)
273     
274                          dvdx_at_N =  (V_G(IPJK) - V_G(IJK)) * ONEoDX_E_V(IJK)
275     
276                          IF(NOC_UG) dvdx_at_N = dvdx_at_N - ((Vi-VW_g)*ONEoDX_E_V(IJK)/DEL_H*(Sy*Ny+Sz*Nz))
277     
278                       ELSE
279                          dvdx_at_N =  ZERO
280                       ENDIF
281     
282                       IF(V_NODE_AT_SE.AND.V_NODE_AT_SW) THEN
283     
284                          Vi = HALF * (V_G(IPJMK) + V_G(IJMK))
285                          Xi = HALF * (X_V(IPJMK) + X_V(IJMK))
286                          Yi = HALF * (Y_V(IPJMK) + Y_V(IJMK))
287                          Zi = HALF * (Z_V(IPJMK) + Z_V(IJMK))
288                          Sx = X_V(IPJMK) - X_V(IJMK)
289                          Sy = Y_V(IPJMK) - Y_V(IJMK)
290                          Sz = Z_V(IPJMK) - Z_V(IJMK)
291     !
292                          CALL GET_DEL_H(IJK,'U_MOMENTUM',Xi,Yi,Zi,Del_H,Nx,Ny,Nz)
293     
294                          dvdx_at_S =  (V_G(IPJMK) - V_G(IJMK)) * ONEoDX_E_V(IJMK)
295     
296                          IF(NOC_UG) dvdx_at_S = dvdx_at_S - ((Vi-VW_g)*ONEoDX_E_V(IJMK)/DEL_H*(Sy*Ny+Sz*Nz))
297     
298                       ELSE
299                          dvdx_at_S =  ZERO
300                       ENDIF
301     
302                         IF(V_NODE_AT_NW) THEN
303                            CALL GET_DEL_H(IJK,'U_MOMENTUM',X_V(IJK),Y_V(IJK),Z_V(IJK),Del_H,Nx,Ny,Nz)
304                            SSY_CUT = - MU_GT_CUT * (V_G(IJK) - VW_g) / DEL_H * (Nx*Ny) * Area_U_CUT(IJK)
305                         ELSE
306                            SSY_CUT =  ZERO
307                         ENDIF
308     
309     
310                       SSY = AVG_X_H(AVG_Y_H(MU_GT(IJK),MU_GT(IJKN),J),AVG_Y_H(MU_GT(IJKE)&
311                          ,MU_GT(IJKNE),J),I)*dvdx_at_N*AXZ_U(IJK) - &
312                          AVG_X_H(AVG_Y_H(MU_GT(IJKS),MU_GT(IJK),JM),AVG_Y_H(MU_GT(IJKSE),&
313                          MU_GT(IJKE),JM),I)*dvdx_at_S*AXZ_U(IJMK) &
314                         + SSY_CUT
315     
316     !           SSZ:
317     
318                       IF(DO_K) THEN
319     
320                          W_NODE_AT_TE = ((.NOT.BLOCKED_W_CELL_AT(IPJK)).AND.(.NOT.WALL_W_AT(IPJK)))
321                          W_NODE_AT_TW = ((.NOT.BLOCKED_W_CELL_AT(IJK)).AND.(.NOT.WALL_W_AT(IJK)))
322                          W_NODE_AT_BE = ((.NOT.BLOCKED_W_CELL_AT(IPJKM)).AND.(.NOT.WALL_W_AT(IPJKM)))
323                          W_NODE_AT_BW = ((.NOT.BLOCKED_W_CELL_AT(IJKM)).AND.(.NOT.WALL_W_AT(IJKM)))
324     
325                          IF(W_NODE_AT_TE.AND.W_NODE_AT_TW) THEN
326     
327                             Wi = HALF * (W_G(IPJK) + W_G(IJK))
328                             Xi = HALF * (X_W(IPJK) + X_W(IJK))
329                             Yi = HALF * (Y_W(IPJK) + Y_W(IJK))
330                             Zi = HALF * (Z_W(IPJK) + Z_W(IJK))
331                             Sx = X_W(IPJK) - X_W(IJK)
332                             Sy = Y_W(IPJK) - Y_W(IJK)
333                             Sz = Z_W(IPJK) - Z_W(IJK)
334     
335                             CALL GET_DEL_H(IJK,'U_MOMENTUM',Xi,Yi,Zi,Del_H,Nx,Ny,Nz)
336     
337                             dwdx_at_T =  (W_G(IPJK) - W_G(IJK)) * ONEoDX_E_W(IJK)
338     
339                             IF(NOC_UG) dwdx_at_T = dwdx_at_T - ((Wi-WW_g) * ONEoDX_E_W(IJK)/DEL_H*(Sy*Ny+Sz*Nz))
340     
341                          ELSE
342                             dwdx_at_T =  ZERO
343                          ENDIF
344     
345     
346                          IF(W_NODE_AT_BE.AND.W_NODE_AT_BW) THEN
347     
348                             Wi = HALF * (W_G(IPJKM) + W_G(IJKM))
349                             Xi = HALF * (X_W(IPJKM) + X_W(IJKM))
350                             Yi = HALF * (Y_W(IPJKM) + Y_W(IJKM))
351                             Zi = HALF * (Z_W(IPJKM) + Z_W(IJKM))
352                             Sx = X_W(IPJKM) - X_W(IJKM)
353                             Sy = Y_W(IPJKM) - Y_W(IJKM)
354                             Sz = Z_W(IPJKM) - Z_W(IJKM)
355     
356                             CALL GET_DEL_H(IJK,'U_MOMENTUM',Xi,Yi,Zi,Del_H,Nx,Ny,Nz)
357     
358                             dwdx_at_B =  (W_G(IPJKM) - W_G(IJKM)) * ONEoDX_E_W(IJKM)
359     
360                             IF(NOC_UG) dwdx_at_B = dwdx_at_B - ((Wi-WW_g) * ONEoDX_E_W(IJKM)/DEL_H*(Sy*Ny+Sz*Nz))
361     
362                          ELSE
363                             dwdx_at_B =  ZERO
364                          ENDIF
365     
366                          IF(W_NODE_AT_TW) THEN
367                             CALL GET_DEL_H(IJK,'U_MOMENTUM',X_W(IJK),Y_W(IJK),Z_W(IJK),Del_H,Nx,Ny,Nz)
368                             SSZ_CUT = - MU_GT_CUT * (W_G(IJK) - WW_g) / DEL_H * (Nx*Nz) * Area_U_CUT(IJK)
369                          ELSE
370                             SSZ_CUT =  ZERO
371                          ENDIF
372     
373                          EPMU_GTE = AVG_X_H(AVG_Z_H(MU_GT(IJK),MU_GT(IJKT),K),AVG_Z_H(MU_GT(&
374                             IJKE),MU_GT(IJKTE),K),I)
375                          EPMU_GBE = AVG_X_H(AVG_Z_H(MU_GT(IJKB),MU_GT(IJK),KM),AVG_Z_H(MU_GT&
376                             (IJKBE),MU_GT(IJKE),KM),I)
377                          SSZ =   EPMU_GTE*dwdx_at_T*AXY_U(IJK)  &
378                                - EPMU_GBE*dwdx_at_B*AXY_U(IJKM) &
379                                + SSZ_CUT
380     
381                       ELSE
382     
383                         SSZ = ZERO
384     
385                       ENDIF ! DO_K
386     
387                    ENDIF  ! CUT_CELL
388     
389     !  Original terms
390     !            SSX = MU_GT(IJKE)*(U_G(IPJK)-U_G(IJK))*ODX(IP)*AYZ_U(IJK) - MU_GT(&
391     !               IJK)*(U_G(IJK)-U_G(IMJK))*ODX(I)*AYZ_U(IMJK)
392     
393     !            SSY = AVG_X_H(AVG_Y_H(MU_GT(IJK),MU_GT(IJKN),J),AVG_Y_H(MU_GT(IJKE)&
394     !               ,MU_GT(IJKNE),J),I)*(V_G(IPJK)-V_G(IJK))*ODX_E(I)*AXZ_U(IJK) - &
395     !               AVG_X_H(AVG_Y_H(MU_GT(IJKS),MU_GT(IJK),JM),AVG_Y_H(MU_GT(IJKSE),&
396     !               MU_GT(IJKE),JM),I)*(V_G(IPJMK)-V_G(IJMK))*ODX_E(I)*AXZ_U(IJMK)
397     
398     
399     !            EPMU_GTE = AVG_X_H(AVG_Z_H(MU_GT(IJK),MU_GT(IJKT),K),AVG_Z_H(MU_GT(&
400     !               IJKE),MU_GT(IJKTE),K),I)
401     !            EPMU_GBE = AVG_X_H(AVG_Z_H(MU_GT(IJKB),MU_GT(IJK),KM),AVG_Z_H(MU_GT&
402     !               (IJKBE),MU_GT(IJKE),KM),I)
403     !            SSZ = EPMU_GTE*(W_G(IPJK)-W_G(IJK))*ODX_E(I)*AXY_U(IJK) - EPMU_GBE*&
404     !               (W_G(IPJKM)-W_G(IJKM))*ODX_E(I)*AXY_U(IJKM)
405     
406                 ENDIF ! CARTESIAN GRID
407     !=======================================================================
408     ! JFD: END MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
409     !=======================================================================
410     !
411     !         Special terms for cylindrical coordinates
412                 IF (CYLINDRICAL) THEN
413     !
414     !           modify Ssz: integral of (1/x) (d/dz) (mu*(-w/x))
415                    SSZ = SSZ - (EPMU_GTE*(HALF*(W_G(IPJK)+W_G(IJK))*OX_E(I))*AXY_U(&
416                       IJK)-EPMU_GBE*(HALF*(W_G(IPJKM)+W_G(IJKM))*OX_E(I))*AXY_U(&
417                       IJKM))
418     !
419     !           -(2mu/x)*(1/x)*dw/dz part of Tau_zz/X
420                    EPMUGA = AVG_X(MU_G(IJK),MU_G(IJKE),I)
421                    DWOXDZ = HALF*((W_G(IJK)-W_G(IJKM))*OX(I)*ODZ(K)+(W_G(IPJK)-W_G(&
422                       IPJKM))*OX(IP)*ODZ(K))
423                    VTZB = -2.d0*EPMUGA*OX_E(I)*DWOXDZ
424                 ELSE
425                    VTZB = ZERO
426                 ENDIF
427     !
428     !         Add the terms
429                 TAU_U_G(IJK) = SBV + SSX + SSY + SSZ + VTZB*VOL_U(IJK)
430              ELSE
431                 TAU_U_G(IJK) = ZERO
432              ENDIF
433           END DO
434     
435           call send_recv(tau_u_g,2)
436     
437           RETURN
438     
439           END SUBROUTINE CALC_TAU_U_G
440     
441     !// Comments on the modifications for DMP version implementation
442     !// 001 Include header file and common declarations for parallelization
443     !// 350 Changed do loop limits: 1,ijkmax2-> ijkstart3, ijkend3
444