File: N:\mfix\model\cartesian_grid\CG_set_outflow.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Module name: SET_OUTFLOW(BCV, I1, I2, J1, J2, K1, K2)               C
4     !  Purpose: Set specified pressure outflow bc for a specified range of C
5     !           cells                                                      C
6     !                                                                      C
7     !  Author: M. Syamlal                                 Date: 21-JAN-92  C
8     !  Reviewer:M. Syamlal, S. Venkatesan, P. Nicoletti,  Date: 29-JAN-92  C
9     !           W. Rogers                                                  C
10     !                                                                      C
11     !  Revision Number: 1                                                  C
12     !  Purpose: To incorporate Cartesian grid modifications                C
13     !  Author: Jeff Dietiker                              Date: 01-Jul-09  C
14     !                                                                      C
15     !  Literature/Document References:                                     C
16     !                                                                      C
17     !  Variables referenced: MMAX                                          C
18     !                                                                      C
19     !  Variables modified: I, J, K, RO_g, ROP_g,                           C
20     !                      EP_g, C
21     !                      T_g, T_s,  M, ROP_s, U_g, U_s, V_g, V_s,  C
22     !                      W_g, W_s,
23     !                                                                      C
24     !  Local variables: IJK, LFLUID                                      C
25     !                                                                      C
26     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
27     !
28           SUBROUTINE CG_SET_OUTFLOW
29     !...Translated by Pacific-Sierra Research VAST-90 2.06G5  12:17:31  12/09/98
30     !...Switches: -xf
31     !
32     !-----------------------------------------------
33     !   M o d u l e s
34     !-----------------------------------------------
35           USE bc
36           USE compar        !//d
37           USE cutcell
38           USE eos, ONLY: EOSG
39           USE fldvar
40           USE geometry
41           USE indices
42           USE mflux
43           USE param
44           USE param1
45           USE physprop
46           USE quadric
47           USE run
48           USE scalars
49           IMPLICIT NONE
50     !-----------------------------------------------
51     !   G l o b a l   P a r a m e t e r s
52     !-----------------------------------------------
53     !-----------------------------------------------
54     !   D u m m y   A r g u m e n t s
55     !-----------------------------------------------
56     !
57     !                      indices
58           INTEGER          I, J, K, M, NN
59     !
60     !                      Local index for boundary cell
61           INTEGER          IJK
62     !
63     !                      Boundary condition number
64           INTEGER          BCV
65     !
66     !                      Locall index for a fluid cell near the boundary cell
67           INTEGER          LFLUID
68     
69           INTEGER          IJKW,IJKWW,IJKS,IJKSS,IJKB
70     
71           INTEGER :: BCT1,BCT2,BCT3,BCT4
72     
73           LOGICAL :: TEST1,TEST2
74     !-----------------------------------------------
75     !
76     !      print*,'top of cg_set_outflow'
77           DO IJK = IJKSTART3, IJKEND3
78           IF(INTERIOR_CELL_AT(IJK)) THEN
79     
80     
81              I = I_OF(IJK)
82              J = J_OF(IJK)
83              K = K_OF(IJK)
84     
85     !//SP Check if current i,j,k resides on this PE
86              IF (.NOT.IS_ON_myPE_plus2layers(I,J,K)) CYCLE
87     !        IJK = FUNIJK(I,J,K)
88     !
89     ! Fluid cell at West
90     !
91     !      print*,'west'
92               BCV = BC_U_ID(IJK)
93     
94     
95               IF(BCV>0)THEN
96     
97                 IF (BC_TYPE_ENUM(BCV) == CG_PO) THEN
98     
99                    IF (FLUID_AT(IM_OF(IJK))) THEN
100                       LFLUID = IM_OF(IJK)
101     !            print*,'west treatment:IJK,LFLUID=',IJK,LFLUID
102     !            read(*,*)
103     !                  IF (U_G(LFLUID)>=ZERO .OR. EP_G(IJK)==UNDEFINED) THEN
104     !                     IF (BC_TYPE(BCV) /= 'P_OUTFLOW') P_G(IJK) = P_G(LFLUID)
105                          T_G(IJK) = T_G(LFLUID)
106                          NN = 1
107                          IF (NMAX(0) > 0) THEN
108                             X_G(IJK,:NMAX(0)) = X_G(LFLUID,:NMAX(0))
109                             NN = NMAX(0) + 1
110                          ENDIF
111                          MW_MIX_G(IJK) = MW_MIX_G(LFLUID)
112                          IF (RO_G0 == UNDEFINED) RO_G(IJK) = EOSG(MW_MIX_G(IJK),P_G&
113                             (IJK),T_G(IJK))
114     !                  ENDIF
115                       P_STAR(IJK) = P_STAR(LFLUID)
116                       IF (BC_EP_G(BCV) == UNDEFINED) EP_G(IJK) = ONE
117     
118                       DO NN = 1, NScalar
119                          M = Phase4Scalar(NN)
120                          IF(M == 0)Then
121     !                      IF (U_G(LFLUID)>=ZERO) THEN
122                               Scalar(IJK, NN) = Scalar(LFLUID, NN)
123     !                      ENDIF
124                          Else
125     !                      IF (U_s(LFLUID, M)>=ZERO) THEN
126                               Scalar(IJK, NN) = Scalar(LFLUID, NN)
127     !                     ENDIF
128                          Endif
129                       END DO
130     
131                       IF(K_Epsilon) THEN
132     !                    IF (U_G(LFLUID) >= ZERO) THEN
133                          K_Turb_G(IJK) = K_Turb_G(LFLUID)
134                          E_Turb_G(IJK) = E_Turb_G(LFLUID)
135     !                   ENDIF
136                       ENDIF
137     
138     
139                       DO M = 1, MMAX
140                          P_S(IJK,M) = P_S(LFLUID,M)
141     !                     IF (U_S(LFLUID,M) >= ZERO) THEN
142                             ROP_S(IJK,M) = ROP_S(LFLUID,M)
143                             T_S(IJK,M) = T_S(LFLUID,M)
144     !                     ELSE
145     !                        ROP_S(IJK,M) = ZERO
146     !                     ENDIF
147     !
148                          IF(BC_ROP_S(BCV,M)/=UNDEFINED)ROP_S(IJK,M)=BC_ROP_S(BCV,M)
149     !
150                          IF(BC_EP_G(BCV)==UNDEFINED)EP_G(IJK)=EP_G(IJK)-EP_S(IJK,M)
151     !
152                          NN = 1
153                          IF (NMAX(M) > 0) THEN
154                             X_S(IJK,M,:NMAX(M)) = X_S(LFLUID,M,:NMAX(M))
155                             NN = NMAX(M) + 1
156                          ENDIF
157                       END DO
158                       ROP_G(IJK) = RO_G(IJK)*EP_G(IJK)
159     !                  IF (ROP_G(IJK) > ZERO) THEN
160     !                     U_G(IJK) = ROP_G(LFLUID)*U_G(LFLUID)/ROP_G(IJK)
161     !                  ELSE
162     !                     U_G(IJK) = ZERO
163     !                  ENDIF
164     !                  V_G(IJK) = V_G(LFLUID)
165     !                  W_G(IJK) = W_G(LFLUID)
166     !                  Flux_gE(IJK) = Flux_gE(LFLUID)
167     !                 Flux_gN(IJK) = Flux_gN(LFLUID)
168     !                 Flux_gT(IJK) = Flux_gT(LFLUID)
169     
170     !                  IF (MMAX > 0) THEN
171     !                     WHERE (ROP_S(IJK,:MMAX) > ZERO)
172     !                        U_S(IJK,:MMAX) = ROP_S(LFLUID,:MMAX)*U_S(LFLUID,:MMAX)/&
173     !                           ROP_S(IJK,:MMAX)
174     !                     ELSEWHERE
175     !                        U_S(IJK,:MMAX) = ZERO
176     !                    END WHERE
177     !                     V_S(IJK,:MMAX) = V_S(LFLUID,:MMAX)
178     !                     W_S(IJK,:MMAX) = W_S(LFLUID,:MMAX)
179     !                     Flux_sE(IJK,:MMAX) = Flux_sE(LFLUID,:MMAX)
180     !                    Flux_sN(IJK,:MMAX) = Flux_sN(LFLUID,:MMAX)
181     !                    Flux_sT(IJK,:MMAX) = Flux_sT(LFLUID,:MMAX)
182     !                  ENDIF
183                    ENDIF
184     
185     
186                 ENDIF
187              ENDIF
188     
189              IJKW = WEST_OF(IJK)
190              IJKWW = WEST_OF(IJKW)
191     
192              BCT1=BLANK
193              BCT2=BLANK
194              BCT3=BLANK
195              BCT4=BLANK
196              BCV = BC_ID(IJK)
197              IF(BCV>0) BCT1 = BC_TYPE_ENUM(BCV)
198              BCV = BC_U_ID(IJK)
199              IF(BCV>0) BCT2 = BC_TYPE_ENUM(BCV)
200              BCV = BC_U_ID(IJKW)
201              IF(BCV>0) BCT3 = BC_TYPE_ENUM(BCV)
202              BCV = BC_ID(IJKW)
203              IF(BCV>0) BCT4 = BC_TYPE_ENUM(BCV)
204     
205              TEST1= ((BCT1 == CG_PO).AND.(BCT2 /= CG_PO).AND.(BCT3 == CG_PO))
206     
207              TEST2 = (BCT2 == CG_PO).AND.(BCT4 /= CG_PO)
208     
209     !         TEST2 = (FLAG(IJK) == 11)
210     
211     
212              IF(TEST1) THEN
213                 BCV = BC_ID(IJK)
214     
215     !         IF(TEST1.OR.TEST2) THEN
216     
217     !            IF(FLUID_AT(IJKW)) THEN
218     
219     !            ELSEIF(FLUID_AT(IJKWW)) THEN
220     !               IJKW = IJKWW
221     !            ENDIF
222     
223     
224     !         IF((BCT1 == 'CG_PO').AND.(BCT2 /= 'CG_PO').AND.(BCT3 == 'CG_PO')) THEN
225     
226     
227     
228     !         IF(BCT2 == 'CG_PO')
229     
230     
231     !            print*,'IJK,IJKW=',IJK,IJKW
232     !            read(*,*)
233     
234                 T_G(IJK) = T_G(IJKW)
235                 NN = 1
236                 IF (NMAX(0) > 0) THEN
237                    X_G(IJK,:NMAX(0)) = X_G(IJKW,:NMAX(0))
238                    NN = NMAX(0) + 1
239                 ENDIF
240                 MW_MIX_G(IJK) = MW_MIX_G(IJKW)
241                 RO_G(IJK) = RO_G(IJKW)
242                 P_STAR(IJK) = P_STAR(IJKW)
243                 EP_G(IJK) = EP_G(IJKW)
244                 ROP_G(IJK) = ROP_G(IJKW)
245                 DO NN = 1, NScalar
246                    M = Phase4Scalar(NN)
247                    IF(M == 0)Then
248                       Scalar(IJK, NN) = Scalar(IJKW, NN)
249                    ELSE
250                       Scalar(IJK, NN) = Scalar(IJKW, NN)
251                    ENDIF
252                 ENDDO
253                 IF(K_Epsilon) THEN
254                    K_Turb_G(IJK) = K_Turb_G(IJKW)
255                    E_Turb_G(IJK) = E_Turb_G(IJKW)
256                 ENDIF
257     !            DO M = 1, MMAX
258     !               P_S(IJK,M) = P_S(IJKW,M)
259     !               ROP_S(IJK,M) = ROP_S(IJKW,M)
260     !               T_S(IJK,M) = T_S(IJKW,M)
261     !               N = 1
262     !               IF (NMAX(M) > 0) THEN
263     !                 X_S(IJK,M,:NMAX(M)) = X_S(IJKW,M,:NMAX(M))
264     !                 N = NMAX(M) + 1
265     !               ENDIF
266     !            ENDDO
267     
268                 DO M = 1, MMAX
269                    P_S(IJK,M) = P_S(IJKW,M)
270     !               IF (U_S(IJKW,M) >= ZERO) THEN
271                       ROP_S(IJK,M) = ROP_S(IJKW,M)
272                       T_S(IJK,M) = T_S(IJKW,M)
273     !               ELSE
274     !                  ROP_S(IJK,M) = ZERO
275     !               ENDIF
276     !
277                    IF(BC_ROP_S(BCV,M)/=UNDEFINED)ROP_S(IJK,M)=BC_ROP_S(BCV,M)
278     !
279                    IF(BC_EP_G(BCV)==UNDEFINED)EP_G(IJK)=EP_G(IJK)-EP_S(IJK,M)
280     !
281                    NN = 1
282                    IF (NMAX(M) > 0) THEN
283                       X_S(IJK,M,:NMAX(M)) = X_S(IJKW,M,:NMAX(M))
284                       NN = NMAX(M) + 1
285                    ENDIF
286                 END DO
287                 ROP_G(IJK) = RO_G(IJK)*EP_G(IJK)
288     !            IF (ROP_G(IJK) > ZERO) THEN
289     !               U_G(IJK) = ROP_G(IJKW)*U_G(IJKW)/ROP_G(IJK)
290     !            ELSE
291     !               U_G(IJK) = ZERO
292     !            ENDIF
293     !            V_G(IJK) = V_G(IJKW)
294     !            W_G(IJK) = W_G(IJKW)
295     !            Flux_gE(IJK) = Flux_gE(IJKW)
296     !            Flux_gN(IJK) = Flux_gN(IJKW)
297     !            Flux_gT(IJK) = Flux_gT(IJKW)
298     
299     !            IF (MMAX > 0) THEN
300     !               WHERE (ROP_S(IJK,:MMAX) > ZERO)
301     !                  U_S(IJK,:MMAX) = ROP_S(IJKW,:MMAX)*U_S(IJKW,:MMAX)/&
302     !                     ROP_S(IJK,:MMAX)
303     !               ELSEWHERE
304     !                  U_S(IJK,:MMAX) = ZERO
305     !               END WHERE
306     !               V_S(IJK,:MMAX) = V_S(IJKW,:MMAX)
307     !               W_S(IJK,:MMAX) = W_S(IJKW,:MMAX)
308     !               Flux_sE(IJK,:MMAX) = Flux_sE(IJKW,:MMAX)
309     !              Flux_sN(IJK,:MMAX) = Flux_sN(IJKW,:MMAX)
310     !              Flux_sT(IJK,:MMAX) = Flux_sT(IJKW,:MMAX)
311     !            ENDIF
312     
313     
314     
315     !            CYCLE
316     
317              ENDIF
318     
319     
320     
321     
322     !
323     ! Fluid cell at East
324     !      print*,'east'
325     
326               BCV = BC_U_ID(IJK)
327     
328               IF(BCV>0)THEN
329     
330                 IF (BC_TYPE_ENUM(BCV) == CG_PO) THEN
331     
332     
333                    IF (FLUID_AT(IP_OF(IJK))) THEN
334                       LFLUID = IP_OF(IJK)
335     !            print*,'east treatment:IJK,LFLUID=',IJK,LFLUID
336     !            read(*,*)
337     !                  IF (U_G(IJK)<=ZERO .OR. EP_G(IJK)==UNDEFINED) THEN
338     !                     IF (BC_TYPE_ENUM(BCV) /= 'P_OUTFLOW') P_G(IJK) = P_G(LFLUID)
339                          T_G(IJK) = T_G(LFLUID)
340                          NN = 1
341                          IF (NMAX(0) > 0) THEN
342                             X_G(IJK,:NMAX(0)) = X_G(LFLUID,:NMAX(0))
343                             NN = NMAX(0) + 1
344                          ENDIF
345                          MW_MIX_G(IJK) = MW_MIX_G(LFLUID)
346                          IF (RO_G0 == UNDEFINED) RO_G(IJK) = EOSG(MW_MIX_G(IJK),P_G&
347                             (IJK),T_G(IJK))
348     !                  ENDIF
349                       P_STAR(IJK) = P_STAR(LFLUID)
350                       IF (BC_EP_G(BCV) == UNDEFINED) EP_G(IJK) = ONE
351     
352                       DO NN = 1, NScalar
353                          M = Phase4Scalar(NN)
354                          IF(M == 0)Then
355     !                      IF (U_G(LFLUID) <= ZERO) THEN
356                               Scalar(IJK, NN) = Scalar(LFLUID, NN)
357     !                      ENDIF
358                          Else
359     !                      IF (U_s(LFLUID, M) <= ZERO) THEN
360                               Scalar(IJK, NN) = Scalar(LFLUID, NN)
361     !                      ENDIF
362                          Endif
363                       END DO
364     
365                       IF(K_Epsilon) THEN
366     !                    IF (U_G(LFLUID) <= ZERO) THEN
367                             K_Turb_G(IJK) = K_Turb_G(LFLUID)
368                             E_Turb_G(IJK) = E_Turb_G(LFLUID)
369     !                    ENDIF
370                       ENDIF
371     
372                       DO M = 1, MMAX
373                          P_S(IJK,M) = P_S(LFLUID,M)
374     !                     IF (U_S(IJK,M) <= ZERO) THEN
375                             ROP_S(IJK,M) = ROP_S(LFLUID,M)
376                             T_S(IJK,M) = T_S(LFLUID,M)
377     !                     ELSE
378     !                        ROP_S(IJK,M) = ZERO
379     !                     ENDIF
380     !
381                          IF(BC_ROP_S(BCV,M)/=UNDEFINED)ROP_S(IJK,M)=BC_ROP_S(BCV,M)
382     !
383                          IF(BC_EP_G(BCV)==UNDEFINED)EP_G(IJK)=EP_G(IJK)-EP_S(IJK,M)
384                          NN = 1
385                          IF (NMAX(M) > 0) THEN
386                             X_S(IJK,M,:NMAX(M)) = X_S(LFLUID,M,:NMAX(M))
387                             NN = NMAX(M) + 1
388                          ENDIF
389                       END DO
390                       ROP_G(IJK) = RO_G(IJK)*EP_G(IJK)
391     !                  IF (U_G(IJK) == UNDEFINED) THEN
392     !                     IF (ROP_G(IJK) > ZERO) THEN
393     !                        U_G(IJK) = ROP_G(LFLUID)*U_G(LFLUID)/ROP_G(IJK)
394     !                     ELSE
395     !                        U_G(IJK) = ZERO
396     !                     ENDIF
397     !                  ENDIF
398     !                  V_G(IJK) = V_G(LFLUID)
399     !                  W_G(IJK) = W_G(LFLUID)
400     !                  Flux_gE(IJK) = Flux_gE(LFLUID)
401     !                  Flux_gN(IJK) = Flux_gN(LFLUID)
402     !                  Flux_gT(IJK) = Flux_gT(LFLUID)
403     
404     !                  DO M = 1, MMAX
405     !                     IF (U_S(IJK,M) == UNDEFINED) THEN
406     !                        IF (ROP_S(IJK,M) > ZERO) THEN
407     !                           U_S(IJK,M) = ROP_S(LFLUID,M)*U_S(LFLUID,M)/ROP_S(IJK&
408     !                              ,M)
409     !                        ELSE
410     !                           U_S(IJK,M) = ZERO
411     !                        ENDIF
412     !                     ENDIF
413     !                     V_S(IJK,M) = V_S(LFLUID,M)
414     !                     W_S(IJK,M) = W_S(LFLUID,M)
415     !                     Flux_sE(IJK,M) = Flux_sE(LFLUID,M)
416     !                    Flux_sN(IJK,M) = Flux_sN(LFLUID,M)
417     !                    Flux_sT(IJK,M) = Flux_sT(LFLUID,M)
418     !                  END DO
419                    ENDIF
420     
421                 ENDIF
422              ENDIF
423     
424     !
425     ! Fluid cell at South
426     !
427     !      print*,'south'
428               BCV = BC_V_ID(IJK)
429     !      print*,ijk,I,J,K,bcv,BC_TYPE_ENUM(BCV)
430               IF(BCV>0)THEN
431     
432                 IF (BC_TYPE_ENUM(BCV) == CG_PO) THEN
433     
434     
435                    IF (FLUID_AT(JM_OF(IJK))) THEN
436                       LFLUID = JM_OF(IJK)
437     !            print*,'south treatment:IJK,LFLUID=',IJK,LFLUID
438     !            read(*,*)
439     
440     
441     !                  IF (V_G(LFLUID)>=ZERO .OR. EP_G(IJK)==UNDEFINED) THEN
442     !                     IF (BC_TYPE_ENUM(BCV) /= 'P_OUTFLOW') P_G(IJK) = P_G(LFLUID)
443                          T_G(IJK) = T_G(LFLUID)
444                          NN = 1
445                          IF (NMAX(0) > 0) THEN
446                             X_G(IJK,:NMAX(0)) = X_G(LFLUID,:NMAX(0))
447                             NN = NMAX(0) + 1
448                          ENDIF
449                          MW_MIX_G(IJK) = MW_MIX_G(LFLUID)
450                          IF (RO_G0 == UNDEFINED) RO_G(IJK) = EOSG(MW_MIX_G(IJK),P_G&
451                             (IJK),T_G(IJK))
452     !                  ENDIF
453                       P_STAR(IJK) = P_STAR(LFLUID)
454                       IF (BC_EP_G(BCV) == UNDEFINED) EP_G(IJK) = ONE
455     
456                       DO NN = 1, NScalar
457                          M = Phase4Scalar(NN)
458                          IF(M == 0)Then
459     !                      IF (V_G(LFLUID) >= ZERO) THEN
460                               Scalar(IJK, NN) = Scalar(LFLUID, NN)
461     !                      ENDIF
462                          Else
463     !                      IF (V_s(LFLUID, M) >= ZERO) THEN
464                            Scalar(IJK, NN) = Scalar(LFLUID, NN)
465     !                     ENDIF
466                          Endif
467                       END DO
468     
469                       IF(K_Epsilon) THEN
470     !                    IF (V_G(LFLUID) >= ZERO) THEN
471                             K_Turb_G(IJK) = K_Turb_G(LFLUID)
472                             E_Turb_G(IJK) = E_Turb_G(LFLUID)
473     !                   ENDIF
474                       ENDIF
475     
476                       DO M = 1, MMAX
477                          P_S(IJK,M) = P_S(LFLUID,M)
478     !                     IF (V_S(LFLUID,M) >= 0.) THEN
479                             ROP_S(IJK,M) = ROP_S(LFLUID,M)
480                             T_S(IJK,M) = T_S(LFLUID,M)
481     !                     ELSE
482     !                        ROP_S(IJK,M) = ZERO
483     !                     ENDIF
484     !
485                          IF(BC_ROP_S(BCV,M)/=UNDEFINED)ROP_S(IJK,M)=BC_ROP_S(BCV,M)
486     !
487                          IF(BC_EP_G(BCV)==UNDEFINED)EP_G(IJK)=EP_G(IJK)-EP_S(IJK,M)
488                          NN = 1
489                          IF (NMAX(M) > 0) THEN
490                             X_S(IJK,M,:NMAX(M)) = X_S(LFLUID,M,:NMAX(M))
491                             NN = NMAX(M) + 1
492                          ENDIF
493                       END DO
494                       ROP_G(IJK) = RO_G(IJK)*EP_G(IJK)
495     !                  U_G(IJK) = U_G(LFLUID)
496     !                  IF (ROP_G(IJK) > ZERO) THEN
497     !                     V_G(IJK) = ROP_G(LFLUID)*V_G(LFLUID)/ROP_G(IJK)
498     !                  ELSE
499     !                     V_G(IJK) = ZERO
500     !                  ENDIF
501     !                  W_G(IJK) = W_G(LFLUID)
502     !                  Flux_gE(IJK) = Flux_gE(LFLUID)
503     !                  Flux_gN(IJK) = Flux_gN(LFLUID)
504     !                  Flux_gT(IJK) = Flux_gT(LFLUID)
505     
506     !                  IF (MMAX > 0) THEN
507     !                     U_S(IJK,:MMAX) = U_S(LFLUID,:MMAX)
508     !                     V_S(IJK,:MMAX) = V_S(LFLUID,:MMAX)
509     !                     W_S(IJK,:MMAX) = W_S(LFLUID,:MMAX)
510     !                     Flux_sE(IJK,:MMAX) = Flux_sE(LFLUID,:MMAX)
511     !                    Flux_sN(IJK,:MMAX) = Flux_sN(LFLUID,:MMAX)
512     !                    Flux_sT(IJK,:MMAX) = Flux_sT(LFLUID,:MMAX)
513     !                  ENDIF
514                    ENDIF
515     
516                 ENDIF
517              ENDIF
518     
519              IJKS = SOUTH_OF(IJK)
520              IJKSS = SOUTH_OF(IJKS)
521     !         print*,'ijk,ijks=',ijk,ijks,JM_OF(IJK)
522              BCT1=BLANK
523              BCT2=BLANK
524              BCT3=BLANK
525              BCV = BC_ID(IJK)
526     !         print*,'bcv=',bcv
527              IF(BCV>0) BCT1 = BC_TYPE_ENUM(BCV)
528              BCV = BC_V_ID(IJK)
529     !         print*,'bcv=',bcv
530              IF(BCV>0) BCT2 = BC_TYPE_ENUM(BCV)
531              BCV = BC_V_ID(IJKS)
532     !         print*,'bcv=',bcv
533              IF(BCV>0) BCT3 = BC_TYPE_ENUM(BCV)
534     
535     !         IF((BCT1 == 'CG_PO').AND.(BCT2 /= 'CG_PO').AND.(BCT3 == 'CG_PO')) THEN
536     
537     
538     
539              TEST1 = (BCT1 == CG_PO).AND.(BCT2 /= CG_PO).AND.(BCT3 == CG_PO)
540              TEST2 = (FLAG(IJK) == 11).AND.(BCT2 == CG_PO)
541     
542              TEST2 = (FLAG(IJK) == 11)
543     
544              IF(TEST1) THEN
545     
546                 BCV = BC_ID(IJK)
547     
548     !         IF(TEST1.OR.TEST2) THEN
549     
550     !            IF(FLUID_AT(IJKS)) THEN
551     
552     !            ELSEIF(FLUID_AT(IJKSS)) THEN
553     !               IJKS = IJKSS
554     !            ENDIF
555     
556     
557     !            print*,'IJK,IJKS=',IJK,IJKS
558     !            read(*,*)
559     
560                 T_G(IJK) = T_G(IJKS)
561                 NN = 1
562                 IF (NMAX(0) > 0) THEN
563                    X_G(IJK,:NMAX(0)) = X_G(IJKS,:NMAX(0))
564                    NN = NMAX(0) + 1
565                 ENDIF
566                 MW_MIX_G(IJK) = MW_MIX_G(IJKS)
567                 RO_G(IJK) = RO_G(IJKS)
568                 P_STAR(IJK) = P_STAR(IJKS)
569                 EP_G(IJK) = EP_G(IJKS)
570                 ROP_G(IJK) = ROP_G(IJKS)
571                 DO NN = 1, NScalar
572                    M = Phase4Scalar(NN)
573                    IF(M == 0)Then
574                       Scalar(IJK, NN) = Scalar(IJKS, NN)
575                    ELSE
576                       Scalar(IJK, NN) = Scalar(IJKS, NN)
577                    ENDIF
578                 ENDDO
579                 IF(K_Epsilon) THEN
580                    K_Turb_G(IJK) = K_Turb_G(IJKS)
581                    E_Turb_G(IJK) = E_Turb_G(IJKS)
582                 ENDIF
583     !            DO M = 1, MMAX
584     !               P_S(IJK,M) = P_S(IJKS,M)
585     !               ROP_S(IJK,M) = ROP_S(IJKS,M)
586     !               T_S(IJK,M) = T_S(IJKS,M)
587     !               N = 1
588     !               IF (NMAX(M) > 0) THEN
589     !                 X_S(IJK,M,:NMAX(M)) = X_S(IJKS,M,:NMAX(M))
590     !                 N = NMAX(M) + 1
591     !               ENDIF
592     !            ENDDO
593     
594                 DO M = 1, MMAX
595                    P_S(IJK,M) = P_S(IJKS,M)
596     !               IF (V_S(IJKS,M) >= 0.) THEN
597                       ROP_S(IJK,M) = ROP_S(IJKS,M)
598                       T_S(IJK,M) = T_S(IJKS,M)
599     !               ELSE
600     !                  ROP_S(IJK,M) = ZERO
601     !               ENDIF
602     !
603                    IF(BC_ROP_S(BCV,M)/=UNDEFINED)ROP_S(IJK,M)=BC_ROP_S(BCV,M)
604     !
605                    IF(BC_EP_G(BCV)==UNDEFINED)EP_G(IJK)=EP_G(IJK)-EP_S(IJK,M)
606                    NN = 1
607                    IF (NMAX(M) > 0) THEN
608                       X_S(IJK,M,:NMAX(M)) = X_S(IJKS,M,:NMAX(M))
609                       NN = NMAX(M) + 1
610                    ENDIF
611                 END DO
612                 ROP_G(IJK) = RO_G(IJK)*EP_G(IJK)
613     !            U_G(IJK) = U_G(IJKS)
614     !            IF (ROP_G(IJK) > ZERO) THEN
615     !               V_G(IJK) = ROP_G(IJKS)*V_G(IJKS)/ROP_G(IJK)
616     !            ELSE
617     !               V_G(IJK) = ZERO
618     !            ENDIF
619     !            W_G(IJK) = W_G(IJKS)
620     !            Flux_gE(IJK) = Flux_gE(IJKS)
621     !           Flux_gN(IJK) = Flux_gN(IJKS)
622     !           Flux_gT(IJK) = Flux_gT(IJKS)
623     
624     !            IF (MMAX > 0) THEN
625     !               U_S(IJK,:MMAX) = U_S(IJKS,:MMAX)
626     !               V_S(IJK,:MMAX) = V_S(IJKS,:MMAX)
627     !               W_S(IJK,:MMAX) = W_S(IJKS,:MMAX)
628     !               Flux_sE(IJK,:MMAX) = Flux_sE(IJKS,:MMAX)
629     !              Flux_sN(IJK,:MMAX) = Flux_sN(IJKS,:MMAX)
630     !              Flux_sT(IJK,:MMAX) = Flux_sT(IJKS,:MMAX)
631     !            ENDIF
632     
633     
634     
635     
636     !            CYCLE
637     
638              ENDIF
639     
640     
641     
642     !
643     ! Fluid cell at North
644     !
645     
646     !      print*,'north'
647               BCV = BC_V_ID(IJK)
648     
649               IF(BCV>0)THEN
650     
651                 IF (BC_TYPE_ENUM(BCV) == CG_PO) THEN
652     
653     
654                    IF (FLUID_AT(JP_OF(IJK))) THEN
655                       LFLUID = JP_OF(IJK)
656     
657     !            print*,'north treatment:IJK,LFLUID=',IJK,LFLUID
658     !            read(*,*)
659     
660     !                  IF (V_G(IJK)<=ZERO .OR. EP_G(IJK)==UNDEFINED) THEN
661     !                     IF (BC_TYPE_ENUM(BCV) /= 'P_OUTFLOW') P_G(IJK) = P_G(LFLUID)
662                          T_G(IJK) = T_G(LFLUID)
663                          NN = 1
664                          IF (NMAX(0) > 0) THEN
665                             X_G(IJK,:NMAX(0)) = X_G(LFLUID,:NMAX(0))
666                             NN = NMAX(0) + 1
667                          ENDIF
668                          MW_MIX_G(IJK) = MW_MIX_G(LFLUID)
669                          IF (RO_G0 == UNDEFINED) RO_G(IJK) = EOSG(MW_MIX_G(IJK),P_G&
670                             (IJK),T_G(IJK))
671     !                  ENDIF
672                       P_STAR(IJK) = P_STAR(LFLUID)
673                       IF (BC_EP_G(BCV) == UNDEFINED) EP_G(IJK) = ONE
674     
675                       DO NN = 1, NScalar
676                          M = Phase4Scalar(NN)
677                          IF(M == 0)Then
678     !                      IF (V_G(LFLUID) <= ZERO) THEN
679                              Scalar(IJK, NN) = Scalar(LFLUID, NN)
680     !                     ENDIF
681                          Else
682     !                      IF (V_s(LFLUID, M) <= ZERO) THEN
683                              Scalar(IJK, NN) = Scalar(LFLUID, NN)
684     !                     ENDIF
685                          Endif
686                       END DO
687     
688                       IF(K_Epsilon) THEN
689     !                    IF (V_G(LFLUID) <= ZERO) THEN
690                            K_Turb_G(IJK) = K_Turb_G(LFLUID)
691                            E_Turb_G(IJK) = E_Turb_G(LFLUID)
692     !                   ENDIF
693                       ENDIF
694     
695                       DO M = 1, MMAX
696                          P_S(IJK,M) = P_S(LFLUID,M)
697     !                     IF (V_S(IJK,M) <= ZERO) THEN
698                             ROP_S(IJK,M) = ROP_S(LFLUID,M)
699                             T_S(IJK,M) = T_S(LFLUID,M)
700     !                     ELSE
701     !                        ROP_S(IJK,M) = ZERO
702     !                     ENDIF
703     !
704                          IF(BC_ROP_S(BCV,M)/=UNDEFINED)ROP_S(IJK,M)=BC_ROP_S(BCV,M)
705     !
706                          IF(BC_EP_G(BCV)==UNDEFINED)EP_G(IJK)=EP_G(IJK)-EP_S(IJK,M)
707                          NN = 1
708                          IF (NMAX(M) > 0) THEN
709                             X_S(IJK,M,:NMAX(M)) = X_S(LFLUID,M,:NMAX(M))
710                             NN = NMAX(M) + 1
711                          ENDIF
712                       END DO
713                       ROP_G(IJK) = RO_G(IJK)*EP_G(IJK)
714     !                  U_G(IJK) = U_G(LFLUID)
715     !                  IF (V_G(IJK) == UNDEFINED) THEN
716     !                     IF (ROP_G(IJK) > ZERO) THEN
717     !                        V_G(IJK) = ROP_G(LFLUID)*V_G(LFLUID)/ROP_G(IJK)
718     !                     ELSE
719     !                        V_G(IJK) = ZERO
720     !                     ENDIF
721     !                  ENDIF
722     !                  W_G(IJK) = W_G(LFLUID)
723     !                  Flux_gE(IJK) = Flux_gE(LFLUID)
724     !                 Flux_gN(IJK) = Flux_gN(LFLUID)
725     !                 Flux_gT(IJK) = Flux_gT(LFLUID)
726     
727     !                  IF (MMAX > 0) THEN
728     !                     U_S(IJK,:MMAX) = U_S(LFLUID,:MMAX)
729     !                     WHERE (V_S(IJK,:MMAX) == UNDEFINED) V_S(IJK,:MMAX) = V_S(&
730     !                        LFLUID,:MMAX)
731     !                     W_S(IJK,:MMAX) = W_S(LFLUID,:MMAX)
732     !                     Flux_sE(IJK,:MMAX) = Flux_sE(LFLUID,:MMAX)
733     !                    Flux_sN(IJK,:MMAX) = Flux_sN(LFLUID,:MMAX)
734     !                    Flux_sT(IJK,:MMAX) = Flux_sT(LFLUID,:MMAX)
735     !                  ENDIF
736                    ENDIF
737     
738                 ENDIF
739              ENDIF
740     
741     !
742     ! Fluid cell at Bottom
743     !
744     !      print*,'bottom'
745     
746               BCV = BC_W_ID(IJK)
747     
748               IF(BCV>0)THEN
749     
750                 IF (BC_TYPE_ENUM(BCV) == CG_PO) THEN
751     
752                    IF (FLUID_AT(KM_OF(IJK))) THEN
753                       LFLUID = KM_OF(IJK)
754     !                  IF (W_G(LFLUID)>=ZERO .OR. EP_G(IJK)==UNDEFINED) THEN
755     !                     IF (BC_TYPE_ENUM(BCV) /= 'P_OUTFLOW') P_G(IJK) = P_G(LFLUID)
756                          T_G(IJK) = T_G(LFLUID)
757                          NN = 1
758                          IF (NMAX(0) > 0) THEN
759                             X_G(IJK,:NMAX(0)) = X_G(LFLUID,:NMAX(0))
760                             NN = NMAX(0) + 1
761                          ENDIF
762                          MW_MIX_G(IJK) = MW_MIX_G(LFLUID)
763                          IF (RO_G0 == UNDEFINED) RO_G(IJK) = EOSG(MW_MIX_G(IJK),P_G&
764                             (IJK),T_G(IJK))
765     !                  ENDIF
766                       P_STAR(IJK) = P_STAR(LFLUID)
767                       IF (BC_EP_G(BCV) == UNDEFINED) EP_G(IJK) = ONE
768     
769                       DO NN = 1, NScalar
770                          M = Phase4Scalar(NN)
771                          IF(M == 0)Then
772     !                      IF (W_G(LFLUID) >= ZERO) THEN
773                              Scalar(IJK, NN) = Scalar(LFLUID, NN)
774     !                      ENDIF
775                          Else
776     !                      IF (W_s(LFLUID, M) >= ZERO) THEN
777                              Scalar(IJK, NN) = Scalar(LFLUID, NN)
778     !                     ENDIF
779                          Endif
780                       END DO
781     
782                       IF(K_Epsilon) THEN
783     !                    IF (W_G(LFLUID) >= ZERO) THEN
784                            K_Turb_G(IJK) = K_Turb_G(LFLUID)
785                            E_Turb_G(IJK) = E_Turb_G(LFLUID)
786     !                   ENDIF
787                       ENDIF
788     
789                       DO M = 1, MMAX
790                          P_S(IJK,M) = P_S(LFLUID,M)
791     !                     IF (W_S(LFLUID,M) >= 0.) THEN
792                             ROP_S(IJK,M) = ROP_S(LFLUID,M)
793                             T_S(IJK,M) = T_S(LFLUID,M)
794     !                     ELSE
795                             ROP_S(IJK,M) = ZERO
796     !                     ENDIF
797     !
798                          IF(BC_ROP_S(BCV,M)/=UNDEFINED)ROP_S(IJK,M)=BC_ROP_S(BCV,M)
799     !
800                          IF(BC_EP_G(BCV)==UNDEFINED)EP_G(IJK)=EP_G(IJK)-EP_S(IJK,M)
801                          NN = 1
802                          IF (NMAX(M) > 0) THEN
803                             X_S(IJK,M,:NMAX(M)) = X_S(LFLUID,M,:NMAX(M))
804                             NN = NMAX(M) + 1
805                          ENDIF
806                       END DO
807                       ROP_G(IJK) = RO_G(IJK)*EP_G(IJK)
808     !                  U_G(IJK) = U_G(LFLUID)
809     !                  V_G(IJK) = V_G(LFLUID)
810     !                  IF (ROP_G(IJK) > ZERO) THEN
811     !                     W_G(IJK) = ROP_G(LFLUID)*W_G(LFLUID)/ROP_G(IJK)
812     !                  ELSE
813     !                     W_G(IJK) = ZERO
814     !                  ENDIF
815     !                  Flux_gE(IJK) = Flux_gE(LFLUID)
816     !                 Flux_gN(IJK) = Flux_gN(LFLUID)
817     !                 Flux_gT(IJK) = Flux_gT(LFLUID)
818     
819     !                  IF (MMAX > 0) THEN
820     !                     U_S(IJK,:MMAX) = U_S(LFLUID,:MMAX)
821     !                     V_S(IJK,:MMAX) = V_S(LFLUID,:MMAX)
822     !                     W_S(IJK,:MMAX) = W_S(LFLUID,:MMAX)
823     !                     Flux_sE(IJK,:MMAX) = Flux_sE(LFLUID,:MMAX)
824     !                    Flux_sN(IJK,:MMAX) = Flux_sN(LFLUID,:MMAX)
825     !                    Flux_sT(IJK,:MMAX) = Flux_sT(LFLUID,:MMAX)
826     !                  ENDIF
827                    ENDIF
828     
829     
830                 ENDIF
831              ENDIF
832     
833              IJKB = BOTTOM_OF(IJK)
834              BCT1=BLANK
835              BCT2=BLANK
836              BCT3=BLANK
837              BCV = BC_ID(IJK)
838              IF(BCV>0) BCT1 = BC_TYPE_ENUM(BCV)
839              BCV = BC_W_ID(IJK)
840              IF(BCV>0) BCT2 = BC_TYPE_ENUM(BCV)
841              BCV = BC_W_ID(IJKB)
842              IF(BCV>0) BCT3 = BC_TYPE_ENUM(BCV)
843     
844              IF((BCT1 == CG_PO).AND.(BCT2 /= CG_PO).AND.(BCT3 == CG_PO)) THEN
845     !            print*,'IJK,IJKB=',IJK,IJKB
846     !            read(*,*)
847     
848                 BCV = BC_ID(IJK)
849     
850                 T_G(IJK) = T_G(IJKB)
851                 NN = 1
852                 IF (NMAX(0) > 0) THEN
853                    X_G(IJK,:NMAX(0)) = X_G(IJKB,:NMAX(0))
854                    NN = NMAX(0) + 1
855                 ENDIF
856                 MW_MIX_G(IJK) = MW_MIX_G(IJKB)
857                 RO_G(IJK) = RO_G(IJKB)
858                 P_STAR(IJK) = P_STAR(IJKB)
859                 EP_G(IJK) = EP_G(IJKB)
860                 ROP_G(IJK) = ROP_G(IJKB)
861                 DO NN = 1, NScalar
862                    M = Phase4Scalar(NN)
863                    IF(M == 0)Then
864                       Scalar(IJK, NN) = Scalar(IJKB, NN)
865                    ELSE
866                       Scalar(IJK, NN) = Scalar(IJKB, NN)
867                    ENDIF
868                 ENDDO
869                 IF(K_Epsilon) THEN
870                    K_Turb_G(IJK) = K_Turb_G(IJKB)
871                    E_Turb_G(IJK) = E_Turb_G(IJKB)
872                 ENDIF
873     !            DO M = 1, MMAX
874     !               P_S(IJK,M) = P_S(IJKB,M)
875     !               ROP_S(IJK,M) = ROP_S(IJKB,M)
876     !               T_S(IJK,M) = T_S(IJKB,M)
877     !               N = 1
878     !               IF (NMAX(M) > 0) THEN
879     !                 X_S(IJK,M,:NMAX(M)) = X_S(IJKB,M,:NMAX(M))
880     !                 N = NMAX(M) + 1
881     !               ENDIF
882     !            ENDDO
883     
884                 DO M = 1, MMAX
885                    P_S(IJK,M) = P_S(IJKB,M)
886     !               IF (W_S(IJKB,M) >= 0.) THEN
887                       ROP_S(IJK,M) = ROP_S(IJKB,M)
888                       T_S(IJK,M) = T_S(IJKB,M)
889     !               ELSE
890     !                  ROP_S(IJK,M) = ZERO
891     !               ENDIF
892     !
893                    IF(BC_ROP_S(BCV,M)/=UNDEFINED)ROP_S(IJK,M)=BC_ROP_S(BCV,M)
894     !
895                    IF(BC_EP_G(BCV)==UNDEFINED)EP_G(IJK)=EP_G(IJK)-EP_S(IJK,M)
896                    NN = 1
897                    IF (NMAX(M) > 0) THEN
898                       X_S(IJK,M,:NMAX(M)) = X_S(IJKB,M,:NMAX(M))
899                       NN = NMAX(M) + 1
900                    ENDIF
901                 END DO
902                 ROP_G(IJK) = RO_G(IJK)*EP_G(IJK)
903     !            U_G(IJK) = U_G(IJKB)
904     !            V_G(IJK) = V_G(IJKB)
905     !            IF (ROP_G(IJK) > ZERO) THEN
906     !               W_G(IJK) = ROP_G(IJKB)*W_G(IJKB)/ROP_G(IJK)
907     !            ELSE
908     !               W_G(IJK) = ZERO
909     !            ENDIF
910     !            Flux_gE(IJK) = Flux_gE(IJKB)
911     !           Flux_gN(IJK) = Flux_gN(IJKB)
912     !           Flux_gT(IJK) = Flux_gT(IJKB)
913     
914     !            IF (MMAX > 0) THEN
915     !               U_S(IJK,:MMAX) = U_S(IJKB,:MMAX)
916     !               V_S(IJK,:MMAX) = V_S(IJKB,:MMAX)
917     !               W_S(IJK,:MMAX) = W_S(IJKB,:MMAX)
918     !               Flux_sE(IJK,:MMAX) = Flux_sE(IJKB,:MMAX)
919     !              Flux_sN(IJK,:MMAX) = Flux_sN(IJKB,:MMAX)
920     !              Flux_sT(IJK,:MMAX) = Flux_sT(IJKB,:MMAX)
921     !            ENDIF
922     
923     !            CYCLE
924     
925              ENDIF
926     
927     
928     !
929     ! Fluid cell at Top
930     !
931     !      print*,'top'
932               BCV = BC_W_ID(IJK)
933     
934               IF(BCV>0)THEN
935     
936                 IF (BC_TYPE_ENUM(BCV) == CG_PO) THEN
937     
938     
939                    IF (FLUID_AT(KP_OF(IJK))) THEN
940                       LFLUID = KP_OF(IJK)
941     !                  IF (W_G(IJK)<=ZERO .OR. EP_G(IJK)==UNDEFINED) THEN
942     !                     IF (BC_TYPE_ENUM(BCV) /= 'P_OUTFLOW') P_G(IJK) = P_G(LFLUID)
943                          T_G(IJK) = T_G(LFLUID)
944                          NN = 1
945                          IF (NMAX(0) > 0) THEN
946                             X_G(IJK,:NMAX(0)) = X_G(LFLUID,:NMAX(0))
947                             NN = NMAX(0) + 1
948                          ENDIF
949                          MW_MIX_G(IJK) = MW_MIX_G(LFLUID)
950                          IF (RO_G0 == UNDEFINED) RO_G(IJK) = EOSG(MW_MIX_G(IJK),P_G&
951                             (IJK),T_G(IJK))
952     !                  ENDIF
953                       P_STAR(IJK) = P_STAR(LFLUID)
954                       IF (BC_EP_G(BCV) == UNDEFINED) EP_G(IJK) = ONE
955     
956                       DO NN = 1, NScalar
957                          M = Phase4Scalar(NN)
958                          IF(M == 0)Then
959     !                      IF (W_G(LFLUID) <= ZERO) THEN
960                              Scalar(IJK, NN) = Scalar(LFLUID, NN)
961     !                     ENDIF
962                          Else
963     !                      IF (W_s(LFLUID, M) <= ZERO) THEN
964                              Scalar(IJK, NN) = Scalar(LFLUID, NN)
965     !                      ENDIF
966                          Endif
967                       END DO
968     
969                       IF(K_Epsilon) THEN
970     !                    IF (W_G(LFLUID) <= ZERO) THEN
971                            K_Turb_G(IJK) = K_Turb_G(LFLUID)
972                            E_Turb_G(IJK) = E_Turb_G(LFLUID)
973     !                   ENDIF
974                       ENDIF
975     
976                       DO M = 1, MMAX
977                          P_S(IJK,M) = P_S(LFLUID,M)
978     !                     IF (W_S(IJK,M) <= ZERO) THEN
979                             ROP_S(IJK,M) = ROP_S(LFLUID,M)
980                             T_S(IJK,M) = T_S(LFLUID,M)
981     !                     ELSE
982     !                        ROP_S(IJK,M) = ZERO
983     !                     ENDIF
984     !
985                          IF(BC_ROP_S(BCV,M)/=UNDEFINED)ROP_S(IJK,M)=BC_ROP_S(BCV,M)
986     !
987                          IF(BC_EP_G(BCV)==UNDEFINED)EP_G(IJK)=EP_G(IJK)-EP_S(IJK,M)
988                          NN = 1
989                          IF (NMAX(M) > 0) THEN
990                             X_S(IJK,M,:NMAX(M)) = X_S(LFLUID,M,:NMAX(M))
991                             NN = NMAX(M) + 1
992                          ENDIF
993                       END DO
994                       ROP_G(IJK) = RO_G(IJK)*EP_G(IJK)
995     !                  U_G(IJK) = U_G(LFLUID)
996     !                  V_G(IJK) = V_G(LFLUID)
997     !                  IF (W_G(IJK) == UNDEFINED) THEN
998     !                     IF (ROP_G(IJK) > ZERO) THEN
999     !                        W_G(IJK) = ROP_G(LFLUID)*W_G(LFLUID)/ROP_G(IJK)
1000     !                     ELSE
1001     !                        W_G(IJK) = ZERO
1002     !                     ENDIF
1003     !                  ENDIF
1004     !                  Flux_gE(IJK) = Flux_gE(LFLUID)
1005     !                 Flux_gN(IJK) = Flux_gN(LFLUID)
1006     !                 Flux_gT(IJK) = Flux_gT(LFLUID)
1007     
1008     !                  IF (MMAX > 0) THEN
1009     !                     U_S(IJK,:MMAX) = U_S(LFLUID,:MMAX)
1010     !                     V_S(IJK,:MMAX) = V_S(LFLUID,:MMAX)
1011     !                     WHERE (W_S(IJK,:MMAX) == UNDEFINED) W_S(IJK,:MMAX) = W_S(&
1012     !                        LFLUID,:MMAX)
1013     !                     Flux_sE(IJK,:MMAX) = Flux_sE(LFLUID,:MMAX)
1014     !                    Flux_sN(IJK,:MMAX) = Flux_sN(LFLUID,:MMAX)
1015     !                    Flux_sT(IJK,:MMAX) = Flux_sT(LFLUID,:MMAX)
1016     !                  ENDIF
1017                    ENDIF
1018     
1019                 ENDIF
1020              ENDIF
1021     
1022           ENDIF
1023           END DO
1024     
1025     !      print*,'bottom of cg_set_outflow'
1026           RETURN
1027     
1028         CONTAINS
1029     
1030           INCLUDE 'functions.inc'
1031     
1032           END SUBROUTINE CG_SET_OUTFLOW
1033