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