MFIX  2016-1
bc_phi.f
Go to the documentation of this file.
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
bc::bc_k_b
integer, dimension(dimension_bc) bc_k_b
Definition: bc_mod.f:70
geometry::imax2
integer imax2
Definition: geometry_mod.f:61
indices::i_of
integer, dimension(:), allocatable i_of
Definition: indices_mod.f:45
compar::ijkend3
integer ijkend3
Definition: compar_mod.f:80
geometry::imax3
integer imax3
Definition: geometry_mod.f:91
param1::one
double precision, parameter one
Definition: param1_mod.f:29
param::dimension_3
integer dimension_3
Definition: param_mod.f:11
bc::bc_i_w
integer, dimension(dimension_bc) bc_i_w
Definition: bc_mod.f:54
bc::bc_j_n
integer, dimension(dimension_bc) bc_j_n
Definition: bc_mod.f:66
indices::im1
integer, dimension(:), allocatable im1
Definition: indices_mod.f:50
cutcell::cg_safe_mode
integer, dimension(10) cg_safe_mode
Definition: cutcell_mod.f:415
param::dimension_bc
integer, parameter dimension_bc
Definition: param_mod.f:61
bc::bc_type_enum
integer, dimension(dimension_bc) bc_type_enum
Definition: bc_mod.f:146
param1::undefined
double precision, parameter undefined
Definition: param1_mod.f:18
param::east
integer east
Definition: param_mod.f:29
geometry::imin3
integer imin3
Definition: geometry_mod.f:90
bc::bc_plane
character, dimension(dimension_bc) bc_plane
Definition: bc_mod.f:217
indices::k_of
integer, dimension(:), allocatable k_of
Definition: indices_mod.f:47
geometry::ody_n
double precision, dimension(:), allocatable ody_n
Definition: geometry_mod.f:123
fun_avg
Definition: fun_avg_mod.f:1
geometry::kmax1
integer kmax1
Definition: geometry_mod.f:58
bc::bc_k_t
integer, dimension(dimension_bc) bc_k_t
Definition: bc_mod.f:74
compar::dead_cell_at
logical, dimension(:,:,:), allocatable dead_cell_at
Definition: compar_mod.f:127
indices::j_of
integer, dimension(:), allocatable j_of
Definition: indices_mod.f:46
geometry::odx_e
double precision, dimension(:), allocatable odx_e
Definition: geometry_mod.f:121
indices::jm1
integer, dimension(:), allocatable jm1
Definition: indices_mod.f:51
geometry::imax1
integer imax1
Definition: geometry_mod.f:54
geometry::jmax2
integer jmax2
Definition: geometry_mod.f:63
bc::bc_j_s
integer, dimension(dimension_bc) bc_j_s
Definition: bc_mod.f:62
geometry::ox
double precision, dimension(:), allocatable ox
Definition: geometry_mod.f:140
bc_phi
subroutine bc_phi(VAR, BC_PHIF, BC_PHIW, BC_HW_PHI, BC_C_PHI, M, A_M, B_M)
Definition: bc_phi.f:13
geometry::jmax3
integer jmax3
Definition: geometry_mod.f:91
param::north
integer north
Definition: param_mod.f:37
bc::bc_defined
logical, dimension(dimension_bc) bc_defined
Definition: bc_mod.f:207
param::south
integer south
Definition: param_mod.f:41
geometry::kmax2
integer kmax2
Definition: geometry_mod.f:65
param1::half
double precision, parameter half
Definition: param1_mod.f:28
geometry::jmax1
integer jmax1
Definition: geometry_mod.f:56
param
Definition: param_mod.f:2
geometry::jmin3
integer jmin3
Definition: geometry_mod.f:90
cutcell::cartesian_grid
logical cartesian_grid
Definition: cutcell_mod.f:13
geometry::jmin1
integer jmin1
Definition: geometry_mod.f:42
geometry::kmax3
integer kmax3
Definition: geometry_mod.f:91
geometry::do_k
logical do_k
Definition: geometry_mod.f:30
cutcell::cut_cell_at
logical, dimension(:), allocatable cut_cell_at
Definition: cutcell_mod.f:355
indices::km1
integer, dimension(:), allocatable km1
Definition: indices_mod.f:52
compar::ijkstart3
integer ijkstart3
Definition: compar_mod.f:80
param::west
integer west
Definition: param_mod.f:33
geometry::kmin3
integer kmin3
Definition: geometry_mod.f:90
cutcell::bc_id
integer, dimension(:), allocatable bc_id
Definition: cutcell_mod.f:433
param::top
integer top
Definition: param_mod.f:45
geometry::imin2
integer imin2
Definition: geometry_mod.f:89
bc_phi_cg
subroutine bc_phi_cg(VAR, BC_PHIF, BC_PHIW, BC_HW_PHI, BC_C_PHI, M, A_M, B_M)
Definition: bc_phi.f:457
cutcell::blocked_cell_at
logical, dimension(:), allocatable blocked_cell_at
Definition: cutcell_mod.f:361
geometry::imin1
integer imin1
Definition: geometry_mod.f:40
param::dimension_m
integer dimension_m
Definition: param_mod.f:18
geometry::odz_t
double precision, dimension(:), allocatable odz_t
Definition: geometry_mod.f:125
geometry::kmin1
integer kmin1
Definition: geometry_mod.f:44
param::bottom
integer bottom
Definition: param_mod.f:49
bc::bc_i_e
integer, dimension(dimension_bc) bc_i_e
Definition: bc_mod.f:58
param1::zero
double precision, parameter zero
Definition: param1_mod.f:27
bc
Definition: bc_mod.f:23