File: N:\mfix\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 param
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,east) = D_FE
378           DF_GU(IJK,west) = D_FW
379           DF_GU(IJK,north) = D_FN
380           DF_GU(IJK,south) = D_FS
381           DF_GU(IJK,top) = D_FT
382           DF_GU(IJK,bottom) = 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     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
403     !                                                                      C
404     !  Subroutine: STORE_A_U_g0                                            C
405     !  Purpose: Determine convection diffusion terms for U_g momentum eqs. C
406     !  The off-diagonal coefficients calculated here must be positive.     C
407     !  The center coefficient and the source vector are negative. See      C
408     !  source_u_g.                                                         C
409     !  Implement FOUP discretization                                       C
410     !                                                                      C
411     !  Author: M. Syamlal                                 Date: 29-APR-96  C
412     !                                                                      C
413     !  Revision Number: 1                                                  C
414     !  Purpose: To incorporate Cartesian grid modifications                C
415     !  Author: Jeff Dietiker                              Date: 01-Jul-09  C
416     !                                                                      C
417     !                                                                      C
418     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
419           SUBROUTINE STORE_A_U_G0(A_U_G)
420     
421     ! Modules
422     !---------------------------------------------------------------------//
423           USE compar, only: ijkstart3, ijkend3
424     
425           USE functions, only: flow_at_e
426           USE functions, only: ip_of, jp_of, kp_of
427           USE functions, only: im_of, jm_of, km_of
428     
429           USE geometry, only: do_k
430     
431           USE param
432           USE param1, only: zero
433           IMPLICIT NONE
434     
435     ! Dummy arguments
436     !---------------------------------------------------------------------//
437     ! Septadiagonal matrix A_U_g
438           DOUBLE PRECISION, INTENT(INOUT) :: A_U_g(DIMENSION_3, -3:3, 0:DIMENSION_M)
439     
440     ! Local variables
441     !---------------------------------------------------------------------//
442     ! Indices
443           INTEGER :: IJK
444           INTEGER :: IMJK, IPJK, IJMK, IJPK, IJKM, IJKP
445     ! Face mass flux
446           DOUBLE PRECISION :: flux_e, flux_w, flux_n, flux_s
447           DOUBLE PRECISION :: flux_t, flux_b
448     ! Diffusion parameter
449           DOUBLE PRECISION :: D_fe, d_fw, d_fn, d_fs, d_ft, d_fb
450     
451     !---------------------------------------------------------------------//
452     
453     !$omp     parallel do default(none)                                &
454     !$omp     private(IJK, IPJK, IJPK, IJKP, IMJK, IJMK, IJKM,         &
455     !$omp             D_fe, d_fw, d_fn, d_fs, d_ft, d_fb,              &
456     !$omp             flux_e, flux_w, flux_n, flux_s, flux_t, flux_b)  &
457     !$omp     shared(ijkstart3, ijkend3, do_k, a_u_g)
458     
459           DO IJK = ijkstart3, ijkend3
460     
461              IF (FLOW_AT_E(IJK)) THEN
462     
463     ! Calculate convection-diffusion fluxes through each of the faces
464                 CALL GET_UCELL_GCFLUX_TERMS(flux_e, flux_w, flux_n, &
465                    flux_s, flux_t, flux_b, ijk)
466                 CALL GET_UCELL_GDIFF_TERMS(d_fe, d_fw, d_fn, d_fs, &
467                    d_ft, d_fb, ijk)
468     
469                 IPJK = IP_OF(IJK)
470                 IJPK = JP_OF(IJK)
471                 IMJK = IM_OF(IJK)
472                 IJMK = JM_OF(IJK)
473     
474     ! East face (i+1, j, k)
475                 IF (Flux_e >= ZERO) THEN
476                    A_U_G(IJK,east,0) = D_Fe
477                    A_U_G(IPJK,west,0) = D_Fe + Flux_e
478                 ELSE
479                    A_U_G(IJK,east,0) = D_Fe - Flux_e
480                    A_U_G(IPJK,west,0) = D_Fe
481                 ENDIF
482     ! West face (i, j, k)
483                 IF (.NOT.FLOW_AT_E(IMJK)) THEN
484                    IF (Flux_w >= ZERO) THEN
485                       A_U_G(IJK,west,0) = D_Fw + Flux_w
486                    ELSE
487                       A_U_G(IJK,west,0) = D_Fw
488                    ENDIF
489                 ENDIF
490     
491     
492     ! North face (i+1/2, j+1/2, k)
493                 IF (Flux_n >= ZERO) THEN
494                    A_U_G(IJK,north,0) = D_Fn
495                    A_U_G(IJPK,south,0) = D_Fn + Flux_n
496                 ELSE
497                    A_U_G(IJK,north,0) = D_Fn - Flux_n
498                    A_U_G(IJPK,south,0) = D_Fn
499                 ENDIF
500     ! South face (i+1/2, j-1/2, k)
501                 IF (.NOT.FLOW_AT_E(IJMK)) THEN
502                    IF (Flux_s >= ZERO) THEN
503                       A_U_G(IJK,south,0) = D_Fs + Flux_s
504                    ELSE
505                       A_U_G(IJK,south,0) = D_Fs
506                    ENDIF
507                 ENDIF
508     
509     
510                 IF (DO_K) THEN
511                    IJKP = KP_OF(IJK)
512                    IJKM = KM_OF(IJK)
513     
514     ! Top face (i+1/2, j, k+1/2)
515                    IF (Flux_t >= ZERO) THEN
516                       A_U_G(IJK,top,0) = D_Ft
517                       A_U_G(IJKP,bottom,0) = D_Ft + Flux_t
518                    ELSE
519                       A_U_G(IJK,top,0) = D_Ft - Flux_t
520                       A_U_G(IJKP,bottom,0) = D_Ft
521                    ENDIF
522     ! Bottom face (i+1/2, j, k-1/2)
523                    IF (.NOT.FLOW_AT_E(IJKM)) THEN
524                       IF (Flux_b >= ZERO) THEN
525                          A_U_G(IJK,bottom,0) = D_Fb + Flux_b
526                       ELSE
527                          A_U_G(IJK,bottom,0) = D_Fb
528                       ENDIF
529                    ENDIF
530                 ENDIF   ! end if (do_k)
531     
532              ENDIF   ! end if (flow_at_e)
533           ENDDO   !end do (ijk)
534     !$omp end parallel do
535     
536           RETURN
537           END SUBROUTINE STORE_A_U_G0
538     
539     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
540     !                                                                      C
541     !  Subroutine: STORE_A_U_GDC                                           C
542     !  Purpose: Use deferred correction method to solve the u-momentum     C
543     !  equation. This method combines first order upwind and a user        C
544     !  specified higher order method                                       C
545     !                                                                      C
546     !  Author: C. GUENTHER                                Date: 8-APR-99   C
547     !                                                                      C
548     !  Revision Number: 1                                                  C
549     !  Purpose: To incorporate Cartesian grid modifications                C
550     !  Author: Jeff Dietiker                              Date: 01-Jul-09  C
551     !                                                                      C
552     !                                                                      C
553     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
554           SUBROUTINE STORE_A_U_GDC(B_M)
555     
556     ! Modules
557     !---------------------------------------------------------------------//
558           USE compar, only: ijkstart3, ijkend3
559     
560           USE discretization, only: fpfoi_of
561     
562           USE fldvar, only: u_g
563     
564           USE function3, only: funijk3
565           USE functions, only: flow_at_e
566           USE functions, only: ip_of, jp_of, kp_of
567           USE functions, only: im_of, jm_of, km_of
568     
569           USE geometry, only: do_k
570     
571           USE indices, only: i_of, j_of, k_of
572     
573           USE param, only: dimension_3, dimension_4
574           USE param1, only: zero
575     
576           USE run, only: discretize, fpfoi
577           USE sendrecv3, only: send_recv3
578     
579           USE xsi, only: calc_xsi
580     
581           IMPLICIT NONE
582     
583     ! Dummy arguments
584     !---------------------------------------------------------------------//
585     ! Vector b_m
586           DOUBLE PRECISION, INTENT(INOUT) :: B_m(DIMENSION_3)
587     
588     ! Local variables
589     !---------------------------------------------------------------------//
590     ! Indices
591           INTEGER :: I, J, K, IJK
592           INTEGER :: IMJK, IPJK, IJMK, IJPK, IJKM, IJKP
593           INTEGER :: IJK4, IPPP, IPPP4, JPPP, JPPP4, KPPP, KPPP4
594           INTEGER :: IMMM, IMMM4, JMMM, JMMM4, KMMM, KMMM4
595     ! indication for shear
596           INTEGER :: incr
597     ! Deferred correction contribution from high order method
598           DOUBLE PRECISION :: MOM_HO
599     ! low order approximation
600           DOUBLE PRECISION :: MOM_LO
601     ! convection factor at each face
602           DOUBLE PRECISION :: flux_e, flux_w, flux_n, flux_s
603           DOUBLE PRECISION :: flux_t, flux_b
604     ! deferred correction contribution from each face
605           DOUBLE PRECISION :: EAST_DC
606           DOUBLE PRECISION :: WEST_DC
607           DOUBLE PRECISION :: NORTH_DC
608           DOUBLE PRECISION :: SOUTH_DC
609           DOUBLE PRECISION :: TOP_DC
610           DOUBLE PRECISION :: BOTTOM_DC
611     
612     ! temporary use of global arrays:
613     ! array1 (locally u)  - the x directional velocity
614           DOUBLE PRECISION :: U(DIMENSION_3)
615     ! array2 (locally v)  - the y directional velocity
616           DOUBLE PRECISION :: V(DIMENSION_3)
617     ! array3 (locally ww) - the z directional velocity
618           DOUBLE PRECISION :: WW(DIMENSION_3)
619     !---------------------------------------------------------------------
620           DOUBLE PRECISION :: TMP4(DIMENSION_4)
621           DOUBLE PRECISION, DIMENSION(DIMENSION_3) :: XSI_e, XSI_n, XSI_t
622     
623           CALL GET_UCELL_GVTERMS(U, V, WW)
624     
625     ! Send recv the third ghost layer
626           IF (FPFOI) THEN
627              Do IJK = ijkstart3, ijkend3
628                 I = I_OF(IJK)
629                 J = J_OF(IJK)
630                 K = K_OF(IJK)
631                 IJK4 = funijk3(I,J,K)
632                 TMP4(IJK4) = U_G(IJK)
633              ENDDO
634              CALL send_recv3(tmp4)
635           ENDIF
636     
637     ! shear indicator:
638           incr=1
639           CALL CALC_XSI (DISCRETIZE(3), U_G, U, V, WW, XSI_E, XSI_N, &
640                          XSI_T, incr)
641     
642     
643     !!!$omp      parallel do                                              &
644     !!!$omp&     private(IJK, IPJK, IJPK, IJKP, IMJK, IJMK, IJKM,         &
645     !!!$omp&             flux_e, flux_w, flux_n, flux_s, flux_t, flux_b,  &
646     !!!$omp&             MOM_HO, MOM_LO, EAST_DC, WEST_DC, NORTH_DC,      &
647     !!!$omp&             SOUTH_DC, TOP_DC, BOTTOM_DC)
648           DO IJK = ijkstart3, ijkend3
649     
650              IF (FLOW_AT_E(IJK)) THEN
651     
652     ! Calculate convection fluxes through each of the faces
653                 CALL GET_UCELL_GCFLUX_TERMS(flux_e, flux_w, flux_n, &
654                    flux_s, flux_t, flux_b, ijk)
655     
656                 IPJK = IP_OF(IJK)
657                 IMJK = IM_OF(IJK)
658                 IJMK = JM_OF(IJK)
659                 IJKM = KM_OF(IJK)
660                 IJPK = JP_OF(IJK)
661                 IJKP = KP_OF(IJK)
662     
663     ! Third Ghost layer information
664                 IPPP  = IP_OF(IP_OF(IPJK))
665                 IPPP4 = funijk3(I_OF(IPPP), J_OF(IPPP), K_OF(IPPP))
666                 IMMM  = IM_OF(IM_OF(IMJK))
667                 IMMM4 = funijk3(I_OF(IMMM), J_OF(IMMM), K_OF(IMMM))
668                 JPPP  = JP_OF(JP_OF(IJPK))
669                 JPPP4 = funijk3(I_OF(JPPP), J_OF(JPPP), K_OF(JPPP))
670                 JMMM  = JM_OF(JM_OF(IJMK))
671                 JMMM4 = funijk3(I_OF(JMMM), J_OF(JMMM), K_OF(JMMM))
672                 KPPP  = KP_OF(KP_OF(IJKP))
673                 KPPP4 = funijk3(I_OF(KPPP), J_OF(KPPP), K_OF(KPPP))
674                 KMMM  = KM_OF(KM_OF(IJKM))
675                 KMMM4 = funijk3(I_OF(KMMM), J_OF(KMMM), K_OF(KMMM))
676     
677     
678     ! East face (i+1, j, k)
679                 IF(U(IJK) >= ZERO)THEN
680                    MOM_LO = U_G(IJK)
681                    IF (FPFOI) MOM_HO = FPFOI_OF(U_G(IPJK), U_G(IJK), &
682                                        U_G(IMJK), U_G(IM_OF(IMJK)))
683                 ELSE
684                    MOM_LO = U_G(IPJK)
685                    IF (FPFOI) MOM_HO = FPFOI_OF(U_G(IJK), U_G(IPJK), &
686                                        U_G(IP_OF(IPJK)), TMP4(IPPP4))
687                 ENDIF
688                 IF (.NOT. FPFOI) MOM_HO = XSI_E(IJK)*U_G(IPJK)+ &
689                                           (1.0-XSI_E(IJK))*U_G(IJK)
690                 EAST_DC = Flux_e *(MOM_LO - MOM_HO)
691     
692     ! West face (i, j, k)
693                 IF(U(IMJK) >= ZERO)THEN
694                    MOM_LO = U_G(IMJK)
695                    IF (FPFOI) MOM_HO = FPFOI_OF(U_G(IJK), U_G(IMJK), &
696                                        U_G(IM_OF(IMJK)), TMP4(IMMM4))
697                 ELSE
698                    MOM_LO = U_G(IJK)
699                    IF (FPFOI) MOM_HO = FPFOI_OF(U_G(IMJK), U_G(IJK), &
700                                        U_G(IPJK), U_G(IP_OF(IPJK)))
701                 ENDIF
702                 IF (.NOT. FPFOI) MOM_HO = XSI_E(IMJK)*U_G(IJK)+ &
703                                           (1.0-XSI_E(IMJK))*U_G(IMJK)
704                 WEST_DC = Flux_w * (MOM_LO - MOM_HO)
705     
706     
707     ! North face (i+1/2, j+1/2, k)
708                 IF(V(IJK) >= ZERO)THEN
709                    MOM_LO = U_G(IJK)
710                    IF (FPFOI) MOM_HO = FPFOI_OF(U_G(IJPK), U_G(IJK), &
711                                        U_G(IJMK), U_G(JM_OF(IJMK)))
712                 ELSE
713                    MOM_LO = U_G(IJPK)
714                    IF (FPFOI) MOM_HO = FPFOI_OF(U_G(IJK), U_G(IJPK), &
715                                        U_G(JP_OF(IJPK)), TMP4(JPPP4))
716                 ENDIF
717                 IF (.NOT. FPFOI) MOM_HO = XSI_N(IJK)*U_G(IJPK)+ &
718                                           (1.0-XSI_N(IJK))*U_G(IJK)
719                 NORTH_DC = Flux_n *(MOM_LO - MOM_HO)
720     
721     ! South face (i+1/2, j-1/2, k)
722                IF(V(IJMK) >= ZERO)THEN
723                    MOM_LO = U_G(IJMK)
724                    IF (FPFOI) MOM_HO = FPFOI_OF(U_G(IJK), U_G(IJMK), &
725                                        U_G(JM_OF(IJMK)), TMP4(JMMM4))
726                 ELSE
727                    MOM_LO = U_G(IJK)
728                    IF (FPFOI) MOM_HO = FPFOI_OF(U_G(IJMK), U_G(IJK), &
729                                        U_G(IJPK), U_G(JP_OF(IJPK)))
730                 ENDIF
731                 IF (.NOT. FPFOI) MOM_HO = XSI_N(IJMK)*U_G(IJK)+ &
732                                           (1.0-XSI_N(IJMK))*U_G(IJMK)
733                 SOUTH_DC = Flux_s *(MOM_LO - MOM_HO)
734     
735     
736                 IF (DO_K) THEN
737     ! Top face (i+1/2, j, k+1/2)
738                    IF(WW(IJK) >= ZERO)THEN
739                       MOM_LO = U_G(IJK)
740                       IF (FPFOI) MOM_HO = FPFOI_OF(U_G(IJKP), U_G(IJK), &
741                                           U_G(IJKM), U_G(KM_OF(IJKM)))
742                    ELSE
743                       MOM_LO = U_G(IJKP)
744                       IF (FPFOI) MOM_HO = FPFOI_OF(U_G(IJK), U_G(IJKP), &
745                                           U_G(KP_OF(IJKP)), TMP4(KPPP4))
746                    ENDIF
747                    IF (.NOT. FPFOI) MOM_HO = XSI_T(IJK)*U_G(IJKP)+ &
748                                              (1.0-XSI_T(IJK))*U_G(IJK)
749                    TOP_DC = Flux_t *(MOM_LO - MOM_HO)
750     
751     ! Bottom face (i+1/2, j, k-1/2)
752                    IF(WW(IJK) >= ZERO)THEN
753                       MOM_LO = U_G(IJKM)
754                       IF (FPFOI) MOM_HO = FPFOI_OF(U_G(IJK), U_G(IJKM), &
755                                           U_G(KM_OF(IJKM)), TMP4(KMMM4))
756                    ELSE
757                       MOM_LO = U_G(IJK)
758                       IF (FPFOI) MOM_HO = FPFOI_OF(U_G(IJKM), U_G(IJK), &
759                                           U_G(IJKP), U_G(KP_OF(IJKP)))
760                    ENDIF
761                    IF (.NOT. FPFOI) MOM_HO = XSI_T(IJKM)*U_G(IJK)+ &
762                                             (1.0-XSI_T(IJKM))*U_G(IJKM)
763                    BOTTOM_DC = Flux_b * (MOM_LO - MOM_HO)
764                 ELSE
765                    TOP_DC = ZERO
766                    BOTTOM_DC = ZERO
767                 ENDIF   ! end if (do_k)
768     
769     
770     ! CONTRIBUTION DUE TO DEFERRED CORRECTION
771                 B_M(IJK) = B_M(IJK)+WEST_DC-EAST_DC+SOUTH_DC-NORTH_DC+&
772                            BOTTOM_DC-TOP_DC
773     
774              ENDIF ! end if flow_at_e
775           ENDDO   ! end do ijk
776     
777           RETURN
778           END SUBROUTINE STORE_A_U_GDC
779     
780     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
781     !                                                                      C
782     !  Subroutine: STORE_A_U_g1                                            C
783     !  Purpose: Determine convection diffusion terms for U_g momentum eqs  C
784     !  The off-diagonal coefficients calculated here must be positive.     C
785     !  The center coefficient and the source vector are negative.          C
786     !  Implements higher order discretization.                             C
787     !  See source_u_g                                                      C
788     !                                                                      C
789     !  Author: M. Syamlal                                 Date: 20-MAR-97  C
790     !                                                                      C
791     !  Revision Number: 1                                                  C
792     !  Purpose: To incorporate Cartesian grid modifications                C
793     !  Author: Jeff Dietiker                              Date: 01-Jul-09  C
794     !                                                                      C
795     !                                                                      C
796     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
797           SUBROUTINE STORE_A_U_G1(A_U_G)
798     
799     ! Modules
800     !---------------------------------------------------------------------//
801           USE compar, only: ijkstart3, ijkend3
802           USE fldvar, only: u_g
803     
804           USE functions, only: flow_at_e
805           USE functions, only: ip_of, jp_of, kp_of
806           USE functions, only: im_of, jm_of, km_of
807     
808           USE geometry, only: do_k
809     
810           USE param
811           USE param1, only: one
812     
813           USE run, only: discretize
814     
815           USE xsi, only: calc_xsi
816     
817           IMPLICIT NONE
818     
819     ! Dummy arguments
820     !---------------------------------------------------------------------//
821     ! Septadiagonal matrix A_U_g
822           DOUBLE PRECISION, INTENT(INOUT) :: A_U_g(DIMENSION_3, -3:3, 0:DIMENSION_M)
823     
824     ! Local variables
825     !---------------------------------------------------------------------//
826     ! Indices
827           INTEGER :: IJK, IPJK, IMJK, IJPK, IJMK, IJKP, IJKM
828     ! indicator for shear
829           INTEGER :: incr
830     ! Diffusion parameter
831           DOUBLE PRECISION :: d_fe, d_fw, d_fn, d_fs, d_ft, d_fb
832     ! Face mass flux
833           DOUBLE PRECISION :: Flux_e, flux_w, flux_n, flux_s
834           DOUBLE PRECISION :: flux_t, flux_b
835     
836     ! temporary use of global arrays:
837     ! array1 (locally u)  - the x directional velocity
838           DOUBLE PRECISION :: U(DIMENSION_3)
839     ! array2 (locally v)  - the y directional velocity
840           DOUBLE PRECISION :: V(DIMENSION_3)
841     ! array3 (locally ww) - the z directional velocity
842           DOUBLE PRECISION :: WW(DIMENSION_3)
843     !---------------------------------------------------------------------//
844           DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: TMP4
845           DOUBLE PRECISION, DIMENSION(DIMENSION_3) :: XSI_e, XSI_n, XSI_t
846     
847           CALL GET_UCELL_GVTERMS(U, V, WW)
848     
849     ! shear indicator:
850           incr=1
851           CALL CALC_XSI (DISCRETIZE(3), U_G, U, V, WW, XSI_E, XSI_N, &
852                          XSI_T, incr)
853     
854     !!!$omp      parallel do                                                 &
855     !!!$omp&     private(IJK, IPJK, IJPK, IJKP, IMJK, IJMK, IJKM,            &
856     !!!$omp&             d_fe, d_fw, d_fn, d_fs, d_ft, d_fb,                 &
857     !!!$omp&             flux_e, flux_w, flux_n, flux_s, flux_t, flux_b)
858           DO IJK = ijkstart3, ijkend3
859     
860              IF (FLOW_AT_E(IJK)) THEN
861     
862     ! Calculate convection-diffusion fluxes through each of the faces
863                 CALL GET_UCELL_GCFLUX_TERMS(flux_e, flux_w, flux_n, &
864                    flux_s, flux_t, flux_b, ijk)
865                 CALL GET_UCELL_GDIFF_TERMS(d_fe, d_fw, d_fn, d_fs, &
866                    d_ft, d_fb, ijk)
867     
868                 IPJK = IP_OF(IJK)
869                 IMJK = IM_OF(IJK)
870                 IJPK = JP_OF(IJK)
871                 IJMK = JM_OF(IJK)
872     
873     ! East face (i+1, j, k)
874                 A_U_G(IJK,east,0) = D_Fe - XSI_E(IJK) * Flux_e
875                 A_U_G(IPJK,west,0) = D_Fe + (ONE - XSI_E(IJK)) * Flux_e
876     ! West face (i, j, k)
877                 IF (.NOT.FLOW_AT_E(IMJK)) THEN
878                    A_U_G(IJK,west,0) = D_Fw + (ONE - XSI_E(IMJK)) * Flux_w
879                 ENDIF
880     
881     
882     ! North face (i+1/2, j+1/2, k)
883                 A_U_G(IJK,north,0) = D_Fn - XSI_N(IJK) * Flux_n
884                 A_U_G(IJPK,south,0) = D_Fn + (ONE - XSI_N(IJK)) * Flux_n
885     ! South face (i+1/2, j-1/2, k)
886                 IF (.NOT.FLOW_AT_E(IJMK)) THEN
887                    A_U_G(IJK,south,0) = D_Fs + (ONE - XSI_N(IJMK)) * Flux_s
888                 ENDIF
889     
890     
891     ! Top face (i+1/2, j, k+1/2)
892                 IF (DO_K) THEN
893                    IJKP = KP_OF(IJK)
894                    IJKM = KM_OF(IJK)
895                    A_U_G(IJK,top,0) = D_Ft - XSI_T(IJK) * Flux_t
896                    A_U_G(IJKP,bottom,0) = D_Ft + (ONE - XSI_T(IJK)) * Flux_t
897     ! Bottom face (i+1/2, j, k-1/2)
898                    IF (.NOT.FLOW_AT_E(IJKM)) THEN
899                       A_U_G(IJK,bottom,0) = D_Fb + (ONE - XSI_T(IJKM)) * Flux_b
900                    ENDIF
901                 ENDIF   ! end if (do_k)
902     
903              ENDIF   ! end if flow_at_e
904           ENDDO   ! end do ijk
905     
906           RETURN
907           END SUBROUTINE STORE_A_U_G1
908