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