File: RELATIVE:/../../../mfix.git/model/bc_phi.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Subroutine: BC_phi                                                  C
4     !  Purpose: Set up the phi boundary conditions                         C
5     !                                                                      C
6     !  Author: M. Syamlal                                 Date: 30-APR-97  C
7     !                                                                      C
8     !                                                                      C
9     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
10     
11           SUBROUTINE BC_PHI(VAR, BC_PHIF, BC_PHIW, BC_HW_PHI, &
12                             BC_C_PHI, M, A_M, B_M)
13     
14     
15     ! Modules
16     !--------------------------------------------------------------------//
17           USE param
18           USE param1
19           USE matrix
20           USE geometry
21           USE indices
22           USE bc
23           USE compar
24           USE cutcell, only : CARTESIAN_GRID, CG_SAFE_MODE
25           USE fun_avg
26           USE functions
27           IMPLICIT NONE
28     
29     ! Dummy arguments
30     !--------------------------------------------------------------------//
31     ! The field variable being solved for:
32     !     e.g., T_g, T_s, x_g, x_s, Theta_m, scalar, K_Turb_G,
33     !     e_Turb_G
34           DOUBLE PRECISION, INTENT(IN) :: VAR(DIMENSION_3)
35     ! Boundary conditions specifications
36     ! bc_phif = flow boundary value
37     ! bc_phiw = wall boundary value
38     ! bc_hw_phi = transfer coefficient
39     !      = 0 value means specified flux (neumann type)
40     !      = undefined value means specified wall value (dirichlet type)
41     !      = other value means mixed type
42     ! bc_C_phi = transfer flux
43           DOUBLE PRECISION, INTENT(IN) :: BC_phif(DIMENSION_BC), &
44                                           BC_Phiw(DIMENSION_BC), &
45                                           BC_hw_Phi(DIMENSION_BC), &
46                                           BC_C_Phi(DIMENSION_BC)
47     ! Phase index
48           INTEGER, INTENT(IN) :: M
49     ! Septadiagonal matrix A_m
50           DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
51     ! Vector b_m
52           DOUBLE PRECISION, INTENT(INOUT) :: B_m(DIMENSION_3, 0:DIMENSION_M)
53     
54     ! Local variables
55     !--------------------------------------------------------------------//
56     ! Boundary condition index
57           INTEGER :: L
58     ! Indices
59           INTEGER :: I, J, K, I1, I2, J1, J2, K1, K2, IJK, &
60                      IM, JM, KM
61     !--------------------------------------------------------------------//
62     
63     ! Set up the default walls (i.e., bc_type='dummy' or undefined/default
64     ! boundaries) as non-conducting...
65     ! ---------------------------------------------------------------->>>
66           IF(.NOT.CARTESIAN_GRID) THEN
67     ! when setting up default walls do not use cutcells to avoid conflict
68     
69           IF (DO_K) THEN
70     ! bottom xy plane
71              K1 = 1
72     !!$omp    parallel do private(IJK, J1, I1)
73              DO J1 = jmin3, jmax3
74                 DO I1 = imin3, imax3
75                    IF (.NOT.IS_ON_myPE_plus2layers(I1,J1,K1)) CYCLE
76                    IF (DEAD_CELL_AT(I1,J1,K1)) CYCLE  ! skip dead cells
77                    IJK = FUNIJK(I1,J1,K1)
78     
79                    IF (DEFAULT_WALL_AT(IJK)) THEN
80     ! Cutting the neighbor link between fluid cell and wall cell
81                       A_M(KP_OF(IJK),B,M) = ZERO
82     ! Setting the wall value equal to the adjacent fluid cell value (set
83     ! the boundary cell value equal to adjacent fluid cell value)
84                       A_M(IJK,E,M) = ZERO
85                       A_M(IJK,W,M) = ZERO
86                       A_M(IJK,N,M) = ZERO
87                       A_M(IJK,S,M) = ZERO
88                       A_M(IJK,T,M) = ONE
89                       A_M(IJK,B,M) = ZERO
90                       A_M(IJK,0,M) = -ONE
91                       B_M(IJK,M) = ZERO
92                    ENDIF
93                 ENDDO
94              ENDDO
95     
96     ! top xy plane
97              K1 = KMAX2
98     !!$omp    parallel do private(IJK, J1, I1)
99              DO J1 = jmin3, jmax3
100                 DO I1 = imin3, imax3
101                    IF (.NOT.IS_ON_myPE_plus2layers(I1,J1,K1)) CYCLE
102                    IF (DEAD_CELL_AT(I1,J1,K1)) CYCLE  ! skip dead cells
103                    IJK = FUNIJK(I1,J1,K1)
104                    IF (DEFAULT_WALL_AT(IJK)) THEN
105     ! Cutting the neighbor link between fluid cell and wall cell
106                       A_M(KM_OF(IJK),T,M) = ZERO
107     ! Setting the wall value equal to the adjacent fluid cell value
108                       A_M(IJK,E,M) = ZERO
109                       A_M(IJK,W,M) = ZERO
110                       A_M(IJK,N,M) = ZERO
111                       A_M(IJK,S,M) = ZERO
112                       A_M(IJK,T,M) = ZERO
113                       A_M(IJK,B,M) = ONE
114                       A_M(IJK,0,M) = -ONE
115                       B_M(IJK,M) = ZERO
116                    ENDIF
117                 ENDDO
118              ENDDO
119           ENDIF
120     
121     ! south xz plane
122           J1 = 1
123     !!$omp    parallel do private(IJK, K1, I1)
124           DO K1 = kmin3, kmax3
125              DO I1 = imin3, imax3
126                 IF (.NOT.IS_ON_myPE_plus2layers(I1,J1,K1)) CYCLE
127                 IF (DEAD_CELL_AT(I1,J1,K1)) CYCLE  ! skip dead cells
128                 IJK = FUNIJK(I1,J1,K1)
129                 IF (DEFAULT_WALL_AT(IJK)) THEN
130     ! Cutting the neighbor link between fluid cell and wall cell
131                    A_M(JP_OF(IJK),S,M) = ZERO
132     ! Setting the wall value equal to the adjacent fluid cell value
133                    A_M(IJK,E,M) = ZERO
134                    A_M(IJK,W,M) = ZERO
135                    A_M(IJK,N,M) = ONE
136                    A_M(IJK,S,M) = ZERO
137                    A_M(IJK,T,M) = ZERO
138                    A_M(IJK,B,M) = ZERO
139                    A_M(IJK,0,M) = -ONE
140                    B_M(IJK,M) = ZERO
141                 ENDIF
142              ENDDO
143           ENDDO
144     
145     ! north xz plane
146           J1 = JMAX2
147     !!$omp    parallel do private(IJK, K1, I1)
148           DO K1 = kmin3, kmax3
149              DO I1 = imin3, imax3
150                 IF (.NOT.IS_ON_myPE_plus2layers(I1,J1,K1)) CYCLE
151                 IF (DEAD_CELL_AT(I1,J1,K1)) CYCLE  ! skip dead cells
152                 IJK = FUNIJK(I1,J1,K1)
153                 IF (DEFAULT_WALL_AT(IJK)) THEN
154     ! Cutting the neighbor link between fluid cell and wall cell
155                    A_M(JM_OF(IJK),N,M) = ZERO
156     ! Setting the wall value equal to the adjacent fluid cell value
157                    A_M(IJK,E,M) = ZERO
158                    A_M(IJK,W,M) = ZERO
159                    A_M(IJK,N,M) = ZERO
160                    A_M(IJK,S,M) = ONE
161                    A_M(IJK,T,M) = ZERO
162                    A_M(IJK,B,M) = ZERO
163                    A_M(IJK,0,M) = -ONE
164                    B_M(IJK,M) = ZERO
165                 ENDIF
166              ENDDO
167           ENDDO
168     
169     ! west yz plane
170           I1 = imin2
171     !!$omp    parallel do private(IJK, K1, J1)
172           DO K1 = kmin3, kmax3
173              DO J1 = jmin3, jmax3
174                 IF (.NOT.IS_ON_myPE_plus2layers(I1,J1,K1)) CYCLE
175                 IF (DEAD_CELL_AT(I1,J1,K1)) CYCLE  ! skip dead cells
176                 IJK = FUNIJK(I1,J1,K1)
177                 IF (DEFAULT_WALL_AT(IJK)) THEN
178     ! Cutting the neighbor link between fluid cell and wall cell
179                    A_M(IP_OF(IJK),W,M) = ZERO
180     ! Setting the wall value equal to the adjacent fluid cell value
181                    A_M(IJK,E,M) = ONE
182                    A_M(IJK,W,M) = ZERO
183                    A_M(IJK,N,M) = ZERO
184                    A_M(IJK,S,M) = ZERO
185                    A_M(IJK,T,M) = ZERO
186                    A_M(IJK,B,M) = ZERO
187                    A_M(IJK,0,M) = -ONE
188                    B_M(IJK,M) = ZERO
189                 ENDIF
190              ENDDO
191           ENDDO
192     
193     ! east yz plane
194           I1 = IMAX2
195     !!$omp    parallel do private(IJK, K1, J1)
196           DO K1 = kmin3, kmax3
197              DO J1 = jmin3, jmax3
198                 IF (.NOT.IS_ON_myPE_plus2layers(I1,J1,K1)) CYCLE
199                 IF (DEAD_CELL_AT(I1,J1,K1)) CYCLE  ! skip dead cells
200                 IJK = FUNIJK(I1,J1,K1)
201                 IF (DEFAULT_WALL_AT(IJK)) THEN
202     ! Cutting the neighbor link between fluid cell and wall cell
203                    A_M(IM_OF(IJK),E,M) = ZERO
204     ! Setting the wall value equal to the adjacent fluid cell value
205                    A_M(IJK,E,M) = ZERO
206                    A_M(IJK,W,M) = ONE
207                    A_M(IJK,N,M) = ZERO
208                    A_M(IJK,S,M) = ZERO
209                    A_M(IJK,T,M) = ZERO
210                    A_M(IJK,B,M) = ZERO
211                    A_M(IJK,0,M) = -ONE
212                    B_M(IJK,M) = ZERO
213                 ENDIF
214              ENDDO
215           ENDDO
216     
217           ENDIF   !(.NOT.CARTESIAN_GRID)
218     
219     ! End setting the default boundary conditions
220     ! ----------------------------------------------------------------<<<
221     
222     
223     ! Set user defined wall boundary conditions
224     ! ---------------------------------------------------------------->>>
225           DO L = 1, DIMENSION_BC
226              IF (BC_DEFINED(L)) THEN
227                 IF (BC_TYPE(L)=='NO_SLIP_WALL' .OR. &
228                     BC_TYPE(L)=='FREE_SLIP_WALL' .OR. &
229                     BC_TYPE(L)=='PAR_SLIP_WALL') THEN
230                    I1 = BC_I_W(L)
231                    I2 = BC_I_E(L)
232                    J1 = BC_J_S(L)
233                    J2 = BC_J_N(L)
234                    K1 = BC_K_B(L)
235                    K2 = BC_K_T(L)
236     !!$omp    parallel do private(IJK, K, J, I, IM, JM, KM)
237                    DO K = K1, K2
238                       DO J = J1, J2
239                          DO I = I1, I2
240                             IF (.NOT.IS_ON_myPE_plus2layers(I,J,K)) CYCLE
241                             IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
242                             IJK = FUNIJK(I,J,K)
243                             IM = IM1(I)
244                             JM = JM1(J)
245                             KM = KM1(K)
246     ! first set the boundary cell value equal to the known value in that
247     ! cell
248                             A_M(IJK,E,M) = ZERO
249                             A_M(IJK,W,M) = ZERO
250                             A_M(IJK,N,M) = ZERO
251                             A_M(IJK,S,M) = ZERO
252                             A_M(IJK,T,M) = ZERO
253                             A_M(IJK,B,M) = ZERO
254                             A_M(IJK,0,M) = -ONE
255                             B_M(IJK,M) = VAR(IJK)
256     ! second modify the matrix equation according to the user specified
257     ! boundary condition
258                             IF (FLUID_AT(EAST_OF(IJK))) THEN
259                                IF (BC_HW_PHI(L) == UNDEFINED) THEN
260     ! specified wall value (i.e., dirichlet type boundary)
261                                   A_M(IJK,E,M) = -HALF
262                                   A_M(IJK,0,M) = -HALF
263                                   B_M(IJK,M) = -BC_PHIW(L)
264                                ELSE
265     ! if bc_hw__phi=0 then specified flux boundary (i.e., neumann type
266     ! boundary) otherwise a mixed type boundary
267                                   A_M(IJK,0,M) = -(HALF*BC_HW_PHI(L)+ODX_E(I))
268                                   A_M(IJK,E,M) = -(HALF*BC_HW_PHI(L)-ODX_E(I))
269                                   B_M(IJK,M) = -(BC_HW_PHI(L)*BC_PHIW(L)+&
270                                                  BC_C_PHI(L))
271                                ENDIF
272                             ELSEIF (FLUID_AT(WEST_OF(IJK))) THEN
273                                IF (BC_HW_PHI(L) == UNDEFINED) THEN
274                                   A_M(IJK,W,M) = -HALF
275                                   A_M(IJK,0,M) = -HALF
276                                   B_M(IJK,M) = -BC_PHIW(L)
277                                ELSE
278                                   A_M(IJK,W,M) = -(HALF*BC_HW_PHI(L)-ODX_E(IM))
279                                   A_M(IJK,0,M) = -(HALF*BC_HW_PHI(L)+ODX_E(IM))
280                                   B_M(IJK,M) = -(BC_HW_PHI(L)*BC_PHIW(L)+&
281                                                  BC_C_PHI(L))
282                                ENDIF
283                             ELSEIF (FLUID_AT(NORTH_OF(IJK))) THEN
284                                IF (BC_HW_PHI(L) == UNDEFINED) THEN
285                                   A_M(IJK,N,M) = -HALF
286                                   A_M(IJK,0,M) = -HALF
287                                   B_M(IJK,M) = -BC_PHIW(L)
288                                ELSE
289                                   A_M(IJK,0,M) = -(HALF*BC_HW_PHI(L)+ODY_N(J))
290                                   A_M(IJK,N,M) = -(HALF*BC_HW_PHI(L)-ODY_N(J))
291                                   B_M(IJK,M) = -(BC_HW_PHI(L)*BC_PHIW(L)+&
292                                                  BC_C_PHI(L))
293                                ENDIF
294                             ELSEIF (FLUID_AT(SOUTH_OF(IJK))) THEN
295                                IF (BC_HW_PHI(L) == UNDEFINED) THEN
296                                   A_M(IJK,S,M) = -HALF
297                                   A_M(IJK,0,M) = -HALF
298                                   B_M(IJK,M) = -BC_PHIW(L)
299                                ELSE
300                                   A_M(IJK,S,M) = -(HALF*BC_HW_PHI(L)-ODY_N(JM))
301                                   A_M(IJK,0,M) = -(HALF*BC_HW_PHI(L)+ODY_N(JM))
302                                   B_M(IJK,M) = -(BC_HW_PHI(L)*BC_PHIW(L)+&
303                                                  BC_C_PHI(L))
304                                ENDIF
305                             ELSEIF (FLUID_AT(TOP_OF(IJK))) THEN
306                                IF (BC_HW_PHI(L) == UNDEFINED) THEN
307                                   A_M(IJK,T,M) = -HALF
308                                   A_M(IJK,0,M) = -HALF
309                                   B_M(IJK,M) = -BC_PHIW(L)
310                                ELSE
311                                   A_M(IJK,0,M)=-(HALF*BC_HW_PHI(L)+OX(I)*ODZ_T(K))
312                                   A_M(IJK,T,M)=-(HALF*BC_HW_PHI(L)-OX(I)*ODZ_T(K))
313                                   B_M(IJK,M) = -(BC_HW_PHI(L)*BC_PHIW(L)+&
314                                                  BC_C_PHI(L))
315                                ENDIF
316                             ELSEIF (FLUID_AT(BOTTOM_OF(IJK))) THEN
317                                IF (BC_HW_PHI(L) == UNDEFINED) THEN
318                                   A_M(IJK,B,M) = -HALF
319                                   A_M(IJK,0,M) = -HALF
320                                   B_M(IJK,M) = -BC_PHIW(L)
321                                ELSE
322                                   A_M(IJK,B,M) = -(HALF*BC_HW_PHI(L)-&
323                                                    OX(I)*ODZ_T(KM))
324                                   A_M(IJK,0,M) = -(HALF*BC_HW_PHI(L)+&
325                                                    OX(I)*ODZ_T(KM))
326                                   B_M(IJK,M) = -(BC_HW_PHI(L)*BC_PHIW(L)+&
327                                                  BC_C_PHI(L))
328                                ENDIF
329                             ENDIF
330                          ENDDO
331                       ENDDO
332                    ENDDO
333                ENDIF   ! end if (ns, fs, psw)
334              ENDIF   ! end if (bc_defined)
335           ENDDO   ! end L do loop over dimension_bc
336     ! end setting of wall boundary conditions
337     ! ----------------------------------------------------------------<<<
338     
339     
340     ! Set user defined boundary conditions for non-wall cells
341     ! Setting p_inflow, p_outflow, mass_outflow or outflow flow boundary
342     ! conditions
343     ! ---------------------------------------------------------------->>>
344           DO L = 1, DIMENSION_BC
345              IF (BC_DEFINED(L)) THEN
346                 IF (BC_TYPE(L)=='P_OUTFLOW' .OR. &
347                     BC_TYPE(L)=='MASS_OUTFLOW' .OR. &
348                     BC_TYPE(L)=='OUTFLOW') THEN
349                    I1 = BC_I_W(L)
350                    I2 = BC_I_E(L)
351                    J1 = BC_J_S(L)
352                    J2 = BC_J_N(L)
353                    K1 = BC_K_B(L)
354                    K2 = BC_K_T(L)
355     !!$omp    parallel do private(IJK, K, J, I)
356                    DO K = K1, K2
357                       DO J = J1, J2
358                          DO I = I1, I2
359                            IF (.NOT.IS_ON_myPE_plus2layers(I,J,K)) CYCLE
360                            IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
361                             IJK = FUNIJK(I,J,K)
362     ! first set the flow boundary cell value equal to zero
363                             A_M(IJK,E,M) = ZERO
364                             A_M(IJK,W,M) = ZERO
365                             A_M(IJK,N,M) = ZERO
366                             A_M(IJK,S,M) = ZERO
367                             A_M(IJK,T,M) = ZERO
368                             A_M(IJK,B,M) = ZERO
369                             A_M(IJK,0,M) = -ONE
370                             B_M(IJK,M) = ZERO
371     ! now set the flow boundary cell value equal to the adjacent fluid
372     ! cell value
373                             SELECT CASE (TRIM(BC_PLANE(L)))
374                             CASE ('E')
375     ! fluid cell on the east side
376                                A_M(IJK,E,M) = ONE
377                             CASE ('W')
378     ! fluid cell on the west side
379                                A_M(IJK,W,M) = ONE
380                             CASE ('N')
381                                A_M(IJK,N,M) = ONE
382                             CASE ('S')
383                                A_M(IJK,S,M) = ONE
384                             CASE ('T')
385                                A_M(IJK,T,M) = ONE
386                             CASE ('B')
387                                A_M(IJK,B,M) = ONE
388                             END SELECT
389                          ENDDO
390                       ENDDO
391                    ENDDO
392     ! end setting p_outflow, mass_outflow or outflow flow boundary
393     ! conditions
394     ! ----------------------------------------------------------------<<<
395     
396                 ELSEIF(BC_TYPE(L)=='P_INFLOW' .OR. &
397                        BC_TYPE(L)=='MASS_INFLOW') THEN
398     
399     ! Setting bc that are defined but not nsw, fsw, psw, p_outflow,
400     ! mass_outflow or outflow (at this time, this section addresses
401     ! p_inflow and mass_inflow type boundaries)
402     ! ----------------------------------------------------------------<<<
403                    I1 = BC_I_W(L)
404                    I2 = BC_I_E(L)
405                    J1 = BC_J_S(L)
406                    J2 = BC_J_N(L)
407                    K1 = BC_K_B(L)
408                    K2 = BC_K_T(L)
409     !!$omp    parallel do private(IJK, K, J, I)
410                    DO K = K1, K2
411                       DO J = J1, J2
412                          DO I = I1, I2
413                             IF (.NOT.IS_ON_myPE_plus2layers(I,J,K)) CYCLE
414                             IF (DEAD_CELL_AT(I,J,K)) CYCLE  ! skip dead cells
415                             IJK = FUNIJK(I,J,K)
416     ! setting the value in the boundary cell equal to what is known
417                             A_M(IJK,E,M) = ZERO
418                             A_M(IJK,W,M) = ZERO
419                             A_M(IJK,N,M) = ZERO
420                             A_M(IJK,S,M) = ZERO
421                             A_M(IJK,T,M) = ZERO
422                             A_M(IJK,B,M) = ZERO
423                             A_M(IJK,0,M) = -ONE
424     !                        B_M(IJK,M) = -BC_PHIF(L)  !does not allow the profile to be changed, e.g., from usr1
425                             B_M(IJK,M) = -VAR(IJK)
426                          ENDDO
427                       ENDDO
428                    ENDDO
429                 ENDIF   ! end if/else (bc_type)
430     ! end setting of p_inflow or mass_inflow boundary conditions
431     ! ----------------------------------------------------------------<<<
432     
433              ENDIF   ! end if (bc_defined)
434           ENDDO   ! end L do loop over dimension_bc
435     
436     
437     
438     ! modifications for cartesian grid implementation
439           IF(CARTESIAN_GRID .AND. .NOT.(CG_SAFE_MODE(1)==1)) &
440              CALL BC_PHI_CG(VAR, BC_PHIF, BC_PHIW, BC_HW_PHI, &
441                             BC_C_PHI, M, A_M, B_M)
442     
443           RETURN
444           END SUBROUTINE BC_PHI
445     
446     
447     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
448     !                                                                      C
449     !  Subroutine: BC_PHI_CG                                               C
450     !  Purpose: Modify boundary conditions for cartesian grid cut-cell     C
451     !           implementation                                             C
452     !                                                                      C
453     !  Author: Jeff Dietiker                                               C
454     !                                                                      C
455     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
456           SUBROUTINE BC_PHI_CG(VAR, BC_PHIF, BC_PHIW, BC_HW_PHI, &
457                             BC_C_PHI, M, A_M, B_M)
458     
459     !-----------------------------------------------
460     ! Modules
461     !-----------------------------------------------
462           USE param
463           USE param1
464           USE matrix
465           USE geometry
466           USE indices
467           USE bc
468           USE compar
469           USE cutcell
470           USE fun_avg
471           USE functions
472           IMPLICIT NONE
473     !-----------------------------------------------
474     ! Dummy arguments
475     !-----------------------------------------------
476     ! The field variable being solved for:
477     !     e.g., T_g, T_s, x_g, x_s, Theta_m, scalar, K_Turb_G,
478     !     e_Turb_G
479           DOUBLE PRECISION, INTENT(IN) :: VAR(DIMENSION_3)
480     ! Boundary conditions specifications
481     ! bc_phif = flow boundary value
482     ! bc_phiw = wall boundary value
483     ! bc_hw_phi = transfer coefficient
484     !      = 0 value means specified flux (neumann type)
485     !      = undefined value means specified wall value (dirichlet type)
486     !      = other value means mixed type
487     ! bc_C_phi = transfer flux
488           DOUBLE PRECISION, INTENT(IN) :: BC_phif(DIMENSION_BC), &
489                                           BC_Phiw(DIMENSION_BC), &
490                                           BC_hw_Phi(DIMENSION_BC), &
491                                           BC_C_Phi(DIMENSION_BC)
492     ! Phase index
493           INTEGER, INTENT(IN) :: M
494     ! Septadiagonal matrix A_m
495           DOUBLE PRECISION, INTENT(INOUT) :: A_m(DIMENSION_3, -3:3, 0:DIMENSION_M)
496     ! Vector b_m
497           DOUBLE PRECISION, INTENT(INOUT) :: B_m(DIMENSION_3, 0:DIMENSION_M)
498     !-----------------------------------------------
499     ! Local variables
500     !-----------------------------------------------
501     ! Indices
502           INTEGER :: I, J, K, IJK, IM, JM, KM
503     ! Boundary condition index
504           INTEGER :: L
505     ! Boundary identifiers
506           INTEGER :: BCV
507           CHARACTER(LEN=9) :: BCT
508     
509           LOGICAL :: ALONG_GLOBAL_GHOST_LAYER
510     
511     !-----------------------------------------------
512     
513           DO IJK = ijkstart3, ijkend3
514              I = I_OF(IJK)
515              J = J_OF(IJK)
516              K = K_OF(IJK)
517     
518              ALONG_GLOBAL_GHOST_LAYER = (I<IMIN1).OR.(I>IMAX1).OR.(J<JMIN1).OR.(J>JMAX1)
519     
520              IF(DO_K) ALONG_GLOBAL_GHOST_LAYER = ALONG_GLOBAL_GHOST_LAYER.OR.(K<KMIN1).OR.(K>KMAX1)
521     
522              IF(BLOCKED_CELL_AT(IJK)) THEN
523     ! setting the value in the boundary cell equal to what is known
524                 A_M(IJK,E,M) = ZERO
525                 A_M(IJK,W,M) = ZERO
526                 A_M(IJK,N,M) = ZERO
527                 A_M(IJK,S,M) = ZERO
528                 A_M(IJK,T,M) = ZERO
529                 A_M(IJK,B,M) = ZERO
530                 A_M(IJK,0,M) = -ONE
531                 B_M(IJK,M) = -VAR(IJK)
532              ENDIF
533     
534     
535              IF(BLOCKED_CELL_AT(IJK).OR.ALONG_GLOBAL_GHOST_LAYER) THEN
536                 IM = IM1(I)
537                 JM = JM1(J)
538                 KM = KM1(K)
539     
540                 IF (CUT_CELL_AT(IP_OF(IJK))) THEN
541                    BCV = BC_ID(IP_OF(IJK))
542                    IF(BCV > 0 ) THEN
543                       BCT = BC_TYPE(BCV)
544                    ELSE
545                       BCT = 'NONE'
546                    ENDIF
547     
548                    IF (BCT=='CG_NSW'.OR.BCT=='CG_FSW'.OR.BCT=='CG_PSW') THEN
549                       L = BCV
550                       IF (BC_HW_PHI(L) == UNDEFINED) THEN
551                          A_M(IJK,E,M) = -HALF
552                          A_M(IJK,0,M) = -HALF
553                          B_M(IJK,M) = -BC_PHIW(L)
554                       ELSE
555                          A_M(IJK,0,M) = -(HALF*BC_HW_PHI(L)+ODX_E(I))
556                          A_M(IJK,E,M) = -(HALF*BC_HW_PHI(L)-ODX_E(I))
557                          B_M(IJK,M) = -(BC_HW_PHI(L)*BC_PHIW(L)+BC_C_PHI(L))
558                       ENDIF
559                    ENDIF
560     
561                 ELSEIF (CUT_CELL_AT(IM_OF(IJK))) THEN
562                    BCV = BC_ID(IM_OF(IJK))
563                    IF(BCV > 0 ) THEN
564                       BCT = BC_TYPE(BCV)
565                    ELSE
566                       BCT = 'NONE'
567                    ENDIF
568     
569                    IF (BCT=='CG_NSW'.OR.BCT=='CG_FSW'.OR.BCT=='CG_PSW') THEN
570                       L = BCV
571                       IF (BC_HW_PHI(L) == UNDEFINED) THEN
572                          A_M(IJK,W,M) = -HALF
573                          A_M(IJK,0,M) = -HALF
574                          B_M(IJK,M) = -BC_PHIW(L)
575                       ELSE
576                          A_M(IJK,W,M) = -(HALF*BC_HW_PHI(L)-ODX_E(IM))
577                          A_M(IJK,0,M) = -(HALF*BC_HW_PHI(L)+ODX_E(IM))
578                          B_M(IJK,M) = -(BC_HW_PHI(L)*BC_PHIW(L)+BC_C_PHI(L))
579                       ENDIF
580                    ENDIF
581     
582                 ELSEIF (CUT_CELL_AT(JP_OF(IJK))) THEN
583                    BCV = BC_ID(JP_OF(IJK))
584                    IF(BCV > 0 ) THEN
585                       BCT = BC_TYPE(BCV)
586                    ELSE
587                       BCT = 'NONE'
588                    ENDIF
589     
590                    IF (BCT=='CG_NSW'.OR.BCT=='CG_FSW'.OR.BCT=='CG_PSW') THEN
591                       L = BCV
592                       IF (BC_HW_PHI(L) == UNDEFINED) THEN
593                          A_M(IJK,N,M) = -HALF
594                          A_M(IJK,0,M) = -HALF
595                          B_M(IJK,M) = -BC_PHIW(L)
596                       ELSE
597                          A_M(IJK,0,M) = -(HALF*BC_HW_PHI(L)+ODY_N(J))
598                          A_M(IJK,N,M) = -(HALF*BC_HW_PHI(L)-ODY_N(J))
599                          B_M(IJK,M) = -(BC_HW_PHI(L)*BC_PHIW(L)+BC_C_PHI(L))
600                       ENDIF
601                    ENDIF
602     
603                 ELSEIF (CUT_CELL_AT(JM_OF(IJK))) THEN
604                    BCV = BC_ID(JM_OF(IJK))
605                    IF(BCV > 0 ) THEN
606                       BCT = BC_TYPE(BCV)
607                    ELSE
608                       BCT = 'NONE'
609                    ENDIF
610     
611                    IF (BCT=='CG_NSW'.OR.BCT=='CG_FSW'.OR.BCT=='CG_PSW') THEN
612                       L = BCV
613                       IF (BC_HW_PHI(L) == UNDEFINED) THEN
614                          A_M(IJK,S,M) = -HALF
615                          A_M(IJK,0,M) = -HALF
616                          B_M(IJK,M) = -BC_PHIW(L)
617                       ELSE
618                          A_M(IJK,S,M) = -(HALF*BC_HW_PHI(L)-ODY_N(JM))
619                          A_M(IJK,0,M) = -(HALF*BC_HW_PHI(L)+ODY_N(JM))
620                          B_M(IJK,M) = -(BC_HW_PHI(L)*BC_PHIW(L)+BC_C_PHI(L))
621                       ENDIF
622                    ENDIF
623     
624                 ELSEIF (CUT_CELL_AT(KP_OF(IJK))) THEN
625                    BCV = BC_ID(KP_OF(IJK))
626                    IF(BCV > 0 ) THEN
627                       BCT = BC_TYPE(BCV)
628                    ELSE
629                       BCT = 'NONE'
630                    ENDIF
631     
632                    IF (BCT=='CG_NSW'.OR.BCT=='CG_FSW'.OR.BCT=='CG_PSW') THEN
633                       L = BCV
634                       IF (BC_HW_PHI(L) == UNDEFINED) THEN
635                          A_M(IJK,T,M) = -HALF
636                          A_M(IJK,0,M) = -HALF
637                          B_M(IJK,M) = -BC_PHIW(L)
638                       ELSE
639                          A_M(IJK,0,M)=-(HALF*BC_HW_PHI(L)+OX(I)*ODZ_T(K))
640                          A_M(IJK,T,M)=-(HALF*BC_HW_PHI(L)-OX(I)*ODZ_T(K))
641                          B_M(IJK,M) = -(BC_HW_PHI(L)*BC_PHIW(L)+BC_C_PHI(L))
642                       ENDIF
643                    ENDIF
644     
645                 ELSEIF (CUT_CELL_AT(KM_OF(IJK))) THEN
646                    BCV = BC_ID(KM_OF(IJK))
647                    IF(BCV > 0 ) THEN
648                       BCT = BC_TYPE(BCV)
649                    ELSE
650                       BCT = 'NONE'
651                    ENDIF
652     
653                    IF (BCT=='CG_NSW'.OR.BCT=='CG_FSW'.OR.BCT=='CG_PSW') THEN
654                       L = BCV
655                       IF (BC_HW_PHI(L) == UNDEFINED) THEN
656                          A_M(IJK,B,M) = -HALF
657                          A_M(IJK,0,M) = -HALF
658                          B_M(IJK,M) = -BC_PHIW(L)
659                       ELSE
660                          A_M(IJK,B,M) = -(HALF*BC_HW_PHI(L)-OX(I)*ODZ_T(KM))
661                          A_M(IJK,0,M) = -(HALF*BC_HW_PHI(L)+OX(I)*ODZ_T(KM))
662                          B_M(IJK,M) = -(BC_HW_PHI(L)*BC_PHIW(L)+BC_C_PHI(L))
663                       ENDIF
664                    ENDIF  ! BCT
665     
666                 ENDIF   ! end if/else cut cell next to blocked cell
667     
668              ENDIF   ! end if (blocked_cell_at(ijk))
669     
670              IF(CUT_CELL_AT(IJK)) THEN
671                 BCV = BC_ID(IJK)
672                 IF(BCV > 0 ) THEN
673                    BCT = BC_TYPE(BCV)
674                 ELSE
675                    BCT = 'NONE'
676                 ENDIF
677     
678                 IF (BCT=='CG_MI'.OR.BCT=='CG_PO') THEN
679                    L = BCV
680                    A_M(IJK,E,M) = ZERO
681                    A_M(IJK,W,M) = ZERO
682                    A_M(IJK,N,M) = ZERO
683                    A_M(IJK,S,M) = ZERO
684                    A_M(IJK,T,M) = ZERO
685                    A_M(IJK,B,M) = ZERO
686                    A_M(IJK,0,M) = -ONE
687                    B_M(IJK,M) = -BC_PHIF(L)
688                 ENDIF
689              ENDIF
690     
691           ENDDO   ! end do (ijk=ijkstart3,ijkend3)
692     
693     
694           RETURN
695           END SUBROUTINE BC_PHI_CG
696