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