File: N:\mfix\model\conv_dif_v_g.f

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