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

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