File: RELATIVE:/../../../mfix.git/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, only: bc_type
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(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(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(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(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(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(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, only: bc_type
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(BCV) /= 'P_OUTFLOW' .AND. &
383               BC_TYPE(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, only: bc_type, bc_ep_g, bc_rop_s
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           use discretelement, only: des_rop_s, des_ro_s
424     ! Global parameters
425     !---------------------------------------------------------------------//
426           use param, only: dimension_m
427           use param1, only: undefined, zero, one
428           IMPLICIT NONE
429     
430     ! Dummy arguments
431     !---------------------------------------------------------------------//
432     ! Boundary condition number
433           INTEGER, INTENT(IN) :: BCV
434     ! ijk index for boundary cell
435           INTEGER, INTENT(IN) :: IJK
436     ! ijk index for adjacent fluid cell
437           INTEGER, INTENT(IN) :: FIJK
438     ! the gas or solids velocity in the fluid cell adjacent to the boundary
439     ! cell dot with the outward normal of that bc plane; defines the gas or
440     ! solids velocity component normal to the bc plane as positive when it
441     ! is flowing into the bc cell from the fluid cell. so, for example, for
442     ! an outflow on the eastern boundary this is the u component of velocity
443     ! while for an outflow on the western boundary this is the -u component,
444     ! etc.
445           DOUBLE PRECISION, INTENT(IN) :: RVEL_G
446           DOUBLE PRECISION, INTENT(IN), DIMENSION(DIMENSION_M) :: RVEL_S
447     
448     ! Local variables
449     !---------------------------------------------------------------------//
450     ! indices
451           INTEGER :: M
452     ! solids volume fraction
453           DOUBLE PRECISION :: EPs
454     ! sum of solids phases volume fractions
455           DOUBLE PRECISION :: SUM_EPs
456     ! sum of solids phases bulk densities
457           DOUBLE PRECISION :: SUM_ROPS
458     !---------------------------------------------------------------------//
459     
460     ! initializing summation quantities
461           SUM_ROPS = ZERO
462           SUM_EPS = ZERO
463     
464           DO M = 1, SMAX
465     
466              IF(BC_TYPE(BCV) == 'P_INFLOW') THEN
467                 ROP_S(IJK,M) = ROP_S(FIJK,M)
468              ELSE
469     ! the outflow type bc do not permit 're-entering' solids, in which
470     ! case the solids are removed
471                 IF (RVEL_S(M) >= ZERO) THEN
472     ! solids are leaving the domain
473                    ROP_S(IJK,M) = ROP_S(FIJK,M)
474                 ELSE
475     ! solids cannot enter the domain at an outflow cell
476                    ROP_S(IJK,M) = ZERO
477                 ENDIF
478              ENDIF
479     
480     ! if bc_rop_s is defined, set value of rop_s in the ijk boundary cell
481     ! according to user definition
482              IF(BC_ROP_S(BCV,M)/=UNDEFINED) ROP_S(IJK,M)=BC_ROP_S(BCV,M)
483     
484     ! add to total solids phase bulk density and solids volume fraction
485              SUM_ROPS = SUM_ROPS + ROP_S(IJK,M)
486              SUM_EPS = SUM_EPS + EP_S(IJK,M)
487           ENDDO   ! end do (m=1,smax)
488     
489           IF (KT_TYPE_ENUM == GHD_2007) ROP_S(IJK,MMAX) = SUM_ROPS
490     
491     ! this section must be skipped until after the initial setup of the
492     ! discrete element portion of the simulation (set_bc1 is called once
493     ! before the initial setup).
494           IF (DISCRETE_ELEMENT .AND. ALLOCATED(DES_ROP_S)) THEN
495              DO M = 1, DES_MMAX
496     ! unlike in the two fluid model, in the discrete element model it is
497     ! possible to actually calculate the bulk density in a flow boundary
498     ! cell. Currently, however, such calculations are not strictly enforced.
499     ! therefore use the bulk density of the adjacent fluid cell
500                 DES_ROP_S(IJK,M) = DES_ROP_S(FIJK,M)
501                 SUM_ROPS = SUM_ROPS + DES_ROP_S(IJK,M)
502                 EPS = DES_ROP_S(IJK,M)/DES_RO_S(M)
503                 SUM_EPS = SUM_EPS + EPS
504              ENDDO
505           ENDIF
506     
507     ! if bc_ep_g undefined, set ep_g accordingly (based on flow condition
508     ! or based on bc_rop_s). if bc_ep_g is defined its set value will be
509     ! maintained (from set_bc0).
510           IF (BC_EP_G(BCV) == UNDEFINED) EP_G(IJK) = ONE - SUM_EPS
511     
512     ! now that ep_g in the boundary cell is known, define the bulk density
513     ! of the gas phase in the boundary cell
514           ROP_G(IJK) = RO_G(IJK)*EP_G(IJK)
515     
516           RETURN
517           END SUBROUTINE SET_OUTFLOW_EP
518     
519     
520     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
521     !                                                                      C
522     !  Purpose: Update convective fluxes....                               C
523     !  Set the value of the convective fluxes in the specified boundary    C
524     !  cell according to their value in the adjacent fluid cell.           C
525     !                                                                      C
526     !  Comment/concern:                                                    C
527     !  Should these be assigned in the same method as the velocity? Note   C
528     !  if bc_plane is W, S, B then the normal component of velocity may be C
529     !  assigned a zero value as opposed to value of its neighboring fluid  C
530     !  cell. This routine would seem to introduce some inconsistency       C
531     !  between velocity and flux at boundary.                              C
532     !                                                                      C
533     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
534     
535           SUBROUTINE SET_OUTFLOW_FLUXES(IJK,FIJK)
536     
537     ! Global variables
538     !---------------------------------------------------------------------//
539           use physprop, only: mmax
540           use run, only: kt_type_enum, ghd_2007
541           use run, only: added_mass
542           use mflux, only: flux_ge, flux_gn, flux_gt
543           use mflux, only: flux_se, flux_sn, flux_st
544           use mflux, only: flux_ne, flux_nn, flux_nt
545           use mflux, only: flux_gse, flux_gsn, flux_gst
546           use mflux, only: flux_sse, flux_ssn, flux_sst
547     
548     ! Dummy arguments
549     !---------------------------------------------------------------------//
550     ! ijk index for boundary cell
551           INTEGER, INTENT(IN) :: IJK
552     ! ijk index for adjacent fluid cell
553           INTEGER, INTENT(IN) :: FIJK
554     !---------------------------------------------------------------------//
555     
556           Flux_gE(IJK) = Flux_gE(FIJK)
557           Flux_gN(IJK) = Flux_gN(FIJK)
558           Flux_gT(IJK) = Flux_gT(FIJK)
559           IF (MMAX >0) THEN
560              Flux_sE(IJK,:MMAX) = Flux_sE(FIJK,:MMAX)
561              Flux_sN(IJK,:MMAX) = Flux_sN(FIJK,:MMAX)
562              Flux_sT(IJK,:MMAX) = Flux_sT(FIJK,:MMAX)
563           ENDIF
564     
565           IF(ADDED_MASS) THEN
566             Flux_gSE(IJK) = Flux_gSE(FIJK)
567             Flux_gSN(IJK) = Flux_gSN(FIJK)
568             Flux_gST(IJK) = Flux_gST(FIJK)
569             IF (MMAX >0) THEN
570                Flux_sSE(IJK) = Flux_sSE(FIJK)
571                Flux_sSN(IJK) = Flux_sSN(FIJK)
572                Flux_sST(IJK) = Flux_sST(FIJK)
573              ENDIF
574           ENDIF
575     
576           IF (KT_TYPE_ENUM == GHD_2007) THEN
577              Flux_nE(IJK) = Flux_nE(FIJK)
578              Flux_nN(IJK) = Flux_nN(FIJK)
579              Flux_nT(IJK) = Flux_nT(FIJK)
580           ENDIF
581     
582           RETURN
583           END SUBROUTINE SET_OUTFLOW_FLUXES
584     
585     
586     !vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
587     !                                                                      C
588     ! Purpose: Specify the field variable in the bc cell depending on the  C
589     ! flow direction. If flow is into domain apply user specified BC. If   C
590     ! the flow is out of the domain follow outflow type approach wherein   C
591     ! the value of the adjacent fluid cell is applied to the bc cell.      C
592     !                                                                      C
593     ! WARNING: this routine only serves to specify the field variables     C
594     ! whose governing solver routine invokes bc_phi! do not insert any     C
595     ! other settings here (such as derived quantities) as it is likely     C
596     ! not appropriate!!                                                    C
597     !                                                                      C
598     !^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
599     
600           SUBROUTINE SET_PINOUTFLOW(BCV,IJK,FIJK,RVEL_G,RVEL_S)
601     
602     ! Global variables
603     !---------------------------------------------------------------------//
604           use bc, only: bc_t_g, bc_x_g
605           use bc, only: bc_scalar, bc_k_turb_g, bc_e_turb_g
606           use bc, only: bc_t_s, bc_x_s, bc_theta_m
607           use fldvar, only: t_g, t_s, theta_m
608           use fldvar, only: x_g, x_s, scalar
609           use fldvar, only: k_turb_g, e_turb_g
610     ! needed because of ghd theory
611           use fldvar, only: rop_s, ro_s, d_p
612           use run, only: kt_type_enum, ghd_2007
613           use run, only: k_epsilon
614           use physprop, only: nmax, smax, mmax
615           use scalars, only: nscalar, phase4scalar
616     
617     ! Global parameters
618     !---------------------------------------------------------------------//
619           use param, only: dimension_m
620           use param1, only: zero
621           use constant, only: pi
622           IMPLICIT NONE
623     
624     ! Dummy arguments
625     !---------------------------------------------------------------------//
626     ! Boundary condition number
627           INTEGER, INTENT(IN) :: BCV
628     ! index for boundary cell
629           INTEGER, INTENT(IN) :: IJK
630     ! index for a fluid cell adjacent to the boundary cell
631           INTEGER, INTENT(IN) :: FIJK
632     ! gas and solids velocity defined to be positive when flow is leaving
633     ! the domain (outflow!)
634           DOUBLE PRECISION, INTENT(IN) :: RVEL_G, RVEL_S(DIMENSION_M)
635     
636     ! Local variables
637     !---------------------------------------------------------------------//
638     ! number density
639           DOUBLE PRECISION :: Nm, NTot
640     ! indices
641           INTEGER :: M, N, Ms
642     !---------------------------------------------------------------------//
643     
644     ! address gas phase
645           IF (RVEL_G < 0) THEN   ! inflow (apply BC)
646              T_G(IJK) = BC_T_G(BCV)
647              IF (NMAX(0) > 0) &
648                X_G(IJK,:NMAX(0)) = BC_X_G(BCV,:NMAX(0))
649              IF (K_Epsilon) THEN
650                K_Turb_G(IJK) = BC_K_Turb_G(BCV)
651                E_Turb_G(IJK) = BC_E_Turb_G(BCV)
652              ENDIF
653              IF (NScalar > 0) THEN
654                 DO N = 1,NScalar
655                    Ms = Phase4Scalar(N)
656                    IF (Ms == 0) &
657                       Scalar(IJK, N) = BC_Scalar(BCV,N)
658                 ENDDO
659              ENDIF
660           ELSE   ! outflow (apply neighbor fluid)
661              T_G(IJK) = T_G(FIJK)
662              IF (NMAX(0) > 0) &
663                 X_G(IJK,:NMAX(0)) = X_G(FIJK,:NMAX(0))
664              IF(K_Epsilon) THEN
665                 K_Turb_G(IJK) = K_Turb_G(FIJK)
666                 E_Turb_G(IJK) = E_Turb_G(FIJK)
667              ENDIF
668              IF (NScalar >0) THEN
669                 DO N = 1, NScalar
670                    Ms = Phase4Scalar(N)
671                    IF (Ms == 0) &
672                       Scalar(IJK, N) = Scalar(FIJK, N)
673                 ENDDO
674              ENDIF
675           ENDIF   ! end gas phase
676     
677     ! address solids phases
678           DO M = 1, SMAX
679              IF (RVEL_S(M) < 0) THEN   ! inflow (apply BC)
680                 T_S(IJK,M) = BC_T_S(BCV,M)
681                 THETA_M(IJK,M) = BC_THETA_M(BCV,M)
682                 IF (NMAX(M) > 0) &
683                    X_S(IJK,M,:NMAX(M)) = BC_X_S(BCV,M,:NMAX(M))
684                 IF (NScalar > 0) THEN
685                    DO N = 1,NScalar
686                       Ms = Phase4Scalar(N)
687                       IF (Ms == M) &
688                          Scalar(IJK, N) = BC_Scalar(BCV,N)
689                    ENDDO
690                 ENDIF
691              ELSE   ! outflow (apply neighbor fluid)
692                 T_S(IJK,M) = T_S(FIJK,M)
693                 THETA_M(IJK,M) =  THETA_M(FIJK,M)
694                 IF (NMAX(M) > 0) &
695                    X_S(IJK,M,:NMAX(M)) = X_S(FIJK,M,:NMAX(M))
696                 IF (NScalar > 0) THEN
697                    DO N = 1,NScalar
698                       Ms = Phase4Scalar(N)
699                       IF (Ms == M) &
700                          Scalar(IJK, N) = Scalar(FIJK,N)
701                    ENDDO
702                 ENDIF
703              ENDIF
704           ENDDO
705     
706     
707     ! derived BC field quantity (since it has no explicit user defined BC)
708           IF(KT_TYPE_ENUM == GHD_2007) THEN
709              nTOT = zero
710              DO M = 1, SMAX
711                  nM = ROP_S(IJK,M)*6d0/ &
712                      (PI*D_p(IJK,M)**3*RO_S(IJK,M))
713                  nTOT = nTOT + nM
714                  THETA_M(IJK,MMAX) = THETA_M(IJK,MMAX) + &
715                     nM*THETA_M(IJK,M)
716              ENDDO
717              IF (NTOT > zero) &
718                 THETA_M(IJK,MMAX) = THETA_M(IJK,MMAX) / nTOT
719           ENDIF
720     
721           RETURN
722           END SUBROUTINE SET_PINOUTFLOW
723