File: RELATIVE:/../../../mfix.git/model/conv_dif_u_g.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Subroutine: CONV_DIF_U_g                                            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     !                                                                      C
11     !                                                                      C
12     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
13           SUBROUTINE CONV_DIF_U_G(A_M, B_M)
14     
15     ! Modules
16     !---------------------------------------------------------------------//
17           USE param, only: dimension_3, dimension_m
18           USE run, only: momentum_x_eq
19           USE run, only: discretize
20           USE run, only: def_cor
21           USE visc_g, only: epmu_gt
22           IMPLICIT NONE
23     
24     ! Dummy arguments
25     !---------------------------------------------------------------------//
26     ! Septadiagonal matrix A_m
27           DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
28     ! Vector b_m
29           DOUBLE PRECISION, INTENT(INOUT) :: B_m(DIMENSION_3, 0:DIMENSION_M)
30     !---------------------------------------------------------------------//
31     
32           IF (.NOT.MOMENTUM_X_EQ(0)) RETURN
33     
34           IF(DEF_COR)THEN
35     ! USE DEFERRED CORRECTION TO SOLVE U_G
36              CALL STORE_A_U_G0(A_M)
37              IF (DISCRETIZE(3) > 1) CALL STORE_A_U_GDC(B_M(1,0))
38     
39           ELSE
40     ! DO NOT USE DEFERRED CORRECTION TO SOLVE FOR U_G
41              IF (DISCRETIZE(3) == 0) THEN               ! 0 & 1 => FOUP
42                 CALL STORE_A_U_G0(A_M)
43              ELSE
44                 CALL STORE_A_U_G1(A_M)
45              ENDIF
46           ENDIF
47     
48           CALL DIF_U_IS(EPMU_GT, A_M, 0)
49     
50           RETURN
51           END SUBROUTINE CONV_DIF_U_G
52     
53     
54     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
55     !                                                                      C
56     !  Purpose: Calculate the components of velocity on the east, north,   C
57     !  and top face of a u-momentum cell                                   C
58     !                                                                      C
59     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
60           SUBROUTINE GET_UCELL_GVTERMS(U, V, WW)
61     
62     ! Modules
63     !---------------------------------------------------------------------//
64           USE compar, only: ijkstart3, ijkend3
65     
66           USE cutcell, only: cut_u_treatment_at
67           USE cutcell, only: theta_ue, theta_ue_bar
68           USE cutcell, only: theta_u_ne, theta_u_nw
69           USE cutcell, only: theta_u_te, theta_u_tw
70           USE cutcell, only: alpha_ue_c, alpha_un_c, alpha_ut_c
71     
72           USE fldvar, only: u_g, v_g, w_g
73     
74           USE fun_avg, only: avg_x_e, avg_x
75           USE functions, only: ip_of
76           USE geometry, only: do_k
77           USE indices, only: i_of, ip1
78     
79           USE param, only: dimension_3
80           IMPLICIT NONE
81     
82     ! Dummy arguments
83     !---------------------------------------------------------------------//
84           DOUBLE PRECISION, INTENT(OUT) :: U(DIMENSION_3)
85           DOUBLE PRECISION, INTENT(OUT) :: V(DIMENSION_3)
86           DOUBLE PRECISION, INTENT(OUT) :: WW(DIMENSION_3)
87     
88     ! Local variables
89     !---------------------------------------------------------------------//
90     ! indices
91           INTEGER :: IJK, I, IP, IPJK
92     ! for cartesian grid
93           DOUBLE PRECISION :: AW, HW, VELW
94     !---------------------------------------------------------------------//
95     
96     !!!$omp parallel do private(IJK,I,IP,IPJK)
97           DO IJK = ijkstart3, ijkend3
98              I = I_OF(IJK)
99              IP = IP1(I)
100              IPJK = IP_OF(IJK)
101     
102              IF(CUT_U_TREATMENT_AT(IJK)) THEN
103     
104     ! East face (i+1, j, k)
105                 U(IJK) = (Theta_Ue_bar(IJK) * U_g(IJK) + &
106                           Theta_Ue(IJK) * U_g(IPJK))
107                 CALL GET_INTERPOLATION_TERMS_G(IJK,'U_MOMENTUM',&
108                    alpha_Ue_c(IJK), AW, HW, VELW)
109                 U(IJK) = U(IJK) * AW
110     
111     ! North face (i+1/2, j+1/2, k)
112                 V(IJK) = (Theta_U_nw(IJK) * V_g(IJK) + &
113                           Theta_U_ne(IJK) * V_g(IPJK))
114                 CALL GET_INTERPOLATION_TERMS_G(IJK,'U_MOMENTUM',&
115                    ALPHA_Un_c(IJK), AW, HW, VELW)
116                 V(IJK) = V(IJK) * AW
117     
118     ! Top face (i+1/2, j, k+1/2)
119                 IF (DO_K) THEN
120                    WW(IJK) = (Theta_U_tw(IJK) * W_g(IJK) + &
121                               Theta_U_te(IJK) * W_g(IPJK))
122                    CALL GET_INTERPOLATION_TERMS_G(IJK,'U_MOMENTUM',&
123                       ALPHA_Ut_c(IJK), AW, HW, VELW)
124                    WW(IJK) = WW(IJK) * AW
125                 ENDIF
126     
127              ELSE
128                 U(IJK) = AVG_X_E(U_G(IJK),U_G(IPJK),IP)
129                 V(IJK) = AVG_X(V_G(IJK),V_G(IPJK),I)
130                 IF (DO_K) WW(IJK) = AVG_X(W_G(IJK),W_G(IPJK),I)
131              ENDIF
132           ENDDO   ! end do ijk
133     
134           RETURN
135           END SUBROUTINE GET_UCELL_GVTERMS
136     
137     
138     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
139     !                                                                      C
140     !  Purpose: Calculate the convective fluxes through the faces of a     C
141     !  u-momentum cell. Note the fluxes are calculated at all faces of     C
142     !  regardless of flow_at_e of condition of the west, south, or         C
143     !  bottom face.                                                        C
144     !                                                                      C
145     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
146           SUBROUTINE GET_UCELL_GCFLUX_TERMS(FLUX_E, FLUX_W, FLUX_N, &
147              FLUX_S, FLUX_T, FLUX_B, IJK)
148     
149     ! Modules
150     !---------------------------------------------------------------------//
151           USE cutcell, only: cut_u_treatment_at
152           USE cutcell, only: theta_ue, theta_ue_bar
153           USE cutcell, only: theta_u_ne, theta_u_nw
154           USE cutcell, only: theta_u_te, theta_u_tw
155           USE cutcell, only: alpha_ue_c, alpha_un_c, alpha_ut_c
156           USE functions, only: ip_of, im_of, jm_of, km_of
157           USE geometry, only: do_k
158           USE mflux, only: flux_ge, flux_gn, flux_gt
159     
160           USE param1, only: half
161           IMPLICIT NONE
162     
163     ! Dummy arguments
164     !---------------------------------------------------------------------//
165     ! fluxes through faces of given ijk u-momentum cell
166           DOUBLE PRECISION, INTENT(OUT) :: flux_e, flux_w
167           DOUBLE PRECISION, INTENT(OUT) :: flux_n, flux_s
168           DOUBLE PRECISION, INTENT(OUT) :: flux_t, flux_b
169     ! ijk index
170           INTEGER, INTENT(IN) :: ijk
171     
172     ! Local variables
173     !---------------------------------------------------------------------//
174     ! indices
175           INTEGER :: imjk, ijmk, ijkm
176           INTEGER :: ipjk, ipjmk, ipjkm
177     
178     ! for cartesian grid
179           DOUBLE PRECISION :: AW, HW, VELW
180     !---------------------------------------------------------------------//
181     
182     ! indices
183           IPJK = IP_OF(IJK)
184           IMJK = IM_OF(IJK)
185           IJMK = JM_OF(IJK)
186           IJKM = KM_OF(IJK)
187           IPJMK = IP_OF(IJMK)
188           IPJKM = IP_OF(IJKM)
189     
190     ! First calculate the fluxes at the faces
191           IF(CUT_U_TREATMENT_AT(IJK)) THEN
192     ! east face: i+1, j, k
193              Flux_e = (Theta_Ue_bar(IJK) * Flux_gE(IJK) + &
194                        Theta_Ue(IJK) * Flux_gE(IPJK))
195              CALL GET_INTERPOLATION_TERMS_G(IJK,'U_MOMENTUM',&
196                 alpha_Ue_c(IJK), AW, HW, VELW)
197              Flux_e = Flux_e * AW
198     ! west face: i, j, k
199              Flux_w = (Theta_Ue_bar(IMJK) * Flux_gE(IMJK) + &
200                        Theta_Ue(IMJK) * Flux_gE(IJK))
201              CALL GET_INTERPOLATION_TERMS_G(IJK,'U_MOMENTUM',&
202                 alpha_Ue_c(IMJK), AW, HW, VELW)
203              Flux_w = Flux_w * AW
204     
205     
206     ! north face: i+1/2, j+1/2, k
207              Flux_n = (Theta_U_nw(IJK) * Flux_gN(IJK) + &
208                        Theta_U_ne(IJK) * Flux_gN(IPJK))
209              CALL GET_INTERPOLATION_TERMS_G(IJK,'U_MOMENTUM', &
210                 ALPHA_Un_c(IJK), AW, HW, VELW)
211              Flux_n = Flux_n * AW
212     ! south face: i+1/2, j-1/2, k
213              Flux_s = (Theta_U_nw(IJMK) * Flux_gN(IJMK) + &
214                        Theta_U_ne(IJMK) * Flux_gN(IPJMK))
215              CALL GET_INTERPOLATION_TERMS_G(IJK,'U_MOMENTUM',&
216                 ALPHA_Un_c(IJMK), AW, HW, VELW)
217              Flux_s = Flux_s * AW
218     
219              IF (DO_K) THEN
220     ! top face: i+1/2, j, k+1/2
221                 Flux_t = (Theta_U_tw(IJK) * Flux_gT(IJK) + &
222                          Theta_U_te(IJK) * Flux_gT(IPJK))
223                 CALL GET_INTERPOLATION_TERMS_G(IJK,'U_MOMENTUM', &
224                    ALPHA_Ut_c(IJK), AW, HW, VELW)
225                 Flux_t = Flux_t * AW
226     ! bottom face: i+1/2, j, k-1/2
227                 Flux_b = (Theta_U_tw(IJKM) * Flux_gT(IJKM) + &
228                           Theta_U_te(IJKM) * Flux_gT(IPJKM))
229                 CALL GET_INTERPOLATION_TERMS_G(IJK,'U_MOMENTUM',&
230                    ALPHA_Ut_c(IJKM), AW, HW, VELW)
231                 Flux_b = Flux_b * AW
232              ENDIF
233           ELSE
234              Flux_e = HALF * (Flux_gE(IJK) + Flux_gE(IPJK))
235              Flux_w = HALF * (Flux_gE(IMJK) + Flux_gE(IJK))
236              Flux_n = HALF * (Flux_gN(IJK) + Flux_gN(IPJK))
237              Flux_s = HALF * (Flux_gN(IJMK) + Flux_gN(IPJMK))
238     
239              IF (DO_K) THEN
240                 Flux_t = HALF * (Flux_gT(IJK) + Flux_gT(IPJK))
241                 Flux_b = HALF * (Flux_gT(IJKM) + Flux_gT(IPJKM))
242              ENDIF
243           ENDIF   ! end if/else cut_u_treatment_at
244           RETURN
245           END SUBROUTINE GET_UCELL_GCFLUX_TERMS
246     
247     
248     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
249     !                                                                      C
250     !  Purpose: Calculate the components of diffusive flux through the     C
251     !  faces of a u-momentum cell. Note the fluxes are calculated at       C
252     !  all faces regardless of flow_at_e condition of the west, south      C
253     !  or bottom face.                                                     C
254     !                                                                      C
255     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
256           SUBROUTINE GET_UCELL_GDIFF_TERMS(D_FE, D_FW, D_FN, D_FS, &
257              D_FT, D_FB, IJK)
258     
259     ! Modules
260     !---------------------------------------------------------------------//
261           USE cutcell, only: cut_u_treatment_at
262           USE cutcell, only: oneodx_e_u, oneody_n_u, oneodz_t_u
263     
264           USE fldvar, only: epg_jfac
265     
266           USE functions, only: wall_at
267           USE functions, only: east_of, north_of, top_of
268           USE functions, only: south_of, bottom_of
269           USE functions, only: im_of, jm_of, km_of
270     
271           USE geometry, only: odx, ody_n, odz_t
272           USE geometry, only: do_k
273           USE geometry, only: ox_e
274           USE geometry, only: ayz_u, axz_u, axy_u
275     
276           USE indices, only: i_of, j_of, k_of
277           USE indices, only: ip1, jm1, km1
278     
279           USE matrix, only: e, w, s, n, t, b
280           USE param1, only: zero
281           USE visc_g, only: epmu_gt, DF_GU
282           IMPLICIT NONE
283     
284     ! Dummy arguments
285     !---------------------------------------------------------------------//
286     ! diffusion through faces of given ijk u-momentum cell
287           DOUBLE PRECISION, INTENT(OUT) :: d_fe, d_fw
288           DOUBLE PRECISION, INTENT(OUT) :: d_fn, d_fs
289           DOUBLE PRECISION, INTENT(OUT) :: d_ft, d_fb
290     ! ijk index
291           INTEGER, INTENT(IN) :: ijk
292     
293     ! Local variables
294     !---------------------------------------------------------------------//
295     ! indices
296           INTEGER :: imjk, ijmk, ijkm
297           INTEGER :: i, j, k, ip, jm, km
298           INTEGER :: ijkc, ijke, ijkn, ijkne, ijks, ijkse
299           INTEGER :: ijkt, ijkte, ijkb, ijkbe
300     ! length terms
301           DOUBLE PRECISION :: C_AE, C_AW, C_AN, C_AS, C_AT, C_AB
302     ! avg voidage
303           DOUBLE PRECISION :: EPGA
304     !---------------------------------------------------------------------//
305     
306           IMJK = IM_OF(IJK)
307           IJMK = JM_OF(IJK)
308           IJKM = KM_OF(IJK)
309     
310           I = I_OF(IJK)
311           J = J_OF(IJK)
312           K = K_OF(IJK)
313           IP = IP1(I)
314           JM = JM1(J)
315           KM = KM1(K)
316     
317           IJKE = EAST_OF(IJK)
318           IF (WALL_AT(IJK)) THEN
319              IJKC = IJKE
320           ELSE
321              IJKC = IJK
322           ENDIF
323           IJKN = NORTH_OF(IJK)
324           IJKNE = EAST_OF(IJKN)
325           IJKS = SOUTH_OF(IJK)
326           IJKSE = EAST_OF(IJKS)
327     
328           IF(CUT_U_TREATMENT_AT(IJK)) THEN
329              C_AE = ONEoDX_E_U(IJK)
330              C_AW = ONEoDX_E_U(IMJK)
331              C_AN = ONEoDY_N_U(IJK)
332              C_AS = ONEoDY_N_U(IJMK)
333              C_AT = ONEoDZ_T_U(IJK)
334              C_AB = ONEoDZ_T_U(IJKM)
335           ELSE
336              C_AE = ODX(IP)
337              C_AW = ODX(I)
338              C_AN = ODY_N(J)
339              C_AS = ODY_N(JM)
340              C_AT = ODZ_T(K)
341              C_AB = ODZ_T(KM)
342           ENDIF
343     
344     ! East face (i+1, j, k)
345           D_FE = EPMU_GT(IJKE)*C_AE*AYZ_U(IJK)
346     ! West face (i, j, k)
347           D_FW = EPMU_GT(IJKC)*C_AW*AYZ_U(IMJK)
348     
349     
350     ! North face (i+1/2, j+1/2, k)
351           D_FN = AVG_X_H(AVG_Y_H(EPMU_GT(IJKC),EPMU_GT(IJKN),J),&
352                          AVG_Y_H(EPMU_GT(IJKE),EPMU_GT(IJKNE),J),I)*&
353                  C_AN*AXZ_U(IJK)
354     ! South face (i+1/2, j-1/2, k)
355           D_FS = AVG_X_H(AVG_Y_H(EPMU_GT(IJKS),EPMU_GT(IJKC),JM),&
356                          AVG_Y_H(EPMU_GT(IJKSE),EPMU_GT(IJKE),JM),I)*&
357                  C_AS*AXZ_U(IJMK)
358     
359           D_FT = ZERO
360           D_FB = ZERO
361           IF (DO_K) THEN
362              IJKT = TOP_OF(IJK)
363              IJKTE = EAST_OF(IJKT)
364              IJKB = BOTTOM_OF(IJK)
365              IJKBE = EAST_OF(IJKB)
366     
367     ! Top face (i+1/2, j, k+1/2)
368              D_FT = AVG_X_H(AVG_Z_H(EPMU_GT(IJKC),EPMU_GT(IJKT),K),&
369                             AVG_Z_H(EPMU_GT(IJKE),EPMU_GT(IJKTE),K),I)*&
370                     OX_E(I)*C_AT*AXY_U(IJK)
371     ! Bottom face (i+1/2, j, k-1/2)
372              D_FB = AVG_X_H(AVG_Z_H(EPMU_GT(IJKB),EPMU_GT(IJKC),KM),&
373                             AVG_Z_H(EPMU_GT(IJKBE),EPMU_GT(IJKE),KM),I)*&
374                     OX_E(I)*C_AB*AXY_U(IJKM)
375           ENDIF
376     
377           DF_GU(IJK,E) = D_FE
378           DF_GU(IJK,W) = D_FW
379           DF_GU(IJK,N) = D_FN
380           DF_GU(IJK,S) = D_FS
381           DF_GU(IJK,T) = D_FT
382           DF_GU(IJK,B) = D_FB
383     
384     ! if jackson, implement jackson style governing equations: multiply by
385     ! the void fraction otherwise multiply by 1
386           EPGA = AVG_X(EPG_jfac(IJKC), EPG_jfac(IJKE), I)
387           D_FE = EPGA*D_FE
388           D_FW = EPGA*D_FW
389           D_FN = EPGA*D_FN
390           D_FS = EPGA*D_FS
391           D_FT = EPGA*D_FT
392           D_FB = EPGA*D_FB
393     
394           RETURN
395     
396         CONTAINS
397     
398           INCLUDE 'fun_avg.inc'
399     
400         END SUBROUTINE GET_UCELL_GDIFF_TERMS
401     
402     
403     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
404     !                                                                      C
405     !  Subroutine: STORE_A_U_g0                                            C
406     !  Purpose: Determine convection diffusion terms for U_g momentum eqs. C
407     !  The off-diagonal coefficients calculated here must be positive.     C
408     !  The center coefficient and the source vector are negative. See      C
409     !  source_u_g.                                                         C
410     !  Implement FOUP discretization                                       C
411     !                                                                      C
412     !  Author: M. Syamlal                                 Date: 29-APR-96  C
413     !                                                                      C
414     !  Revision Number: 1                                                  C
415     !  Purpose: To incorporate Cartesian grid modifications                C
416     !  Author: Jeff Dietiker                              Date: 01-Jul-09  C
417     !                                                                      C
418     !                                                                      C
419     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
420           SUBROUTINE STORE_A_U_G0(A_U_G)
421     
422     ! Modules
423     !---------------------------------------------------------------------//
424           USE compar, only: ijkstart3, ijkend3
425     
426           USE functions, only: flow_at_e
427           USE functions, only: ip_of, jp_of, kp_of
428           USE functions, only: im_of, jm_of, km_of
429     
430           USE geometry, only: do_k
431     
432           USE param, only: dimension_3, dimension_m
433           USE param1, only: zero
434           USE matrix, only: e, w, n, s, t, b
435           IMPLICIT NONE
436     
437     ! Dummy arguments
438     !---------------------------------------------------------------------//
439     ! Septadiagonal matrix A_U_g
440           DOUBLE PRECISION, INTENT(INOUT) :: A_U_g(DIMENSION_3, -3:3, 0:DIMENSION_M)
441     
442     ! Local variables
443     !---------------------------------------------------------------------//
444     ! Indices
445           INTEGER :: IJK
446           INTEGER :: IMJK, IPJK, IJMK, IJPK, IJKM, IJKP
447     ! Face mass flux
448           DOUBLE PRECISION :: flux_e, flux_w, flux_n, flux_s
449           DOUBLE PRECISION :: flux_t, flux_b
450     ! Diffusion parameter
451           DOUBLE PRECISION :: D_fe, d_fw, d_fn, d_fs, d_ft, d_fb
452     
453     !---------------------------------------------------------------------//
454     
455     !$omp     parallel do default(none)                                &
456     !$omp     private(IJK, IPJK, IJPK, IJKP, IMJK, IJMK, IJKM,         &
457     !$omp             D_fe, d_fw, d_fn, d_fs, d_ft, d_fb,              &
458     !$omp             flux_e, flux_w, flux_n, flux_s, flux_t, flux_b)  &
459     !$omp     shared(ijkstart3, ijkend3, do_k, a_u_g)
460     
461           DO IJK = ijkstart3, ijkend3
462     
463              IF (FLOW_AT_E(IJK)) THEN
464     
465     ! Calculate convection-diffusion fluxes through each of the faces
466                 CALL GET_UCELL_GCFLUX_TERMS(flux_e, flux_w, flux_n, &
467                    flux_s, flux_t, flux_b, ijk)
468                 CALL GET_UCELL_GDIFF_TERMS(d_fe, d_fw, d_fn, d_fs, &
469                    d_ft, d_fb, ijk)
470     
471                 IPJK = IP_OF(IJK)
472                 IJPK = JP_OF(IJK)
473                 IMJK = IM_OF(IJK)
474                 IJMK = JM_OF(IJK)
475     
476     ! East face (i+1, j, k)
477                 IF (Flux_e >= ZERO) THEN
478                    A_U_G(IJK,E,0) = D_Fe
479                    A_U_G(IPJK,W,0) = D_Fe + Flux_e
480                 ELSE
481                    A_U_G(IJK,E,0) = D_Fe - Flux_e
482                    A_U_G(IPJK,W,0) = D_Fe
483                 ENDIF
484     ! West face (i, j, k)
485                 IF (.NOT.FLOW_AT_E(IMJK)) THEN
486                    IF (Flux_w >= ZERO) THEN
487                       A_U_G(IJK,W,0) = D_Fw + Flux_w
488                    ELSE
489                       A_U_G(IJK,W,0) = D_Fw
490                    ENDIF
491                 ENDIF
492     
493     
494     ! North face (i+1/2, j+1/2, k)
495                 IF (Flux_n >= ZERO) THEN
496                    A_U_G(IJK,N,0) = D_Fn
497                    A_U_G(IJPK,S,0) = D_Fn + Flux_n
498                 ELSE
499                    A_U_G(IJK,N,0) = D_Fn - Flux_n
500                    A_U_G(IJPK,S,0) = D_Fn
501                 ENDIF
502     ! South face (i+1/2, j-1/2, k)
503                 IF (.NOT.FLOW_AT_E(IJMK)) THEN
504                    IF (Flux_s >= ZERO) THEN
505                       A_U_G(IJK,S,0) = D_Fs + Flux_s
506                    ELSE
507                       A_U_G(IJK,S,0) = D_Fs
508                    ENDIF
509                 ENDIF
510     
511     
512                 IF (DO_K) THEN
513                    IJKP = KP_OF(IJK)
514                    IJKM = KM_OF(IJK)
515     
516     ! Top face (i+1/2, j, k+1/2)
517                    IF (Flux_t >= ZERO) THEN
518                       A_U_G(IJK,T,0) = D_Ft
519                       A_U_G(IJKP,B,0) = D_Ft + Flux_t
520                    ELSE
521                       A_U_G(IJK,T,0) = D_Ft - Flux_t
522                       A_U_G(IJKP,B,0) = D_Ft
523                    ENDIF
524     ! Bottom face (i+1/2, j, k-1/2)
525                    IF (.NOT.FLOW_AT_E(IJKM)) THEN
526                       IF (Flux_b >= ZERO) THEN
527                          A_U_G(IJK,B,0) = D_Fb + Flux_b
528                       ELSE
529                          A_U_G(IJK,B,0) = D_Fb
530                       ENDIF
531                    ENDIF
532                 ENDIF   ! end if (do_k)
533     
534              ENDIF   ! end if (flow_at_e)
535           ENDDO   !end do (ijk)
536     !$omp end parallel do
537     
538           RETURN
539           END SUBROUTINE STORE_A_U_G0
540     
541     
542     
543     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
544     !                                                                      C
545     !  Subroutine: STORE_A_U_GDC                                           C
546     !  Purpose: Use deferred correction method to solve the u-momentum     C
547     !  equation. This method combines first order upwind and a user        C
548     !  specified higher order method                                       C
549     !                                                                      C
550     !  Author: C. GUENTHER                                Date: 8-APR-99   C
551     !                                                                      C
552     !  Revision Number: 1                                                  C
553     !  Purpose: To incorporate Cartesian grid modifications                C
554     !  Author: Jeff Dietiker                              Date: 01-Jul-09  C
555     !                                                                      C
556     !                                                                      C
557     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
558           SUBROUTINE STORE_A_U_GDC(B_M)
559     
560     ! Modules
561     !---------------------------------------------------------------------//
562           USE compar, only: ijkstart3, ijkend3
563     
564           USE discretization, only: fpfoi_of
565     
566           USE fldvar, only: u_g
567     
568           USE function3, only: funijk3
569           USE functions, only: flow_at_e
570           USE functions, only: ip_of, jp_of, kp_of
571           USE functions, only: im_of, jm_of, km_of
572     
573           USE geometry, only: do_k
574     
575           USE indices, only: i_of, j_of, k_of
576     
577           USE param, only: dimension_3
578           USE param1, only: zero
579     
580           USE run, only: discretize, fpfoi
581           USE sendrecv3, only: send_recv3
582     
583           USE tmp_array, only: U => Array1, V => Array2, WW => Array3
584           USE tmp_array, only: tmp4
585           USE tmp_array, only: lock_tmp_array, unlock_tmp_array
586           USE tmp_array, only: lock_tmp4_array, unlock_tmp4_array
587     
588           USE xsi, only: calc_xsi
589           USE xsi_array, only: xsi_e, xsi_n, xsi_t
590           USE xsi_array, only: lock_xsi_array, unlock_xsi_array
591     
592           IMPLICIT NONE
593     
594     ! Dummy arguments
595     !---------------------------------------------------------------------//
596     ! Vector b_m
597           DOUBLE PRECISION, INTENT(INOUT) :: B_m(DIMENSION_3)
598     
599     ! Local variables
600     !---------------------------------------------------------------------//
601     ! Indices
602           INTEGER :: I, J, K, IJK
603           INTEGER :: IMJK, IPJK, IJMK, IJPK, IJKM, IJKP
604           INTEGER :: IJK4, IPPP, IPPP4, JPPP, JPPP4, KPPP, KPPP4
605           INTEGER :: IMMM, IMMM4, JMMM, JMMM4, KMMM, KMMM4
606     ! indication for shear
607           INTEGER :: incr
608     ! Deferred correction contribution from high order method
609           DOUBLE PRECISION :: MOM_HO
610     ! low order approximation
611           DOUBLE PRECISION :: MOM_LO
612     ! convection factor at each face
613           DOUBLE PRECISION :: flux_e, flux_w, flux_n, flux_s
614           DOUBLE PRECISION :: flux_t, flux_b
615     ! deferred correction contribution from each face
616           DOUBLE PRECISION :: EAST_DC
617           DOUBLE PRECISION :: WEST_DC
618           DOUBLE PRECISION :: NORTH_DC
619           DOUBLE PRECISION :: SOUTH_DC
620           DOUBLE PRECISION :: TOP_DC
621           DOUBLE PRECISION :: BOTTOM_DC
622     
623     ! temporary use of global arrays:
624     ! array1 (locally u)  - the x directional velocity
625     !      DOUBLE PRECISION :: U(DIMENSION_3)
626     ! array2 (locally v)  - the y directional velocity
627     !      DOUBLE PRECISION :: V(DIMENSION_3)
628     ! array3 (locally ww) - the z directional velocity
629     !      DOUBLE PRECISION :: WW(DIMENSION_3)
630     !---------------------------------------------------------------------//
631     
632           call lock_tmp4_array
633           call lock_tmp_array
634           call lock_xsi_array
635     
636           CALL GET_UCELL_GVTERMS(U, V, WW)
637     
638     ! Send recv the third ghost layer
639           IF (FPFOI) THEN
640              Do IJK = ijkstart3, ijkend3
641                 I = I_OF(IJK)
642                 J = J_OF(IJK)
643                 K = K_OF(IJK)
644                 IJK4 = funijk3(I,J,K)
645                 TMP4(IJK4) = U_G(IJK)
646              ENDDO
647              CALL send_recv3(tmp4)
648           ENDIF
649     
650     ! shear indicator:
651           incr=1
652           CALL CALC_XSI (DISCRETIZE(3), U_G, U, V, WW, XSI_E, XSI_N, &
653                          XSI_T, incr)
654     
655     
656     !!!$omp      parallel do                                              &
657     !!!$omp&     private(IJK, IPJK, IJPK, IJKP, IMJK, IJMK, IJKM,         &
658     !!!$omp&             flux_e, flux_w, flux_n, flux_s, flux_t, flux_b,  &
659     !!!$omp&             MOM_HO, MOM_LO, EAST_DC, WEST_DC, NORTH_DC,      &
660     !!!$omp&             SOUTH_DC, TOP_DC, BOTTOM_DC)
661           DO IJK = ijkstart3, ijkend3
662     
663              IF (FLOW_AT_E(IJK)) THEN
664     
665     ! Calculate convection fluxes through each of the faces
666                 CALL GET_UCELL_GCFLUX_TERMS(flux_e, flux_w, flux_n, &
667                    flux_s, flux_t, flux_b, ijk)
668     
669                 IPJK = IP_OF(IJK)
670                 IMJK = IM_OF(IJK)
671                 IJMK = JM_OF(IJK)
672                 IJKM = KM_OF(IJK)
673                 IJPK = JP_OF(IJK)
674                 IJKP = KP_OF(IJK)
675     
676     ! Third Ghost layer information
677                 IPPP  = IP_OF(IP_OF(IPJK))
678                 IPPP4 = funijk3(I_OF(IPPP), J_OF(IPPP), K_OF(IPPP))
679                 IMMM  = IM_OF(IM_OF(IMJK))
680                 IMMM4 = funijk3(I_OF(IMMM), J_OF(IMMM), K_OF(IMMM))
681                 JPPP  = JP_OF(JP_OF(IJPK))
682                 JPPP4 = funijk3(I_OF(JPPP), J_OF(JPPP), K_OF(JPPP))
683                 JMMM  = JM_OF(JM_OF(IJMK))
684                 JMMM4 = funijk3(I_OF(JMMM), J_OF(JMMM), K_OF(JMMM))
685                 KPPP  = KP_OF(KP_OF(IJKP))
686                 KPPP4 = funijk3(I_OF(KPPP), J_OF(KPPP), K_OF(KPPP))
687                 KMMM  = KM_OF(KM_OF(IJKM))
688                 KMMM4 = funijk3(I_OF(KMMM), J_OF(KMMM), K_OF(KMMM))
689     
690     
691     ! East face (i+1, j, k)
692                 IF(U(IJK) >= ZERO)THEN
693                    MOM_LO = U_G(IJK)
694                    IF (FPFOI) MOM_HO = FPFOI_OF(U_G(IPJK), U_G(IJK), &
695                                        U_G(IMJK), U_G(IM_OF(IMJK)))
696                 ELSE
697                    MOM_LO = U_G(IPJK)
698                    IF (FPFOI) MOM_HO = FPFOI_OF(U_G(IJK), U_G(IPJK), &
699                                        U_G(IP_OF(IPJK)), TMP4(IPPP4))
700                 ENDIF
701                 IF (.NOT. FPFOI) MOM_HO = XSI_E(IJK)*U_G(IPJK)+ &
702                                           (1.0-XSI_E(IJK))*U_G(IJK)
703                 EAST_DC = Flux_e *(MOM_LO - MOM_HO)
704     
705     ! West face (i, j, k)
706                 IF(U(IMJK) >= ZERO)THEN
707                    MOM_LO = U_G(IMJK)
708                    IF (FPFOI) MOM_HO = FPFOI_OF(U_G(IJK), U_G(IMJK), &
709                                        U_G(IM_OF(IMJK)), TMP4(IMMM4))
710                 ELSE
711                    MOM_LO = U_G(IJK)
712                    IF (FPFOI) MOM_HO = FPFOI_OF(U_G(IMJK), U_G(IJK), &
713                                        U_G(IPJK), U_G(IP_OF(IPJK)))
714                 ENDIF
715                 IF (.NOT. FPFOI) MOM_HO = XSI_E(IMJK)*U_G(IJK)+ &
716                                           (1.0-XSI_E(IMJK))*U_G(IMJK)
717                 WEST_DC = Flux_w * (MOM_LO - MOM_HO)
718     
719     
720     ! North face (i+1/2, j+1/2, k)
721                 IF(V(IJK) >= ZERO)THEN
722                    MOM_LO = U_G(IJK)
723                    IF (FPFOI) MOM_HO = FPFOI_OF(U_G(IJPK), U_G(IJK), &
724                                        U_G(IJMK), U_G(JM_OF(IJMK)))
725                 ELSE
726                    MOM_LO = U_G(IJPK)
727                    IF (FPFOI) MOM_HO = FPFOI_OF(U_G(IJK), U_G(IJPK), &
728                                        U_G(JP_OF(IJPK)), TMP4(JPPP4))
729                 ENDIF
730                 IF (.NOT. FPFOI) MOM_HO = XSI_N(IJK)*U_G(IJPK)+ &
731                                           (1.0-XSI_N(IJK))*U_G(IJK)
732                 NORTH_DC = Flux_n *(MOM_LO - MOM_HO)
733     
734     ! South face (i+1/2, j-1/2, k)
735                IF(V(IJMK) >= ZERO)THEN
736                    MOM_LO = U_G(IJMK)
737                    IF (FPFOI) MOM_HO = FPFOI_OF(U_G(IJK), U_G(IJMK), &
738                                        U_G(JM_OF(IJMK)), TMP4(JMMM4))
739                 ELSE
740                    MOM_LO = U_G(IJK)
741                    IF (FPFOI) MOM_HO = FPFOI_OF(U_G(IJMK), U_G(IJK), &
742                                        U_G(IJPK), U_G(JP_OF(IJPK)))
743                 ENDIF
744                 IF (.NOT. FPFOI) MOM_HO = XSI_N(IJMK)*U_G(IJK)+ &
745                                           (1.0-XSI_N(IJMK))*U_G(IJMK)
746                 SOUTH_DC = Flux_s *(MOM_LO - MOM_HO)
747     
748     
749                 IF (DO_K) THEN
750     ! Top face (i+1/2, j, k+1/2)
751                    IF(WW(IJK) >= ZERO)THEN
752                       MOM_LO = U_G(IJK)
753                       IF (FPFOI) MOM_HO = FPFOI_OF(U_G(IJKP), U_G(IJK), &
754                                           U_G(IJKM), U_G(KM_OF(IJKM)))
755                    ELSE
756                       MOM_LO = U_G(IJKP)
757                       IF (FPFOI) MOM_HO = FPFOI_OF(U_G(IJK), U_G(IJKP), &
758                                           U_G(KP_OF(IJKP)), TMP4(KPPP4))
759                    ENDIF
760                    IF (.NOT. FPFOI) MOM_HO = XSI_T(IJK)*U_G(IJKP)+ &
761                                              (1.0-XSI_T(IJK))*U_G(IJK)
762                    TOP_DC = Flux_t *(MOM_LO - MOM_HO)
763     
764     ! Bottom face (i+1/2, j, k-1/2)
765                    IF(WW(IJK) >= ZERO)THEN
766                       MOM_LO = U_G(IJKM)
767                       IF (FPFOI) MOM_HO = FPFOI_OF(U_G(IJK), U_G(IJKM), &
768                                           U_G(KM_OF(IJKM)), TMP4(KMMM4))
769                    ELSE
770                       MOM_LO = U_G(IJK)
771                       IF (FPFOI) MOM_HO = FPFOI_OF(U_G(IJKM), U_G(IJK), &
772                                           U_G(IJKP), U_G(KP_OF(IJKP)))
773                    ENDIF
774                    IF (.NOT. FPFOI) MOM_HO = XSI_T(IJKM)*U_G(IJK)+ &
775                                             (1.0-XSI_T(IJKM))*U_G(IJKM)
776                    BOTTOM_DC = Flux_b * (MOM_LO - MOM_HO)
777                 ELSE
778                    TOP_DC = ZERO
779                    BOTTOM_DC = ZERO
780                 ENDIF   ! end if (do_k)
781     
782     
783     ! CONTRIBUTION DUE TO DEFERRED CORRECTION
784                 B_M(IJK) = B_M(IJK)+WEST_DC-EAST_DC+SOUTH_DC-NORTH_DC+&
785                            BOTTOM_DC-TOP_DC
786     
787              ENDIF ! end if flow_at_e
788           ENDDO   ! end do ijk
789     
790           call unlock_tmp4_array
791           call unlock_tmp_array
792           call unlock_xsi_array
793     
794           RETURN
795           END SUBROUTINE STORE_A_U_GDC
796     
797     
798     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
799     !                                                                      C
800     !  Subroutine: STORE_A_U_g1                                            C
801     !  Purpose: Determine convection diffusion terms for U_g momentum eqs  C
802     !  The off-diagonal coefficients calculated here must be positive.     C
803     !  The center coefficient and the source vector are negative.          C
804     !  Implements higher order discretization.                             C
805     !  See source_u_g                                                      C
806     !                                                                      C
807     !  Author: M. Syamlal                                 Date: 20-MAR-97  C
808     !                                                                      C
809     !  Revision Number: 1                                                  C
810     !  Purpose: To incorporate Cartesian grid modifications                C
811     !  Author: Jeff Dietiker                              Date: 01-Jul-09  C
812     !                                                                      C
813     !                                                                      C
814     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
815           SUBROUTINE STORE_A_U_G1(A_U_G)
816     
817     ! Modules
818     !---------------------------------------------------------------------//
819           USE compar, only: ijkstart3, ijkend3
820           USE fldvar, only: u_g
821     
822           USE functions, only: flow_at_e
823           USE functions, only: ip_of, jp_of, kp_of
824           USE functions, only: im_of, jm_of, km_of
825     
826           USE geometry, only: do_k
827     
828           USE param, only: dimension_3, dimension_m
829           USE param1, only: one
830     
831           USE matrix, only: e, w, n, s, t, b
832     
833           USE run, only: discretize
834     
835           USE tmp_array, only: U => Array1, V => Array2, WW => Array3
836           USE tmp_array, only: lock_tmp_array, unlock_tmp_array
837     
838           USE xsi, only: calc_xsi
839           USE xsi_array, only: xsi_e, xsi_n, xsi_t
840           USE xsi_array, only: lock_xsi_array, unlock_xsi_array
841           IMPLICIT NONE
842     
843     ! Dummy arguments
844     !---------------------------------------------------------------------//
845     ! Septadiagonal matrix A_U_g
846           DOUBLE PRECISION, INTENT(INOUT) :: A_U_g(DIMENSION_3, -3:3, 0:DIMENSION_M)
847     
848     ! Local variables
849     !---------------------------------------------------------------------//
850     ! Indices
851           INTEGER :: IJK, IPJK, IMJK, IJPK, IJMK, IJKP, IJKM
852     ! indicator for shear
853           INTEGER :: incr
854     ! Diffusion parameter
855           DOUBLE PRECISION :: d_fe, d_fw, d_fn, d_fs, d_ft, d_fb
856     ! Face mass flux
857           DOUBLE PRECISION :: Flux_e, flux_w, flux_n, flux_s
858           DOUBLE PRECISION :: flux_t, flux_b
859     
860     ! temporary use of global arrays:
861     ! array1 (locally u)  - the x directional velocity
862     !      DOUBLE PRECISION :: U(DIMENSION_3)
863     ! array2 (locally v)  - the y directional velocity
864     !      DOUBLE PRECISION :: V(DIMENSION_3)
865     ! array3 (locally ww) - the z directional velocity
866     !      DOUBLE PRECISION :: WW(DIMENSION_3)
867     !---------------------------------------------------------------------//
868     
869           call lock_tmp_array
870           call lock_xsi_array
871     
872           CALL GET_UCELL_GVTERMS(U, V, WW)
873     
874     ! shear indicator:
875           incr=1
876           CALL CALC_XSI (DISCRETIZE(3), U_G, U, V, WW, XSI_E, XSI_N, &
877                          XSI_T, incr)
878     
879     
880     !!!$omp      parallel do                                                 &
881     !!!$omp&     private(IJK, IPJK, IJPK, IJKP, IMJK, IJMK, IJKM,            &
882     !!!$omp&             d_fe, d_fw, d_fn, d_fs, d_ft, d_fb,                 &
883     !!!$omp&             flux_e, flux_w, flux_n, flux_s, flux_t, flux_b)
884           DO IJK = ijkstart3, ijkend3
885     
886              IF (FLOW_AT_E(IJK)) THEN
887     
888     ! Calculate convection-diffusion fluxes through each of the faces
889                 CALL GET_UCELL_GCFLUX_TERMS(flux_e, flux_w, flux_n, &
890                    flux_s, flux_t, flux_b, ijk)
891                 CALL GET_UCELL_GDIFF_TERMS(d_fe, d_fw, d_fn, d_fs, &
892                    d_ft, d_fb, ijk)
893     
894                 IPJK = IP_OF(IJK)
895                 IMJK = IM_OF(IJK)
896                 IJPK = JP_OF(IJK)
897                 IJMK = JM_OF(IJK)
898     
899     ! East face (i+1, j, k)
900                 A_U_G(IJK,E,0) = D_Fe - XSI_E(IJK) * Flux_e
901                 A_U_G(IPJK,W,0) = D_Fe + (ONE - XSI_E(IJK)) * Flux_e
902     ! West face (i, j, k)
903                 IF (.NOT.FLOW_AT_E(IMJK)) THEN
904                    A_U_G(IJK,W,0) = D_Fw + (ONE - XSI_E(IMJK)) * Flux_w
905                 ENDIF
906     
907     
908     ! North face (i+1/2, j+1/2, k)
909                 A_U_G(IJK,N,0) = D_Fn - XSI_N(IJK) * Flux_n
910                 A_U_G(IJPK,S,0) = D_Fn + (ONE - XSI_N(IJK)) * Flux_n
911     ! South face (i+1/2, j-1/2, k)
912                 IF (.NOT.FLOW_AT_E(IJMK)) THEN
913                    A_U_G(IJK,S,0) = D_Fs + (ONE - XSI_N(IJMK)) * Flux_s
914                 ENDIF
915     
916     
917     ! Top face (i+1/2, j, k+1/2)
918                 IF (DO_K) THEN
919                    IJKP = KP_OF(IJK)
920                    IJKM = KM_OF(IJK)
921                    A_U_G(IJK,T,0) = D_Ft - XSI_T(IJK) * Flux_t
922                    A_U_G(IJKP,B,0) = D_Ft + (ONE - XSI_T(IJK)) * Flux_t
923     ! Bottom face (i+1/2, j, k-1/2)
924                    IF (.NOT.FLOW_AT_E(IJKM)) THEN
925                       A_U_G(IJK,B,0) = D_Fb + (ONE - XSI_T(IJKM)) * Flux_b
926                    ENDIF
927                 ENDIF   ! end if (do_k)
928     
929              ENDIF   ! end if flow_at_e
930           ENDDO   ! end do ijk
931     
932           call unlock_tmp_array
933           call unlock_xsi_array
934     
935           RETURN
936           END SUBROUTINE STORE_A_U_G1
937     
938