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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Module name: CONV_DIF_V_g(A_m, B_m, IER)                            C
4     !  Purpose: Determine convection diffusion terms for V_g momentum eqs  C
5     !  The off-diagonal coefficients calculated here must be positive. The C
6     !  center coefficient and the source vector are negative;              C
7     !  See source_v_g                                                      C
8     !                                                                      C
9     !  Author: M. Syamlal                                 Date: 24-DEC-96  C
10     !                                                                      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 matrix, only: e, w, n, s, t, b
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,E) = D_FE
380           DF_GV(IJK,W) = D_FW
381           DF_GV(IJK,N) = D_FN
382           DF_GV(IJK,S) = D_FS
383           DF_GV(IJK,T) = D_FT
384           DF_GV(IJK,B) = 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, only: dimension_3, dimension_m
434           USE param1, only: zero
435           USE matrix, only: e, w, n, s, t, b
436     
437           IMPLICIT NONE
438     
439     ! Dummy arguments
440     !---------------------------------------------------------------------//
441     ! Septadiagonal matrix A_V_g
442           DOUBLE PRECISION, INTENT(INOUT) :: A_V_g(DIMENSION_3, -3:3, 0:DIMENSION_M)
443     
444     ! Local variables
445     !---------------------------------------------------------------------//
446     ! Indices
447           INTEGER :: IJK
448           INTEGER :: IMJK, IPJK, IJMK, IJPK, IJKM, IJKP
449     ! Face mass flux
450           DOUBLE PRECISION :: flux_e, flux_w, flux_n, flux_s
451           DOUBLE PRECISION :: flux_t, flux_b
452     ! Diffusion parameter
453           DOUBLE PRECISION :: D_fe, d_fw, d_fn, d_fs, d_ft, d_fb
454     
455     !---------------------------------------------------------------------//
456     
457     !$omp     parallel do default(none)                              &
458     !$omp     private(IJK, IMJK, IJMK, IPJK, IJPK, IJKM, IJKP,       &
459     !$omp             D_fe, d_fw, d_fn, d_fs, d_ft, d_fb,            &
460     !$omp             flux_e, flux_w, flux_n, flux_s, flux_t,        &
461     !$omp             flux_b)                                        &
462     !$omp     shared(ijkstart3, ijkend3, do_k, a_v_g)
463           DO IJK = ijkstart3, ijkend3
464     
465              IF (FLOW_AT_N(IJK)) THEN
466     
467     ! Calculate convection-diffusion fluxes through each of the faces
468                 CALL GET_VCELL_GCFLUX_TERMS(flux_e, flux_w, flux_n, &
469                    flux_s, flux_t, flux_b, ijk)
470                 CALL GET_VCELL_GDIFF_TERMS(d_fe, d_fw, d_fn, d_fs, &
471                    d_ft, d_fb, ijk)
472     
473                 IPJK = IP_OF(IJK)
474                 IJPK = JP_OF(IJK)
475                 IMJK = IM_OF(IJK)
476                 IJMK = JM_OF(IJK)
477     
478     ! East face (i+1/2, j+1/2, k)
479                 IF (Flux_e >= ZERO) THEN
480                    A_V_G(IJK,E,0) = D_Fe
481                    A_V_G(IPJK,W,0) = D_Fe + Flux_e
482                 ELSE
483                    A_V_G(IJK,E,0) = D_Fe - Flux_e
484                    A_V_G(IPJK,W,0) = D_Fe
485                 ENDIF
486     ! West face (i-1/2, j+1/2, k)
487                 IF (.NOT.FLOW_AT_N(IMJK)) THEN
488                    IF (Flux_w >= ZERO) THEN
489                       A_V_G(IJK,W,0) = D_Fw + Flux_w
490                    ELSE
491                       A_V_G(IJK,W,0) = D_Fw
492                    ENDIF
493                 ENDIF
494     
495     
496     ! North face (i, j+1, k)
497                 IF (Flux_n >= ZERO) THEN
498                    A_V_G(IJK,N,0) = D_Fn
499                    A_V_G(IJPK,S,0) = D_Fn + Flux_n
500                 ELSE
501                    A_V_G(IJK,N,0) = D_Fn - Flux_n
502                    A_V_G(IJPK,S,0) = D_Fn
503                 ENDIF
504     ! South face (i, j, k)
505                 IF (.NOT.FLOW_AT_N(IJMK)) THEN
506                    IF (Flux_s >= ZERO) THEN
507                       A_V_G(IJK,S,0) = D_Fs + Flux_s
508                    ELSE
509                       A_V_G(IJK,S,0) = D_Fs
510                    ENDIF
511                 ENDIF
512     
513     
514                 IF (DO_K) THEN
515                    IJKP = KP_OF(IJK)
516                    IJKM = KM_OF(IJK)
517     
518     ! Top face (i, j+1/2, k+1/2)
519                    IF (Flux_t >= ZERO) THEN
520                       A_V_G(IJK,T,0) = D_Ft
521                       A_V_G(IJKP,B,0) = D_Ft + Flux_t
522                    ELSE
523                       A_V_G(IJK,T,0) = D_Ft - Flux_t
524                       A_V_G(IJKP,B,0) = D_Ft
525                    ENDIF
526     ! Bottom face (i, j+1/2, k-1/2)
527                    IF (.NOT.FLOW_AT_N(IJKM)) THEN
528                       IF (Flux_b >= ZERO) THEN
529                          A_V_G(IJK,B,0) = D_Fb + Flux_b
530                       ELSE
531                          A_V_G(IJK,B,0) = D_Fb
532                       ENDIF
533                    ENDIF
534                 ENDIF   ! end if (do_k)
535     
536              ENDIF   ! end if flow_at_n
537           END DO   ! end do ijk
538     !$omp end parallel do
539     
540           RETURN
541           END SUBROUTINE STORE_A_V_G0
542     
543     
544     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
545     !                                                                      C
546     !  Subroutine: STORE_A_V_GDC                                           C
547     !  Purpose: Use deferred correction method to solve the v-momentum     C
548     !  equation. This method combines first order upwind and a user        C
549     !  specified higher order method                                       C
550     !                                                                      C
551     !  Author: C. GUENTHER                                Date: 8-APR-99   C
552     !                                                                      C
553     !  Revision Number: 1                                                  C
554     !  Purpose: To incorporate Cartesian grid modifications                C
555     !  Author: Jeff Dietiker                              Date: 01-Jul-09  C
556     !                                                                      C
557     !                                                                      C
558     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
559           SUBROUTINE STORE_A_V_GDC(B_M)
560     
561     ! Modules
562     !---------------------------------------------------------------------//
563           USE compar, only: ijkstart3, ijkend3
564     
565           USE discretization, only: fpfoi_of
566     
567           USE fldvar, only: v_g
568     
569           USE function3, only: funijk3
570           USE functions, only: flow_at_n
571           USE functions, only: ip_of, jp_of, kp_of
572           USE functions, only: im_of, jm_of, km_of
573     
574           USE geometry, only: do_k
575     
576           USE indices, only: i_of, j_of, k_of
577     
578           USE param, only: dimension_3
579           USE param1, only: zero
580     
581           USE run, only: discretize, fpfoi
582           USE sendrecv3, only: send_recv3
583     
584           USE tmp_array, only: U => Array1, V => Array2, WW => Array3
585           USE tmp_array, only: tmp4
586           USE tmp_array, only: lock_tmp_array, unlock_tmp_array
587           USE tmp_array, only: lock_tmp4_array, unlock_tmp4_array
588     
589           USE xsi, only: calc_xsi
590           USE xsi_array, only: xsi_e, xsi_n, xsi_t
591           USE xsi_array, only: lock_xsi_array, unlock_xsi_array
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 :: IPJK, IMJK, 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 corrction contribution form 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 contributions 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_VCELL_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) = V_G(IJK)
646              ENDDO
647              CALL send_recv3(tmp4)
648           ENDIF
649     
650     ! shear indicator: v-momentum
651           incr=2
652           CALL CALC_XSI (DISCRETIZE(4), V_G, U, V, WW, XSI_E, XSI_N, &
653                          XSI_T, incr)
654     
655     
656     !!!$omp      parallel do                                             &
657     !!!$omp&     private( I, J, K, IJK, IPJK, IMJK, IJPK, IJMK,          &
658     !!!$omp&             flux_e, flux_w, flux_n, flux_s, flux_b, flux_t, &
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_N(IJK)) THEN
664     
665     ! Calculate convection fluxes through each of the faces
666                 CALL GET_VCELL_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                 IJPK = JP_OF(IJK)
672                 IJMK = JM_OF(IJK)
673                 IJKP = KP_OF(IJK)
674                 IJKM = KM_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/2, j+1/2, k)
692                 IF(U(IJK) >= ZERO)THEN
693                    MOM_LO = V_G(IJK)
694                    IF (FPFOI) MOM_HO = FPFOI_OF(V_G(IPJK), V_G(IJK), &
695                                        V_G(IMJK), V_G(IM_OF(IMJK)))
696                 ELSE
697                    MOM_LO = V_G(IPJK)
698                     IF (FPFOI) MOM_HO = FPFOI_OF(V_G(IJK), V_G(IPJK), &
699                                         V_G(IP_OF(IPJK)), TMP4(IPPP4))
700                 ENDIF
701                 IF (.NOT. FPFOI)  MOM_HO = XSI_E(IJK)*V_G(IPJK)+&
702                                            (1.0-XSI_E(IJK))*V_G(IJK)
703                 EAST_DC = Flux_e*(MOM_LO-MOM_HO)
704     
705     
706     ! North face (i, j+1, k)
707                 IF(V(IJK) >= ZERO)THEN
708                    MOM_LO = V_G(IJK)
709                    IF (FPFOI) MOM_HO = FPFOI_OF(V_G(IJPK), V_G(IJK), &
710                                        V_G(IJMK), V_G(JM_OF(IJMK)))
711                 ELSE
712                    MOM_LO = V_G(IJPK)
713                    IF (FPFOI) MOM_HO = FPFOI_OF(V_G(IJK), V_G(IJPK), &
714                                        V_G(JP_OF(IJPK)), TMP4(JPPP4))
715                 ENDIF
716                 IF (.NOT. FPFOI) MOM_HO = XSI_N(IJK)*V_G(IJPK)+&
717                                           (1.0-XSI_N(IJK))*V_G(IJK)
718                 NORTH_DC = Flux_n*(MOM_LO-MOM_HO)
719     
720     
721     ! Top face (i, j+1/2, k+1/2)
722                 IF (DO_K) THEN
723                    IF(WW(IJK) >= ZERO)THEN
724                       MOM_LO = V_G(IJK)
725                       IF (FPFOI) MOM_HO = FPFOI_OF(V_G(IJKP), V_G(IJK), &
726                                           V_G(IJKM), V_G(KM_OF(IJKM)))
727                    ELSE
728                       MOM_LO = V_G(IJKP)
729                       IF (FPFOI) MOM_HO = FPFOI_OF(V_G(IJK), V_G(IJKP), &
730                                           V_G(KP_OF(IJKP)), TMP4(KPPP4))
731                    ENDIF
732                    IF (.NOT. FPFOI) MOM_HO = XSI_T(IJK)*V_G(IJKP)+&
733                                              (1.0-XSI_T(IJK))*V_G(IJK)
734                    TOP_DC = Flux_t*(MOM_LO-MOM_HO)
735                 ELSE
736                    TOP_DC = ZERO
737                 ENDIF
738     
739     
740     ! West face (i-1/2, j+1/2, k)
741                 IF(U(IMJK) >= ZERO)THEN
742                    MOM_LO = V_G(IMJK)
743                    IF (FPFOI) MOM_HO = FPFOI_OF(V_G(IJK), V_G(IMJK), &
744                                        V_G(IM_OF(IMJK)), TMP4(IMMM4))
745                 ELSE
746                    MOM_LO = V_G(IJK)
747                    IF (FPFOI) MOM_HO = FPFOI_OF(V_G(IMJK), V_G(IJK), &
748                                        V_G(IPJK), V_G(IP_OF(IPJK)))
749                 ENDIF
750                 IF (.NOT. FPFOI)MOM_HO = XSI_E(IMJK)*V_G(IJK)+&
751                                          (1.0-XSI_E(IMJK))*V_G(IMJK)
752                 WEST_DC = Flux_w*(MOM_LO-MOM_HO)
753     
754     
755     ! South face (i, j-1, k)
756                 IF(V(IJMK) >= ZERO)THEN
757                    MOM_LO = V_G(IJMK)
758                    IF (FPFOI) MOM_HO = FPFOI_OF(V_G(IJK), V_G(IJMK), &
759                                        V_G(JM_OF(IJMK)), TMP4(JMMM4))
760                 ELSE
761                    MOM_LO = V_G(IJK)
762                    IF (FPFOI) MOM_HO = FPFOI_OF(V_G(IJMK), V_G(IJK), &
763                                        V_G(IJPK), V_G(JP_OF(IJPK)))
764                 ENDIF
765                 IF (.NOT. FPFOI) MOM_HO = XSI_N(IJMK)*V_G(IJK)+&
766                                           (1.0-XSI_N(IJMK))*V_G(IJMK)
767                 SOUTH_DC = Flux_s*(MOM_LO-MOM_HO)
768     
769     
770     ! Bottom face (i, j+1/2, k-1/2)
771                 IF (DO_K) THEN
772                    IF(WW(IJK) >= ZERO)THEN
773                       MOM_LO = V_G(IJKM)
774                       IF (FPFOI) MOM_HO = FPFOI_OF(V_G(IJK), V_G(IJKM), &
775                                           V_G(KM_OF(IJKM)), TMP4(KMMM4))
776                    ELSE
777                       MOM_LO = V_G(IJK)
778                       IF (FPFOI) MOM_HO = FPFOI_OF(V_G(IJKM), V_G(IJK), &
779                                           V_G(IJKP), V_G(KP_OF(IJKP)))
780                    ENDIF
781                    IF (.NOT. FPFOI) MOM_HO = XSI_T(IJKM)*V_G(IJK)+&
782                                              (1.0-XSI_T(IJKM))*V_G(IJKM)
783                    BOTTOM_DC = Flux_b*(MOM_LO-MOM_HO)
784                 ELSE
785                    BOTTOM_DC = ZERO
786                 ENDIF  ! end if (do_k)
787     
788     
789     ! CONTRIBUTION DUE TO DEFERRED CORRECTION
790                 B_M(IJK) = B_M(IJK)+WEST_DC-EAST_DC+SOUTH_DC-NORTH_DC&
791                                     +BOTTOM_DC-TOP_DC
792     
793              ENDIF   !end if flow_at_n
794           ENDDO   ! end do ijk
795     
796           call unlock_tmp4_array
797           call unlock_tmp_array
798           call unlock_xsi_array
799     
800           RETURN
801           END SUBROUTINE STORE_A_V_GDC
802     
803     
804     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
805     !                                                                      C
806     !  Subroutine: STORE_A_V_g1                                            C
807     !  Purpose: Determine convection diffusion terms for V_g momentum eqs  C
808     !  The off-diagonal coefficients calculated here must be positive.     C
809     !  The center coefficient and the source vector are negative.          C
810     !  Implements higher order discretization.                             C
811     !  See source_v_g                                                      C
812     !                                                                      C
813     !  Author: M. Syamlal                                 Date: 20-MAR-97  C
814     !                                                                      C
815     !  Revision Number: 1                                                  C
816     !  Purpose: To incorporate Cartesian grid modifications                C
817     !  Author: Jeff Dietiker                              Date: 01-Jul-09  C
818     !                                                                      C
819     !                                                                      C
820     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
821           SUBROUTINE STORE_A_V_G1(A_V_G)
822     
823     ! Modules
824     !---------------------------------------------------------------------//
825           USE compar, only: ijkstart3, ijkend3
826           USE fldvar, only: v_g
827     
828           USE functions, only: flow_at_n
829           USE functions, only: ip_of, jp_of, kp_of
830           USE functions, only: im_of, jm_of, km_of
831     
832           USE geometry, only: do_k
833     
834           USE param, only: dimension_3
835           USE param1, only: one
836     
837           USE matrix, only: e, w, n, s, t, b
838     
839           USE run, only: discretize
840     
841           USE tmp_array, only: U => Array1, V => Array2, WW => Array3
842           USE tmp_array, only: lock_tmp_array, unlock_tmp_array
843     
844           USE xsi, only: calc_xsi
845           USE xsi_array, only: xsi_e, xsi_n, xsi_t
846           USE xsi_array, only: lock_xsi_array, unlock_xsi_array
847           IMPLICIT NONE
848     
849     ! Dummy arguments
850     !---------------------------------------------------------------------//
851     ! Septadiagonal matrix A_V_g
852           DOUBLE PRECISION, INTENT(INOUT) :: A_V_g(DIMENSION_3, -3:3)
853     
854     ! Local variables
855     !---------------------------------------------------------------------//
856     ! Indices
857           INTEGER :: IJK, IPJK, IMJK, IJPK, IJMK, IJKP, IJKM
858     ! indicator for shear
859           INTEGER :: incr
860     ! Diffusion parameter
861           DOUBLE PRECISION :: d_fe, d_fw, d_fn, d_fs, d_ft, d_fb
862     ! Face mass flux
863           DOUBLE PRECISION :: Flux_e, flux_w, flux_n, flux_s
864           DOUBLE PRECISION :: flux_t, flux_b
865     
866     ! temporary use of global arrays:
867     ! array1 (locally u)  - the x directional velocity
868     !      DOUBLE PRECISION :: U(DIMENSION_3)
869     ! array2 (locally v)  - the y directional velocity
870     !      DOUBLE PRECISION :: V(DIMENSION_3)
871     ! array3 (locally ww) - the z directional velocity
872     !      DOUBLE PRECISION :: WW(DIMENSION_3)
873     !---------------------------------------------------------------------//
874     
875           call lock_tmp_array
876           call lock_xsi_array
877     
878           CALL GET_VCELL_GVTERMS(U, V, WW)
879     
880     ! shear indicator: y-momentum
881           incr=2
882           CALL CALC_XSI (DISCRETIZE(4), V_G, U, V, WW, XSI_E, XSI_N, &
883              XSI_T, incr)
884     
885     
886     !!!$omp      parallel do                                             &
887     !!!$omp&     private(IJK, IMJK, IJMK, IPJK, IJPK, IJKM, IJKP,        &
888     !!!$omp&             d_fe, d_fw, d_fn, d_fs, d_ft, d_fb,             &
889     !!!$omp&             flux_e, flux_w, flux_n, flux_s, flux_t, flux_b)
890           DO IJK = ijkstart3, ijkend3
891     
892              IF (FLOW_AT_N(IJK)) THEN
893     
894     ! Calculate convection-diffusion fluxes through each of the faces
895                 CALL GET_VCELL_GCFLUX_TERMS(flux_e, flux_w, flux_n, &
896                    flux_s, flux_t, flux_b, ijk)
897                 CALL GET_VCELL_GDIFF_TERMS(d_fe, d_fw, d_fn, d_fs, &
898                    d_ft, d_fb, ijk)
899     
900                 IMJK = IM_OF(IJK)
901                 IJMK = JM_OF(IJK)
902                 IPJK = IP_OF(IJK)
903                 IJPK = JP_OF(IJK)
904     
905     ! East face (i+1/2, j+1/2, k)
906                 A_V_G(IJK,E) = D_Fe - XSI_E(IJK)*Flux_e
907                 A_V_G(IPJK,W) = D_Fe + (ONE - XSI_E(IJK))*Flux_e
908     ! West face (i-1/2, j+1/2, k)
909                 IF (.NOT.FLOW_AT_N(IMJK)) THEN
910                    A_V_G(IJK,W) = D_Fw + (ONE - XSI_E(IMJK))*Flux_w
911                 ENDIF
912     
913     
914     ! North face (i, j+1, k)
915                 A_V_G(IJK,N) = D_Fn - XSI_N(IJK)*Flux_n
916                 A_V_G(IJPK,S) = D_Fn + (ONE - XSI_N(IJK))*Flux_n
917     ! South face (i, j, k)
918                 IF (.NOT.FLOW_AT_N(IJMK)) THEN
919                    A_V_G(IJK,S) = D_Fs + (ONE - XSI_N(IJMK))*Flux_s
920                 ENDIF
921     
922     
923                 IF (DO_K) THEN
924                    IJKP = KP_OF(IJK)
925                    IJKM = KM_OF(IJK)
926     ! Top face (i, j+1/2, k+1/2)
927                    A_V_G(IJK,T) = D_Ft - XSI_T(IJK)*Flux_t
928                    A_V_G(IJKP,B) = D_Ft + (ONE - XSI_T(IJK))*Flux_t
929     ! Bottom face (i, j+1/2, k-1/2)
930                    IF (.NOT.FLOW_AT_N(IJKM)) THEN
931                       A_V_G(IJK,B) = D_Fb + (ONE - XSI_T(IJKM))*Flux_b
932                    ENDIF
933                 ENDIF   ! end if (do_k)
934     
935              ENDIF   ! end if (flow_at_n)
936           ENDDO   ! end do (ijk)
937     
938           call unlock_tmp_array
939           call unlock_xsi_array
940     
941           RETURN
942           END SUBROUTINE STORE_A_V_G1
943     
944