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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Module name: CONV_DIF_W_g(A_m, B_m, IER)                            C
4     !  Purpose: Determine convection diffusion terms for W_g momentum eqs  C
5     !  The off-diagonal coefficients calculated here must be positive. The C
6     !  center coefficient and the source vector are negative;              C
7     !  See source_w_g                                                      C
8     !                                                                      C
9     !  Author: M. Syamlal                                 Date: 24-DEC-96  C
10     !  Reviewer:                                          Date:            C
11     !                                                                      C
12     !                                                                      C
13     !  Literature/Document References:                                     C
14     !                                                                      C
15     !  Variables referenced:                                               C
16     !  Variables modified:                                                 C
17     !                                                                      C
18     !  Local variables:                                                    C
19     !                                                                      C
20     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
21     !
22           SUBROUTINE CONV_DIF_W_G(A_M, B_M, IER)
23     !...Translated by Pacific-Sierra Research VAST-90 2.06G5  12:17:31  12/09/98
24     !...Switches: -xf
25     !
26     !  Include param.inc file to specify parameter values
27     !
28     !-----------------------------------------------
29     !   M o d u l e s
30     !-----------------------------------------------
31           USE param
32           USE param1
33           USE parallel
34           USE matrix
35           USE geometry
36           USE indices
37           USE run
38           USE visc_g
39           USE compar
40     
41           IMPLICIT NONE
42     !-----------------------------------------------
43     !   G l o b a l   P a r a m e t e r s
44     !-----------------------------------------------
45     !-----------------------------------------------
46     !   D u m m y   A r g u m e n t s
47     !-----------------------------------------------
48     !
49     !
50     !
51     !                      Error index
52           INTEGER          IER
53     !
54     !                      Septadiagonal matrix A_m
55           DOUBLE PRECISION A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
56     !
57     !                      Vector b_m
58           DOUBLE PRECISION B_m(DIMENSION_3, 0:DIMENSION_M)
59     
60     !
61           IF (.NOT.MOMENTUM_Z_EQ(0)) RETURN
62     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
63     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
64     !       IF DEFERRED CORRECTION IS USED TO SOLVE W_G
65           IF (DEF_COR) THEN
66             CALL STORE_A_W_G0 (A_M(1,-3,0), IER)
67             IF (DISCRETIZE(5) > 1)CALL STORE_A_W_GDC (A_M(1,-3,0), B_M(1,0))
68           ELSE
69     !
70             IF (DISCRETIZE(5) == 0) THEN               ! 0 & 1 => FOUP
71               CALL STORE_A_W_G0 (A_M(1,-3,0), IER)
72             ELSE
73               CALL STORE_A_W_G1 (A_M(1,-3,0))
74             ENDIF
75           ENDIF
76     !
77           CALL DIF_W_IS (MU_GT, A_M, B_M, 0, IER)
78     !
79     
80           RETURN
81           END SUBROUTINE CONV_DIF_W_G
82     !
83     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
84     !                                                                      C
85     !  Module name: STORE_A_W_g0(A_W_g, IER)                               C
86     !  Purpose: Determine convection diffusion terms for W_g momentum eqs  C
87     !  The off-diagonal coefficients calculated here must be positive. The C
88     !  center coefficient and the source vector are negative; FOUP         C
89     !  See source_w_g                                                      C
90     !                                                                      C
91     !  Author: M. Syamlal                                 Date: 7-JUN-96   C
92     !  Reviewer:                                          Date:            C
93     !                                                                      C
94     !  Revision Number: 1                                                  C
95     !  Purpose: To incorporate Cartesian grid modifications                C
96     !  Author: Jeff Dietiker                              Date: 01-Jul-09  C
97     !                                                                      C
98     !  Literature/Document References:                                     C
99     !                                                                      C
100     !  Variables referenced:                                               C
101     !  Variables modified:                                                 C
102     !                                                                      C
103     !  Local variables:                                                    C
104     !                                                                      C
105     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
106     !
107           SUBROUTINE STORE_A_W_G0(A_W_G, IER)
108     !...Translated by Pacific-Sierra Research VAST-90 2.06G5  12:17:31  12/09/98
109     !...Switches: -xf
110     !
111     !  Include param.inc file to specify parameter values
112     !
113     !-----------------------------------------------
114     !   M o d u l e s
115     !-----------------------------------------------
116           USE param
117           USE param1
118           USE parallel
119           USE matrix
120           USE geometry
121           USE indices
122           USE run
123           USE visc_g
124           USE toleranc
125           USE physprop
126           USE fldvar
127           USE output
128           USE compar
129           USE mflux
130           USE fun_avg
131           USE functions
132     !=======================================================================
133     ! JFD: START MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
134     !=======================================================================
135           USE cutcell
136     !=======================================================================
137     ! JFD: END MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
138     !=======================================================================
139           IMPLICIT NONE
140     !-----------------------------------------------
141     !   G l o b a l   P a r a m e t e r s
142     !-----------------------------------------------
143     !-----------------------------------------------
144     !   D u m m y   A r g u m e n t s
145     !-----------------------------------------------
146     !
147     !
148     !
149     !                      Error index
150           INTEGER          IER
151     !
152     !                      Indices
153           INTEGER          I,  J, K, IPJK, IJPK, IJKN, IJKC, KP, IJKE,&
154                            IJKTE, IJKP, IJKT, IJKTN, IJK
155           INTEGER          IMJK, IM, IJKW, IJKWT, IMJKP
156           INTEGER          IJMK, JM, IJMKP, IJKS, IJKST
157           INTEGER          IJKM, KM, IJKB
158     !
159     !                      Solids phase
160           INTEGER          M
161     !
162     !                      Face mass flux
163           DOUBLE PRECISION Flux
164     !
165     !                      Diffusion parameter
166           DOUBLE PRECISION D_f
167     !
168     !                      Septadiagonal matrix A_W_g
169           DOUBLE PRECISION A_W_g(DIMENSION_3, -3:3)
170     !=======================================================================
171     ! JFD: START MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
172     !=======================================================================
173           DOUBLE PRECISION :: AW,HW,VELW
174     !=======================================================================
175     ! JFD: END MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
176     !=======================================================================
177     
178     !  Calculate convection-diffusion fluxes through each of the faces
179     !
180     !     Fluid phase
181           M = 0
182     
183     !$omp     parallel do default(none)                                 &
184     !$omp     private( I,  J, K, IPJK, IJPK, IJKN, IJKC, KP,     &
185     !$omp             IJKE, IJKTE, IJKP, IJKT, IJKTN, IJK, D_f,    &
186     !$omp             IMJK, IM, IJKW, IJKWT, IMJKP,                     &
187     !$omp             IJMK, JM, IJMKP, IJKS, IJKST,                     &
188     !$omp             IJKM, KM, IJKB, FLUX,VELW,HW,AW) &
189     !$omp     shared (ijkstart3,ijkend3,i_of,j_of,k_of,kp1,cut_w_treatment_at,theta_vn_bar,flux_ge,theta_vn,alpha_ve_c,mu_gt, &
190     !$omp            oneody_n_v,ayz_v,ody,theta_v_st,flux_gn,theta_v_nt,oneody_n_u,axz_v,axy_v,ody_n,do_k,theta_u_tw,flux_gt, &
191     !$omp            theta_u_te,alpha_vt_c,oneodz_t_v,odz_t,im1,jm1,km1,a_w_g,theta_v_ne,oneodx_e_v,alpha_vn_c,theta_v_se,    &
192     !$omp            odx_e,ox,theta_w_be,theta_w_te,alpha_we_c,ONEODX_E_W,AYZ_W,theta_w_bn,theta_w_tn,alpha_wn_c,oneody_n_w,  &
193     !$omp            axz_w,theta_wt_bar,theta_wt,alpha_wt_c,oneodz_t_w,axy_w,odz)
194           DO IJK = ijkstart3, ijkend3
195     !
196              IF (FLOW_AT_T(IJK)) THEN
197                 I = I_OF(IJK)
198                 J = J_OF(IJK)
199                 K = K_OF(IJK)
200                 IPJK = IP_OF(IJK)
201                 IJPK = JP_OF(IJK)
202                 IJKN = NORTH_OF(IJK )
203                 IJKT = TOP_OF(IJK)
204                 IF (WALL_AT(IJK)) THEN
205                    IJKC = IJKT
206                 ELSE
207                    IJKC = IJK
208                 ENDIF
209                 KP = KP1(K)
210                 IJKE = EAST_OF(IJK)
211                 IJKP = KP_OF(IJK)
212                 IJKTN = NORTH_OF(IJKT)
213                 IJKTE = EAST_OF(IJKT)
214     !
215     !           East face (i+1/2, j, k+1/2)
216     !=======================================================================
217     ! JFD: START MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
218     !=======================================================================
219                 IF(CUT_W_TREATMENT_AT(IJK)) THEN
220                    Flux = (Theta_W_be(IJK) * Flux_gE(IJK) + Theta_W_te(IJK) * Flux_gE(IJKP))
221                    CALL GET_INTERPOLATION_TERMS_G(IJK,'W_MOMENTUM',ALPHA_We_c(IJK),AW,HW,VELW)
222                    Flux = Flux * AW
223                    D_F = AVG_Z_H(AVG_X_H(MU_GT(IJKC),MU_GT(IJKE),I),AVG_X_H(MU_GT(IJKT&
224                       ),MU_GT(IJKTE),I),K)*ONEoDX_E_W(IJK)*AYZ_W(IJK)
225                 ELSE   ! Original terms
226                    Flux = HALF * (Flux_gE(IJK) + Flux_gE(IJKP))
227                    D_F = AVG_Z_H(AVG_X_H(MU_GT(IJKC),MU_GT(IJKE),I),AVG_X_H(MU_GT(IJKT&
228                       ),MU_GT(IJKTE),I),K)*ODX_E(I)*AYZ_W(IJK)
229                 ENDIF
230     !=======================================================================
231     ! JFD: END MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
232     !=======================================================================
233                 IF (Flux >= ZERO) THEN
234                    A_W_G(IJK,E) = D_F
235                    A_W_G(IPJK,W) = D_F + Flux
236                 ELSE
237                    A_W_G(IJK,E) = D_F - Flux
238                    A_W_G(IPJK,W) = D_F
239                 ENDIF
240     !
241     !           North face (i, j+1/2, k+1/2)
242     !=======================================================================
243     ! JFD: START MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
244     !=======================================================================
245                 IF(CUT_W_TREATMENT_AT(IJK)) THEN
246                    Flux = (Theta_W_bn(IJK) * Flux_gN(IJK) + Theta_W_tn(IJK) * Flux_gN(IJKP))
247                    CALL GET_INTERPOLATION_TERMS_G(IJK,'W_MOMENTUM',ALPHA_Wn_c(IJK),AW,HW,VELW)
248                    Flux = Flux * AW
249                    D_F = AVG_Z_H(AVG_Y_H(MU_GT(IJKC),MU_GT(IJKN),J),AVG_Y_H(MU_GT(IJKT&
250                       ),MU_GT(IJKTN),J),K)*ONEoDY_N_W(IJK)*AXZ_W(IJK)
251                 ELSE   ! Original terms
252                    Flux = HALF * (Flux_gN(IJK) + Flux_gN(IJKP))
253                    D_F = AVG_Z_H(AVG_Y_H(MU_GT(IJKC),MU_GT(IJKN),J),AVG_Y_H(MU_GT(IJKT&
254                       ),MU_GT(IJKTN),J),K)*ODY_N(J)*AXZ_W(IJK)
255                 ENDIF
256     !=======================================================================
257     ! JFD: END MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
258     !=======================================================================
259                 IF (Flux >= ZERO) THEN
260                    A_W_G(IJK,N) = D_F
261                    A_W_G(IJPK,S) = D_F + Flux
262                 ELSE
263                    A_W_G(IJK,N) = D_F - Flux
264                    A_W_G(IJPK,S) = D_F
265                 ENDIF
266     !
267     !           Top face (i, j, k+1)
268     !=======================================================================
269     ! JFD: START MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
270     !=======================================================================
271                 IF(CUT_W_TREATMENT_AT(IJK)) THEN
272                    Flux = (Theta_Wt_bar(IJK) * Flux_gT(IJK) + Theta_Wt(IJK) * Flux_gT(IJKP))
273                    CALL GET_INTERPOLATION_TERMS_G(IJK,'W_MOMENTUM',alpha_Wt_c(IJK),AW,HW,VELW)
274                    Flux = Flux * AW
275                    D_F = MU_GT(IJKT)*ONEoDZ_T_W(IJK)*AXY_W(IJK)
276                 ELSE   ! Original terms
277                    Flux = HALF * (Flux_gT(IJK) + Flux_gT(IJKP))
278                    D_F = MU_GT(IJKT)*OX(I)*ODZ(KP)*AXY_W(IJK)
279                 ENDIF
280     !=======================================================================
281     ! JFD: END MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
282     !=======================================================================
283                 IF (Flux >= ZERO) THEN
284                    A_W_G(IJK,T) = D_F
285                    A_W_G(IJKP,B) = D_F + Flux
286                 ELSE
287                    A_W_G(IJK,T) = D_F - Flux
288                    A_W_G(IJKP,B) = D_F
289                 ENDIF
290     !
291     !           West face (i-1/2, j, k+1/2)
292                 IMJK = IM_OF(IJK)
293                 IF (.NOT.FLOW_AT_T(IMJK)) THEN
294                    IM = IM1(I)
295                    IJKW = WEST_OF(IJK)
296                    IJKWT = TOP_OF(IJKW)
297                    IMJKP = KP_OF(IMJK)
298     !=======================================================================
299     ! JFD: START MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
300     !=======================================================================
301                    IF(CUT_W_TREATMENT_AT(IJK)) THEN
302                       Flux = (Theta_W_be(IMJK) * Flux_gE(IMJK) + Theta_W_te(IMJK) * Flux_gE(IMJKP))
303                       CALL GET_INTERPOLATION_TERMS_G(IJK,'W_MOMENTUM',ALPHA_We_c(IMJK),AW,HW,VELW)
304                       Flux = Flux * AW
305                       D_F = AVG_Z_H(AVG_X_H(MU_GT(IJKW),MU_GT(IJKC),IM),AVG_X_H(MU_GT(&
306                           IJKWT),MU_GT(IJKT),IM),K)*ONEoDX_E_W(IMJK)*AYZ_W(IMJK)
307                     ELSE   ! Original terms
308                        Flux = HALF * (Flux_gE(IMJK) + Flux_gE(IMJKP))
309                        D_F = AVG_Z_H(AVG_X_H(MU_GT(IJKW),MU_GT(IJKC),IM),AVG_X_H(MU_GT(&
310                        IJKWT),MU_GT(IJKT),IM),K)*ODX_E(IM)*AYZ_W(IMJK)
311                    ENDIF
312     !=======================================================================
313     ! JFD: END MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
314     !=======================================================================
315                    IF (Flux >= ZERO) THEN
316                       A_W_G(IJK,W) = D_F + Flux
317                    ELSE
318                       A_W_G(IJK,W) = D_F
319                    ENDIF
320                 ENDIF
321     !
322     !           South face (i, j-1/2, k+1/2)
323                 IJMK = JM_OF(IJK)
324                 IF (.NOT.FLOW_AT_T(IJMK)) THEN
325                    JM = JM1(J)
326                    IJMKP = KP_OF(IJMK)
327                    IJKS = SOUTH_OF(IJK)
328                    IJKST = TOP_OF(IJKS)
329     !=======================================================================
330     ! JFD: START MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
331     !=======================================================================
332                    IF(CUT_W_TREATMENT_AT(IJK)) THEN
333                       Flux = (Theta_W_bn(IJMK) * Flux_gN(IJMK) + Theta_W_tn(IJMK) * Flux_gN(IJMKP))
334                       CALL GET_INTERPOLATION_TERMS_G(IJK,'W_MOMENTUM',ALPHA_Wn_c(IJMK),AW,HW,VELW)
335                       Flux = Flux * AW
336                       D_F = AVG_Z_H(AVG_Y_H(MU_GT(IJKS),MU_GT(IJKC),JM),AVG_Y_H(MU_GT(&
337                          IJKST),MU_GT(IJKT),JM),K)*ONEoDY_N_W(IJMK)*AXZ_W(IJMK)
338                    ELSE   ! Original terms
339                       Flux = HALF * (Flux_gN(IJMK) + Flux_gN(IJMKP))
340                       D_F = AVG_Z_H(AVG_Y_H(MU_GT(IJKS),MU_GT(IJKC),JM),AVG_Y_H(MU_GT(&
341                          IJKST),MU_GT(IJKT),JM),K)*ODY_N(JM)*AXZ_W(IJMK)
342                    ENDIF
343     !=======================================================================
344     ! JFD: END MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
345     !=======================================================================
346                    IF (Flux >= ZERO) THEN
347                       A_W_G(IJK,S) = D_F + Flux
348                    ELSE
349                       A_W_G(IJK,S) = D_F
350                    ENDIF
351                 ENDIF
352     !
353     !           Bottom face (i, j, k)
354                 IJKM = KM_OF(IJK)
355                 IF (.NOT.FLOW_AT_T(IJKM)) THEN
356                    KM = KM1(K)
357                    IJKB = BOTTOM_OF(IJK)
358     !=======================================================================
359     ! JFD: START MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
360     !=======================================================================
361                    IF(CUT_W_TREATMENT_AT(IJK)) THEN
362                       Flux = (Theta_Wt_bar(IJKM) * Flux_gT(IJKM) + Theta_Wt(IJKM) * Flux_gT(IJK))
363                       CALL GET_INTERPOLATION_TERMS_G(IJK,'W_MOMENTUM',alpha_Wt_c(IJKM),AW,HW,VELW)
364                       Flux = Flux * AW
365                       D_F = MU_GT(IJK)*ONEoDZ_T_W(IJKM)*AXY_W(IJKM)
366                    ELSE   ! Original terms
367                       Flux = HALF * (Flux_gT(IJKM) + Flux_gT(IJK))
368                       D_F = MU_GT(IJK)*OX(I)*ODZ(K)*AXY_W(IJKM)
369                    ENDIF
370     !=======================================================================
371     ! JFD: END MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
372     !=======================================================================
373                    IF (Flux >= ZERO) THEN
374                       A_W_G(IJK,B) = D_F + Flux
375                    ELSE
376                       A_W_G(IJK,B) = D_F
377                    ENDIF
378                 ENDIF
379              ENDIF
380           END DO
381     !$omp end parallel do
382     
383           RETURN
384           END SUBROUTINE STORE_A_W_G0
385     
386     !
387     !
388     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
389     !                                                                      C
390     !  Module name: STORE_A_W_gdc(A_W_g, B_M, IER)                         C
391     !  Purpose: TO USE DEFERRED CORRECTION METHOD TO SOLVE THE W-MOMENTUM  C
392     !  EQUATION. THIS METHOD COMBINES FIRST ORDER UPWIND AND A USER        C
393     !  SPECIFIED HIGH ORDER METHOD                                         C
394     !                                                                      C
395     !  Author: C. GUENTHER                                 Date:8-APR-99   C
396     !  Reviewer:                                          Date:            C
397     !                                                                      C
398     !  Revision Number: 1                                                  C
399     !  Purpose: To incorporate Cartesian grid modifications                C
400     !  Author: Jeff Dietiker                              Date: 01-Jul-09  C
401     !                                                                      C
402     !  Literature/Document References:                                     C
403     !                                                                      C
404     !  Variables referenced:                                               C
405     !  Variables modified:                                                 C
406     !                                                                      C
407     !  Local variables:                                                    C
408     !                                                                      C
409     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
410     !
411           SUBROUTINE STORE_A_W_GDC(A_W_G, B_M)
412     !...Translated by Pacific-Sierra Research VAST-90 2.06G5  12:17:31  12/09/98
413     !...Switches: -xf
414     !
415     !  Include param.inc file to specify parameter values
416     !
417     !-----------------------------------------------
418     !   M o d u l e s
419     !-----------------------------------------------
420           USE compar
421           USE discretization, ONLY: fpfoi_of
422           USE fldvar
423           USE fun_avg
424           USE function3
425           USE functions
426           USE geometry
427           USE indices
428           USE matrix
429           USE mflux
430           USE output
431           USE parallel
432           USE param
433           USE param1
434           USE physprop
435           USE run
436           USE sendrecv
437           USE sendrecv3
438           USE toleranc
439           USE visc_g
440           USE xsi
441           USE xsi_array
442           Use tmp_array,  U => Array1, V => Array2, WW => Array3
443     !=======================================================================
444     ! JFD: START MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
445     !=======================================================================
446           USE cutcell
447     !=======================================================================
448     ! JFD: END MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
449     !=======================================================================
450           IMPLICIT NONE
451     !-----------------------------------------------
452     !   G l o b a l   P a r a m e t e r s
453     !-----------------------------------------------
454     !-----------------------------------------------
455     !   D u m m y   A r g u m e n t s
456     !-----------------------------------------------
457     !
458     !                      Indices
459           INTEGER          I,  J, K, IPJK, IJPK, IJKN, IJKC, KP, IJKE,&
460                            IJKTE, IJKP, IJKT, IJKTN, IJK
461           INTEGER          IMJK, IM, IJKW, IJKWT, IMJKP
462           INTEGER          IJMK, JM, IJMKP, IJKS, IJKST
463           INTEGER          IJKM, KM, IJKB
464           INTEGER          IJK4, IPPP, IPPP4, JPPP, JPPP4, KPPP, KPPP4
465           INTEGER          IMMM, IMMM4, JMMM, JMMM4, KMMM, KMMM4
466     !
467     ! loezos
468             INTEGER incr
469     ! loezos
470     
471     !                      Septadiagonal matrix A_W_g
472           DOUBLE PRECISION A_W_g(DIMENSION_3, -3:3)
473     !
474     !                      Vector b_m
475           DOUBLE PRECISION B_m(DIMENSION_3)
476     !
477     !       DEFERRED CORRCTION CONTRIBUTION FORM HIGH ORDER METHOD
478             DOUBLE PRECISION MOM_HO
479     !
480     !       LOW ORDER APPROXIMATION
481             DOUBLE PRECISION MOM_LO
482     !
483     !       CONVECTION FACTOR AT THE FACE
484             DOUBLE PRECISION Flux
485     !
486     !       DEFERRED CORRECTION CONTRIBUTIONS FROM EACH FACE
487             DOUBLE PRECISION        EAST_DC
488             DOUBLE PRECISION        WEST_DC
489             DOUBLE PRECISION        NORTH_DC
490             DOUBLE PRECISION        SOUTH_DC
491             DOUBLE PRECISION  TOP_DC
492             DOUBLE PRECISION  BOTTOM_DC
493     !=======================================================================
494     ! JFD: START MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
495     !=======================================================================
496           DOUBLE PRECISION :: AW,HW,VELW
497     !=======================================================================
498     ! JFD: END MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
499     !=======================================================================
500     !
501           call lock_tmp4_array
502           call lock_tmp_array
503           call lock_xsi_array
504     !
505     !  Calculate convection factors
506     !
507     !
508     ! Send recv the third ghost layer
509           IF ( FPFOI ) THEN
510              Do IJK = ijkstart3, ijkend3
511                 I = I_OF(IJK)
512                 J = J_OF(IJK)
513                 K = K_OF(IJK)
514                 IJK4 = funijk3(I,J,K)
515                 TMP4(IJK4) = W_G(IJK)
516              ENDDO
517              CALL send_recv3(tmp4)
518           ENDIF
519     
520     
521     !!!$omp parallel do private(IJK,K,IJKT,IJKP)
522           DO IJK = ijkstart3, ijkend3
523              K = K_OF(IJK)
524              IJKT = TOP_OF(IJK)
525              IJKP = KP_OF(IJK)
526     !
527     !=======================================================================
528     ! JFD: START MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
529     !=======================================================================
530     !           East face (i+1/2, j, k+1/2)
531              IF(CUT_W_TREATMENT_AT(IJK)) THEN
532                 U(IJK) = (Theta_W_be(IJK) * U_G(IJK) + Theta_W_te(IJK) * U_G(IJKP))
533                 CALL GET_INTERPOLATION_TERMS_G(IJK,'W_MOMENTUM',ALPHA_We_c(IJK),AW,HW,VELW)
534                 U(IJK) = U(IJK) * AW
535              ELSE   ! Original terms
536                 U(IJK) = AVG_Z(U_G(IJK),U_G(IJKP),K)
537              ENDIF
538     !
539     !
540     !           North face (i, j+1/2, k+1/2)
541              IF(CUT_W_TREATMENT_AT(IJK)) THEN
542                 V(IJK) = (Theta_W_bn(IJK) * V_G(IJK) + Theta_W_tn(IJK) * V_G(IJKP))
543                 CALL GET_INTERPOLATION_TERMS_G(IJK,'W_MOMENTUM',ALPHA_Wn_c(IJK),AW,HW,VELW)
544                 V(IJK) = V(IJK) * AW
545              ELSE   ! Original terms
546                 V(IJK) = AVG_Z(V_G(IJK),V_G(IJKP),K)
547              ENDIF
548     !
549     !
550     !           Top face (i, j, k+1)
551              IF(CUT_W_TREATMENT_AT(IJK)) THEN
552                 WW(IJK) = (Theta_Wt_bar(IJK) * W_G(IJK) + Theta_Wt(IJK) * W_G(IJKP))
553                 CALL GET_INTERPOLATION_TERMS_G(IJK,'W_MOMENTUM',alpha_Wt_c(IJK),AW,HW,VELW)
554                 WW(IJK) = WW(IJK) * AW
555              ELSE   ! Original terms
556                 WW(IJK) = AVG_Z_T(W_G(IJK),W_G(IJKP))
557              ENDIF
558     !=======================================================================
559     ! JFD: END MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
560     !=======================================================================
561           END DO
562     
563     ! loezos
564             incr=0
565     ! loezos
566     
567           CALL CALC_XSI (DISCRETIZE(5), W_G, U, V, WW, XSI_E, XSI_N,&
568              XSI_T,incr)
569     !
570     !
571     !
572     !  Calculate convection-diffusion fluxes through each of the faces
573     !
574     !
575     
576     !!!$omp      parallel do                                               &
577     !!!$omp&     private( I,  J, K, IPJK, IJPK, IJKN, IJKC, KP,     &
578     !!!$omp&             IJKE, IJKTE, IJKP, IJKT, IJKTN, IJK,  D_f,        &
579     !!!$omp&             IMJK, IM, IJKW, IJKWT, IMJKP,                     &
580     !!!$omp&             IJMK, JM, IJMKP, IJKS, IJKST,                     &
581     !!!$omp&             IJKM, KM, IJKB, &
582     !!!$omp&              MOM_HO, MOM_LO, EAST_DC,WEST_DC,NORTH_DC,&
583     !!!$omp&              SOUTH_DC, TOP_DC,BOTTOM_DC)
584           DO IJK = ijkstart3, ijkend3
585     !
586              IF (FLOW_AT_T(IJK)) THEN
587                 I = I_OF(IJK)
588                 J = J_OF(IJK)
589                 K = K_OF(IJK)
590                 IPJK = IP_OF(IJK)
591                 IMJK = IM_OF(IJK)
592                 IJPK = JP_OF(IJK)
593                 IJMK = JM_OF(IJK)
594                 IJKP = KP_OF(IJK)
595                 IJKM = KM_OF(IJK)
596                 IJKN = NORTH_OF(IJK)
597                 IJKT = TOP_OF(IJK)
598                 IF (WALL_AT(IJK)) THEN
599                    IJKC = IJKT
600                 ELSE
601                    IJKC = IJK
602                 ENDIF
603                 KP = KP1(K)
604                 IJKE = EAST_OF(IJK)
605                 IJKP = KP_OF(IJK)
606                 IJKTN = NORTH_OF(IJKT)
607                 IJKTE = EAST_OF(IJKT)
608     !
609     !           Third Ghost layer information
610                 IPPP  = IP_OF(IP_OF(IPJK))
611                 IPPP4 = funijk3(I_OF(IPPP), J_OF(IPPP), K_OF(IPPP))
612                 IMMM  = IM_OF(IM_OF(IMJK))
613                 IMMM4 = funijk3(I_OF(IMMM), J_OF(IMMM), K_OF(IMMM))
614                 JPPP  = JP_OF(JP_OF(IJPK))
615                 JPPP4 = funijk3(I_OF(JPPP), J_OF(JPPP), K_OF(JPPP))
616                 JMMM  = JM_OF(JM_OF(IJMK))
617                 JMMM4 = funijk3(I_OF(JMMM), J_OF(JMMM), K_OF(JMMM))
618                 KPPP  = KP_OF(KP_OF(IJKP))
619                 KPPP4 = funijk3(K_OF(IPPP), J_OF(KPPP), K_OF(KPPP))
620                 KMMM  = KM_OF(KM_OF(IJKM))
621                 KMMM4 = funijk3(I_OF(KMMM), J_OF(KMMM), K_OF(KMMM))
622     !
623     !
624     
625     !
626     !           DEFERRED CORRECTION CONTRIBUTION AT THE East face (i+1/2, j, k+1/2)
627     !
628                 IF(U(IJK) >= ZERO)THEN
629                   MOM_LO = W_G(IJK)
630                   IF ( FPFOI ) &
631                           MOM_HO = FPFOI_OF(W_G(IPJK), W_G(IJK), &
632                                 W_G(IMJK), W_G(IM_OF(IMJK)))
633                 ELSE
634                   MOM_LO = W_G(IPJK)
635                   IF ( FPFOI ) &
636                           MOM_HO = FPFOI_OF(W_G(IJK), W_G(IPJK), &
637                                 W_G(IP_OF(IPJK)), TMP4(IPPP4))
638                 ENDIF
639                 IF (.NOT. FPFOI ) &
640                           MOM_HO = XSI_E(IJK)*W_G(IPJK)+ &
641                                    (1.0-XSI_E(IJK))*W_G(IJK)
642     !=======================================================================
643     ! JFD: START MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
644     !=======================================================================
645                 IF(CUT_W_TREATMENT_AT(IJK)) THEN
646                    Flux = (Theta_W_be(IJK) * Flux_gE(IJK) + Theta_W_te(IJK) * Flux_gE(IJKP))
647                    CALL GET_INTERPOLATION_TERMS_G(IJK,'W_MOMENTUM',ALPHA_We_c(IJK),AW,HW,VELW)
648                    Flux = Flux * AW
649                 ELSE   ! Original terms
650                    Flux = HALF * (Flux_gE(IJK) + Flux_gE(IJKP))
651                 ENDIF
652     !=======================================================================
653     ! JFD: END MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
654     !=======================================================================
655                 EAST_DC = Flux*(MOM_LO-MOM_HO)
656     !
657     !           DEFERRED CORRECTION CONTRIBUTION AT THE North face (i, j+1/2, k+1/2)
658     !
659                 IF(V(IJK) >= ZERO)THEN
660                   MOM_LO = W_G(IJK)
661                   IF ( FPFOI ) &
662                           MOM_HO = FPFOI_OF(W_G(IJPK), W_G(IJK), &
663                                 W_G(IJMK), W_G(JM_OF(IJMK)))
664                 ELSE
665                   MOM_LO = W_G(IJPK)
666                   IF ( FPFOI ) &
667                           MOM_HO = FPFOI_OF(W_G(IJK), W_G(IJPK), &
668                                 W_G(JP_OF(IJPK)), TMP4(JPPP4))
669                 ENDIF
670                 IF (.NOT. FPFOI ) &
671                            MOM_HO = XSI_N(IJK)*W_G(IJPK)+ &
672                                     (1.0-XSI_N(IJK))*W_G(IJK)
673     !=======================================================================
674     ! JFD: START MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
675     !=======================================================================
676                 IF(CUT_W_TREATMENT_AT(IJK)) THEN
677                    Flux = (Theta_W_bn(IJK) * Flux_gN(IJK) + Theta_W_tn(IJK) * Flux_gN(IJKP))
678                    CALL GET_INTERPOLATION_TERMS_G(IJK,'W_MOMENTUM',ALPHA_Wn_c(IJK),AW,HW,VELW)
679                    Flux = Flux * AW
680                 ELSE   ! Original terms
681                    Flux = HALF * (Flux_gN(IJK) + Flux_gN(IJKP))
682                 ENDIF
683     !=======================================================================
684     ! JFD: END MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
685     !=======================================================================
686                 NORTH_DC = Flux*(MOM_LO-MOM_HO)
687     !
688     !           DEFERRED CORRECTION CONTRIBUTION AT THE Top face (i, j, k+1)
689     !
690                 IF(WW(IJK) >= ZERO)THEN
691                   MOM_LO = W_G(IJK)
692                   IF ( FPFOI ) &
693                           MOM_HO = FPFOI_OF(W_G(IJKP), W_G(IJK), &
694                                 W_G(IJKM), W_G(KM_OF(IJKM)))
695                 ELSE
696                   MOM_LO = W_G(IJKP)
697                   IF ( FPFOI ) &
698                           MOM_HO = FPFOI_OF(W_G(IJK), W_G(IJKP), &
699                                 W_G(KP_OF(IJKP)), TMP4(KPPP4))
700                 ENDIF
701                 IF (.NOT. FPFOI ) &
702                            MOM_HO = XSI_T(IJK)*W_G(IJKP)+ &
703                                     (1.0-XSI_T(IJK))*W_G(IJK)
704     !=======================================================================
705     ! JFD: START MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
706     !=======================================================================
707                 IF(CUT_W_TREATMENT_AT(IJK)) THEN
708                    Flux = (Theta_Wt_bar(IJK) * Flux_gT(IJK) + Theta_Wt(IJK) * Flux_gT(IJKP))
709                    CALL GET_INTERPOLATION_TERMS_G(IJK,'W_MOMENTUM',alpha_Wt_c(IJK),AW,HW,VELW)
710                    Flux = Flux * AW
711                 ELSE   ! Original terms
712                    Flux = HALF * (Flux_gT(IJK) + Flux_gT(IJKP))
713                 ENDIF
714     !=======================================================================
715     ! JFD: END MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
716     !=======================================================================
717                 TOP_DC = Flux*(MOM_LO-MOM_HO)
718     !
719     !           DEFERRED CORRECTION CONTRIBUTION AT THE West face (i-1/2, j, k+1/2)
720     !
721                 IMJK = IM_OF(IJK)
722                 IM = IM1(I)
723                 IJKW = WEST_OF(IJK)
724                 IJKWT = TOP_OF(IJKW)
725                 IMJKP = KP_OF(IMJK)
726                 IF(U(IMJK) >= ZERO)THEN
727                   MOM_LO = W_G(IMJK)
728                   IF ( FPFOI ) &
729                           MOM_HO = FPFOI_OF(W_G(IJK), W_G(IMJK), &
730                                 W_G(IM_OF(IMJK)), TMP4(IMMM4))
731                 ELSE
732                   MOM_LO = W_G(IJK)
733                   IF ( FPFOI ) &
734                           MOM_HO = FPFOI_OF(W_G(IMJK), W_G(IJK), &
735                                 W_G(IPJK), W_G(IP_OF(IPJK)))
736                 ENDIF
737                 IF (.NOT. FPFOI ) &
738                            MOM_HO = XSI_E(IMJK)*W_G(IJK)+ &
739                                     (1.0-XSI_E(IMJK))*W_G(IMJK)
740     !=======================================================================
741     ! JFD: START MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
742     !=======================================================================
743                 IF(CUT_W_TREATMENT_AT(IJK)) THEN
744                    Flux = (Theta_W_be(IMJK) * Flux_gE(IMJK) + Theta_W_te(IMJK) * Flux_gE(IMJKP))
745                    CALL GET_INTERPOLATION_TERMS_G(IJK,'W_MOMENTUM',ALPHA_We_c(IMJK),AW,HW,VELW)
746                    Flux = Flux * AW
747                 ELSE   ! Original terms
748                    Flux = HALF * (Flux_gE(IMJK) + Flux_gE(IMJKP))
749                 ENDIF
750     !=======================================================================
751     ! JFD: END MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
752     !=======================================================================
753                 WEST_DC = Flux*(MOM_LO-MOM_HO)
754     !
755     !           DEFERRED CORRECTION CONTRIBUTION AT THE South face (i, j-1/2, k+1/2)
756     !
757                 IJMK = JM_OF(IJK)
758                 JM = JM1(J)
759                 IJMKP = KP_OF(IJMK)
760                 IJKS = SOUTH_OF(IJK)
761                 IJKST = TOP_OF(IJKS)
762                 IF(V(IJMK) >= ZERO)THEN
763                   MOM_LO = W_G(IJMK)
764                   IF ( FPFOI ) &
765                           MOM_HO = FPFOI_OF(W_G(IJK), W_G(IJMK), &
766                                 W_G(JM_OF(IJMK)), TMP4(JMMM4))
767                 ELSE
768                   MOM_LO = W_G(IJK)
769                   IF ( FPFOI ) &
770                           MOM_HO = FPFOI_OF(W_G(IJMK), W_G(IJK), &
771                                 W_G(IJPK), W_G(JP_OF(IJPK)))
772                 ENDIF
773                 IF (.NOT. FPFOI ) &
774                             MOM_HO = XSI_N(IJMK)*W_G(IJK)+ &
775                                      (1.0-XSI_N(IJMK))*W_G(IJMK)
776     !=======================================================================
777     ! JFD: START MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
778     !=======================================================================
779                 IF(CUT_W_TREATMENT_AT(IJK)) THEN
780                    Flux = (Theta_W_bn(IJMK) * Flux_gN(IJMK) + Theta_W_tn(IJMK) * Flux_gN(IJMKP))
781                    CALL GET_INTERPOLATION_TERMS_G(IJK,'W_MOMENTUM',ALPHA_Wn_c(IJMK),AW,HW,VELW)
782                    Flux = Flux * AW
783                 ELSE   ! Original terms
784                    Flux = HALF * (Flux_gN(IJMK) + Flux_gN(IJMKP))
785                 ENDIF
786     !=======================================================================
787     ! JFD: END MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
788     !=======================================================================
789                 SOUTH_DC = Flux*(MOM_LO-MOM_HO)
790     !
791     !           DEFERRED CORRECTION CONTRIBUTION AT THE Bottom face (i, j, k)
792     !
793                 IJKM = KM_OF(IJK)
794                 KM = KM1(K)
795                 IJKB = BOTTOM_OF(IJK)
796                 IF(WW(IJK) >= ZERO)THEN
797                   MOM_LO = W_G(IJKM)
798                   IF ( FPFOI ) &
799                           MOM_HO = FPFOI_OF(W_G(IJK), W_G(IJKM), &
800                                 W_G(KM_OF(IJKM)), TMP4(KMMM4))
801                 ELSE
802                   MOM_LO = W_G(IJK)
803                   IF ( FPFOI ) &
804                           MOM_HO = FPFOI_OF(W_G(IJKM), W_G(IJK), &
805                                 W_G(IJKP), W_G(KP_OF(IJKP)))
806                 ENDIF
807                 IF (.NOT. FPFOI ) &
808                            MOM_HO = XSI_T(IJKM)*W_G(IJK)+ &
809                                     (1.0-XSI_T(IJKM))*W_G(IJKM)
810     !=======================================================================
811     ! JFD: START MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
812     !=======================================================================
813                 IF(CUT_W_TREATMENT_AT(IJK)) THEN
814                    Flux = (Theta_Wt_bar(IJKM) * Flux_gT(IJKM) + Theta_Wt(IJKM) * Flux_gT(IJK))
815                    CALL GET_INTERPOLATION_TERMS_G(IJK,'W_MOMENTUM',alpha_Wt_c(IJKM),AW,HW,VELW)
816                    Flux = Flux * AW
817                 ELSE   ! Original terms
818                    Flux = HALF * (Flux_gT(IJKM) + Flux_gT(IJK))
819                 ENDIF
820     !=======================================================================
821     ! JFD: END MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
822     !=======================================================================
823                 BOTTOM_DC = Flux*(MOM_LO-MOM_HO)
824     !
825     !               CONTRIBUTION DUE TO DEFERRED CORRECTION
826     !
827                 B_M(IJK) = B_M(IJK)+WEST_DC-EAST_DC+SOUTH_DC-NORTH_DC&
828                                     +BOTTOM_DC-TOP_DC
829     !
830              ENDIF
831           END DO
832     
833           call unlock_tmp4_array
834           call unlock_tmp_array
835           call unlock_xsi_array
836     
837           RETURN
838           END SUBROUTINE STORE_A_W_GDC
839     
840     
841     !
842     !
843     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
844     !                                                                      C
845     !  Module name: STORE_A_W_g1(A_W_g, IER)                               C
846     !  Purpose: Determine convection diffusion terms for W_g momentum eqs  C
847     !  The off-diagonal coefficients calculated here must be positive. The C
848     !  center coefficient and the source vector are negative. Higher order C
849     !  discretization.                                                     C
850     !  See source_w_g                                                      C
851     !                                                                      C
852     !  Author: M. Syamlal                                 Date:19-MAR-97   C
853     !  Reviewer:                                          Date:            C
854     !                                                                      C
855     !  Revision Number: 1                                                  C
856     !  Purpose: To incorporate Cartesian grid modifications                C
857     !  Author: Jeff Dietiker                              Date: 01-Jul-09  C
858     !                                                                      C
859     !  Literature/Document References:                                     C
860     !                                                                      C
861     !  Variables referenced:                                               C
862     !  Variables modified:                                                 C
863     !                                                                      C
864     !  Local variables:                                                    C
865     !                                                                      C
866     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
867     !
868           SUBROUTINE STORE_A_W_G1(A_W_G)
869     !...Translated by Pacific-Sierra Research VAST-90 2.06G5  12:17:31  12/09/98
870     !...Switches: -xf
871     !
872     !  Include param.inc file to specify parameter values
873     !
874     !-----------------------------------------------
875     !   M o d u l e s
876     !-----------------------------------------------
877           USE param
878           USE param1
879           USE parallel
880           USE matrix
881           USE geometry
882           USE indices
883           USE run
884           USE visc_g
885           USE toleranc
886           USE physprop
887           USE fldvar
888           USE output
889           USE vshear
890           USE xsi
891           USE xsi_array
892           Use tmp_array,  U => Array1, V => Array2, WW => Array3
893           USE compar
894           USE mflux
895           USE fun_avg
896           USE functions
897     !=======================================================================
898     ! JFD: START MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
899     !=======================================================================
900           USE cutcell
901     !=======================================================================
902     ! JFD: END MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
903     !=======================================================================
904     
905           IMPLICIT NONE
906     !-----------------------------------------------
907     !   G l o b a l   P a r a m e t e r s
908     !-----------------------------------------------
909     !-----------------------------------------------
910     !   D u m m y   A r g u m e n t s
911     !-----------------------------------------------
912     !
913     !                      Indices
914           INTEGER          I,  J, K, IPJK, IJPK, IJKN, IJKC, KP, IJKE,&
915                            IJKTE, IJKP, IJKT, IJKTN, IJK
916           INTEGER          IMJK, IM, IJKW, IJKWT, IMJKP
917           INTEGER          IJMK, JM, IJMKP, IJKS, IJKST
918           INTEGER          IJKM, KM, IJKB
919     !
920     ! loezos
921             INTEGER incr
922     ! loezos
923     
924     !                      Face mass flux
925           DOUBLE PRECISION Flux
926     
927     !                      Diffusion parameter
928           DOUBLE PRECISION D_f
929     !
930     !                      Septadiagonal matrix A_W_g
931           DOUBLE PRECISION A_W_g(DIMENSION_3, -3:3)
932     !
933     !                      Convection weighting factors
934     !      DOUBLE PRECISION XSI_e(DIMENSION_3), XSI_n(DIMENSION_3),&
935     !                       XSI_t(DIMENSION_3)
936     !      DOUBLE PRECISION U(DIMENSION_3),&
937     !                       V(DIMENSION_3), WW(DIMENSION_3)
938     !=======================================================================
939     ! JFD: START MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
940     !=======================================================================
941           DOUBLE PRECISION :: AW,HW,VELW
942     !=======================================================================
943     ! JFD: END MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
944     !=======================================================================
945     
946           call lock_tmp_array
947           call lock_xsi_array
948     
949     !
950     !  Calculate convection factors
951     !
952     
953     !!!$omp parallel do private(IJK,K,IJKT,IJKP)
954           DO IJK = ijkstart3, ijkend3
955              K = K_OF(IJK)
956              IJKT = TOP_OF(IJK)
957              IJKP = KP_OF(IJK)
958     !
959     !
960     !=======================================================================
961     ! JFD: START MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
962     !=======================================================================
963              IF(CUT_W_TREATMENT_AT(IJK)) THEN
964     !           East face (i+1/2, j, k+1/2)
965                 U(IJK) = (Theta_W_be(IJK) * U_G(IJK) + Theta_W_te(IJK) * U_G(IJKP))
966                 CALL GET_INTERPOLATION_TERMS_G(IJK,'W_MOMENTUM',ALPHA_We_c(IJK),AW,HW,VELW)
967                 U(IJK) = U(IJK) * AW
968     !           North face (i, j+1/2, k+1/2)
969                 V(IJK) = (Theta_W_bn(IJK) * V_G(IJK) + Theta_W_tn(IJK) * V_G(IJKP))
970                 CALL GET_INTERPOLATION_TERMS_G(IJK,'W_MOMENTUM',ALPHA_Wn_c(IJK),AW,HW,VELW)
971                 V(IJK) = V(IJK) * AW
972     !           Top face (i, j, k+1)
973                 WW(IJK) = (Theta_Wt_bar(IJK) * W_G(IJK) + Theta_Wt(IJK) * W_G(IJKP))
974                 CALL GET_INTERPOLATION_TERMS_G(IJK,'W_MOMENTUM',alpha_Wt_c(IJK),AW,HW,VELW)
975                 WW(IJK) = WW(IJK) * AW
976              ELSE   ! Original terms
977                 U(IJK) = AVG_Z(U_G(IJK),U_G(IJKP),K)
978                 V(IJK) = AVG_Z(V_G(IJK),V_G(IJKP),K)
979                 WW(IJK) = AVG_Z_T(W_G(IJK),W_G(IJKP))
980              ENDIF
981     
982     !=======================================================================
983     ! JFD: END MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
984     !=======================================================================
985           END DO
986     
987     ! loezos
988             incr=0
989     ! loezos
990     
991           CALL CALC_XSI (DISCRETIZE(5), W_G, U, V, WW, XSI_E, XSI_N,&
992                    XSI_T,incr)
993     
994     ! loezos
995     ! update to true velocity
996           IF (SHEAR) THEN
997     !!!$omp parallel do private(IJK)
998              DO IJK = ijkstart3, ijkend3
999              IF (FLUID_AT(IJK)) THEN
1000                V(IJK)=V(IJK)+VSH(IJK)
1001               END IF
1002             END DO
1003           END IF
1004     ! loezos
1005     
1006     !
1007     !
1008     !
1009     !  Calculate convection-diffusion fluxes through each of the faces
1010     !
1011     !
1012     
1013     !!!!$omp      parallel do                                               &
1014     !!!!$omp&     private( I,  J, K, IPJK, IJPK, IJKN, IJKC, KP,     &
1015     !!!!$omp&             IJKE, IJKTE, IJKP, IJKT, IJKTN, IJK,  D_f,        &
1016     !!!!$omp&             IMJK, IM, IJKW, IJKWT, IMJKP,                     &
1017     !!!!$omp&             IJMK, JM, IJMKP, IJKS, IJKST,                     &
1018     !!!!$omp&             IJKM, KM, IJKB)
1019           DO IJK = ijkstart3, ijkend3
1020     !
1021              IF (FLOW_AT_T(IJK)) THEN
1022                 I = I_OF(IJK)
1023                 J = J_OF(IJK)
1024                 K = K_OF(IJK)
1025                 IPJK = IP_OF(IJK)
1026                 IJPK = JP_OF(IJK)
1027                 IJKN = NORTH_OF(IJK)
1028                 IJKT = TOP_OF(IJK)
1029                 IF (WALL_AT(IJK)) THEN
1030                    IJKC = IJKT
1031                 ELSE
1032                    IJKC = IJK
1033                 ENDIF
1034                 KP = KP1(K)
1035                 IJKE = EAST_OF(IJK)
1036                 IJKP = KP_OF(IJK)
1037                 IJKTN = NORTH_OF(IJKT)
1038                 IJKTE = EAST_OF(IJKT)
1039     !
1040     !           East face (i+1/2, j, k+1/2)
1041     !=======================================================================
1042     ! JFD: START MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
1043     !=======================================================================
1044                 IF(CUT_W_TREATMENT_AT(IJK)) THEN
1045                    Flux = (Theta_W_be(IJK) * Flux_gE(IJK) + Theta_W_te(IJK) * Flux_gE(IJKP))
1046                    CALL GET_INTERPOLATION_TERMS_G(IJK,'W_MOMENTUM',ALPHA_We_c(IJK),AW,HW,VELW)
1047                    Flux = Flux * AW
1048                    D_F = AVG_Z_H(AVG_X_H(MU_GT(IJKC),MU_GT(IJKE),I),AVG_X_H(MU_GT(IJKT&
1049                       ),MU_GT(IJKTE),I),K)*ONEoDX_E_W(IJK)*AYZ_W(IJK)
1050                 ELSE   ! Original terms
1051                    Flux = HALF * (Flux_gE(IJK) + Flux_gE(IJKP))
1052                    D_F = AVG_Z_H(AVG_X_H(MU_GT(IJKC),MU_GT(IJKE),I),AVG_X_H(MU_GT(IJKT&
1053                       ),MU_GT(IJKTE),I),K)*ODX_E(I)*AYZ_W(IJK)
1054                 ENDIF
1055     !=======================================================================
1056     ! JFD: END MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
1057     !=======================================================================
1058     !
1059                 A_W_G(IJK,E) = D_F - XSI_E(IJK)*Flux
1060     !
1061                 A_W_G(IPJK,W) = D_F + (ONE - XSI_E(IJK))*Flux
1062     !
1063     !           North face (i, j+1/2, k+1/2)
1064     !=======================================================================
1065     ! JFD: START MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
1066     !=======================================================================
1067                 IF(CUT_W_TREATMENT_AT(IJK)) THEN
1068                    Flux = (Theta_W_bn(IJK) * Flux_gN(IJK) + Theta_W_tn(IJK) * Flux_gN(IJKP))
1069                    CALL GET_INTERPOLATION_TERMS_G(IJK,'W_MOMENTUM',ALPHA_Wn_c(IJK),AW,HW,VELW)
1070                    Flux = Flux * AW
1071                    D_F = AVG_Z_H(AVG_Y_H(MU_GT(IJKC),MU_GT(IJKN),J),AVG_Y_H(MU_GT(IJKT&
1072                       ),MU_GT(IJKTN),J),K)*ONEoDY_N_W(IJK)*AXZ_W(IJK)
1073                 ELSE   ! Original terms
1074                    Flux = HALF * (Flux_gN(IJK) + Flux_gN(IJKP))
1075                    D_F = AVG_Z_H(AVG_Y_H(MU_GT(IJKC),MU_GT(IJKN),J),AVG_Y_H(MU_GT(IJKT&
1076                       ),MU_GT(IJKTN),J),K)*ODY_N(J)*AXZ_W(IJK)
1077                 ENDIF
1078     !=======================================================================
1079     ! JFD: END MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
1080     !=======================================================================
1081     !
1082                 A_W_G(IJK,N) = D_F - XSI_N(IJK)*Flux
1083     !
1084                 A_W_G(IJPK,S) = D_F + (ONE - XSI_N(IJK))*Flux
1085     !
1086     !           Top face (i, j, k+1)
1087     !=======================================================================
1088     ! JFD: START MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
1089     !=======================================================================
1090                 IF(CUT_W_TREATMENT_AT(IJK)) THEN
1091                    Flux = (Theta_Wt_bar(IJK) * Flux_gT(IJK) + Theta_Wt(IJK) * Flux_gT(IJKP))
1092                    CALL GET_INTERPOLATION_TERMS_G(IJK,'W_MOMENTUM',alpha_Wt_c(IJK),AW,HW,VELW)
1093                    Flux = Flux * AW
1094                    D_F = MU_GT(IJKT)*ONEoDZ_T_W(IJK)*AXY_W(IJK)
1095                 ELSE   ! Original terms
1096                    Flux = HALF * (Flux_gT(IJK) + Flux_gT(IJKP))
1097                    D_F = MU_GT(IJKT)*OX(I)*ODZ(KP)*AXY_W(IJK)
1098                 ENDIF
1099     !=======================================================================
1100     ! JFD: END MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
1101     !=======================================================================
1102                 A_W_G(IJK,T) = D_F - XSI_T(IJK)*Flux
1103                 A_W_G(IJKP,B) = D_F + (ONE - XSI_T(IJK))*Flux
1104     !
1105     !           West face (i-1/2, j, k+1/2)
1106                 IMJK = IM_OF(IJK)
1107                 IF (.NOT.FLOW_AT_T(IMJK)) THEN
1108                    IM = IM1(I)
1109                    IJKW = WEST_OF(IJK)
1110                    IJKWT = TOP_OF(IJKW)
1111                    IMJKP = KP_OF(IMJK)
1112     !=======================================================================
1113     ! JFD: START MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
1114     !=======================================================================
1115                    IF(CUT_W_TREATMENT_AT(IJK)) THEN
1116                       Flux = (Theta_W_be(IMJK) * Flux_gE(IMJK) + Theta_W_te(IMJK) * Flux_gE(IMJKP))
1117                       CALL GET_INTERPOLATION_TERMS_G(IJK,'W_MOMENTUM',ALPHA_We_c(IMJK),AW,HW,VELW)
1118                       Flux = Flux * AW
1119                       D_F = AVG_Z_H(AVG_X_H(MU_GT(IJKW),MU_GT(IJKC),IM),AVG_X_H(MU_GT(&
1120                           IJKWT),MU_GT(IJKT),IM),K)*ONEoDX_E_W(IMJK)*AYZ_W(IMJK)
1121                    ELSE   ! Original terms
1122                       Flux = HALF * (Flux_gE(IMJK) + Flux_gE(IMJKP))
1123                       D_F = AVG_Z_H(AVG_X_H(MU_GT(IJKW),MU_GT(IJKC),IM),AVG_X_H(MU_GT(&
1124                          IJKWT),MU_GT(IJKT),IM),K)*ODX_E(IM)*AYZ_W(IMJK)
1125                    ENDIF
1126     !=======================================================================
1127     ! JFD: END MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
1128     !=======================================================================
1129                    A_W_G(IJK,W) = D_F + (ONE - XSI_E(IMJK))*Flux
1130                 ENDIF
1131     !
1132     !           South face (i, j-1/2, k+1/2)
1133                 IJMK = JM_OF(IJK)
1134                 IF (.NOT.FLOW_AT_T(IJMK)) THEN
1135                    JM = JM1(J)
1136                    IJMKP = KP_OF(IJMK)
1137                    IJKS = SOUTH_OF(IJK)
1138                    IJKST = TOP_OF(IJKS)
1139     !=======================================================================
1140     ! JFD: START MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
1141     !=======================================================================
1142                    IF(CUT_W_TREATMENT_AT(IJK)) THEN
1143                       Flux = (Theta_W_bn(IJMK) * Flux_gN(IJMK) + Theta_W_tn(IJMK) * Flux_gN(IJMKP))
1144                       CALL GET_INTERPOLATION_TERMS_G(IJK,'W_MOMENTUM',ALPHA_Wn_c(IJMK),AW,HW,VELW)
1145                       Flux = Flux * AW
1146                       D_F = AVG_Z_H(AVG_Y_H(MU_GT(IJKS),MU_GT(IJKC),JM),AVG_Y_H(MU_GT(&
1147                          IJKST),MU_GT(IJKT),JM),K)*ONEoDY_N_W(IJMK)*AXZ_W(IJMK)
1148                    ELSE   ! Original terms
1149                       Flux = HALF * (Flux_gN(IJMK) + Flux_gN(IJMKP))
1150                       D_F = AVG_Z_H(AVG_Y_H(MU_GT(IJKS),MU_GT(IJKC),JM),AVG_Y_H(MU_GT(&
1151                          IJKST),MU_GT(IJKT),JM),K)*ODY_N(JM)*AXZ_W(IJMK)
1152                    ENDIF
1153     !=======================================================================
1154     ! JFD: END MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
1155     !=======================================================================
1156                    A_W_G(IJK,S) = D_F + (ONE - XSI_N(IJMK))*Flux
1157                 ENDIF
1158     !
1159     !           Bottom face (i, j, k)
1160                 IJKM = KM_OF(IJK)
1161                 IF (.NOT.FLOW_AT_T(IJKM)) THEN
1162                    KM = KM1(K)
1163                    IJKB = BOTTOM_OF(IJK)
1164     !=======================================================================
1165     ! JFD: START MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
1166     !=======================================================================
1167                    IF(CUT_W_TREATMENT_AT(IJK)) THEN
1168                       Flux = (Theta_Wt_bar(IJKM) * Flux_gT(IJKM) + Theta_Wt(IJKM) * Flux_gT(IJK))
1169                       CALL GET_INTERPOLATION_TERMS_G(IJK,'W_MOMENTUM',alpha_Wt_c(IJKM),AW,HW,VELW)
1170                       Flux = Flux * AW
1171                       D_F = MU_GT(IJK)*ONEoDZ_T_W(IJKM)*AXY_W(IJKM)
1172                    ELSE   ! Original terms
1173                       Flux = HALF * (Flux_gT(IJKM) + Flux_gT(IJK))
1174                       D_F = MU_GT(IJK)*OX(I)*ODZ(K)*AXY_W(IJKM)
1175                    ENDIF
1176     !=======================================================================
1177     ! JFD: END MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
1178     !=======================================================================
1179                    A_W_G(IJK,B) = D_F + (ONE - XSI_T(IJKM))*Flux
1180                 ENDIF
1181              ENDIF
1182           END DO
1183     
1184           call unlock_tmp_array
1185           call unlock_xsi_array
1186     
1187           RETURN
1188           END SUBROUTINE STORE_A_W_G1
1189     
1190     !// Comments on the modifications for DMP version implementation
1191     !// 001 Include header file and common declarations for parallelization
1192     !// 350 Changed do loop limits: 1,ijkmax2-> ijkstart3, ijkend3
1193     
1194