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

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