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