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

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