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