File: RELATIVE:/../../../mfix.git/model/conv_dif_w_s.f

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