File: /nfs/home/0/users/jenkins/mfix.git/model/conv_dif_phi.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Module name:
4     !  CONV_DIF_Phi(Phi, Dif, Disc, Uf, Vf, Wf, Flux_E, Flux_N, Flux_T, M, A_m, B_m, IER)    C
5     !  Purpose: Determine convection diffusion terms for a sclar phi       C
6     !  The off-diagonal coefficients calculated here must be positive. The C
7     !  center coefficient and the source vector are negative;              C
8     
9     !  The diffusion at the flow boundaries is prevented by setting the
10     !  diffusion coefficients at boundary cells to zero and then using a
11     !  harmonic average to calculate the boundary diffusivity.  The value
12     !  diffusivities at the boundaries are checked in check_data_30.  Ensure
13     !  that harmonic avergaing is used in this routine.
14     !  See source_phi                                                      C
15     !                                                                      C
16     !  Author: M. Syamlal                                 Date: 21-APR-97  C
17     !  Reviewer:                                          Date:            C
18     !                                                                      C
19     !                                                                      C
20     !  Literature/Document References:                                     C
21     !                                                                      C
22     !  Variables referenced:                                               C
23     !  Variables modified:                                                 C
24     !                                                                      C
25     !  Local variables:                                                    C
26     !                                                                      C
27     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
28     !
29           SUBROUTINE CONV_DIF_PHI(PHI,DIF,DISC,UF,VF,WF,Flux_E,Flux_N,Flux_T,M,A_M,B_M,IER)
30     !...Translated by Pacific-Sierra Research VAST-90 2.06G5  12:17:31  12/09/98
31     !...Switches: -xf
32     !
33     !  Include param.inc file to specify parameter values
34     !
35     !-----------------------------------------------
36     !   M o d u l e s
37     !-----------------------------------------------
38           USE param
39           USE param1
40           USE run
41           USE geometry
42           USE compar
43           USE sendrecv
44           Use xsi_array
45           USE mpi_utility
46           USE indices
47           IMPLICIT NONE
48     !-----------------------------------------------
49     !   G l o b a l   P a r a m e t e r s
50     !-----------------------------------------------
51     !-----------------------------------------------
52     !   D u m m y   A r g u m e n t s
53     !-----------------------------------------------
54     !
55     !
56     !                      Scalar
57           DOUBLE PRECISION Phi(DIMENSION_3)
58     !
59     !                      Gamma -- diffusion coefficient
60           DOUBLE PRECISION Dif(DIMENSION_3)
61     !
62     !                      Discretizationindex
63           INTEGER          Disc
64     !
65     !                      Velocity components
66           DOUBLE PRECISION Uf(DIMENSION_3), Vf(DIMENSION_3), Wf(DIMENSION_3)
67     !
68     !                      Mass flux components
69           DOUBLE PRECISION Flux_E(DIMENSION_3), Flux_N(DIMENSION_3), Flux_T(DIMENSION_3)
70     !
71     !                      Phase index
72           INTEGER          M
73     !
74     !                      Septadiagonal matrix A_m
75           DOUBLE PRECISION A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
76     !
77     !                      Vector b_m
78           DOUBLE PRECISION B_m(DIMENSION_3, 0:DIMENSION_M)
79     !
80     !                      Error index
81           INTEGER          IER
82     
83     !
84     
85     !
86     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
87     !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
88     !       IF DEFERRED CORRECTION IS USED WITH THE SCALAR TRANSPORT EQN.
89     !
90             IF(DEF_COR)THEN
91               CALL CONV_DIF_PHI0(PHI,DIF,DISC,UF,VF,WF,Flux_E,Flux_N,Flux_T,M,A_M,B_M)
92               if (DISC > 1) CALL CONV_DIF_PHI_DC(PHI,DIF,DISC,UF,VF,WF,Flux_E,Flux_N,Flux_T,M,A_M,B_M,IER)
93             ELSE
94     !
95     !       NO DEFERRED CORRECTION IS USED WITH THE SCALAR TRANSPORT EQN.
96     !
97               IF (DISC == 0) THEN
98                 CALL CONV_DIF_PHI0(PHI,DIF,DISC,UF,VF,WF,Flux_E,Flux_N,Flux_T,M,A_M,B_M)
99               ELSE
100                 CALL CONV_DIF_PHI1(PHI,DIF,DISC,UF,VF,WF,Flux_E,Flux_N,Flux_T,M,A_M,B_M,IER)
101               ENDIF
102             ENDIF
103     
104             CALL DIF_PHI_IS (DIF, A_M, B_M, M)
105     
106             RETURN
107           END SUBROUTINE CONV_DIF_PHI
108     !
109     !
110     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
111     !                                                                      C
112     !  Module name:
113     !   CONV_DIF_Phi0(Phi, Dif, Disc, Uf, Vf, Wf, Flux_E,Flux_N,Flux_T, M, A_m, B_m, IER)
114     !  Purpose: Determine convection diffusion terms for Phi balance       C
115     !  The off-diagonal coefficients calculated here must be positive. The C
116     !  center coefficient and the source vector are negative;              C
117     !  See source_phi                                                      C
118     !                                                                      C
119     !  Author: M. Syamlal                                 Date: 21-APR-97  C
120     !  Reviewer:                                          Date:            C
121     !                                                                      C
122     !                                                                      C
123     !  Literature/Document References:                                     C
124     !                                                                      C
125     !  Variables referenced:                                               C
126     !  Variables modified:                                                 C
127     !                                                                      C
128     !  Local variables:                                                    C
129     !                                                                      C
130     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
131     !
132           SUBROUTINE CONV_DIF_PHI0(PHI,DIF,DISC,UF,VF,WF,Flux_E,Flux_N,Flux_T,M,A_M,B_M)
133     !...Translated by Pacific-Sierra Research VAST-90 2.06G5  12:17:31  12/09/98
134     !...Switches: -xf
135     !
136     !  Include param.inc file to specify parameter values
137     !
138     !-----------------------------------------------
139     !   M o d u l e s
140     !-----------------------------------------------
141           USE param
142           USE param1
143           USE parallel
144           USE matrix
145           USE toleranc
146           USE run
147           USE geometry
148           USE compar
149           USE sendrecv
150           USE indices
151           USE fun_avg
152           USE functions
153     !=======================================================================
154     ! JFD: START MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
155     !=======================================================================
156           USE cutcell
157     !=======================================================================
158     ! JFD: END MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
159     !=======================================================================
160           IMPLICIT NONE
161     !-----------------------------------------------
162     !   G l o b a l   P a r a m e t e r s
163     !-----------------------------------------------
164     !-----------------------------------------------
165     !   D u m m y   A r g u m e n t s
166     !-----------------------------------------------
167     !
168     !                      Scalar
169           DOUBLE PRECISION Phi(DIMENSION_3)
170     !
171     !                      Gamma -- diffusion coefficient
172           DOUBLE PRECISION Dif(DIMENSION_3)
173     !
174     !                      Discretizationindex
175           INTEGER          Disc
176     !
177     !                      Velocity components
178           DOUBLE PRECISION Uf(DIMENSION_3), Vf(DIMENSION_3), Wf(DIMENSION_3)
179     !
180     !                      Mass flux components
181           DOUBLE PRECISION Flux_E(DIMENSION_3), Flux_N(DIMENSION_3), Flux_T(DIMENSION_3)
182     !
183     !                      Phase index
184           INTEGER          M
185     !
186     !                      Septadiagonal matrix A_m
187           DOUBLE PRECISION A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
188     !
189     !                      Vector b_m
190           DOUBLE PRECISION B_m(DIMENSION_3, 0:DIMENSION_M)
191     !
192     !                      Indices
193           INTEGER          I,  J, K, IJK, IPJK, IJPK, IJKE, IJKN,&
194                            IJKP, IJKT
195     
196           INTEGER          IMJK, IM, IJKW
197           INTEGER          IJMK, JM, IJKS
198           INTEGER          IJKM, KM, IJKB
199     !
200     !                      Face velocity
201           DOUBLE PRECISION V_f
202     !
203     !                      Difusion parameter
204           DOUBLE PRECISION D_f
205     !
206     !-----------------------------------------------
207     !
208     !  Calculate convection-diffusion fluxes through each of the faces
209     !
210     !
211     !!!$omp      parallel do                                              &
212     !!!$omp&     private(I,  J, K,  IJK,  IPJK, IJPK, IJKE, IJKN,         &
213     !!!$omp&             IJKP, IJKT,  V_f, D_f,                    &
214     !!!$omp&             IMJK, IM, IJKW,                                  &
215     !!!$omp&             IJMK, JM, IJKS,                                  &
216     !!!$omp&             IJKM, KM,  IJKB)
217           DO IJK = ijkstart3, ijkend3
218     !
219            I = I_OF(IJK)
220            J = J_OF(IJK)
221            K = K_OF(IJK)
222     !
223              IF (FLUID_AT(IJK)) THEN
224     !
225                 IPJK = IP_OF(IJK)
226                 IJPK = JP_OF(IJK)
227                 IJKE = EAST_OF(IJK)
228                 IJKN = NORTH_OF(IJK)
229     !
230     !
231     !           East face (i+1/2, j, k)
232                 V_F = UF(IJK)
233                 D_F = AVG_X_H(DIF(IJK),DIF(IJKE),I)*ODX_E(I)*AYZ(IJK)
234     !=======================================================================
235     ! JFD: START MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
236     !=======================================================================
237                 IF(CUT_TREATMENT_AT(IJK)) THEN
238                    IF(CUT_CELL_AT(IJK).AND.(.NOT.FLUID_AT(IPJK))) THEN
239                       D_F = AVG_X_H(DIF(IJK),DIF(IJKE),I)*ODX_E(I)*DY(J)*DZ(K)
240                    ENDIF
241                 ENDIF
242     !=======================================================================
243     ! JFD: END MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
244     !=======================================================================
245     
246                 IF (V_F >= ZERO) THEN
247                    A_M(IJK,E,M) = D_F
248                    A_M(IPJK,W,M) = D_F + FLUX_E(IJK)
249                 ELSE
250                    A_M(IJK,E,M) = D_F - FLUX_E(IJK)
251                    A_M(IPJK,W,M) = D_F
252                 ENDIF
253     !
254     !
255     !           North face (i, j+1/2, k)
256                 V_F = VF(IJK)
257                 D_F = AVG_Y_H(DIF(IJK),DIF(IJKN),J)*ODY_N(J)*AXZ(IJK)
258     !=======================================================================
259     ! JFD: START MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
260     !=======================================================================
261                 IF(CUT_TREATMENT_AT(IJK)) THEN
262                    IF(CUT_CELL_AT(IJK).AND.(.NOT.FLUID_AT(IJPK))) THEN
263                       D_F = AVG_Y_H(DIF(IJK),DIF(IJKN),J)*ODY_N(J)*DX(I)*DZ(K)
264                    ENDIF
265                 ENDIF
266     !=======================================================================
267     ! JFD: END MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
268     !=======================================================================
269                 IF (V_F >= ZERO) THEN
270                    A_M(IJK,N,M) = D_F
271                    A_M(IJPK,S,M) = D_F + FLUX_N(IJK)
272                 ELSE
273                    A_M(IJK,N,M) = D_F - FLUX_N(IJK)
274                    A_M(IJPK,S,M) = D_F
275                 ENDIF
276     !
277     !           Top face (i, j, k+1/2)
278                 IF (DO_K) THEN
279                    IJKP = KP_OF(IJK)
280                    IJKT = TOP_OF(IJK)
281                    V_F = WF(IJK)
282                    D_F = AVG_Z_H(DIF(IJK),DIF(IJKT),K)*OX(I)*ODZ_T(K)*AXY(IJK)
283     !=======================================================================
284     ! JFD: START MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
285     !=======================================================================
286                    IF(CUT_TREATMENT_AT(IJK)) THEN
287                       IF(CUT_CELL_AT(IJK).AND.(.NOT.FLUID_AT(IJKP))) THEN
288                          D_F = AVG_Z_H(DIF(IJK),DIF(IJKT),K)*OX(I)*ODZ_T(K)*DX(I)*DY(J)
289                       ENDIF
290                    ENDIF
291     !=======================================================================
292     ! JFD: END MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
293     !=======================================================================
294                    IF (V_F >= ZERO) THEN
295                       A_M(IJK,T,M) = D_F
296                       A_M(IJKP,B,M) = D_F + FLUX_T(IJK)
297                    ELSE
298                       A_M(IJK,T,M) = D_F - FLUX_T(IJK)
299                       A_M(IJKP,B,M) = D_F
300                    ENDIF
301                 ENDIF
302     !
303     !
304     !           West face (i-1/2, j, k)
305                 IMJK = IM_OF(IJK)
306                 IF (.NOT.FLUID_AT(IMJK)) THEN
307                    IM = IM1(I)
308                    IJKW = WEST_OF(IJK)
309                    V_F = UF(IMJK)
310                    D_F = AVG_X_H(DIF(IJKW),DIF(IJK),IM)*ODX_E(IM)*AYZ(IMJK)
311     !=======================================================================
312     ! JFD: START MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
313     !=======================================================================
314                    IF(CUT_TREATMENT_AT(IJK)) THEN
315                       IF(CUT_CELL_AT(IJK).AND.(.NOT.FLUID_AT(IMJK))) THEN
316                          D_F = AVG_X_H(DIF(IJKW),DIF(IJK),IM)*ODX_E(IM)*DY(J)*DZ(K)
317                       ENDIF
318                    ENDIF
319     !=======================================================================
320     ! JFD: END MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
321     !=======================================================================
322                    IF (V_F >= ZERO) THEN
323                       A_M(IJK,W,M) = D_F + FLUX_E(IMJK)
324                    ELSE
325                       A_M(IJK,W,M) = D_F
326                    ENDIF
327                 ENDIF
328     !
329     !           South face (i, j-1/2, k)
330                 IJMK = JM_OF(IJK)
331                 IF (.NOT.FLUID_AT(IJMK)) THEN
332                    JM = JM1(J)
333                    IJKS = SOUTH_OF(IJK)
334                    V_F = VF(IJMK)
335                    D_F = AVG_Y_H(DIF(IJKS),DIF(IJK),JM)*ODY_N(JM)*AXZ(IJMK)
336     !=======================================================================
337     ! JFD: START MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
338     !=======================================================================
339                    IF(CUT_TREATMENT_AT(IJK)) THEN
340                       IF(CUT_CELL_AT(IJK).AND.(.NOT.FLUID_AT(IJMK))) THEN
341                          D_F = AVG_Y_H(DIF(IJKS),DIF(IJK),JM)*ODY_N(JM)*DX(I)*DZ(K)
342                       ENDIF
343                    ENDIF
344     !=======================================================================
345     ! JFD: END MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
346     !=======================================================================
347                    IF (V_F >= ZERO) THEN
348                       A_M(IJK,S,M) = D_F + FLUX_N(IJMK)
349                    ELSE
350                       A_M(IJK,S,M) = D_F
351                    ENDIF
352                 ENDIF
353     !
354     !           Bottom face (i, j, k-1/2)
355                 IF (DO_K) THEN
356                    IJKM = KM_OF(IJK)
357                    IF (.NOT.FLUID_AT(IJKM)) THEN
358                       KM = KM1(K)
359                       IJKB = BOTTOM_OF(IJK)
360                       V_F = WF(IJKM)
361                       D_F = AVG_Z_H(DIF(IJKB),DIF(IJK),KM)*OX(I)*ODZ_T(KM)*AXY(&
362                          IJKM)
363     !=======================================================================
364     ! JFD: START MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
365     !=======================================================================
366                       IF(CUT_TREATMENT_AT(IJK)) THEN
367                          IF(CUT_CELL_AT(IJK).AND.(.NOT.FLUID_AT(IJKM))) THEN
368                             D_F = AVG_Z_H(DIF(IJKB),DIF(IJK),KM)*OX(I)*ODZ_T(KM)*DX(I)*DY(J)
369                          ENDIF
370                       ENDIF
371     !=======================================================================
372     ! JFD: END MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
373     !=======================================================================
374                       IF (V_F >= ZERO) THEN
375                          A_M(IJK,B,M) = D_F + FLUX_T(IJKM)
376                       ELSE
377                          A_M(IJK,B,M) = D_F
378                       ENDIF
379                    ENDIF
380                 ENDIF
381     !
382              ENDIF
383           END DO
384     !
385           RETURN
386           END SUBROUTINE CONV_DIF_PHI0
387     
388     !
389     !
390     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
391     !                                                                      C
392     !  Module name:
393     !    CONV_DIF_Phi_DC(Phi, Dif, Disc, Uf, Vf, Wf, Flux_E,Flux_N,Flux_T, M, A_m, B_m, IER)
394     !  Purpose: TO USE DEFERRED CORRECTION IN SOLVING THE SCALAR TRANSPORT C
395     !  EQN. THIS METHOD COMBINES FIRST ORDER UPWIND AND A USER SPECIFIED   C
396     !  HIGH ORDER METHOD TO SOLVE FOR THE SCALAR PHI.
397     !  See source_Phi                                                  C
398     !                                                                      C
399     !  Author: C. GUENTHER                                Date: 1-ARP-99   C
400     !  Reviewer:                                          Date:            C
401     !                                                                      C
402     !                                                                      C
403     !  Literature/Document References:                                     C
404     !                                                                      C
405     !  Variables referenced:                                               C
406     !  Variables modified:                                                 C
407     !                                                                      C
408     !  Local variables:                                                    C
409     !                                                                      C
410     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
411     !
412           SUBROUTINE CONV_DIF_PHI_DC(PHI,DIF,DISC,UF,VF,WF,Flux_E,Flux_N,Flux_T,M,A_M,B_M,IER)
413     !...Translated by Pacific-Sierra Research VAST-90 2.06G5  12:17:31  12/09/98
414     !...Switches: -xf
415     !
416     !  Include param.inc file to specify parameter values
417     !
418     !-----------------------------------------------
419     !   M o d u l e s
420     !-----------------------------------------------
421           USE compar
422           USE discretization, ONLY: fpfoi_of
423           USE fun_avg
424           USE function3
425           USE functions
426           USE geometry
427           USE indices
428           USE matrix
429           USE parallel
430           USE param
431           USE param1
432           USE run
433           USE sendrecv
434           USE sendrecv3
435           USE toleranc
436           USE xsi
437           Use tmp_array
438           Use xsi_array
439           IMPLICIT NONE
440     !-----------------------------------------------
441     !   G l o b a l   P a r a m e t e r s
442     !-----------------------------------------------
443     !-----------------------------------------------
444     !   D u m m y   A r g u m e n t s
445     !-----------------------------------------------
446     !
447     !                      Scalar
448           DOUBLE PRECISION Phi(DIMENSION_3)
449     !
450     !                      Gamma -- diffusion coefficient
451           DOUBLE PRECISION Dif(DIMENSION_3)
452     !
453     !                      Discretizationindex
454           INTEGER          Disc
455     !
456     !                      Velocity components
457           DOUBLE PRECISION Uf(DIMENSION_3), Vf(DIMENSION_3), Wf(DIMENSION_3)
458     !
459     !                      Mass flux components
460           DOUBLE PRECISION Flux_E(DIMENSION_3), Flux_N(DIMENSION_3), Flux_T(DIMENSION_3)
461     !
462     !                      Phase index
463           INTEGER          M
464     !
465     !                      Septadiagonal matrix A_m
466           DOUBLE PRECISION A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
467     !
468     !                      Vector b_m
469           DOUBLE PRECISION B_m(DIMENSION_3, 0:DIMENSION_M)
470     !
471     !                      Error index
472           INTEGER          IER
473     !
474     !                      Indices
475           INTEGER          I,  J, K, IJK, IPJK, IJPK, IJKE, IJKN,&
476                            IJKP, IJKT
477     
478           INTEGER          IMJK, IJKW
479           INTEGER          IJMK, IJKS
480           INTEGER          IJKM, IJKB
481           INTEGER          IJK4, IPPP, IPPP4, JPPP, JPPP4, KPPP, KPPP4
482           INTEGER          IMMM, IMMM4, JMMM, JMMM4, KMMM, KMMM4
483     
484     ! loezos
485           INTEGER  incr
486     ! loezos
487     
488     !       FACE VELOCITY
489             DOUBLE PRECISION V_F
490     !
491     !       DEFERRED CORRCTION CONTRIBUTION FORM HIGH ORDER METHOD
492             DOUBLE PRECISION PHI_HO
493     !
494     !       LOW ORDER APPROXIMATION
495             DOUBLE PRECISION PHI_LO
496     !
497     !       DEFERRED CORRECTION CONTRIBUTIONS FROM EACH FACE
498             DOUBLE PRECISION  EAST_DC
499             DOUBLE PRECISION  WEST_DC
500             DOUBLE PRECISION  NORTH_DC
501             DOUBLE PRECISION  SOUTH_DC
502             DOUBLE PRECISION  TOP_DC
503             DOUBLE PRECISION  BOTTOM_DC
504     !
505     !---------------------------------------------------------------
506     
507           call lock_xsi_array
508           call lock_tmp4_array
509     !
510     !  Calculate convection factors
511     !
512     !
513     ! Send recv the third ghost layer
514           IF ( FPFOI ) THEN
515              Do IJK = ijkstart3, ijkend3
516                 I = I_OF(IJK)
517                 J = J_OF(IJK)
518                 K = K_OF(IJK)
519                 IJK4 = funijk3(I,J,K)
520                 TMP4(IJK4) = PHI(IJK)
521              ENDDO
522              CALL send_recv3(tmp4)
523           ENDIF
524     
525     ! loezos
526             incr=0
527     ! loezos
528     
529           CALL CALC_XSI (DISC, PHI, UF, VF, WF, XSI_E, XSI_N, XSI_T,incr)
530     !
531     !
532     !  Calculate convection-diffusion fluxes through each of the faces
533     !
534     !
535     !!!$omp      parallel do                                              &
536     !!!$omp&     private(I,  J, K,  IJK,  IPJK, IJPK, IJKE, IJKN,         &
537     !!!$omp&             IJKP, IJKT,  V_f, D_f,                    &
538     !!!$omp&             IMJK, IJKW,                                  &
539     !!!$omp&             IJMK, IJKS,                                  &
540     !!!$omp&             IJKM, IJKB, PHI_HO, PHI_LO,        &
541     !!!$omp&             EAST_DC, WEST_DC, NORTH_DC, SOUTH_DC, TOP_DC, BOTTOM_DC)
542     !
543           DO IJK = ijkstart3, ijkend3
544     !
545     ! Determine whether IJK falls within 1 ghost layer........
546            I = I_OF(IJK)
547            J = J_OF(IJK)
548            K = K_OF(IJK)
549     !
550              IF (FLUID_AT(IJK)) THEN
551     !
552     !
553                 IPJK = IP_OF(IJK)
554                 IMJK = IM_OF(IJK)
555                 IJPK = JP_OF(IJK)
556                 IJMK = JM_OF(IJK)
557                 IJKP = KP_OF(IJK)
558                 IJKM = KM_OF(IJK)
559                 IJKE = EAST_OF(IJK)
560                 IJKN = NORTH_OF(IJK)
561     !
562     !           Third Ghost layer information
563                 IPPP  = IP_OF(IP_OF(IPJK))
564                 IPPP4 = funijk3(I_OF(IPPP), J_OF(IPPP), K_OF(IPPP))
565                 IMMM  = IM_OF(IM_OF(IMJK))
566                 IMMM4 = funijk3(I_OF(IMMM), J_OF(IMMM), K_OF(IMMM))
567                 JPPP  = JP_OF(JP_OF(IJPK))
568                 JPPP4 = funijk3(I_OF(JPPP), J_OF(JPPP), K_OF(JPPP))
569                 JMMM  = JM_OF(JM_OF(IJMK))
570                 JMMM4 = funijk3(I_OF(JMMM), J_OF(JMMM), K_OF(JMMM))
571                 KPPP  = KP_OF(KP_OF(IJKP))
572                 KPPP4 = funijk3(K_OF(IPPP), J_OF(KPPP), K_OF(KPPP))
573                 KMMM  = KM_OF(KM_OF(IJKM))
574                 KMMM4 = funijk3(I_OF(KMMM), J_OF(KMMM), K_OF(KMMM))
575     !
576     !
577     !           DEFERRED CORRECTION CONTRIBUTION AT THE East face (i+1/2, j, k)
578     !
579                     V_F = UF(IJK)
580                     IF(V_F >= ZERO)THEN
581                        PHI_LO = PHI(IJK)
582                        IF ( FPFOI ) &
583                           PHI_HO = FPFOI_OF(PHI(IPJK), PHI(IJK), &
584                                 PHI(IMJK), PHI(IM_OF(IMJK)))
585                     ELSE
586                        PHI_LO = PHI(IPJK)
587                        IF ( FPFOI ) &
588                           PHI_HO = FPFOI_OF(PHI(IJK), PHI(IPJK), &
589                                 PHI(IP_OF(IPJK)), TMP4(IPPP4))
590                     ENDIF
591                     IF (.NOT. FPFOI ) &
592                           PHI_HO = XSI_E(IJK)*PHI(IPJK)+(1.0-XSI_E(IJK))*PHI(IJK)
593                     EAST_DC = FLUX_E(IJK)*(PHI_LO - PHI_HO)
594     !
595     !
596     !           DEFERRED CORRECTION CONTRIBUTION AT THE North face (i, j+1/2, k)
597     !
598                     V_F = VF(IJK)
599                     IF(V_F >= ZERO)THEN
600                        PHI_LO = PHI(IJK)
601                        IF ( FPFOI ) &
602                           PHI_HO = FPFOI_OF(PHI(IJPK), PHI(IJK), &
603                                 PHI(IJMK), PHI(JM_OF(IJMK)))
604                     ELSE
605                        PHI_LO = PHI(IJPK)
606                        IF ( FPFOI ) &
607                           PHI_HO = FPFOI_OF(PHI(IJK), PHI(IJPK), &
608                                 PHI(JP_OF(IJPK)), TMP4(JPPP4))
609                     ENDIF
610                     IF (.NOT. FPFOI ) &
611                          PHI_HO = XSI_N(IJK)*PHI(IJPK)+(1.0-XSI_N(IJK))*PHI(IJK)
612                     NORTH_DC = FLUX_N(IJK)*(PHI_LO - PHI_HO)
613     !
614     !
615     !           DEFERRED CORRECTION CONTRIBUTION AT THE Top face (i, j, k+1/2)
616     !
617                   IF (DO_K) THEN
618                     IJKP = KP_OF(IJK)
619                     IJKT = TOP_OF(IJK)
620                     V_F = WF(IJK)
621                     IF(V_F >= ZERO)THEN
622                        PHI_LO = PHI(IJK)
623                        IF ( FPFOI ) &
624                           PHI_HO = FPFOI_OF(PHI(IJKP),  PHI(IJK), &
625                                 PHI(IJKM), PHI(KM_OF(IJKM)))
626                     ELSE
627                        PHI_LO = PHI(IJKP)
628                        IF ( FPFOI ) &
629                           PHI_HO = FPFOI_OF(PHI(IJK), PHI(IJKP),  &
630                                 PHI(KP_OF(IJKP)), TMP4(KPPP4))
631                     ENDIF
632                     IF (.NOT. FPFOI ) &
633                          PHI_HO = XSI_T(IJK)*PHI(IJKP)+(1.0-XSI_T(IJK))*PHI(IJK)
634                     TOP_DC = FLUX_T(IJK)*(PHI_LO - PHI_HO)
635                   ELSE
636                       TOP_DC = ZERO
637                   ENDIF
638     !
639     !
640     !           DEFERRED CORRECTION CONTRIBUTION AT THE West face (i-1/2, j, k)
641     !
642                     IMJK = IM_OF(IJK)
643                     IJKW = WEST_OF(IJK)
644                     V_F = UF(IMJK)
645                     IF(V_F >= ZERO)THEN
646                        PHI_LO = PHI(IMJK)
647                        IF ( FPFOI ) &
648                           PHI_HO = FPFOI_OF(PHI(IJK), PHI(IMJK), &
649                                 PHI(IM_OF(IMJK)), TMP4(IMMM4))
650                     ELSE
651                        PHI_LO = PHI(IJK)
652                        IF ( FPFOI ) &
653                           PHI_HO = FPFOI_OF(PHI(IMJK), PHI(IJK), &
654                                 PHI(IPJK), PHI(IP_OF(IPJK)))
655                     ENDIF
656                     IF (.NOT. FPFOI ) &
657                           PHI_HO = XSI_E(IMJK)*PHI(IJK)+(ONE-XSI_E(IMJK))*PHI(IMJK)
658                     WEST_DC = FLUX_E(IMJK)*(PHI_LO - PHI_HO)
659     !
660     !
661     !           DEFERRED CORRECTION CONTRIBUTION AT THE South face (i, j-1/2, k)
662     !
663                     IJMK = JM_OF(IJK)
664                     IJKS = SOUTH_OF(IJK)
665                     V_F = VF(IJMK)
666                     IF(V_F >= ZERO)THEN
667                        PHI_LO = PHI(IJMK)
668                        IF ( FPFOI ) &
669                           PHI_HO = FPFOI_OF(PHI(IJK), PHI(IJMK), &
670                                 PHI(JM_OF(IJMK)), TMP4(JMMM4))
671                     ELSE
672                        PHI_LO = PHI(IJK)
673                        IF ( FPFOI ) &
674                           PHI_HO = FPFOI_OF(PHI(IJMK), PHI(IJK), &
675                                 PHI(IJPK), PHI(JP_OF(IJPK)))
676                     ENDIF
677                     IF (.NOT. FPFOI ) &
678                           PHI_HO = XSI_N(IJMK)*PHI(IJK)+(ONE-XSI_N(IJMK))*PHI(IJMK)
679                     SOUTH_DC = FLUX_N(IJMK)*(PHI_LO - PHI_HO)
680     !
681     !
682     !           DEFERRED CORRECTION CONTRIBUTION AT THE Bottom face (i, j, k-1/2)
683                   IF (DO_K) THEN
684                      IJKM = KM_OF(IJK)
685                      IJKB = BOTTOM_OF(IJK)
686                      V_F = WF(IJKM)
687                      IF(V_F >= ZERO)THEN
688                        PHI_LO = PHI(IJKM)
689                        IF ( FPFOI ) &
690                           PHI_HO = FPFOI_OF(PHI(IJK), PHI(IJKM), &
691                                 PHI(KM_OF(IJKM)), TMP4(KMMM4))
692                      ELSE
693                        PHI_LO = PHI(IJK)
694                        IF ( FPFOI ) &
695                           PHI_HO = FPFOI_OF(PHI(IJKM), PHI(IJK), &
696                                 PHI(IJKP), PHI(KP_OF(IJKP)))
697                      ENDIF
698                      IF (.NOT. FPFOI ) &
699                           PHI_HO = XSI_T(IJKM)*PHI(IJK)+(1.0-XSI_T(IJKM))*PHI(IJKM)
700                      BOTTOM_DC = FLUX_T(IJKM)*(PHI_LO - PHI_HO)
701                   ELSE
702                        BOTTOM_DC = ZERO
703                   ENDIF
704     
705     !
706     !           CONTRIBUTION DUE TO DEFERRED CORRECTION
707     !
708                     B_M(IJK,M) = B_M(IJK,M)+WEST_DC-EAST_DC+SOUTH_DC-NORTH_DC&
709                                     +BOTTOM_DC-TOP_DC
710     !
711              ENDIF
712           END DO
713           call unlock_tmp4_array
714           call unlock_xsi_array
715     !
716     !
717           RETURN
718           END SUBROUTINE CONV_DIF_PHI_DC
719     !
720     !
721     
722     !
723     !
724     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
725     !                                                                      C
726     !  Module name:
727     !    CONV_DIF_Phi1(Phi, Dif, Disc, Uf, Vf, Wf, Flux_E,Flux_N,Flux_T, M, A_m, B_m, IER)
728     !  Purpose: Determine convection diffusion terms for gas energy eq Phi C
729     !  The off-diagonal coefficients calculated here must be positive. The C
730     !  center coefficient and the source vector are negative;              C
731     !  See source_Phi                                                  C
732     !                                                                      C
733     !  Author: M. Syamlal                                 Date: 21-APR-97  C
734     !  Reviewer:                                          Date:            C
735     !                                                                      C
736     !                                                                      C
737     !  Literature/Document References:                                     C
738     !                                                                      C
739     !  Variables referenced:                                               C
740     !  Variables modified:                                                 C
741     !                                                                      C
742     !  Local variables:                                                    C
743     !                                                                      C
744     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
745     !
746           SUBROUTINE CONV_DIF_PHI1(PHI,DIF,DISC,UF,VF,WF,Flux_E,Flux_N,Flux_T,M,A_M,B_M,IER)
747     !...Translated by Pacific-Sierra Research VAST-90 2.06G5  12:17:31  12/09/98
748     !...Switches: -xf
749     !
750     !  Include param.inc file to specify parameter values
751     !
752     !-----------------------------------------------
753     !   M o d u l e s
754     !-----------------------------------------------
755           USE param
756           USE param1
757           USE parallel
758           USE matrix
759           USE toleranc
760           USE run
761           USE geometry
762           USE compar
763           USE sendrecv
764           USE indices
765           USE vshear
766           USE xsi
767           USE xsi_array
768           USE fun_avg
769           USE functions
770     !=======================================================================
771     ! JFD: START MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
772     !=======================================================================
773           USE cutcell
774     !=======================================================================
775     ! JFD: END MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
776     !=======================================================================
777           IMPLICIT NONE
778     !-----------------------------------------------
779     !   G l o b a l   P a r a m e t e r s
780     !-----------------------------------------------
781     !-----------------------------------------------
782     !   D u m m y   A r g u m e n t s
783     !-----------------------------------------------
784     !
785     !                      Scalar
786           DOUBLE PRECISION Phi(DIMENSION_3)
787     !
788     !                      Gamma -- diffusion coefficient
789           DOUBLE PRECISION Dif(DIMENSION_3)
790     !
791     !                      Discretizationindex
792           INTEGER          Disc
793     !
794     !                      Velocity components
795           DOUBLE PRECISION Uf(DIMENSION_3), Vf(DIMENSION_3), Wf(DIMENSION_3)
796     !
797     !                      Mass flux components
798           DOUBLE PRECISION Flux_E(DIMENSION_3), Flux_N(DIMENSION_3), Flux_T(DIMENSION_3)
799     !
800     !                      Phase index
801           INTEGER          M
802     !
803     !                      Septadiagonal matrix A_m
804           DOUBLE PRECISION A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
805     !
806     !                      Vector b_m
807           DOUBLE PRECISION B_m(DIMENSION_3, 0:DIMENSION_M)
808     !
809     !                      Error index
810           INTEGER          IER
811     !
812     !                      Indices
813           INTEGER          I,  J, K, IJK, IPJK, IJPK, IJKE, IJKN,&
814                            IJKP, IJKT
815     
816           INTEGER          IMJK, IM, IJKW
817           INTEGER          IJMK, JM, IJKS
818           INTEGER          IJKM, KM, IJKB
819     ! start loezos
820           INTEGER incr
821     ! end loezos
822     
823     !
824     !                      Difusion parameter
825           DOUBLE PRECISION D_f
826     !
827     !                      Convection weighting factors
828     !      DOUBLE PRECISION XSI_e(DIMENSION_3), XSI_n(DIMENSION_3),&
829     !                       XSI_t(DIMENSION_3)
830     !-----------------------------------------------
831           call lock_xsi_array
832     !
833     !  Calculate convection factors
834     !
835     !
836     
837     ! loezos
838             incr=0
839     ! loezos
840     
841           CALL CALC_XSI (DISC, PHI, UF, VF, WF, XSI_E, XSI_N, XSI_T,incr)
842     
843     ! loezos
844     !update V to true velocity
845     
846           IF (SHEAR) THEN
847              DO IJK = ijkstart3, ijkend3
848              IF (FLUID_AT(IJK)) THEN
849                VF(IJK)=VF(IJK)+VSH(IJK)
850               END IF
851             END DO
852           END IF
853     ! loezos
854     
855     !
856     !
857     !  Calculate convection-diffusion fluxes through each of the faces
858     !
859     !!!$omp      parallel do                                               &
860     !!!$omp&     private(I,  J, K,  IJK,  IPJK, IJPK, IJKE, IJKN,          &
861     !!!$omp&             IJKP, IJKT,    D_f,                        &
862     !!!$omp&             IMJK, IM, IJKW,                                   &
863     !!!$omp&             IJMK, JM,  IJKS,                                  &
864     !!!$omp&             IJKM, KM,  IJKB )
865     !
866     !
867     !
868           DO IJK = ijkstart3, ijkend3
869     !
870            I = I_OF(IJK)
871            J = J_OF(IJK)
872            K = K_OF(IJK)
873     !
874              IF (FLUID_AT(IJK)) THEN
875     !
876                 IPJK = IP_OF(IJK)
877                 IJPK = JP_OF(IJK)
878                 IJKE = EAST_OF(IJK)
879                 IJKN = NORTH_OF(IJK)
880     !
881     !
882     !           East face (i+1/2, j, k)
883                 D_F = AVG_X_H(DIF(IJK),DIF(IJKE),I)*ODX_E(I)*AYZ(IJK)
884     !=======================================================================
885     ! JFD: START MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
886     !=======================================================================
887                 IF(CUT_TREATMENT_AT(IJK)) THEN
888                    IF(CUT_CELL_AT(IJK).AND.(.NOT.FLUID_AT(IPJK))) THEN
889                       D_F = AVG_X_H(DIF(IJK),DIF(IJKE),I)*ODX_E(I)*DY(J)*DZ(K)
890                    ENDIF
891                 ENDIF
892     !=======================================================================
893     ! JFD: END MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
894     !=======================================================================
895     !
896                 A_M(IJK,E,M) = D_F - XSI_E(IJK)*FLUX_E(IJK)
897     !
898                 A_M(IPJK,W,M) = D_F + (ONE - XSI_E(IJK))*FLUX_E(IJK)
899     !
900     !
901     !           North face (i, j+1/2, k)
902                 D_F = AVG_Y_H(DIF(IJK),DIF(IJKN),J)*ODY_N(J)*AXZ(IJK)
903     !=======================================================================
904     ! JFD: START MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
905     !=======================================================================
906                 IF(CUT_TREATMENT_AT(IJK)) THEN
907                    IF(CUT_CELL_AT(IJK).AND.(.NOT.FLUID_AT(IJPK))) THEN
908                       D_F = AVG_Y_H(DIF(IJK),DIF(IJKN),J)*ODY_N(J)*DX(I)*DZ(K)
909                    ENDIF
910                 ENDIF
911     !=======================================================================
912     ! JFD: END MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
913     !=======================================================================
914     !
915                 A_M(IJK,N,M) = D_F - XSI_N(IJK)*FLUX_N(IJK)
916     !
917                 A_M(IJPK,S,M) = D_F + (ONE - XSI_N(IJK))*FLUX_N(IJK)
918     !
919     !
920     !           Top face (i, j, k+1/2)
921                 IF (DO_K) THEN
922                    IJKP = KP_OF(IJK)
923                    IJKT = TOP_OF(IJK)
924     !
925                    D_F = AVG_Z_H(DIF(IJK),DIF(IJKT),K)*OX(I)*ODZ_T(K)*AXY(IJK)
926     !=======================================================================
927     ! JFD: START MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
928     !=======================================================================
929                    IF(CUT_TREATMENT_AT(IJK)) THEN
930                       IF(CUT_CELL_AT(IJK).AND.(.NOT.FLUID_AT(IJKP))) THEN
931                          D_F = AVG_Z_H(DIF(IJK),DIF(IJKT),K)*OX(I)*ODZ_T(K)*DX(I)*DY(J)
932                       ENDIF
933                    ENDIF
934     !=======================================================================
935     ! JFD: END MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
936     !=======================================================================
937     !
938                    A_M(IJK,T,M) = D_F - XSI_T(IJK)*FLUX_T(IJK)
939     !
940                    A_M(IJKP,B,M)=D_F+(ONE-XSI_T(IJK))*FLUX_T(IJK)
941                 ENDIF
942     !
943     !           West face (i-1/2, j, k)
944                 IMJK = IM_OF(IJK)
945                 IF (.NOT.FLUID_AT(IMJK)) THEN
946                    IM = IM1(I)
947                    IJKW = WEST_OF(IJK)
948     !
949                    D_F = AVG_X_H(DIF(IJKW),DIF(IJK),IM)*ODX_E(IM)*AYZ(IMJK)
950     !=======================================================================
951     ! JFD: START MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
952     !=======================================================================
953                    IF(CUT_TREATMENT_AT(IJK)) THEN
954                       IF(CUT_CELL_AT(IJK).AND.(.NOT.FLUID_AT(IMJK))) THEN
955                          D_F = AVG_X_H(DIF(IJKW),DIF(IJK),IM)*ODX_E(IM)*DY(J)*DZ(K)
956                       ENDIF
957                    ENDIF
958     !=======================================================================
959     ! JFD: END MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
960     !=======================================================================
961     
962     !
963                    A_M(IJK,W,M) = D_F + (ONE - XSI_E(IMJK))*FLUX_E(IMJK)
964                 ENDIF
965     !
966     !           South face (i, j-1/2, k)
967                 IJMK = JM_OF(IJK)
968                 IF (.NOT.FLUID_AT(IJMK)) THEN
969                    JM = JM1(J)
970                    IJKS = SOUTH_OF(IJK)
971     !
972                    D_F = AVG_Y_H(DIF(IJKS),DIF(IJK),JM)*ODY_N(JM)*AXZ(IJMK)
973     !=======================================================================
974     ! JFD: START MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
975     !=======================================================================
976                    IF(CUT_TREATMENT_AT(IJK)) THEN
977                       IF(CUT_CELL_AT(IJK).AND.(.NOT.FLUID_AT(IJMK))) THEN
978                          D_F = AVG_Y_H(DIF(IJKS),DIF(IJK),JM)*ODY_N(JM)*DX(I)*DZ(K)
979                       ENDIF
980                    ENDIF
981     !=======================================================================
982     ! JFD: END MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
983     !=======================================================================
984     
985     !
986                    A_M(IJK,S,M) = D_F + (ONE - XSI_N(IJMK))*FLUX_N(IJMK)
987                 ENDIF
988     !
989     !           Bottom face (i, j, k-1/2)
990                 IF (DO_K) THEN
991                    IJKM = KM_OF(IJK)
992                    IF (.NOT.FLUID_AT(IJKM)) THEN
993                       KM = KM1(K)
994                       IJKB = BOTTOM_OF(IJK)
995     !
996                       D_F = AVG_Z_H(DIF(IJKB),DIF(IJK),KM)*OX(I)*ODZ_T(KM)*AXY(&
997                          IJKM)
998     !=======================================================================
999     ! JFD: START MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
1000     !=======================================================================
1001                       IF(CUT_TREATMENT_AT(IJK)) THEN
1002                          IF(CUT_CELL_AT(IJK).AND.(.NOT.FLUID_AT(IJKM))) THEN
1003                             D_F = AVG_Z_H(DIF(IJKB),DIF(IJK),KM)*OX(I)*ODZ_T(KM)*DX(I)*DY(J)
1004                          ENDIF
1005                       ENDIF
1006     !=======================================================================
1007     ! JFD: END MODIFICATION FOR CARTESIAN GRID IMPLEMENTATION
1008     !=======================================================================
1009     !
1010                       A_M(IJK,B,M) = D_F + (ONE - XSI_T(IJKM))*FLUX_T(IJKM)
1011                    ENDIF
1012                 ENDIF
1013     !
1014              ENDIF
1015           END DO
1016     
1017     ! loezos
1018            IF (SHEAR) THEN
1019              DO IJK = ijkstart3, ijkend3
1020               IF (FLUID_AT(IJK)) THEN
1021                VF(IJK)=VF(IJK)-VSH(IJK)
1022               END IF
1023              END DO
1024             END IF
1025     ! loezos
1026           call unlock_xsi_array
1027     
1028     !
1029           RETURN
1030           END SUBROUTINE CONV_DIF_PHI1
1031     !
1032     !
1033     !
1034     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
1035     !                                                                      C
1036     !  Module name: DIF_phi_IS(Dif, A_m, B_m, M, IER)                      C
1037     !  Purpose: Remove diffusive fluxes across internal surfaces.          C
1038     !                                                                      C
1039     !  Author: M. Syamlal                                 Date: 30-APR-97  C
1040     !  Reviewer:                                          Date:            C
1041     !                                                                      C
1042     !                                                                      C
1043     !  Literature/Document References:                                     C
1044     !                                                                      C
1045     !  Variables referenced:                                               C
1046     !  Variables modified:                                                 C
1047     !                                                                      C
1048     !  Local variables:                                                    C
1049     !                                                                      C
1050     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
1051     !
1052           SUBROUTINE DIF_PHI_IS(DIF, A_M, B_M, M)
1053     !...Translated by Pacific-Sierra Research VAST-90 2.06G5  12:17:31  12/09/98
1054     !...Switches: -xf
1055     !
1056     !  Include param.inc file to specify parameter values
1057     !
1058     !-----------------------------------------------
1059     !   M o d u l e s
1060     !-----------------------------------------------
1061           USE param
1062           USE param1
1063           USE parallel
1064           USE matrix
1065           USE toleranc
1066           USE run
1067           USE geometry
1068           USE compar
1069           USE sendrecv
1070           USE indices
1071           USE scales
1072           USE constant
1073           USE physprop
1074           USE fldvar
1075           USE visc_s
1076           USE output
1077           USE is
1078           USE fun_avg
1079           USE function3
1080           USE functions
1081           IMPLICIT NONE
1082     !-----------------------------------------------
1083     !   G l o b a l   P a r a m e t e r s
1084     !-----------------------------------------------
1085     !-----------------------------------------------
1086     !   D u m m y   A r g u m e n t s
1087     !-----------------------------------------------
1088     !
1089     !                      Internal surface
1090           INTEGER          L
1091     !
1092     !                      Indices
1093           INTEGER          I,  J, K, I1, I2, J1, J2, K1, K2, IJK,&
1094                            IJKE, IJKN, IJKT, IPJK, IJPK, IJKP
1095     !
1096     !                      Solids phase
1097           INTEGER          M
1098     !
1099     !                      Gamma -- diffusion coefficient
1100           DOUBLE PRECISION Dif(DIMENSION_3)
1101     !
1102     !                      Septadiagonal matrix A_m
1103           DOUBLE PRECISION A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
1104     !
1105     !                      Vector b_m
1106           DOUBLE PRECISION B_m(DIMENSION_3, 0:DIMENSION_M)
1107     !
1108     !                      Difusion parameter
1109           DOUBLE PRECISION D_f
1110     !-----------------------------------------------
1111     !
1112     ! Make user defined internal surfaces non-conducting
1113     !
1114           DO L = 1, DIMENSION_IS
1115              IF (IS_DEFINED(L)) THEN
1116                 I1 = IS_I_W(L)
1117                 I2 = IS_I_E(L)
1118                 J1 = IS_J_S(L)
1119                 J2 = IS_J_N(L)
1120                 K1 = IS_K_B(L)
1121                 K2 = IS_K_T(L)
1122     
1123     !       Limit I1, I2 and all to local processor first ghost layer
1124     
1125                 IF(I1.LE.IEND2)   I1 = MAX(I1, ISTART2)
1126     
1127                 IF(J1.LE.JEND2)   J1 = MAX(J1, JSTART2)
1128     
1129                 IF(K1.LE.KEND2)   K1 = MAX(K1, KSTART2)
1130     
1131                 IF(I2.GE.ISTART2) I2 = MIN(I2, IEND2)
1132     
1133                 IF(J2.GE.JSTART2) J2 = MIN(J2, JEND2)
1134     
1135                 IF(K2.GE.KSTART2) K2 = MIN(K2, KEND2)
1136     
1137     !     End of limiting to the first ghost cells of the processor....
1138                 DO K = K1, K2
1139                    DO J = J1, J2
1140                       DO I = I1, I2
1141                          IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
1142                          IJK = FUNIJK(I,J,K)
1143     !
1144                          SELECT CASE (TRIM(IS_PLANE(L)))
1145                          CASE ('E')
1146                             IJKE = EAST_OF(IJK)
1147                             IPJK = IP_OF(IJK)
1148     !
1149                             D_F = AVG_X_H(DIF(IJK),DIF(IJKE),I)*ODX_E(I)*AYZ(IJK)
1150     !
1151                             A_M(IJK,E,M) = A_M(IJK,E,M) - D_F
1152                             A_M(IPJK,W,M) = A_M(IPJK,W,M) - D_F
1153     !
1154                          CASE ('N')
1155                             IJKN = NORTH_OF(IJK)
1156                             IJPK = JP_OF(IJK)
1157     !
1158                             D_F = AVG_Y_H(DIF(IJK),DIF(IJKN),J)*ODY_N(J)*AXZ(IJK)
1159     !
1160                             A_M(IJK,N,M) = A_M(IJK,N,M) - D_F
1161                             A_M(IJPK,S,M) = A_M(IJPK,S,M) - D_F
1162     !
1163                          CASE ('T')
1164                             IF (DO_K) THEN
1165                                IJKT = TOP_OF(IJK)
1166                                IJKP = KP_OF(IJK)
1167     !
1168                                D_F = AVG_Z_H(DIF(IJK),DIF(IJKT),K)*OX(I)*ODZ_T(K)*&
1169                                   AXY(IJK)
1170     !
1171                                A_M(IJK,T,M) = A_M(IJK,T,M) - D_F
1172                                A_M(IJKP,B,M) = A_M(IJKP,B,M) - D_F
1173     !
1174                             ENDIF
1175                          CASE DEFAULT
1176     !
1177                          END SELECT
1178                       END DO
1179                    END DO
1180                 END DO
1181              ENDIF
1182           END DO
1183     !
1184           RETURN
1185           END SUBROUTINE DIF_PHI_IS
1186