File: N:\mfix\model\set_outflow.f

1     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
2     !                                                                      C
3     !  Subroutine: SET_OUTFLOW                                             C
4     !  Purpose: Set specified outflow bc for pressure outflow,             C
5     !  mass outflow, outflow and now also pressure inflow bc               C
6     !                                                                      C
7     !  Author: M. Syamlal                                 Date: 21-JAN-92  C
8     !                                                                      C
9     !  Comments:                                                           C
10     !  If the outflow boundary is on the W, S or B side of the domain and  C
11     !  the component of velocity through the plane is defined then this    C
12     !  routine will NOT modify it (i.e., its value from the momentum       C
13     !  solver is maintained).                                              C
14     !                                                                      C
15     !  In general it would seem this routine does little more than what    C
16     !  is already done within the momentum routines in terms of velocity.  C
17     !  The normal component is either 1) untouched if the outflow is on    C
18     !  the W, S, B side of the domain or 2) is set to value of the         C
19     !  adjacent fluid cell if it is on the E, N, T side (similar to the    C
20     !  respective momentum bc routines). The primary addition here is      C
21     !  that the tangential components of a bc cell are set to that of      C
22     !  the adjacent fluid cell. Note the tangential components are not     C
23     !  explicitly handled in the momentum _BC_ routines; instead their     C
24     !  values are based on solution of the momentum equation which is      C
25     !  replaced here                                                       C
26     !                                                                      C
27     !  Several routines are called which perform the following tasks:      C
28     !  set_outflow_misc - several derived quantities are set in the        C
29     !      boundary                                                        C
30     !  set_outflow_ep - the void/volume fraction and bulk densities are    C
31     !      set in the boundary                                             C
32     !  set_outflow_fluxes - convective fluxes are set in the boundary      C
33     !                                                                      C
34     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
35           SUBROUTINE SET_OUTFLOW(BCV)
36     
37     ! Modules
38     !---------------------------------------------------------------------//
39           use bc
40           use bc, only: bc_k_b, bc_k_t
41           use bc, only: bc_j_s, bc_j_n
42           use bc, only: bc_i_w, bc_i_e
43     
44           use fldvar, only: rop_g, rop_s
45           use fldvar, only: u_g, v_g, w_g
46           use fldvar, only: u_s, v_s, w_s
47           use physprop, only: mmax
48     
49           use functions, only: im_of, ip_of, jm_of, jp_of, km_of, kp_of
50           use functions, only: fluid_at
51           use functions, only: is_on_mype_plus2layers
52           use functions, only: funijk
53           use compar, only: dead_cell_at
54     
55           use param, only: dimension_m
56           use param1, only: undefined, zero
57           IMPLICIT NONE
58     
59     ! Dummy arguments
60     !---------------------------------------------------------------------//
61     ! Boundary condition number
62           INTEGER, INTENT(IN) :: BCV
63     
64     ! Local variables
65     !---------------------------------------------------------------------//
66     ! indices
67           INTEGER :: I, J, K, M
68     ! index for boundary cell
69           INTEGER :: IJK
70     ! index for a fluid cell adjacent to the boundary cell
71           INTEGER :: FIJK
72     ! local value for normal component of gas and solids velocity defined
73     ! such that
74           DOUBLE PRECISION :: RVEL_G, RVEL_S(DIMENSION_M)
75     !---------------------------------------------------------------------//
76     
77     ! Loop over the range of boundary cells
78           DO K = BC_K_B(BCV), BC_K_T(BCV)
79              DO J = BC_J_S(BCV), BC_J_N(BCV)
80                 DO I = BC_I_W(BCV), BC_I_E(BCV)
81     ! Check if current i,j,k resides on this PE
82                    IF (.NOT.IS_ON_myPE_plus2layers(I,J,K)) CYCLE
83                    IF(DEAD_CELL_AT(I,J,K)) CYCLE
84                    IJK = FUNIJK(I,J,K)
85     
86     ! Fluid cell at West
87     ! --------------------------------------------------------------------//
88                    IF (FLUID_AT(IM_OF(IJK))) THEN
89                       FIJK = IM_OF(IJK)
90                       RVEL_G = U_G(FIJK)
91                       IF (MMAX>0) RVEL_S(:MMAX) = U_S(FIJK,:MMAX)
92     
93                       CALL SET_OUTFLOW_MISC(BCV, IJK, FIJK)
94                       CALL SET_OUTFLOW_EP(BCV, IJK, FIJK, RVEL_G, RVEL_S)
95                       IF (BC_TYPE_ENUM(BCV)==P_INFLOW) &
96                          CALL SET_PINOUTFLOW(BCV, IJK, FIJK, RVEL_G, RVEL_S)
97     
98     ! Set the boundary cell value of the normal component of velocity
99     ! according to the value in the adjacent fluid cell. Note the value
100     ! of the boundary velocity is a scaled version of the value of the
101     ! adjacent fluid velocity based on the concentration ratio of the fluid
102     ! cell to the boundary cell.
103     ! - For the gas phase, this ratio is most likely 1 except for
104     !   compressible cases with a PO/PI boundary where P_g of the boundary
105     !   is set and may differ from the value of the adjacent fluid cell.
106     ! - For the solids phase this seems unnecessary..? Differences may arise
107     !   if bc_rop_s is set...
108                       IF (ROP_G(IJK) > ZERO) THEN
109                         U_G(IJK) = ROP_G(FIJK)*U_G(FIJK)/ROP_G(IJK)
110                       ELSE
111                          U_G(IJK) = ZERO
112                       ENDIF
113     
114     ! the tangential components are not explicitly handled in the boundary
115     ! condition routines of the corresponding momentum equation
116                       V_G(IJK) = V_G(FIJK)
117                       W_G(IJK) = W_G(FIJK)
118     
119                       IF (MMAX > 0) THEN
120                          WHERE (ROP_S(IJK,:MMAX) > ZERO)
121                            U_S(IJK,:MMAX) = ROP_S(FIJK,:MMAX)*&
122                                U_S(FIJK,:MMAX)/ROP_S(IJK,:MMAX)
123                          ELSEWHERE
124                             U_S(IJK,:MMAX) = ZERO
125                          END WHERE
126                          V_S(IJK,:MMAX) = V_S(FIJK,:MMAX)
127                          W_S(IJK,:MMAX) = W_S(FIJK,:MMAX)
128                       ENDIF
129     
130                       CALL SET_OUTFLOW_FLUXES(IJK, FIJK)
131                    ENDIF   ! end if (fluid_at(im_of(ijk)))
132     
133     
134     ! Fluid cell at East
135     ! --------------------------------------------------------------------//
136                    IF (FLUID_AT(IP_OF(IJK))) THEN
137                       FIJK = IP_OF(IJK)
138     ! define normal component such that it is positive when exiting the
139     ! domain
140                       RVEL_G = -U_G(IJK)
141                       IF (MMAX >0) RVEL_S(:MMAX) = -U_S(IJK,:MMAX)
142     
143                       CALL SET_OUTFLOW_MISC(BCV, IJK, FIJK)
144                       CALL SET_OUTFLOW_EP(BCV, IJK, FIJK, RVEL_G, RVEL_S)
145                       IF (BC_TYPE_ENUM(BCV)==P_INFLOW) &
146                          CALL SET_PINOUTFLOW(BCV, IJK, FIJK, RVEL_G, RVEL_S)
147     
148     ! provide an initial value for the velocity component through the domain
149     ! otherwise its present value (from solution of the corresponding
150     ! momentum eqn) is kept. values for the velocity components in the off
151     ! directions are modified (needed for PO or O boundaries but not MO or
152     ! PI as velocities should be fully specified by this point)
153                       IF (U_G(IJK) == UNDEFINED) THEN
154                          IF (ROP_G(IJK) > ZERO) THEN
155                             U_G(IJK) = ROP_G(FIJK)*U_G(FIJK)/ROP_G(IJK)
156                          ELSE
157                             U_G(IJK) = ZERO
158                          ENDIF
159                       ENDIF
160                       V_G(IJK) = V_G(FIJK)
161                       W_G(IJK) = W_G(FIJK)
162     
163                       DO M = 1, MMAX
164                          IF (U_S(IJK,M) == UNDEFINED) THEN
165                             IF (ROP_S(IJK,M) > ZERO) THEN
166                                U_S(IJK,M) = ROP_S(FIJK,M)*&
167                                   U_S(FIJK,M)/ROP_S(IJK,M)
168                             ELSE
169                                U_S(IJK,M) = ZERO
170                             ENDIF
171                          ENDIF
172                          V_S(IJK,M) = V_S(FIJK,M)
173                          W_S(IJK,M) = W_S(FIJK,M)
174                       ENDDO
175     
176                       CALL SET_OUTFLOW_FLUXES(IJK, FIJK)
177                    ENDIF   ! end if (fluid_at(ip_of(ijk)))
178     
179     
180     ! Fluid cell at South
181     ! --------------------------------------------------------------------//
182                    IF (FLUID_AT(JM_OF(IJK))) THEN
183                       FIJK = JM_OF(IJK)
184                       RVEL_G = V_G(FIJK)
185                       IF(MMAX>0) RVEL_S(:MMAX) = V_S(FIJK,:MMAX)
186     
187                       CALL SET_OUTFLOW_MISC(BCV, IJK, FIJK)
188                       CALL SET_OUTFLOW_EP(BCV, IJK, FIJK, RVEL_G, RVEL_S)
189                       IF (BC_TYPE_ENUM(BCV)==P_INFLOW) &
190                          CALL SET_PINOUTFLOW(BCV, IJK, FIJK, RVEL_G, RVEL_S)
191     
192                       IF (ROP_G(IJK) > ZERO) THEN
193                          V_G(IJK) = ROP_G(FIJK)*V_G(FIJK)/ROP_G(IJK)
194                       ELSE
195                          V_G(IJK) = ZERO
196                       ENDIF
197                       U_G(IJK) = U_G(FIJK)
198                       W_G(IJK) = W_G(FIJK)
199     
200                       IF (MMAX > 0) THEN
201                           WHERE (ROP_S(IJK,:MMAX) > ZERO)
202                              V_S(IJK,:MMAX) = ROP_S(FIJK,:MMAX)*&
203                                 V_S(FIJK,:MMAX)/ROP_S(IJK,:MMAX)
204                           ELSEWHERE
205                              V_S(IJK,:MMAX) = ZERO
206                           END WHERE
207                           U_S(IJK,:MMAX) = U_S(FIJK,:MMAX)
208                           W_S(IJK,:MMAX) = W_S(FIJK,:MMAX)
209                       ENDIF
210     
211                       CALL SET_OUTFLOW_FLUXES(IJK, FIJK)
212                    ENDIF   ! end if (fluid_at(jm_of(ijk)))
213     
214     
215     ! Fluid cell at North
216     ! --------------------------------------------------------------------//
217                    IF (FLUID_AT(JP_OF(IJK))) THEN
218                       FIJK = JP_OF(IJK)
219                       RVEL_G = -V_G(IJK)
220                       IF (MMAX>0) RVEL_S(:MMAX) = -V_S(IJK,:MMAX)
221     
222                       CALL SET_OUTFLOW_MISC(BCV, IJK, FIJK)
223                       CALL SET_OUTFLOW_EP(BCV, IJK, FIJK, RVEL_G, RVEL_S)
224                       IF (BC_TYPE_ENUM(BCV)==P_INFLOW) &
225                          CALL SET_PINOUTFLOW(BCV, IJK, FIJK, RVEL_G, RVEL_S)
226     
227                       IF (V_G(IJK) == UNDEFINED) THEN
228                          IF (ROP_G(IJK) > ZERO) THEN
229                             V_G(IJK) = ROP_G(FIJK)*V_G(FIJK)/ROP_G(IJK)
230                          ELSE
231                             V_G(IJK) = ZERO
232                          ENDIF
233                       ENDIF
234                       U_G(IJK) = U_G(FIJK)
235                       W_G(IJK) = W_G(FIJK)
236     
237                       DO M = 1, MMAX
238                          IF (V_S(IJK,M) == UNDEFINED) THEN
239                             IF (ROP_S(IJK,M) > ZERO) THEN
240                                V_S(IJK,M) = ROP_S(FIJK,M)*&
241                                   V_S(FIJK,M)/ROP_S(IJK,M)
242                             ELSE
243                                V_S(IJK,M) = ZERO
244                             ENDIF
245                          ENDIF
246                          U_S(IJK,M) = U_S(FIJK,M)
247                          W_S(IJK,M) = W_S(FIJK,M)
248                       ENDDO
249     
250                       CALL SET_OUTFLOW_FLUXES(IJK, FIJK)
251                    ENDIF   ! if (fluid_at(jp_of(ijk)))
252     
253     
254     ! Fluid cell at Bottom
255     ! --------------------------------------------------------------------//
256                    IF (FLUID_AT(KM_OF(IJK))) THEN
257                       FIJK = KM_OF(IJK)
258                       RVEL_G = W_G(FIJK)
259                       IF (MMAX>0) RVEL_S(:MMAX) = W_S(FIJK,:MMAX)
260     
261                       CALL SET_OUTFLOW_MISC(BCV, IJK, FIJK)
262                       CALL SET_OUTFLOW_EP(BCV, IJK, FIJK, RVEL_G, RVEL_S)
263                       IF (BC_TYPE_ENUM(BCV)==P_INFLOW) &
264                          CALL SET_PINOUTFLOW(BCV, IJK, FIJK, RVEL_G, RVEL_S)
265     
266                       IF (ROP_G(IJK) > ZERO) THEN
267                          W_G(IJK) = ROP_G(FIJK)*W_G(FIJK)/ROP_G(IJK)
268                       ELSE
269                          W_G(IJK) = ZERO
270                       ENDIF
271                       U_G(IJK) = U_G(FIJK)
272                       V_G(IJK) = V_G(FIJK)
273     
274                       IF (MMAX > 0) THEN
275                         WHERE (ROP_S(IJK,:MMAX) > ZERO)
276                             W_S(IJK,:MMAX) = ROP_S(FIJK,:MMAX)*&
277                                W_S(FIJK,:MMAX)/ROP_S(IJK,:MMAX)
278                          ELSEWHERE
279                             W_S(IJK,:MMAX) = ZERO
280                          END WHERE
281                          U_S(IJK,:MMAX) = U_S(FIJK,:MMAX)
282                          V_S(IJK,:MMAX) = V_S(FIJK,:MMAX)
283                       ENDIF
284     
285                       CALL SET_OUTFLOW_FLUXES(IJK, FIJK)
286                    ENDIF   ! if (fluid_at(km_of(ijk)))
287     
288     
289     ! Fluid cell at Top
290     ! --------------------------------------------------------------------//
291                    IF (FLUID_AT(KP_OF(IJK))) THEN
292                       FIJK = KP_OF(IJK)
293                       RVEL_G = -W_G(IJK)
294                       IF (MMAX>0) RVEL_S(:MMAX) = -W_S(IJK,:MMAX)
295     
296                       CALL SET_OUTFLOW_MISC(BCV, IJK, FIJK)
297                       CALL SET_OUTFLOW_EP(BCV, IJK, FIJK, RVEL_G, RVEL_S)
298                       IF (BC_TYPE_ENUM(BCV)==P_INFLOW) &
299                          CALL SET_PINOUTFLOW(BCV, IJK, FIJK, RVEL_G, RVEL_S)
300     
301                       IF (W_G(IJK) == UNDEFINED) THEN
302                          IF (ROP_G(IJK) > ZERO) THEN
303                             W_G(IJK) = ROP_G(FIJK)*W_G(FIJK)/ROP_G(IJK)
304                          ELSE
305                             W_G(IJK) = ZERO
306                          ENDIF
307                       ENDIF
308                       U_G(IJK) = U_G(FIJK)
309                       V_G(IJK) = V_G(FIJK)
310     
311                       DO M = 1, MMAX
312                          IF (W_S(IJK,M) == UNDEFINED) THEN
313                             IF (ROP_S(IJK,M) > ZERO) THEN
314                                W_S(IJK,M) = ROP_S(FIJK,M)*&
315                                   W_S(FIJK,M)/ROP_S(IJK,M)
316                             ELSE
317                                W_S(IJK,M) = ZERO
318                             ENDIF
319                          ENDIF
320                          U_S(IJK,M) = U_S(FIJK,M)
321                          V_S(IJK,M) = V_S(FIJK,M)
322                       ENDDO
323     
324                       CALL SET_OUTFLOW_FLUXES(IJK, FIJK)
325                    ENDIF   ! if (fluid_at(kp_of(ijk)))
326     
327                 ENDDO   ! end do (i=i1,i2)
328              ENDDO   ! end do (j=j1,j2)
329           ENDDO   ! end do (k=k1,k2)
330     
331           RETURN
332           END SUBROUTINE SET_OUTFLOW
333     
334     
335     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
336     !                                                                      C
337     !  Subroutine: SET_OUTFLOW_MISC                                        C
338     !  Purpose: Set the value of certain variables in the specified        C
339     !  outflow boundary cell that would not otherwise be set according     C
340     !  to their value in the adjacent fluid cell.                          C
341     !                                                                      C
342     !  READ BEFORE MODIFYING!!                                             C
343     !  For many of the scalar field variables (e.g., T_g, T_s, X_g, X_s,   C
344     !  k_turb_g, e_turb_g, theta_m and scalar) their corresponding         C
345     !  governing equation solver routine sets its own value in the         C
346     !  boundary cell (see bc_phi). So DO NOT set their value here unless   C
347     !  a clear explanation to the need is also provided; more likely it    C
348     !  should simply be dealt with in the governing equation routine.      C
349     !                                                                      C
350     !  This routine sets the value of quantites that are derived or with   C
351     !  unique solver routines (e.g., P_star, P_s, MW_MIX_g, P_g) since no  C
352     !  other routine will.                                                 C
353     !                                                                      C
354     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
355           SUBROUTINE SET_OUTFLOW_MISC(BCV, IJK, FIJK)
356     
357     ! Global variables
358     !---------------------------------------------------------------------//
359           use bc
360           use run, only: kt_type_enum, ghd_2007
361           use fldvar, only: p_g, ro_g, T_g
362           use fldvar, only: p_s, p_star
363           use physprop, only: smax, mmax
364           use physprop, only: ro_g0, mw_mix_g
365           use eos, only: EOSG
366     
367     ! Global parameters
368     !---------------------------------------------------------------------//
369           use param1, only: undefined
370           IMPLICIT NONE
371     
372     ! Dummy arguments
373     !---------------------------------------------------------------------//
374     ! Boundary condition number
375           INTEGER, INTENT(IN) :: BCV
376     ! ijk index for boundary cell
377           INTEGER, INTENT(IN) :: IJK
378     ! ijk index for adjacent fluid cell
379           INTEGER, INTENT(IN) :: FIJK
380     !---------------------------------------------------------------------//
381     
382           IF (BC_TYPE_ENUM(BCV) /= P_OUTFLOW .AND. &
383               BC_TYPE_ENUM(BCV) /= P_INFLOW) P_G(IJK) = P_G(FIJK)
384     
385           MW_MIX_G(IJK) = MW_MIX_G(FIJK)
386     
387     ! T_g (bc_phi), P_g (above depending on bc), and MW_MIX_G (above) have
388     ! now been defined at IJK
389           IF (RO_G0 == UNDEFINED) RO_G(IJK) = &
390              EOSG(MW_MIX_G(IJK),P_G(IJK),T_G(IJK))
391     
392           IF (SMAX >0) THEN
393              P_STAR(IJK) = P_STAR(FIJK)
394              P_S(IJK,:SMAX) = P_S(FIJK,:SMAX)
395           ENDIF
396     
397           IF (KT_TYPE_ENUM == GHD_2007 .AND. MMAX>0) &
398              P_S(IJK,MMAX) = P_S(FIJK,MMAX)
399     
400           RETURN
401           END SUBROUTINE SET_OUTFLOW_MISC
402     
403     
404     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
405     !                                                                      C
406     !  Subroutine: SET_OUTFLOW_EP                                          C
407     !  Purpose: Set the volume fraction/bulk density (i.e., ROP_s, EP_g    C
408     !  ROP_S) in the specified outflow boundary cell that would not        C
409     !  otherwise be set according to their value in the adjacent fluid     C
410     !  cell.                                                               C
411     !                                                                      C
412     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
413           SUBROUTINE SET_OUTFLOW_EP(BCV, IJK, FIJK, RVEL_G, RVEL_S)
414     
415     ! Global variables
416     !---------------------------------------------------------------------//
417           use bc
418           use run, only: kt_type_enum, ghd_2007
419           use physprop, only: smax, mmax
420           use fldvar, only: rop_g, ro_g, ep_g
421           use fldvar, only: rop_s, ep_s
422           use discretelement, only: discrete_element, des_mmax
423     ! Global parameters
424     !---------------------------------------------------------------------//
425           use param, only: dimension_m
426           use param1, only: undefined, zero, one
427           IMPLICIT NONE
428     
429     ! Dummy arguments
430     !---------------------------------------------------------------------//
431     ! Boundary condition number
432           INTEGER, INTENT(IN) :: BCV
433     ! ijk index for boundary cell
434           INTEGER, INTENT(IN) :: IJK
435     ! ijk index for adjacent fluid cell
436           INTEGER, INTENT(IN) :: FIJK
437     ! the gas or solids velocity in the fluid cell adjacent to the boundary
438     ! cell dot with the outward normal of that bc plane; defines the gas or
439     ! solids velocity component normal to the bc plane as positive when it
440     ! is flowing into the bc cell from the fluid cell. so, for example, for
441     ! an outflow on the eastern boundary this is the u component of velocity
442     ! while for an outflow on the western boundary this is the -u component,
443     ! etc.
444           DOUBLE PRECISION, INTENT(IN) :: RVEL_G
445           DOUBLE PRECISION, INTENT(IN), DIMENSION(DIMENSION_M) :: RVEL_S
446     
447     ! Local variables
448     !---------------------------------------------------------------------//
449     ! indices
450           INTEGER :: M
451     ! sum of solids phases volume fractions
452           DOUBLE PRECISION :: SUM_EPs
453     ! sum of solids phases bulk densities
454           DOUBLE PRECISION :: SUM_ROPS
455           LOGICAL, SAVE :: FIRST_PASS = .TRUE.
456     !---------------------------------------------------------------------//
457     
458     ! initializing summation quantities
459           SUM_ROPS = ZERO
460           SUM_EPS = ZERO
461     
462           DO M = 1, SMAX
463     
464              IF(BC_TYPE_ENUM(BCV) == P_INFLOW) THEN
465                 ROP_S(IJK,M) = ROP_S(FIJK,M)
466              ELSE
467     ! the outflow type bc do not permit 're-entering' solids, in which
468     ! case the solids are removed
469                 IF (RVEL_S(M) >= ZERO) THEN
470     ! solids are leaving the domain
471                    ROP_S(IJK,M) = ROP_S(FIJK,M)
472                 ELSE
473     ! solids cannot enter the domain at an outflow cell
474                    ROP_S(IJK,M) = ZERO
475                 ENDIF
476              ENDIF
477     
478     ! if bc_rop_s is defined, set value of rop_s in the ijk boundary cell
479     ! according to user definition
480              IF(BC_ROP_S(BCV,M)/=UNDEFINED) ROP_S(IJK,M)=BC_ROP_S(BCV,M)
481     
482     ! add to total solids phase bulk density and solids volume fraction
483              SUM_ROPS = SUM_ROPS + ROP_S(IJK,M)
484              SUM_EPS = SUM_EPS + EP_S(IJK,M)
485           ENDDO   ! end do (m=1,smax)
486     
487           IF (KT_TYPE_ENUM == GHD_2007) ROP_S(IJK,MMAX) = SUM_ROPS
488     
489     ! this section must be skipped until after the initial setup of the
490     ! discrete element portion of the simulation (set_bc1 is called once
491     ! before the initial dem setup).
492           IF (DISCRETE_ELEMENT .AND. .NOT.FIRST_PASS) THEN
493              FIRST_PASS = .FALSE.
494              DO M = MMAX+1, DES_MMAX+MMAX
495     ! unlike in the two fluid model, in the discrete element model it is
496     ! possible to actually calculate the bulk density in a flow boundary
497     ! cell. Currently, however, such calculations are not strictly enforced.
498     ! therefore use the bulk density of the adjacent fluid cell
499                 ROP_S(IJK,M) = ROP_S(FIJK,M)
500                 SUM_ROPS = SUM_ROPS + ROP_S(IJK,M)
501                 SUM_EPS = SUM_EPS + EP_S(IJK,M)
502              ENDDO
503           ENDIF
504     
505     ! if bc_ep_g undefined, set ep_g accordingly (based on flow condition
506     ! or based on bc_rop_s). if bc_ep_g is defined its set value will be
507     ! maintained (from set_bc0).
508           IF (BC_EP_G(BCV) == UNDEFINED) EP_G(IJK) = ONE - SUM_EPS
509     
510     ! now that ep_g in the boundary cell is known, define the bulk density
511     ! of the gas phase in the boundary cell
512           ROP_G(IJK) = RO_G(IJK)*EP_G(IJK)
513     
514           RETURN
515           END SUBROUTINE SET_OUTFLOW_EP
516     
517     
518     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
519     !                                                                      C
520     !  Purpose: Update convective fluxes....                               C
521     !  Set the value of the convective fluxes in the specified boundary    C
522     !  cell according to their value in the adjacent fluid cell.           C
523     !                                                                      C
524     !  Comment/concern:                                                    C
525     !  Should these be assigned in the same method as the velocity? Note   C
526     !  if bc_plane is W, S, B then the normal component of velocity may be C
527     !  assigned a zero value as opposed to value of its neighboring fluid  C
528     !  cell. This routine would seem to introduce some inconsistency       C
529     !  between velocity and flux at boundary.                              C
530     !                                                                      C
531     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
532     
533           SUBROUTINE SET_OUTFLOW_FLUXES(IJK,FIJK)
534     
535     ! Global variables
536     !---------------------------------------------------------------------//
537           use physprop, only: mmax
538           use run, only: kt_type_enum, ghd_2007
539           use run, only: added_mass
540           use mflux, only: flux_ge, flux_gn, flux_gt
541           use mflux, only: flux_se, flux_sn, flux_st
542           use mflux, only: flux_ne, flux_nn, flux_nt
543           use mflux, only: flux_gse, flux_gsn, flux_gst
544           use mflux, only: flux_sse, flux_ssn, flux_sst
545     
546     ! Dummy arguments
547     !---------------------------------------------------------------------//
548     ! ijk index for boundary cell
549           INTEGER, INTENT(IN) :: IJK
550     ! ijk index for adjacent fluid cell
551           INTEGER, INTENT(IN) :: FIJK
552     !---------------------------------------------------------------------//
553     
554           Flux_gE(IJK) = Flux_gE(FIJK)
555           Flux_gN(IJK) = Flux_gN(FIJK)
556           Flux_gT(IJK) = Flux_gT(FIJK)
557           IF (MMAX >0) THEN
558              Flux_sE(IJK,:MMAX) = Flux_sE(FIJK,:MMAX)
559              Flux_sN(IJK,:MMAX) = Flux_sN(FIJK,:MMAX)
560              Flux_sT(IJK,:MMAX) = Flux_sT(FIJK,:MMAX)
561           ENDIF
562     
563           IF(ADDED_MASS) THEN
564             Flux_gSE(IJK) = Flux_gSE(FIJK)
565             Flux_gSN(IJK) = Flux_gSN(FIJK)
566             Flux_gST(IJK) = Flux_gST(FIJK)
567             IF (MMAX >0) THEN
568                Flux_sSE(IJK) = Flux_sSE(FIJK)
569                Flux_sSN(IJK) = Flux_sSN(FIJK)
570                Flux_sST(IJK) = Flux_sST(FIJK)
571              ENDIF
572           ENDIF
573     
574           IF (KT_TYPE_ENUM == GHD_2007) THEN
575              Flux_nE(IJK) = Flux_nE(FIJK)
576              Flux_nN(IJK) = Flux_nN(FIJK)
577              Flux_nT(IJK) = Flux_nT(FIJK)
578           ENDIF
579     
580           RETURN
581           END SUBROUTINE SET_OUTFLOW_FLUXES
582     
583     
584     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
585     !                                                                      C
586     ! Purpose: Specify the field variable in the bc cell depending on the  C
587     ! flow direction. If flow is into domain apply user specified BC. If   C
588     ! the flow is out of the domain follow outflow type approach wherein   C
589     ! the value of the adjacent fluid cell is applied to the bc cell.      C
590     !                                                                      C
591     ! WARNING: this routine only serves to specify the field variables     C
592     ! whose governing solver routine invokes bc_phi! do not insert any     C
593     ! other settings here (such as derived quantities) as it is likely     C
594     ! not appropriate!!                                                    C
595     !                                                                      C
596     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
597     
598           SUBROUTINE SET_PINOUTFLOW(BCV,IJK,FIJK,RVEL_G,RVEL_S)
599     
600     ! Global variables
601     !---------------------------------------------------------------------//
602           use bc, only: bc_t_g, bc_x_g
603           use bc, only: bc_scalar, bc_k_turb_g, bc_e_turb_g
604           use bc, only: bc_t_s, bc_x_s, bc_theta_m
605           use fldvar, only: t_g, t_s, theta_m
606           use fldvar, only: x_g, x_s, scalar
607           use fldvar, only: k_turb_g, e_turb_g
608     ! needed because of ghd theory
609           use fldvar, only: rop_s, ro_s, d_p
610           use run, only: kt_type_enum, ghd_2007
611           use run, only: k_epsilon
612           use physprop, only: nmax, smax, mmax
613           use scalars, only: nscalar, phase4scalar
614     
615     ! Global parameters
616     !---------------------------------------------------------------------//
617           use param, only: dimension_m
618           use param1, only: zero
619           use constant, only: pi
620           IMPLICIT NONE
621     
622     ! Dummy arguments
623     !---------------------------------------------------------------------//
624     ! Boundary condition number
625           INTEGER, INTENT(IN) :: BCV
626     ! index for boundary cell
627           INTEGER, INTENT(IN) :: IJK
628     ! index for a fluid cell adjacent to the boundary cell
629           INTEGER, INTENT(IN) :: FIJK
630     ! gas and solids velocity defined to be positive when flow is leaving
631     ! the domain (outflow!)
632           DOUBLE PRECISION, INTENT(IN) :: RVEL_G, RVEL_S(DIMENSION_M)
633     
634     ! Local variables
635     !---------------------------------------------------------------------//
636     ! number density
637           DOUBLE PRECISION :: Nm, NTot
638     ! indices
639           INTEGER :: M, N, Ms
640     !---------------------------------------------------------------------//
641     
642     ! address gas phase
643           IF (RVEL_G < 0) THEN   ! inflow (apply BC)
644              T_G(IJK) = BC_T_G(BCV)
645              IF (NMAX(0) > 0) &
646                X_G(IJK,:NMAX(0)) = BC_X_G(BCV,:NMAX(0))
647              IF (K_Epsilon) THEN
648                K_Turb_G(IJK) = BC_K_Turb_G(BCV)
649                E_Turb_G(IJK) = BC_E_Turb_G(BCV)
650              ENDIF
651              IF (NScalar > 0) THEN
652                 DO N = 1,NScalar
653                    Ms = Phase4Scalar(N)
654                    IF (Ms == 0) &
655                       Scalar(IJK, N) = BC_Scalar(BCV,N)
656                 ENDDO
657              ENDIF
658           ELSE   ! outflow (apply neighbor fluid)
659              T_G(IJK) = T_G(FIJK)
660              IF (NMAX(0) > 0) &
661                 X_G(IJK,:NMAX(0)) = X_G(FIJK,:NMAX(0))
662              IF(K_Epsilon) THEN
663                 K_Turb_G(IJK) = K_Turb_G(FIJK)
664                 E_Turb_G(IJK) = E_Turb_G(FIJK)
665              ENDIF
666              IF (NScalar >0) THEN
667                 DO N = 1, NScalar
668                    Ms = Phase4Scalar(N)
669                    IF (Ms == 0) &
670                       Scalar(IJK, N) = Scalar(FIJK, N)
671                 ENDDO
672              ENDIF
673           ENDIF   ! end gas phase
674     
675     ! address solids phases
676           DO M = 1, SMAX
677              IF (RVEL_S(M) < 0) THEN   ! inflow (apply BC)
678                 T_S(IJK,M) = BC_T_S(BCV,M)
679                 THETA_M(IJK,M) = BC_THETA_M(BCV,M)
680                 IF (NMAX(M) > 0) &
681                    X_S(IJK,M,:NMAX(M)) = BC_X_S(BCV,M,:NMAX(M))
682                 IF (NScalar > 0) THEN
683                    DO N = 1,NScalar
684                       Ms = Phase4Scalar(N)
685                       IF (Ms == M) &
686                          Scalar(IJK, N) = BC_Scalar(BCV,N)
687                    ENDDO
688                 ENDIF
689              ELSE   ! outflow (apply neighbor fluid)
690                 T_S(IJK,M) = T_S(FIJK,M)
691                 THETA_M(IJK,M) =  THETA_M(FIJK,M)
692                 IF (NMAX(M) > 0) &
693                    X_S(IJK,M,:NMAX(M)) = X_S(FIJK,M,:NMAX(M))
694                 IF (NScalar > 0) THEN
695                    DO N = 1,NScalar
696                       Ms = Phase4Scalar(N)
697                       IF (Ms == M) &
698                          Scalar(IJK, N) = Scalar(FIJK,N)
699                    ENDDO
700                 ENDIF
701              ENDIF
702           ENDDO
703     
704     
705     ! derived BC field quantity (since it has no explicit user defined BC)
706           IF(KT_TYPE_ENUM == GHD_2007) THEN
707              nTOT = zero
708              DO M = 1, SMAX
709                  nM = ROP_S(IJK,M)*6d0/ &
710                      (PI*D_p(IJK,M)**3*RO_S(IJK,M))
711                  nTOT = nTOT + nM
712                  THETA_M(IJK,MMAX) = THETA_M(IJK,MMAX) + &
713                     nM*THETA_M(IJK,M)
714              ENDDO
715              IF (NTOT > zero) &
716                 THETA_M(IJK,MMAX) = THETA_M(IJK,MMAX) / nTOT
717           ENDIF
718     
719           RETURN
720           END SUBROUTINE SET_PINOUTFLOW
721