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

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