File: RELATIVE:/../../../mfix.git/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     
225           USE geometry, only: do_k
226     
227           USE param, only: dimension_3, dimension_m
228           USE param1, only: zero
229           USE matrix, only: e, w, n, s, t, b
230           IMPLICIT NONE
231     
232     ! Dummy arguments
233     !---------------------------------------------------------------------//
234     ! Scalar
235           DOUBLE PRECISION, INTENT(IN) :: Phi(DIMENSION_3)
236     ! Gamma -- diffusion coefficient
237           DOUBLE PRECISION, INTENT(IN) :: Dif(DIMENSION_3)
238     ! Velocity components
239           DOUBLE PRECISION, INTENT(IN) :: Uf(DIMENSION_3)
240           DOUBLE PRECISION, INTENT(IN) :: Vf(DIMENSION_3)
241           DOUBLE PRECISION, INTENT(IN) :: Wf(DIMENSION_3)
242     ! Mass flux components
243           DOUBLE PRECISION, INTENT(IN) :: Flux_E(DIMENSION_3)
244           DOUBLE PRECISION, INTENT(IN) :: Flux_N(DIMENSION_3)
245           DOUBLE PRECISION, INTENT(IN) :: Flux_T(DIMENSION_3)
246     ! Phase index
247           INTEGER, INTENT(IN) :: M
248     ! Septadiagonal matrix A_m
249           DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
250     
251     ! Local variables
252     !---------------------------------------------------------------------//
253     ! Indices
254           INTEGER :: IJK
255           INTEGER :: IPJK, IJPK, IJKP
256           INTEGER :: IMJK, IJMK, IJKM
257     ! Diffusion parameter
258           DOUBLE PRECISION :: D_fe, d_fw, d_fn, d_fs, d_ft, d_fb
259     !---------------------------------------------------------------------//
260     
261     
262     !!!$omp      parallel do                                              &
263     !!!$omp&     private(IJK, IPJK, IJPK, IJKM, IMJK, IJMK, IJKM,         &
264     !!!$omp&             d_fe, df_w, df_n, df_s, df_t, df_b)
265           DO IJK = ijkstart3, ijkend3
266     
267              IF (FLUID_AT(IJK)) THEN
268     
269     ! Calculate convection-diffusion fluxes through each of the faces
270                 CALL GET_PHICELL_DIFF_TERMS(dif, d_fe, d_fw, d_fn, d_fs, &
271                    d_ft, d_fb, ijk)
272     
273                 IPJK = IP_OF(IJK)
274                 IJPK = JP_OF(IJK)
275                 IMJK = IM_OF(IJK)
276                 IJMK = JM_OF(IJK)
277     
278     ! East face (i+1/2, j, k)
279                 IF (UF(IJK) >= ZERO) THEN
280                    A_M(IJK,E,M) = D_Fe
281                    A_M(IPJK,W,M) = D_Fe + FLUX_E(IJK)
282                 ELSE
283                    A_M(IJK,E,M) = D_Fe - FLUX_E(IJK)
284                    A_M(IPJK,W,M) = D_Fe
285                 ENDIF
286     ! West face (i-1/2, j, k)
287                 IF (.NOT.FLUID_AT(IMJK)) THEN
288                    IF (UF(IMJK) >= ZERO) THEN
289                       A_M(IJK,W,M) = D_Fw + FLUX_E(IMJK)
290                    ELSE
291                       A_M(IJK,W,M) = D_Fw
292                    ENDIF
293                 ENDIF
294     
295     
296     ! North face (i, j+1/2, k)
297                 IF (VF(IJK) >= ZERO) THEN
298                    A_M(IJK,N,M) = D_Fn
299                    A_M(IJPK,S,M) = D_Fn + FLUX_N(IJK)
300                 ELSE
301                    A_M(IJK,N,M) = D_Fn - FLUX_N(IJK)
302                    A_M(IJPK,S,M) = D_Fn
303                 ENDIF
304     ! South face (i, j-1/2, k)
305                 IF (.NOT.FLUID_AT(IJMK)) THEN
306                    IF (VF(IJMK) >= ZERO) THEN
307                       A_M(IJK,S,M) = D_Fs + FLUX_N(IJMK)
308                    ELSE
309                       A_M(IJK,S,M) = D_Fs
310                    ENDIF
311                 ENDIF
312     
313     
314                 IF (DO_K) THEN
315                    IJKP = KP_OF(IJK)
316                    IJKM = KM_OF(IJK)
317     
318     ! Top face (i, j, k+1/2)
319                    IF (WF(IJK) >= ZERO) THEN
320                       A_M(IJK,T,M) = D_FT
321                       A_M(IJKP,B,M) = D_Ft + FLUX_T(IJK)
322                    ELSE
323                       A_M(IJK,T,M) = D_Ft - FLUX_T(IJK)
324                       A_M(IJKP,B,M) = D_Ft
325                    ENDIF
326     
327     ! Bottom face (i, j, k-1/2)
328     
329                    IF (.NOT.FLUID_AT(IJKM)) THEN
330                       IF (WF(IJKM) >= ZERO) THEN
331                          A_M(IJK,B,M) = D_Fb + FLUX_T(IJKM)
332                       ELSE
333                          A_M(IJK,B,M) = D_Fb
334                       ENDIF
335                    ENDIF
336                 ENDIF
337     
338              ENDIF
339           ENDDO
340     
341           RETURN
342           END SUBROUTINE CONV_DIF_PHI0
343     
344     
345     
346     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
347     !                                                                      C
348     !  Subroutine: CONV_DIF_PHI_DC                                         C
349     !  Purpose: Use deferred correction method to solve the scalar         C
350     !  transport equation. This method combines first order upwind and     C
351     !  a user specified higher order method                                C
352     !                                                                      C
353     !  Author: C. GUENTHER                                Date: 1-ARP-99   C
354     !                                                                      C
355     !                                                                      C
356     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
357     
358           SUBROUTINE CONV_DIF_PHI_DC(PHI, DISC, UF, VF, WF, &
359              Flux_E, Flux_N, Flux_T, M, B_M)
360     
361     ! Modules
362     !---------------------------------------------------------------------//
363           USE compar, only: ijkstart3, ijkend3
364     
365           USE discretization, only: fpfoi_of
366     
367           USE function3, only: funijk3
368           USE functions, only: fluid_at
369           USE functions, only: ip_of, jp_of, kp_of
370           USE functions, only: im_of, jm_of, km_of
371     
372           USE geometry, only: do_k
373     
374           USE indices, only: i_of, j_of, k_of
375     
376           USE param, only: dimension_3, dimension_m
377           USE param1, only: zero, one
378     
379           USE run, only: fpfoi
380           USE sendrecv3, only: send_recv3
381     
382           USE tmp_array, only: tmp4
383           USE tmp_array, only: lock_tmp4_array, unlock_tmp4_array
384     
385           USE xsi, only: calc_xsi
386           USE xsi_array, only: xsi_e, xsi_n, xsi_t
387           USE xsi_array, only: lock_xsi_array, unlock_xsi_array
388           IMPLICIT NONE
389     
390     ! Dummy arguments
391     !---------------------------------------------------------------------//
392     ! Scalar
393           DOUBLE PRECISION, INTENT(IN) :: Phi(DIMENSION_3)
394     ! Discretization index
395           INTEGER, INTENT(IN) :: Disc
396     ! Velocity components
397           DOUBLE PRECISION, INTENT(IN) :: Uf(DIMENSION_3)
398           DOUBLE PRECISION, INTENT(IN) :: Vf(DIMENSION_3)
399           DOUBLE PRECISION, INTENT(IN) :: Wf(DIMENSION_3)
400     ! Mass flux components
401           DOUBLE PRECISION, INTENT(IN) :: Flux_E(DIMENSION_3)
402           DOUBLE PRECISION, INTENT(IN) :: Flux_N(DIMENSION_3)
403           DOUBLE PRECISION, INTENT(IN) :: Flux_T(DIMENSION_3)
404     ! Phase index
405           INTEGER, INTENT(IN) :: M
406     ! Vector b_m
407           DOUBLE PRECISION, INTENT(INOUT) :: B_m(DIMENSION_3, 0:DIMENSION_M)
408     
409     ! Dummy arguments
410     !---------------------------------------------------------------------//
411     ! Indices
412           INTEGER :: I, J, K, IJK
413           INTEGER :: IPJK, IJPK, IJKP, IMJK, IJMK, IJKM
414           INTEGER :: IJK4, IPPP, IPPP4, JPPP, JPPP4, KPPP, KPPP4
415           INTEGER :: IMMM, IMMM4, JMMM, JMMM4, KMMM, KMMM4
416     ! indication for shear
417           INTEGER :: incr
418     ! Deferred correction contribution from high order method
419           DOUBLE PRECISION :: PHI_HO
420     ! low order approximation
421           DOUBLE PRECISION :: PHI_LO
422     ! deferred correction contribution from each face
423           DOUBLE PRECISION :: EAST_DC
424           DOUBLE PRECISION :: WEST_DC
425           DOUBLE PRECISION :: NORTH_DC
426           DOUBLE PRECISION :: SOUTH_DC
427           DOUBLE PRECISION :: TOP_DC
428           DOUBLE PRECISION :: BOTTOM_DC
429     !---------------------------------------------------------------------//
430     
431           call lock_xsi_array
432           call lock_tmp4_array
433     
434     ! Send recv the third ghost layer
435           IF (FPFOI) THEN
436              Do IJK = ijkstart3, ijkend3
437                 I = I_OF(IJK)
438                 J = J_OF(IJK)
439                 K = K_OF(IJK)
440                 IJK4 = funijk3(I,J,K)
441                 TMP4(IJK4) = PHI(IJK)
442              ENDDO
443              CALL send_recv3(tmp4)
444           ENDIF
445     
446     ! shear indicator:
447           incr=0
448           CALL CALC_XSI (DISC, PHI, UF, VF, WF, XSI_E, XSI_N, XSI_T, incr)
449     
450     !!!$omp      parallel do                                             &
451     !!!$omp&     private(I, J, K, IJK,  IPJK, IJPK, IJKP,                &
452     !!!$omp&             IMJK, IJMK, IJKM  V_f,                          &
453     !!!$omp&             PHI_HO, PHI_LO, EAST_DC, WEST_DC, NORTH_DC,     &
454     !!!$omp&             SOUTH_DC, TOP_DC, BOTTOM_DC)
455           DO IJK = ijkstart3, ijkend3
456     
457              IF (FLUID_AT(IJK)) THEN
458     
459                 IPJK = IP_OF(IJK)
460                 IMJK = IM_OF(IJK)
461                 IJPK = JP_OF(IJK)
462                 IJMK = JM_OF(IJK)
463                 IJKP = KP_OF(IJK)
464                 IJKM = KM_OF(IJK)
465     
466     ! Third Ghost layer information
467                 IPPP  = IP_OF(IP_OF(IPJK))
468                 IPPP4 = funijk3(I_OF(IPPP), J_OF(IPPP), K_OF(IPPP))
469                 IMMM  = IM_OF(IM_OF(IMJK))
470                 IMMM4 = funijk3(I_OF(IMMM), J_OF(IMMM), K_OF(IMMM))
471                 JPPP  = JP_OF(JP_OF(IJPK))
472                 JPPP4 = funijk3(I_OF(JPPP), J_OF(JPPP), K_OF(JPPP))
473                 JMMM  = JM_OF(JM_OF(IJMK))
474                 JMMM4 = funijk3(I_OF(JMMM), J_OF(JMMM), K_OF(JMMM))
475                 KPPP  = KP_OF(KP_OF(IJKP))
476                 KPPP4 = funijk3(I_OF(KPPP), J_OF(KPPP), K_OF(KPPP))
477                 KMMM  = KM_OF(KM_OF(IJKM))
478                 KMMM4 = funijk3(I_OF(KMMM), J_OF(KMMM), K_OF(KMMM))
479     
480     
481     ! East face (i+1/2, j, k)
482                 IF(UF(IJK)>= ZERO)THEN
483                    PHI_LO = PHI(IJK)
484                    IF (FPFOI) PHI_HO = FPFOI_OF(PHI(IPJK), PHI(IJK), &
485                                        PHI(IMJK), PHI(IM_OF(IMJK)))
486                 ELSE
487                    PHI_LO = PHI(IPJK)
488                    IF (FPFOI) PHI_HO = FPFOI_OF(PHI(IJK), PHI(IPJK), &
489                                        PHI(IP_OF(IPJK)), TMP4(IPPP4))
490                 ENDIF
491                 IF (.NOT. FPFOI) PHI_HO = XSI_E(IJK)*PHI(IPJK)+&
492                                           (1.0-XSI_E(IJK))*PHI(IJK)
493                 EAST_DC = FLUX_E(IJK)*(PHI_LO - PHI_HO)
494     
495     
496     ! North face (i, j+1/2, k)
497                 IF(VF(IJK) >= ZERO)THEN
498                    PHI_LO = PHI(IJK)
499                    IF (FPFOI) PHI_HO = FPFOI_OF(PHI(IJPK), PHI(IJK), &
500                                        PHI(IJMK), PHI(JM_OF(IJMK)))
501                 ELSE
502                    PHI_LO = PHI(IJPK)
503                    IF (FPFOI) PHI_HO = FPFOI_OF(PHI(IJK), PHI(IJPK), &
504                                        PHI(JP_OF(IJPK)), TMP4(JPPP4))
505                 ENDIF
506                 IF (.NOT. FPFOI) PHI_HO = XSI_N(IJK)*PHI(IJPK)+&
507                                           (1.0-XSI_N(IJK))*PHI(IJK)
508                 NORTH_DC = FLUX_N(IJK)*(PHI_LO - PHI_HO)
509     
510     
511     ! Top face (i, j, k+1/2)
512                 IF (DO_K) THEN
513                    IJKP = KP_OF(IJK)
514                    IF(WF(IJK) >= ZERO)THEN
515                       PHI_LO = PHI(IJK)
516                       IF (FPFOI) PHI_HO = FPFOI_OF(PHI(IJKP), PHI(IJK), &
517                                           PHI(IJKM), PHI(KM_OF(IJKM)))
518                    ELSE
519                       PHI_LO = PHI(IJKP)
520                       IF (FPFOI) PHI_HO = FPFOI_OF(PHI(IJK), PHI(IJKP),  &
521                                           PHI(KP_OF(IJKP)), TMP4(KPPP4))
522                    ENDIF
523                    IF (.NOT. FPFOI) PHI_HO = XSI_T(IJK)*PHI(IJKP)+&
524                                              (1.0-XSI_T(IJK))*PHI(IJK)
525                    TOP_DC = FLUX_T(IJK)*(PHI_LO - PHI_HO)
526                 ELSE
527                    TOP_DC = ZERO
528                 ENDIF
529     
530     
531     ! West face (i-1/2, j, k)
532                 IMJK = IM_OF(IJK)
533                 IF(UF(IMJK) >= ZERO)THEN
534                    PHI_LO = PHI(IMJK)
535                    IF (FPFOI) PHI_HO = FPFOI_OF(PHI(IJK), PHI(IMJK), &
536                                        PHI(IM_OF(IMJK)), TMP4(IMMM4))
537                 ELSE
538                    PHI_LO = PHI(IJK)
539                    IF (FPFOI) PHI_HO = FPFOI_OF(PHI(IMJK), PHI(IJK), &
540                                        PHI(IPJK), PHI(IP_OF(IPJK)))
541                 ENDIF
542                 IF (.NOT. FPFOI) PHI_HO = XSI_E(IMJK)*PHI(IJK)+&
543                                           (ONE-XSI_E(IMJK))*PHI(IMJK)
544                 WEST_DC = FLUX_E(IMJK)*(PHI_LO - PHI_HO)
545     
546     
547     ! South face (i, j-1/2, k)
548                 IJMK = JM_OF(IJK)
549                 IF(VF(IJMK) >= ZERO)THEN
550                    PHI_LO = PHI(IJMK)
551                    IF (FPFOI) PHI_HO = FPFOI_OF(PHI(IJK), PHI(IJMK), &
552                                        PHI(JM_OF(IJMK)), TMP4(JMMM4))
553                 ELSE
554                    PHI_LO = PHI(IJK)
555                    IF (FPFOI) PHI_HO = FPFOI_OF(PHI(IJMK), PHI(IJK), &
556                                        PHI(IJPK), PHI(JP_OF(IJPK)))
557                 ENDIF
558                 IF (.NOT. FPFOI) PHI_HO = XSI_N(IJMK)*PHI(IJK)+&
559                                           (ONE-XSI_N(IJMK))*PHI(IJMK)
560                 SOUTH_DC = FLUX_N(IJMK)*(PHI_LO - PHI_HO)
561     
562     
563     ! Bottom face (i, j, k-1/2)
564                 IF (DO_K) THEN
565                    IJKM = KM_OF(IJK)
566                    IF(WF(IJKM) >= ZERO)THEN
567                       PHI_LO = PHI(IJKM)
568                       IF (FPFOI) PHI_HO = FPFOI_OF(PHI(IJK), PHI(IJKM), &
569                                           PHI(KM_OF(IJKM)), TMP4(KMMM4))
570                    ELSE
571                       PHI_LO = PHI(IJK)
572                       IF (FPFOI) PHI_HO = FPFOI_OF(PHI(IJKM), PHI(IJK), &
573                                           PHI(IJKP), PHI(KP_OF(IJKP)))
574                    ENDIF
575                    IF (.NOT. FPFOI) PHI_HO = XSI_T(IJKM)*PHI(IJK)+&
576                                              (1.0-XSI_T(IJKM))*PHI(IJKM)
577                    BOTTOM_DC = FLUX_T(IJKM)*(PHI_LO - PHI_HO)
578                 ELSE
579                    BOTTOM_DC = ZERO
580                 ENDIF
581     
582     
583     ! CONTRIBUTION DUE TO DEFERRED CORRECTION
584                 B_M(IJK,M) = B_M(IJK,M)+WEST_DC-EAST_DC+SOUTH_DC-&
585                              NORTH_DC+BOTTOM_DC-TOP_DC
586     
587              ENDIF   ! end if fluid_at
588           ENDDO   ! end do ijk
589     
590           call unlock_tmp4_array
591           call unlock_xsi_array
592     
593           RETURN
594           END SUBROUTINE CONV_DIF_PHI_DC
595     
596     
597     
598     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
599     !                                                                      C
600     !  Subroutine: CONV_DIF_Phi1                                           C
601     !  Purpose: Determine convection diffusion terms for scalar transport  C
602     !  equations. The off-diagonal coefficients calculated here must be    C
603     !  positive. The center coefficient and the source vector are          C
604     !  negative;                                                           C
605     !  See source_Phi                                                      C
606     !                                                                      C
607     !  Author: M. Syamlal                                 Date: 21-APR-97  C
608     !                                                                      C
609     !                                                                      C
610     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
611           SUBROUTINE CONV_DIF_PHI1(PHI, DIF, DISC, UF, VF, WF, &
612                                    Flux_E, Flux_N, Flux_T, M, A_M)
613     
614     ! Modules
615     !---------------------------------------------------------------------//
616           USE compar, only: ijkstart3, ijkend3
617     
618           USE functions, only: fluid_at
619           USE functions, only: ip_of, jp_of, kp_of
620           USE functions, only: im_of, jm_of, km_of
621     
622           USE geometry, only: do_k
623     
624           USE param, only: dimension_3, dimension_m
625           USE param1, only: one
626     
627           USE matrix, only: e, w, n, s, t, b
628     
629           USE xsi, only: calc_xsi
630           USE xsi_array, only: xsi_e, xsi_n, xsi_t
631           USE xsi_array, only: lock_xsi_array, unlock_xsi_array
632           IMPLICIT NONE
633     
634     ! Dummy arguments
635     !---------------------------------------------------------------------//
636     ! Scalar
637           DOUBLE PRECISION, INTENT(IN) :: Phi(DIMENSION_3)
638     ! Gamma -- diffusion coefficient
639           DOUBLE PRECISION, INTENT(IN) :: Dif(DIMENSION_3)
640     ! Discretization index
641           INTEGER, INTENT(IN) :: Disc
642     ! Velocity components
643           DOUBLE PRECISION, INTENT(IN) :: Uf(DIMENSION_3)
644           DOUBLE PRECISION, INTENT(IN) :: Vf(DIMENSION_3)
645           DOUBLE PRECISION, INTENT(IN) :: Wf(DIMENSION_3)
646     ! Mass flux components
647           DOUBLE PRECISION, INTENT(IN) :: Flux_E(DIMENSION_3)
648           DOUBLE PRECISION, INTENT(IN) :: Flux_N(DIMENSION_3)
649           DOUBLE PRECISION, INTENT(IN) :: Flux_T(DIMENSION_3)
650     ! Phase index
651           INTEGER, INTENT(IN) :: M
652     ! Septadiagonal matrix A_m
653           DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
654     
655     ! Local variables
656     !---------------------------------------------------------------------//
657     ! Indices
658           INTEGER :: IJK
659           INTEGER :: IPJK, IJPK, IJKP
660           INTEGER :: IMJK, IJMK, IJKM
661     ! indicator for shear
662           INTEGER :: incr
663     ! Diffusion parameter
664           DOUBLE PRECISION :: D_fe, d_fw, d_fn, d_fs, d_ft, d_fb
665     !---------------------------------------------------------------------//
666     
667           call lock_xsi_array
668     
669     ! shear indicator
670           incr=0
671           CALL CALC_XSI (DISC, PHI, UF, VF, WF, XSI_E, XSI_N, XSI_T, incr)
672     
673     !!!$omp      parallel do                                                 &
674     !!!$omp&     private(IJK, IPJK, IJPK, IJKP, IMJK, IJMK, IJKM,            &
675     !!!$omp&             d_fe, d_fw, d_fn, d_fs, d_ft, d_fb)
676           DO IJK = ijkstart3, ijkend3
677     
678              IF (FLUID_AT(IJK)) THEN
679     ! Calculate convection-diffusion fluxes through each of the faces
680                 CALL GET_PHICELL_DIFF_TERMS(dif, d_fe, d_fw, d_fn, d_fs, &
681                    d_ft, d_fb, ijk)
682     
683                 IPJK = IP_OF(IJK)
684                 IJPK = JP_OF(IJK)
685                 IMJK = IM_OF(IJK)
686                 IJMK = JM_OF(IJK)
687     
688     ! East face (i+1/2, j, k)
689                 A_M(IJK,E,M) = D_Fe - XSI_E(IJK)*FLUX_E(IJK)
690                 A_M(IPJK,W,M) = D_Fe + (ONE - XSI_E(IJK))*FLUX_E(IJK)
691     !  West face (i-1/2, j, k)
692                 IF (.NOT.FLUID_AT(IMJK)) THEN
693                   A_M(IJK,W,M) = D_Fw + (ONE - XSI_E(IMJK))*FLUX_E(IMJK)
694                 ENDIF
695     
696     
697     ! North face (i, j+1/2, k)
698                 A_M(IJK,N,M) = D_Fn - XSI_N(IJK)*FLUX_N(IJK)
699                 A_M(IJPK,S,M) = D_Fn + (ONE - XSI_N(IJK))*FLUX_N(IJK)
700     ! South face (i, j-1/2, k)
701                 IF (.NOT.FLUID_AT(IJMK)) THEN
702                    A_M(IJK,S,M) = D_Fs + (ONE - XSI_N(IJMK))*FLUX_N(IJMK)
703                 ENDIF
704     
705     
706                 IF (DO_K) THEN
707                    IJKP = KP_OF(IJK)
708                    IJKM = KM_OF(IJK)
709     ! Top face (i, j, k+1/2)
710                    A_M(IJK,T,M) = D_Ft - XSI_T(IJK)*FLUX_T(IJK)
711                    A_M(IJKP,B,M)=D_Ft + (ONE-XSI_T(IJK))*FLUX_T(IJK)
712     ! Bottom face (i, j, k-1/2)
713                    IF (.NOT.FLUID_AT(IJKM)) THEN
714                       A_M(IJK,B,M) = D_Fb + (ONE - XSI_T(IJKM))*FLUX_T(IJKM)
715                    ENDIF
716                 ENDIF
717     
718              ENDIF
719           ENDDO
720     
721           call unlock_xsi_array
722     
723           RETURN
724           END SUBROUTINE CONV_DIF_PHI1
725     
726     
727     
728     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
729     !                                                                      C
730     !  Subroutine: DIF_phi_IS                                              C
731     !  Purpose: Remove diffusive fluxes across internal surfaces.          C
732     !  (Make user defined internal surfaces non-conducting)                C
733     !                                                                      C
734     !  Author: M. Syamlal                                 Date: 30-APR-97  C
735     !                                                                      C
736     !                                                                      C
737     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
738           SUBROUTINE DIF_PHI_IS(DIF, A_M, M)
739     
740     ! Modules
741     !---------------------------------------------------------------------//
742           USE param, only: dimension_is, dimension_3, dimension_m
743           USE matrix, only: e, w, n, s, t, b
744           USE geometry, only: do_k
745           USE geometry, only: ody_n, odx_e, odz_t, ox
746           USE geometry, only: axz, axy, ayz
747     
748           USE is, only: is_defined, is_plane
749           USE is, only: is_i_w, is_i_e, is_j_s, is_j_n, is_k_t, is_k_b
750     
751           USE fun_avg, only: avg_x_h, avg_y_h, avg_z_h
752     
753           USE functions, only: funijk, ip_of, jp_of, kp_of
754           USE functions, only: east_of, north_of, top_of
755     
756           USE compar, only: dead_cell_at
757           USE compar, only: istart2, jstart2, kstart2
758           USE compar, only: iend2, jend2, kend2
759           IMPLICIT NONE
760     
761     ! Dummy arguments
762     !---------------------------------------------------------------------//
763     ! Gamma -- diffusion coefficient
764           DOUBLE PRECISION, INTENT(IN) :: Dif(DIMENSION_3)
765     ! Septadiagonal matrix A_m
766           DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
767     ! Solids phase
768           INTEGER, INTENT(IN) :: M
769     
770     ! Local variables
771     !---------------------------------------------------------------------//
772     ! Internal surface
773           INTEGER :: L
774     ! Indices
775           INTEGER :: I1, I2, J1, J2, K1, K2
776           INTEGER :: I, J, K, IJK
777           INTEGER :: IJKE, IJKN, IJKT, IPJK, IJPK, IJKP
778     ! Diffusion parameter
779           DOUBLE PRECISION :: D_f
780     !---------------------------------------------------------------------//
781     
782     
783           DO L = 1, DIMENSION_IS
784              IF (IS_DEFINED(L)) THEN
785                 I1 = IS_I_W(L)
786                 I2 = IS_I_E(L)
787                 J1 = IS_J_S(L)
788                 J2 = IS_J_N(L)
789                 K1 = IS_K_B(L)
790                 K2 = IS_K_T(L)
791     
792     ! Limit I1, I2 and all to local processor first ghost layer
793                 IF(I1.LE.IEND2)   I1 = MAX(I1, ISTART2)
794                 IF(J1.LE.JEND2)   J1 = MAX(J1, JSTART2)
795                 IF(K1.LE.KEND2)   K1 = MAX(K1, KSTART2)
796                 IF(I2.GE.ISTART2) I2 = MIN(I2, IEND2)
797                 IF(J2.GE.JSTART2) J2 = MIN(J2, JEND2)
798                 IF(K2.GE.KSTART2) K2 = MIN(K2, KEND2)
799     
800                 IF (IS_PLANE(L) == 'E') THEN
801                    DO K = K1, K2
802                    DO J = J1, J2
803                    DO I = I1, I2
804                       IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
805                       IJK = FUNIJK(I,J,K)
806                       IJKE = EAST_OF(IJK)
807                       IPJK = IP_OF(IJK)
808     
809                       D_F = AVG_X_H(DIF(IJK),DIF(IJKE),I)*ODX_E(I)*AYZ(IJK)
810                       A_M(IJK,E,M) = A_M(IJK,E,M) - D_F
811                       A_M(IPJK,W,M) = A_M(IPJK,W,M) - D_F
812                    ENDDO
813                    ENDDO
814                    ENDDO
815     
816                 ELSEIF(IS_PLANE(L) == 'N') THEN
817                    DO K = K1, K2
818                    DO J = J1, J2
819                    DO I = I1, I2
820                       IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
821                       IJK = FUNIJK(I,J,K)
822                       IJKN = NORTH_OF(IJK)
823                       IJPK = JP_OF(IJK)
824     
825                       D_F = AVG_Y_H(DIF(IJK),DIF(IJKN),J)*ODY_N(J)*AXZ(IJK)
826                       A_M(IJK,N,M) = A_M(IJK,N,M) - D_F
827                       A_M(IJPK,S,M) = A_M(IJPK,S,M) - D_F
828                    ENDDO
829                    ENDDO
830                    ENDDO
831     
832                 ELSEIF(IS_PLANE(L) == 'T') THEN
833                    IF (DO_K) THEN
834                       DO K = K1, K2
835                       DO J = J1, J2
836                       DO I = I1, I2
837                          IJKT = TOP_OF(IJK)
838                          IJKP = KP_OF(IJK)
839     
840                          D_F = AVG_Z_H(DIF(IJK),DIF(IJKT),K)*&
841                             OX(I)*ODZ_T(K)*AXY(IJK)
842                          A_M(IJK,T,M) = A_M(IJK,T,M) - D_F
843                          A_M(IJKP,B,M) = A_M(IJKP,B,M) - D_F
844                       ENDDO
845                       ENDDO
846                       ENDDO
847                    ENDIF   ! end if do_k
848                 ENDIF   ! endif/else is_plane
849              ENDIF   ! endif is_defined
850           ENDDO   ! end do dimension_is
851     
852           RETURN
853           END SUBROUTINE DIF_PHI_IS
854