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