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