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

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Subroutine: CONV_DIF_PHI                                            C
4     !  Purpose: Determine convection diffusion terms for a scalar eqn.     C
5     !  The off-diagonal coefficients calculated here must be positive.     C
6     !  The center coefficient and the source vector are negative;          C
7     !  See source_phi                                                      C
8     !                                                                      C
9     !  Diffusion at the flow boundaries is prevented by setting the        C
10     !  diffusion coefficients at boundary cells to zero and then using a   C
11     !  harmonic average to calculate the boundary diffusivity.  The value  C
12     !  diffusivities at the boundaries are checked in check_data_30.       C
13     !  Ensure that harmonic avergaing is used in this routine.             C
14     !  See source_phi                                                      C
15     !                                                                      C
16     !  Author: M. Syamlal                                 Date: 21-APR-97  C
17     !                                                                      C
18     !                                                                      C
19     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
20     
21           SUBROUTINE CONV_DIF_PHI(PHI, DIF, DISC, UF, VF, WF, &
22                                   Flux_E, Flux_N, Flux_T, M, A_M, B_M)
23     
24     
25     ! Modules
26     !---------------------------------------------------------------------//
27           USE param, only: dimension_3, dimension_m
28           USE run, only: def_cor
29           IMPLICIT NONE
30     
31     ! Dummy arguments
32     !---------------------------------------------------------------------//
33     ! Scalar
34           DOUBLE PRECISION Phi(DIMENSION_3)
35     ! Gamma -- diffusion coefficient
36           DOUBLE PRECISION Dif(DIMENSION_3)
37     ! Discretization index
38           INTEGER, INTENT(IN) :: Disc
39     ! Velocity components
40           DOUBLE PRECISION, INTENT(IN) :: Uf(DIMENSION_3)
41           DOUBLE PRECISION, INTENT(IN) :: Vf(DIMENSION_3)
42           DOUBLE PRECISION, INTENT(IN) :: Wf(DIMENSION_3)
43     ! Mass flux components
44           DOUBLE PRECISION, INTENT(IN) :: Flux_E(DIMENSION_3)
45           DOUBLE PRECISION, INTENT(IN) :: Flux_N(DIMENSION_3)
46           DOUBLE PRECISION, INTENT(IN) :: Flux_T(DIMENSION_3)
47     ! Phase index
48           INTEGER, INTENT(IN) :: M
49     ! Septadiagonal matrix A_m
50           DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
51     ! Vector b_m
52           DOUBLE PRECISION, INTENT(INOUT) :: B_m(DIMENSION_3, 0:DIMENSION_M)
53     !---------------------------------------------------------------------//
54     
55     ! DEFERRED CORRECTION IS USED WITH THE SCALAR TRANSPORT EQN.
56           IF(DEF_COR)THEN
57              CALL CONV_DIF_PHI0(PHI, DIF, UF, VF, WF, &
58                 Flux_E, Flux_N, Flux_T, M, A_M)
59              IF (DISC > 1) CALL CONV_DIF_PHI_DC(PHI, DISC, UF, VF, WF,&
60                              Flux_E, Flux_N, Flux_T, M, B_M)
61           ELSE
62     
63     ! DO NOT USE DEFERRED CORRECTION TO SOLVE THE SCALAR TRANSPORT EQN.
64              IF (DISC == 0) THEN
65                 CALL CONV_DIF_PHI0(PHI, DIF, UF, VF, WF, &
66                    Flux_E, Flux_N, Flux_T, M, A_M)
67              ELSE
68                 CALL CONV_DIF_PHI1(PHI, DIF, DISC, UF, VF, WF, &
69                    Flux_E, Flux_N, Flux_T, M, A_M)
70              ENDIF
71           ENDIF
72     
73           CALL DIF_PHI_IS (DIF, A_M, M)
74     
75           RETURN
76           END SUBROUTINE CONV_DIF_PHI
77     
78     
79     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
80     !                                                                      C
81     !  Purpose: Calculate the components of diffusive flux through the     C
82     !  faces of a scalar cell. Note the fluxes are calculated at           C
83     !  all faces regardless of fuid_at condition of the west, south        C
84     !  or bottom cell.                                                     C
85     !                                                                      C
86     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
87           SUBROUTINE GET_PHICELL_DIFF_TERMS(Dif, D_FE, D_FW, D_FN, D_FS, &
88              D_FT, D_FB, IJK)
89     
90     ! Modules
91     !---------------------------------------------------------------------//
92           USE cutcell, only: cut_treatment_at, cut_cell_at
93     
94           USE functions, only: fluid_at
95           USE functions, only: east_of, north_of, top_of
96           USE functions, only: west_of, south_of, bottom_of
97           USE functions, only: im_of, jm_of, km_of
98           USE functions, only: ip_of, jp_of, kp_of
99     
100           USE fun_avg, only: avg_x_h, avg_y_h, avg_z_h
101     
102           USE geometry, only: odx_e, ody_n, odz_t
103           USE geometry, only: do_k
104           USE geometry, only: ox
105           USE geometry, only: dx, dy, dz
106           USE geometry, only: ayz, axz, axy
107     
108           USE indices, only: i_of, j_of, k_of
109           USE indices, only: im1, jm1, km1
110     
111           USE param, only: dimension_3
112           IMPLICIT NONE
113     
114     ! Dummy arguments
115     !---------------------------------------------------------------------//
116     ! Gamma -- diffusion coefficient
117           DOUBLE PRECISION, INTENT(IN) :: Dif(DIMENSION_3)
118     
119     ! fluxes through faces of given ijk u-momentum cell
120           DOUBLE PRECISION, INTENT(OUT) :: d_fe, d_fw
121           DOUBLE PRECISION, INTENT(OUT) :: d_fn, d_fs
122           DOUBLE PRECISION, INTENT(OUT) :: d_ft, d_fb
123     ! ijk index
124           INTEGER, INTENT(IN) :: ijk
125     
126     ! Local variables
127     !---------------------------------------------------------------------//
128     ! indices
129           INTEGER :: imjk, ijmk, ijkm
130           INTEGER :: ipjk, ijpk, ijkp
131           INTEGER :: i, j, k, im, jm, km
132           INTEGER :: ijke, ijkw, ijkn, ijks, ijkt, ijkb
133     
134     ! area terms
135           DOUBLE PRECISION :: C_AE, C_AW, C_AN, C_AS, C_AT, C_AB
136     !---------------------------------------------------------------------//
137     
138           IMJK = IM_OF(IJK)
139           IJMK = JM_OF(IJK)
140           IJKM = KM_OF(IJK)
141     
142     
143           I = I_OF(IJK)
144           J = J_OF(IJK)
145           K = K_OF(IJK)
146           IM = IM1(I)
147           JM = JM1(J)
148           KM = KM1(K)
149     
150           IJKE = EAST_OF(IJK)
151           IJKN = NORTH_OF(IJK)
152           IJKS = SOUTH_OF(IJK)
153           IJKW = WEST_OF(IJK)
154     
155           C_AE = ODX_E(I)*AYZ(IJK)
156           C_AW = ODX_E(IM)*AYZ(IMJK)
157           C_AN = ODY_N(J)*AXZ(IJK)
158           C_AS = ODY_N(JM)*AXZ(IJMK)
159           C_AT = OX(I)*ODZ_T(K)*AXY(IJK)
160           C_AB = OX(I)*ODZ_T(KM)*AXY(IJKM)
161     
162           IF(CUT_TREATMENT_AT(IJK).AND.CUT_CELL_AT(IJK)) THEN
163              IPJK = IP_OF(IJK)
164              IJPK = JP_OF(IJK)
165              IJKP = KP_OF(IJK)
166     
167              IF (.NOT.FLUID_AT(IPJK)) C_AE = ODX_E(I)*DY(J)*DZ(K)
168              IF (.NOT.FLUID_AT(IMJK)) C_AW = ODX_E(IM)*DY(J)*DZ(K)
169              IF (.NOT.FLUID_AT(IJPK)) C_AN = ODY_N(J)*DX(I)*DZ(K)
170              IF (.NOT.FLUID_AT(IJMK)) C_AS = ODY_N(JM)*DX(I)*DZ(K)
171              IF (.NOT.FLUID_AT(IJKP)) C_AT = OX(I)*ODZ_T(K)*DX(I)*DY(J)
172              IF (.NOT.FLUID_AT(IJKM)) C_AB = OX(I)*ODZ_T(KM)*DX(I)*DY(J)
173           ENDIF
174     
175     ! East face (i+1/2, j, k)
176           D_Fe = AVG_X_H(DIF(IJK),DIF(IJKE),I)*C_AE
177     
178     ! West face (i-1/1, j, k)
179           D_FW = AVG_X_H(DIF(IJKW),DIF(IJK),IM)*C_AW
180     
181     ! North face (i, j+1/2, k)
182           D_FN = AVG_Y_H(DIF(IJK),DIF(IJKN),J)*C_AN
183     
184     ! South face (i, j-1/2, k)
185           D_FS = AVG_Y_H(DIF(IJKS),DIF(IJK),JM)*C_AS
186     
187           IF (DO_K) THEN
188              IJKT = TOP_OF(IJK)
189              IJKB = BOTTOM_OF(IJK)
190     
191     ! Top face (i, j, k+1/2)
192              D_FT = AVG_Z_H(DIF(IJK),DIF(IJKT),K)*C_AT
193     ! Bottom face (i, j, k-1/2)
194              D_FB = AVG_Z_H(DIF(IJKB),DIF(IJK),KM)*C_AB
195           ENDIF
196     
197           RETURN
198           END SUBROUTINE GET_PHICELL_DIFF_TERMS
199     
200     
201     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
202     !                                                                      C
203     !  Subroutine: CONV_DIF_PHI0                                           C
204     !  Purpose: Determine convection diffusion terms for Phi balance       C
205     !  The off-diagonal coefficients calculated here must be positive.     C
206     !  The center coefficient and the source vector are negative;          C
207     !  See source_phi                                                      C
208     !  Implement FOUP discretization                                       C
209     !                                                                      C
210     !  Author: M. Syamlal                                 Date: 21-APR-97  C
211     !                                                                      C
212     !                                                                      C
213     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
214           SUBROUTINE CONV_DIF_PHI0(PHI, DIF, UF, VF, WF, &
215                                    Flux_E, Flux_N, Flux_T, M, A_M)
216     
217     ! Modules
218     !---------------------------------------------------------------------//
219           USE compar, only: ijkstart3, ijkend3
220     
221           USE functions, only: fluid_at
222           USE functions, only: ip_of, jp_of, kp_of
223           USE functions, only: im_of, jm_of, km_of
224           USE geometry, only: do_k
225           USE param
226           USE param1, only: zero
227           IMPLICIT NONE
228     
229     ! Dummy arguments
230     !---------------------------------------------------------------------//
231     ! Scalar
232           DOUBLE PRECISION, INTENT(IN) :: Phi(DIMENSION_3)
233     ! Gamma -- diffusion coefficient
234           DOUBLE PRECISION, INTENT(IN) :: Dif(DIMENSION_3)
235     ! Velocity components
236           DOUBLE PRECISION, INTENT(IN) :: Uf(DIMENSION_3)
237           DOUBLE PRECISION, INTENT(IN) :: Vf(DIMENSION_3)
238           DOUBLE PRECISION, INTENT(IN) :: Wf(DIMENSION_3)
239     ! Mass flux components
240           DOUBLE PRECISION, INTENT(IN) :: Flux_E(DIMENSION_3)
241           DOUBLE PRECISION, INTENT(IN) :: Flux_N(DIMENSION_3)
242           DOUBLE PRECISION, INTENT(IN) :: Flux_T(DIMENSION_3)
243     ! Phase index
244           INTEGER, INTENT(IN) :: M
245     ! Septadiagonal matrix A_m
246           DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
247     
248     ! Local variables
249     !---------------------------------------------------------------------//
250     ! Indices
251           INTEGER :: IJK
252           INTEGER :: IPJK, IJPK, IJKP
253           INTEGER :: IMJK, IJMK, IJKM
254     ! Diffusion parameter
255           DOUBLE PRECISION :: D_fe, d_fw, d_fn, d_fs, d_ft, d_fb
256     !---------------------------------------------------------------------//
257     
258     
259     !!!$omp      parallel do                                              &
260     !!!$omp&     private(IJK, IPJK, IJPK, IJKM, IMJK, IJMK, IJKM,         &
261     !!!$omp&             d_fe, df_w, df_n, df_s, df_t, df_b)
262           DO IJK = ijkstart3, ijkend3
263     
264              IF (FLUID_AT(IJK)) THEN
265     
266     ! Calculate convection-diffusion fluxes through each of the faces
267                 CALL GET_PHICELL_DIFF_TERMS(dif, d_fe, d_fw, d_fn, d_fs, &
268                    d_ft, d_fb, ijk)
269     
270                 IPJK = IP_OF(IJK)
271                 IJPK = JP_OF(IJK)
272                 IMJK = IM_OF(IJK)
273                 IJMK = JM_OF(IJK)
274     
275     ! East face (i+1/2, j, k)
276                 IF (UF(IJK) >= ZERO) THEN
277                    A_M(IJK,east,M) = D_Fe
278                    A_M(IPJK,west,M) = D_Fe + FLUX_E(IJK)
279                 ELSE
280                    A_M(IJK,east,M) = D_Fe - FLUX_E(IJK)
281                    A_M(IPJK,west,M) = D_Fe
282                 ENDIF
283     ! West face (i-1/2, j, k)
284                 IF (.NOT.FLUID_AT(IMJK)) THEN
285                    IF (UF(IMJK) >= ZERO) THEN
286                       A_M(IJK,west,M) = D_Fw + FLUX_E(IMJK)
287                    ELSE
288                       A_M(IJK,west,M) = D_Fw
289                    ENDIF
290                 ENDIF
291     
292     
293     ! North face (i, j+1/2, k)
294                 IF (VF(IJK) >= ZERO) THEN
295                    A_M(IJK,north,M) = D_Fn
296                    A_M(IJPK,south,M) = D_Fn + FLUX_N(IJK)
297                 ELSE
298                    A_M(IJK,north,M) = D_Fn - FLUX_N(IJK)
299                    A_M(IJPK,south,M) = D_Fn
300                 ENDIF
301     ! South face (i, j-1/2, k)
302                 IF (.NOT.FLUID_AT(IJMK)) THEN
303                    IF (VF(IJMK) >= ZERO) THEN
304                       A_M(IJK,south,M) = D_Fs + FLUX_N(IJMK)
305                    ELSE
306                       A_M(IJK,south,M) = D_Fs
307                    ENDIF
308                 ENDIF
309     
310     
311                 IF (DO_K) THEN
312                    IJKP = KP_OF(IJK)
313                    IJKM = KM_OF(IJK)
314     
315     ! Top face (i, j, k+1/2)
316                    IF (WF(IJK) >= ZERO) THEN
317                       A_M(IJK,top,M) = D_FT
318                       A_M(IJKP,bottom,M) = D_Ft + FLUX_T(IJK)
319                    ELSE
320                       A_M(IJK,top,M) = D_Ft - FLUX_T(IJK)
321                       A_M(IJKP,bottom,M) = D_Ft
322                    ENDIF
323     
324     ! Bottom face (i, j, k-1/2)
325     
326                    IF (.NOT.FLUID_AT(IJKM)) THEN
327                       IF (WF(IJKM) >= ZERO) THEN
328                          A_M(IJK,bottom,M) = D_Fb + FLUX_T(IJKM)
329                       ELSE
330                          A_M(IJK,bottom,M) = D_Fb
331                       ENDIF
332                    ENDIF
333                 ENDIF
334     
335              ENDIF
336           ENDDO
337     
338           RETURN
339           END SUBROUTINE CONV_DIF_PHI0
340     
341     
342     
343     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
344     !                                                                      C
345     !  Subroutine: CONV_DIF_PHI_DC                                         C
346     !  Purpose: Use deferred correction method to solve the scalar         C
347     !  transport equation. This method combines first order upwind and     C
348     !  a user specified higher order method                                C
349     !                                                                      C
350     !  Author: C. GUENTHER                                Date: 1-ARP-99   C
351     !                                                                      C
352     !                                                                      C
353     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
354     
355           SUBROUTINE CONV_DIF_PHI_DC(PHI, DISC, UF, VF, WF, &
356              Flux_E, Flux_N, Flux_T, M, B_M)
357     
358     ! Modules
359     !---------------------------------------------------------------------//
360           USE compar, only: ijkstart3, ijkend3
361     
362           USE discretization, only: fpfoi_of
363     
364           USE function3, only: funijk3
365           USE functions, only: fluid_at
366           USE functions, only: ip_of, jp_of, kp_of
367           USE functions, only: im_of, jm_of, km_of
368     
369           USE geometry, only: do_k
370     
371           USE indices, only: i_of, j_of, k_of
372     
373           USE param, only: dimension_3, dimension_m, dimension_4
374           USE param1, only: zero, one
375     
376           USE run, only: fpfoi
377           USE sendrecv3, only: send_recv3
378     
379           USE xsi, only: calc_xsi
380           IMPLICIT NONE
381     
382     ! Dummy arguments
383     !---------------------------------------------------------------------//
384     ! Scalar
385           DOUBLE PRECISION, INTENT(IN) :: Phi(DIMENSION_3)
386     ! Discretization index
387           INTEGER, INTENT(IN) :: Disc
388     ! Velocity components
389           DOUBLE PRECISION, INTENT(IN) :: Uf(DIMENSION_3)
390           DOUBLE PRECISION, INTENT(IN) :: Vf(DIMENSION_3)
391           DOUBLE PRECISION, INTENT(IN) :: Wf(DIMENSION_3)
392     ! Mass flux components
393           DOUBLE PRECISION, INTENT(IN) :: Flux_E(DIMENSION_3)
394           DOUBLE PRECISION, INTENT(IN) :: Flux_N(DIMENSION_3)
395           DOUBLE PRECISION, INTENT(IN) :: Flux_T(DIMENSION_3)
396     ! Phase index
397           INTEGER, INTENT(IN) :: M
398     ! Vector b_m
399           DOUBLE PRECISION, INTENT(INOUT) :: B_m(DIMENSION_3, 0:DIMENSION_M)
400     
401     ! Dummy arguments
402     !---------------------------------------------------------------------//
403     ! Indices
404           INTEGER :: I, J, K, IJK
405           INTEGER :: IPJK, IJPK, IJKP, IMJK, IJMK, IJKM
406           INTEGER :: IJK4, IPPP, IPPP4, JPPP, JPPP4, KPPP, KPPP4
407           INTEGER :: IMMM, IMMM4, JMMM, JMMM4, KMMM, KMMM4
408     ! indication for shear
409           INTEGER :: incr
410     ! Deferred correction contribution from high order method
411           DOUBLE PRECISION :: PHI_HO
412     ! low order approximation
413           DOUBLE PRECISION :: PHI_LO
414     ! deferred correction contribution from each face
415           DOUBLE PRECISION :: EAST_DC
416           DOUBLE PRECISION :: WEST_DC
417           DOUBLE PRECISION :: NORTH_DC
418           DOUBLE PRECISION :: SOUTH_DC
419           DOUBLE PRECISION :: TOP_DC
420           DOUBLE PRECISION :: BOTTOM_DC
421           DOUBLE PRECISION, DIMENSION(:), ALLOCATABLE :: TMP4
422     
423           DOUBLE PRECISION, DIMENSION(DIMENSION_3) :: XSI_e, XSI_n, XSI_t
424     
425     !---------------------------------------------------------------------//
426     
427           allocate(tmp4(DIMENSION_4))
428     
429     ! Send recv the third ghost layer
430           IF (FPFOI) THEN
431              Do IJK = ijkstart3, ijkend3
432                 I = I_OF(IJK)
433                 J = J_OF(IJK)
434                 K = K_OF(IJK)
435                 IJK4 = funijk3(I,J,K)
436                 TMP4(IJK4) = PHI(IJK)
437              ENDDO
438              CALL send_recv3(tmp4)
439           ENDIF
440     
441     ! shear indicator:
442           incr=0
443           CALL CALC_XSI (DISC, PHI, UF, VF, WF, XSI_E, XSI_N, XSI_T, incr)
444     
445     !!!$omp      parallel do                                             &
446     !!!$omp&     private(I, J, K, IJK,  IPJK, IJPK, IJKP,                &
447     !!!$omp&             IMJK, IJMK, IJKM  V_f,                          &
448     !!!$omp&             PHI_HO, PHI_LO, EAST_DC, WEST_DC, NORTH_DC,     &
449     !!!$omp&             SOUTH_DC, TOP_DC, BOTTOM_DC)
450           DO IJK = ijkstart3, ijkend3
451     
452              IF (FLUID_AT(IJK)) THEN
453     
454                 IPJK = IP_OF(IJK)
455                 IMJK = IM_OF(IJK)
456                 IJPK = JP_OF(IJK)
457                 IJMK = JM_OF(IJK)
458                 IJKP = KP_OF(IJK)
459                 IJKM = KM_OF(IJK)
460     
461     ! Third Ghost layer information
462                 IPPP  = IP_OF(IP_OF(IPJK))
463                 IPPP4 = funijk3(I_OF(IPPP), J_OF(IPPP), K_OF(IPPP))
464                 IMMM  = IM_OF(IM_OF(IMJK))
465                 IMMM4 = funijk3(I_OF(IMMM), J_OF(IMMM), K_OF(IMMM))
466                 JPPP  = JP_OF(JP_OF(IJPK))
467                 JPPP4 = funijk3(I_OF(JPPP), J_OF(JPPP), K_OF(JPPP))
468                 JMMM  = JM_OF(JM_OF(IJMK))
469                 JMMM4 = funijk3(I_OF(JMMM), J_OF(JMMM), K_OF(JMMM))
470                 KPPP  = KP_OF(KP_OF(IJKP))
471                 KPPP4 = funijk3(I_OF(KPPP), J_OF(KPPP), K_OF(KPPP))
472                 KMMM  = KM_OF(KM_OF(IJKM))
473                 KMMM4 = funijk3(I_OF(KMMM), J_OF(KMMM), K_OF(KMMM))
474     
475     
476     ! East face (i+1/2, j, k)
477                 IF(UF(IJK)>= ZERO)THEN
478                    PHI_LO = PHI(IJK)
479                    IF (FPFOI) PHI_HO = FPFOI_OF(PHI(IPJK), PHI(IJK), &
480                                        PHI(IMJK), PHI(IM_OF(IMJK)))
481                 ELSE
482                    PHI_LO = PHI(IPJK)
483                    IF (FPFOI) PHI_HO = FPFOI_OF(PHI(IJK), PHI(IPJK), &
484                                        PHI(IP_OF(IPJK)), TMP4(IPPP4))
485                 ENDIF
486                 IF (.NOT. FPFOI) PHI_HO = XSI_E(IJK)*PHI(IPJK)+&
487                                           (1.0-XSI_E(IJK))*PHI(IJK)
488                 EAST_DC = FLUX_E(IJK)*(PHI_LO - PHI_HO)
489     
490     
491     ! North face (i, j+1/2, k)
492                 IF(VF(IJK) >= ZERO)THEN
493                    PHI_LO = PHI(IJK)
494                    IF (FPFOI) PHI_HO = FPFOI_OF(PHI(IJPK), PHI(IJK), &
495                                        PHI(IJMK), PHI(JM_OF(IJMK)))
496                 ELSE
497                    PHI_LO = PHI(IJPK)
498                    IF (FPFOI) PHI_HO = FPFOI_OF(PHI(IJK), PHI(IJPK), &
499                                        PHI(JP_OF(IJPK)), TMP4(JPPP4))
500                 ENDIF
501                 IF (.NOT. FPFOI) PHI_HO = XSI_N(IJK)*PHI(IJPK)+&
502                                           (1.0-XSI_N(IJK))*PHI(IJK)
503                 NORTH_DC = FLUX_N(IJK)*(PHI_LO - PHI_HO)
504     
505     
506     ! Top face (i, j, k+1/2)
507                 IF (DO_K) THEN
508                    IJKP = KP_OF(IJK)
509                    IF(WF(IJK) >= ZERO)THEN
510                       PHI_LO = PHI(IJK)
511                       IF (FPFOI) PHI_HO = FPFOI_OF(PHI(IJKP), PHI(IJK), &
512                                           PHI(IJKM), PHI(KM_OF(IJKM)))
513                    ELSE
514                       PHI_LO = PHI(IJKP)
515                       IF (FPFOI) PHI_HO = FPFOI_OF(PHI(IJK), PHI(IJKP),  &
516                                           PHI(KP_OF(IJKP)), TMP4(KPPP4))
517                    ENDIF
518                    IF (.NOT. FPFOI) PHI_HO = XSI_T(IJK)*PHI(IJKP)+&
519                                              (1.0-XSI_T(IJK))*PHI(IJK)
520                    TOP_DC = FLUX_T(IJK)*(PHI_LO - PHI_HO)
521                 ELSE
522                    TOP_DC = ZERO
523                 ENDIF
524     
525     
526     ! West face (i-1/2, j, k)
527                 IMJK = IM_OF(IJK)
528                 IF(UF(IMJK) >= ZERO)THEN
529                    PHI_LO = PHI(IMJK)
530                    IF (FPFOI) PHI_HO = FPFOI_OF(PHI(IJK), PHI(IMJK), &
531                                        PHI(IM_OF(IMJK)), TMP4(IMMM4))
532                 ELSE
533                    PHI_LO = PHI(IJK)
534                    IF (FPFOI) PHI_HO = FPFOI_OF(PHI(IMJK), PHI(IJK), &
535                                        PHI(IPJK), PHI(IP_OF(IPJK)))
536                 ENDIF
537                 IF (.NOT. FPFOI) PHI_HO = XSI_E(IMJK)*PHI(IJK)+&
538                                           (ONE-XSI_E(IMJK))*PHI(IMJK)
539                 WEST_DC = FLUX_E(IMJK)*(PHI_LO - PHI_HO)
540     
541     
542     ! South face (i, j-1/2, k)
543                 IJMK = JM_OF(IJK)
544                 IF(VF(IJMK) >= ZERO)THEN
545                    PHI_LO = PHI(IJMK)
546                    IF (FPFOI) PHI_HO = FPFOI_OF(PHI(IJK), PHI(IJMK), &
547                                        PHI(JM_OF(IJMK)), TMP4(JMMM4))
548                 ELSE
549                    PHI_LO = PHI(IJK)
550                    IF (FPFOI) PHI_HO = FPFOI_OF(PHI(IJMK), PHI(IJK), &
551                                        PHI(IJPK), PHI(JP_OF(IJPK)))
552                 ENDIF
553                 IF (.NOT. FPFOI) PHI_HO = XSI_N(IJMK)*PHI(IJK)+&
554                                           (ONE-XSI_N(IJMK))*PHI(IJMK)
555                 SOUTH_DC = FLUX_N(IJMK)*(PHI_LO - PHI_HO)
556     
557     
558     ! Bottom face (i, j, k-1/2)
559                 IF (DO_K) THEN
560                    IJKM = KM_OF(IJK)
561                    IF(WF(IJKM) >= ZERO)THEN
562                       PHI_LO = PHI(IJKM)
563                       IF (FPFOI) PHI_HO = FPFOI_OF(PHI(IJK), PHI(IJKM), &
564                                           PHI(KM_OF(IJKM)), TMP4(KMMM4))
565                    ELSE
566                       PHI_LO = PHI(IJK)
567                       IF (FPFOI) PHI_HO = FPFOI_OF(PHI(IJKM), PHI(IJK), &
568                                           PHI(IJKP), PHI(KP_OF(IJKP)))
569                    ENDIF
570                    IF (.NOT. FPFOI) PHI_HO = XSI_T(IJKM)*PHI(IJK)+&
571                                              (1.0-XSI_T(IJKM))*PHI(IJKM)
572                    BOTTOM_DC = FLUX_T(IJKM)*(PHI_LO - PHI_HO)
573                 ELSE
574                    BOTTOM_DC = ZERO
575                 ENDIF
576     
577     
578     ! CONTRIBUTION DUE TO DEFERRED CORRECTION
579                 B_M(IJK,M) = B_M(IJK,M)+WEST_DC-EAST_DC+SOUTH_DC-&
580                              NORTH_DC+BOTTOM_DC-TOP_DC
581     
582              ENDIF   ! end if fluid_at
583           ENDDO   ! end do ijk
584     
585           DEALLOCATE(tmp4)
586     
587           RETURN
588           END SUBROUTINE CONV_DIF_PHI_DC
589     
590     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
591     !                                                                      C
592     !  Subroutine: CONV_DIF_Phi1                                           C
593     !  Purpose: Determine convection diffusion terms for scalar transport  C
594     !  equations. The off-diagonal coefficients calculated here must be    C
595     !  positive. The center coefficient and the source vector are          C
596     !  negative;                                                           C
597     !  See source_Phi                                                      C
598     !                                                                      C
599     !  Author: M. Syamlal                                 Date: 21-APR-97  C
600     !                                                                      C
601     !                                                                      C
602     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
603           SUBROUTINE CONV_DIF_PHI1(PHI, DIF, DISC, UF, VF, WF, &
604                                    Flux_E, Flux_N, Flux_T, M, A_M)
605     
606     ! Modules
607     !---------------------------------------------------------------------//
608           USE compar, only: ijkstart3, ijkend3
609     
610           USE functions, only: fluid_at
611           USE functions, only: ip_of, jp_of, kp_of
612           USE functions, only: im_of, jm_of, km_of
613           USE geometry, only: do_k
614           USE param
615           USE param1, only: one
616           USE xsi, only: calc_xsi
617           IMPLICIT NONE
618     
619     ! Dummy arguments
620     !---------------------------------------------------------------------//
621     ! Scalar
622           DOUBLE PRECISION, INTENT(IN) :: Phi(DIMENSION_3)
623     ! Gamma -- diffusion coefficient
624           DOUBLE PRECISION, INTENT(IN) :: Dif(DIMENSION_3)
625     ! Discretization index
626           INTEGER, INTENT(IN) :: Disc
627     ! Velocity components
628           DOUBLE PRECISION, INTENT(IN) :: Uf(DIMENSION_3)
629           DOUBLE PRECISION, INTENT(IN) :: Vf(DIMENSION_3)
630           DOUBLE PRECISION, INTENT(IN) :: Wf(DIMENSION_3)
631     ! Mass flux components
632           DOUBLE PRECISION, INTENT(IN) :: Flux_E(DIMENSION_3)
633           DOUBLE PRECISION, INTENT(IN) :: Flux_N(DIMENSION_3)
634           DOUBLE PRECISION, INTENT(IN) :: Flux_T(DIMENSION_3)
635     ! Phase index
636           INTEGER, INTENT(IN) :: M
637     ! Septadiagonal matrix A_m
638           DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
639     
640     ! Local variables
641     !---------------------------------------------------------------------//
642     ! Indices
643           INTEGER :: IJK
644           INTEGER :: IPJK, IJPK, IJKP
645           INTEGER :: IMJK, IJMK, IJKM
646     ! indicator for shear
647           INTEGER :: incr
648     ! Diffusion parameter
649           DOUBLE PRECISION :: D_fe, d_fw, d_fn, d_fs, d_ft, d_fb
650     
651           DOUBLE PRECISION, DIMENSION(DIMENSION_3) :: XSI_e, XSI_n, XSI_t
652     
653     !---------------------------------------------------------------------//
654     
655     ! shear indicator
656           incr=0
657           CALL CALC_XSI (DISC, PHI, UF, VF, WF, XSI_E, XSI_N, XSI_T, incr)
658     
659     !!!$omp      parallel do                                                 &
660     !!!$omp&     private(IJK, IPJK, IJPK, IJKP, IMJK, IJMK, IJKM,            &
661     !!!$omp&             d_fe, d_fw, d_fn, d_fs, d_ft, d_fb)
662           DO IJK = ijkstart3, ijkend3
663     
664              IF (FLUID_AT(IJK)) THEN
665     ! Calculate convection-diffusion fluxes through each of the faces
666                 CALL GET_PHICELL_DIFF_TERMS(dif, d_fe, d_fw, d_fn, d_fs, &
667                    d_ft, d_fb, ijk)
668     
669                 IPJK = IP_OF(IJK)
670                 IJPK = JP_OF(IJK)
671                 IMJK = IM_OF(IJK)
672                 IJMK = JM_OF(IJK)
673     
674     ! East face (i+1/2, j, k)
675                 A_M(IJK,east,M) = D_Fe - XSI_E(IJK)*FLUX_E(IJK)
676                 A_M(IPJK,west,M) = D_Fe + (ONE - XSI_E(IJK))*FLUX_E(IJK)
677     !  West face (i-1/2, j, k)
678                 IF (.NOT.FLUID_AT(IMJK)) THEN
679                   A_M(IJK,west,M) = D_Fw + (ONE - XSI_E(IMJK))*FLUX_E(IMJK)
680                 ENDIF
681     
682     
683     ! North face (i, j+1/2, k)
684                 A_M(IJK,north,M) = D_Fn - XSI_N(IJK)*FLUX_N(IJK)
685                 A_M(IJPK,south,M) = D_Fn + (ONE - XSI_N(IJK))*FLUX_N(IJK)
686     ! South face (i, j-1/2, k)
687                 IF (.NOT.FLUID_AT(IJMK)) THEN
688                    A_M(IJK,south,M) = D_Fs + (ONE - XSI_N(IJMK))*FLUX_N(IJMK)
689                 ENDIF
690     
691     
692                 IF (DO_K) THEN
693                    IJKP = KP_OF(IJK)
694                    IJKM = KM_OF(IJK)
695     ! Top face (i, j, k+1/2)
696                    A_M(IJK,top,M) = D_Ft - XSI_T(IJK)*FLUX_T(IJK)
697                    A_M(IJKP,bottom,M)=D_Ft + (ONE-XSI_T(IJK))*FLUX_T(IJK)
698     ! Bottom face (i, j, k-1/2)
699                    IF (.NOT.FLUID_AT(IJKM)) THEN
700                       A_M(IJK,bottom,M) = D_Fb + (ONE - XSI_T(IJKM))*FLUX_T(IJKM)
701                    ENDIF
702                 ENDIF
703     
704              ENDIF
705           ENDDO
706     
707           RETURN
708           END SUBROUTINE CONV_DIF_PHI1
709     
710     
711     
712     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
713     !                                                                      C
714     !  Subroutine: DIF_phi_IS                                              C
715     !  Purpose: Remove diffusive fluxes across internal surfaces.          C
716     !  (Make user defined internal surfaces non-conducting)                C
717     !                                                                      C
718     !  Author: M. Syamlal                                 Date: 30-APR-97  C
719     !                                                                      C
720     !                                                                      C
721     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
722           SUBROUTINE DIF_PHI_IS(DIF, A_M, M)
723     
724     ! Modules
725     !---------------------------------------------------------------------//
726           USE param
727           USE geometry, only: do_k
728           USE geometry, only: ody_n, odx_e, odz_t, ox
729           USE geometry, only: axz, axy, ayz
730     
731           USE is, only: is_defined, is_plane
732           USE is, only: is_i_w, is_i_e, is_j_s, is_j_n, is_k_t, is_k_b
733     
734           USE fun_avg, only: avg_x_h, avg_y_h, avg_z_h
735     
736           USE functions, only: funijk, ip_of, jp_of, kp_of
737           USE functions, only: east_of, north_of, top_of
738     
739           USE compar, only: dead_cell_at
740           USE compar, only: istart2, jstart2, kstart2
741           USE compar, only: iend2, jend2, kend2
742           IMPLICIT NONE
743     
744     ! Dummy arguments
745     !---------------------------------------------------------------------//
746     ! Gamma -- diffusion coefficient
747           DOUBLE PRECISION, INTENT(IN) :: Dif(DIMENSION_3)
748     ! Septadiagonal matrix A_m
749           DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
750     ! Solids phase
751           INTEGER, INTENT(IN) :: M
752     
753     ! Local variables
754     !---------------------------------------------------------------------//
755     ! Internal surface
756           INTEGER :: L
757     ! Indices
758           INTEGER :: I1, I2, J1, J2, K1, K2
759           INTEGER :: I, J, K, IJK
760           INTEGER :: IJKE, IJKN, IJKT, IPJK, IJPK, IJKP
761     ! Diffusion parameter
762           DOUBLE PRECISION :: D_f
763     !---------------------------------------------------------------------//
764     
765     
766           DO L = 1, DIMENSION_IS
767              IF (IS_DEFINED(L)) THEN
768                 I1 = IS_I_W(L)
769                 I2 = IS_I_E(L)
770                 J1 = IS_J_S(L)
771                 J2 = IS_J_N(L)
772                 K1 = IS_K_B(L)
773                 K2 = IS_K_T(L)
774     
775     ! Limit I1, I2 and all to local processor first ghost layer
776                 IF(I1.LE.IEND2)   I1 = MAX(I1, ISTART2)
777                 IF(J1.LE.JEND2)   J1 = MAX(J1, JSTART2)
778                 IF(K1.LE.KEND2)   K1 = MAX(K1, KSTART2)
779                 IF(I2.GE.ISTART2) I2 = MIN(I2, IEND2)
780                 IF(J2.GE.JSTART2) J2 = MIN(J2, JEND2)
781                 IF(K2.GE.KSTART2) K2 = MIN(K2, KEND2)
782     
783                 IF (IS_PLANE(L) == 'E') THEN
784                    DO K = K1, K2
785                    DO J = J1, J2
786                    DO I = I1, I2
787                       IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
788                       IJK = FUNIJK(I,J,K)
789                       IJKE = EAST_OF(IJK)
790                       IPJK = IP_OF(IJK)
791     
792                       D_F = AVG_X_H(DIF(IJK),DIF(IJKE),I)*ODX_E(I)*AYZ(IJK)
793                       A_M(IJK,east,M) = A_M(IJK,east,M) - D_F
794                       A_M(IPJK,west,M) = A_M(IPJK,west,M) - D_F
795                    ENDDO
796                    ENDDO
797                    ENDDO
798     
799                 ELSEIF(IS_PLANE(L) == 'N') THEN
800                    DO K = K1, K2
801                    DO J = J1, J2
802                    DO I = I1, I2
803                       IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
804                       IJK = FUNIJK(I,J,K)
805                       IJKN = NORTH_OF(IJK)
806                       IJPK = JP_OF(IJK)
807     
808                       D_F = AVG_Y_H(DIF(IJK),DIF(IJKN),J)*ODY_N(J)*AXZ(IJK)
809                       A_M(IJK,north,M) = A_M(IJK,north,M) - D_F
810                       A_M(IJPK,south,M) = A_M(IJPK,south,M) - D_F
811                    ENDDO
812                    ENDDO
813                    ENDDO
814     
815                 ELSEIF(IS_PLANE(L) == 'T') THEN
816                    IF (DO_K) THEN
817                       DO K = K1, K2
818                       DO J = J1, J2
819                       DO I = I1, I2
820                          IJKT = TOP_OF(IJK)
821                          IJKP = KP_OF(IJK)
822     
823                          D_F = AVG_Z_H(DIF(IJK),DIF(IJKT),K)*&
824                             OX(I)*ODZ_T(K)*AXY(IJK)
825                          A_M(IJK,top,M) = A_M(IJK,top,M) - D_F
826                          A_M(IJKP,bottom,M) = A_M(IJKP,bottom,M) - D_F
827                       ENDDO
828                       ENDDO
829                       ENDDO
830                    ENDIF   ! end if do_k
831                 ENDIF   ! endif/else is_plane
832              ENDIF   ! endif is_defined
833           ENDDO   ! end do dimension_is
834     
835           RETURN
836           END SUBROUTINE DIF_PHI_IS
837