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

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