set_outflow.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: SET_OUTFLOW                                             C
!  Purpose: Set specified outflow bc for pressure outflow,             C
!  mass outflow, outflow and now also pressure inflow bc               C
!                                                                      C
!  Author: M. Syamlal                                 Date: 21-JAN-92  C
!                                                                      C
!  Comments:                                                           C
!  If the outflow boundary is on the W, S or B side of the domain and  C
!  the component of velocity through the plane is defined then this    C
!  routine will NOT modify it (i.e., its value from the momentum       C
!  solver is maintained).                                              C
!                                                                      C
!  In general it would seem this routine does little more than what    C
!  is already done within the momentum routines in terms of velocity.  C
!  The normal component is either 1) untouched if the outflow is on    C
!  the W, S, B side of the domain or 2) is set to value of the         C
!  adjacent fluid cell if it is on the E, N, T side (similar to the    C
!  respective momentum bc routines). The primary addition here is      C
!  that the tangential components of a bc cell are set to that of      C
!  the adjacent fluid cell. Note the tangential components are not     C
!  explicitly handled in the momentum _BC_ routines; instead their     C
!  values are based on solution of the momentum equation which is      C
!  replaced here                                                       C
!                                                                      C
!  Several routines are called which perform the following tasks:      C
!  set_outflow_misc - several derived quantities are set in the        C
!      boundary                                                        C
!  set_outflow_ep - the void/volume fraction and bulk densities are    C
!      set in the boundary                                             C
!  set_outflow_fluxes - convective fluxes are set in the boundary      C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE set_outflow(BCV)

! Modules
!---------------------------------------------------------------------//
      use bc
      use bc, only: bc_k_b, bc_k_t
      use bc, only: bc_j_s, bc_j_n
      use bc, only: bc_i_w, bc_i_e

      use fldvar, only: rop_g, rop_s
      use fldvar, only: u_g, v_g, w_g
      use fldvar, only: u_s, v_s, w_s
      use physprop, only: mmax

      use functions, only: im_of, ip_of, jm_of, jp_of, km_of, kp_of
      use functions, only: fluid_at
      use functions, only: is_on_mype_plus2layers
      use functions, only: funijk
      use compar, only: dead_cell_at

      use param, only: dimension_m
      use param1, only: undefined, zero
      IMPLICIT NONE

! Dummy arguments
!---------------------------------------------------------------------//
! Boundary condition number
      INTEGER, INTENT(IN) :: BCV

! Local variables
!---------------------------------------------------------------------//
! indices
      INTEGER :: I, J, K, M
! index for boundary cell
      INTEGER :: IJK
! index for a fluid cell adjacent to the boundary cell
      INTEGER :: FIJK
! local value for normal component of gas and solids velocity defined
! such that
      DOUBLE PRECISION :: RVEL_G, RVEL_S(dimension_m)
!---------------------------------------------------------------------//

! Loop over the range of boundary cells
      DO k = bc_k_b(bcv), bc_k_t(bcv)
         DO j = bc_j_s(bcv), bc_j_n(bcv)
            DO i = bc_i_w(bcv), bc_i_e(bcv)
! Check if current i,j,k resides on this PE
               IF (.NOT.is_on_mype_plus2layers(i,j,k)) cycle
               IF(dead_cell_at(i,j,k)) cycle
               ijk = funijk(i,j,k)

! Fluid cell at West
! --------------------------------------------------------------------//
               IF (fluid_at(im_of(ijk))) THEN
                  fijk = im_of(ijk)
                  rvel_g = u_g(fijk)
                  IF (mmax>0) rvel_s(:mmax) = u_s(fijk,:mmax)

                  CALL set_outflow_misc(bcv, ijk, fijk)
                  CALL set_outflow_ep(bcv, ijk, fijk, rvel_g, rvel_s)
                  IF (bc_type_enum(bcv)==p_inflow) &
                     CALL set_pinoutflow(bcv, ijk, fijk, rvel_g, rvel_s)

! Set the boundary cell value of the normal component of velocity
! according to the value in the adjacent fluid cell. Note the value
! of the boundary velocity is a scaled version of the value of the
! adjacent fluid velocity based on the concentration ratio of the fluid
! cell to the boundary cell.
! - For the gas phase, this ratio is most likely 1 except for
!   compressible cases with a PO/PI boundary where P_g of the boundary
!   is set and may differ from the value of the adjacent fluid cell.
! - For the solids phase this seems unnecessary..? Differences may arise
!   if bc_rop_s is set...
                  IF (rop_g(ijk) > zero) THEN
                    u_g(ijk) = rop_g(fijk)*u_g(fijk)/rop_g(ijk)
                  ELSE
                     u_g(ijk) = zero
                  ENDIF

! the tangential components are not explicitly handled in the boundary
! condition routines of the corresponding momentum equation
                  v_g(ijk) = v_g(fijk)
                  w_g(ijk) = w_g(fijk)

                  IF (mmax > 0) THEN
                     WHERE (rop_s(ijk,:mmax) > zero)
                       u_s(ijk,:mmax) = rop_s(fijk,:mmax)*&
                           u_s(fijk,:mmax)/rop_s(ijk,:mmax)
                     ELSEWHERE
                        u_s(ijk,:mmax) = zero
                     END WHERE
                     v_s(ijk,:mmax) = v_s(fijk,:mmax)
                     w_s(ijk,:mmax) = w_s(fijk,:mmax)
                  ENDIF

                  CALL set_outflow_fluxes(ijk, fijk)
               ENDIF   ! end if (fluid_at(im_of(ijk)))


! Fluid cell at East
! --------------------------------------------------------------------//
               IF (fluid_at(ip_of(ijk))) THEN
                  fijk = ip_of(ijk)
! define normal component such that it is positive when exiting the
! domain
                  rvel_g = -u_g(ijk)
                  IF (mmax >0) rvel_s(:mmax) = -u_s(ijk,:mmax)

                  CALL set_outflow_misc(bcv, ijk, fijk)
                  CALL set_outflow_ep(bcv, ijk, fijk, rvel_g, rvel_s)
                  IF (bc_type_enum(bcv)==p_inflow) &
                     CALL set_pinoutflow(bcv, ijk, fijk, rvel_g, rvel_s)

! provide an initial value for the velocity component through the domain
! otherwise its present value (from solution of the corresponding
! momentum eqn) is kept. values for the velocity components in the off
! directions are modified (needed for PO or O boundaries but not MO or
! PI as velocities should be fully specified by this point)
                  IF (u_g(ijk) == undefined) THEN
                     IF (rop_g(ijk) > zero) THEN
                        u_g(ijk) = rop_g(fijk)*u_g(fijk)/rop_g(ijk)
                     ELSE
                        u_g(ijk) = zero
                     ENDIF
                  ENDIF
                  v_g(ijk) = v_g(fijk)
                  w_g(ijk) = w_g(fijk)

                  DO m = 1, mmax
                     IF (u_s(ijk,m) == undefined) THEN
                        IF (rop_s(ijk,m) > zero) THEN
                           u_s(ijk,m) = rop_s(fijk,m)*&
                              u_s(fijk,m)/rop_s(ijk,m)
                        ELSE
                           u_s(ijk,m) = zero
                        ENDIF
                     ENDIF
                     v_s(ijk,m) = v_s(fijk,m)
                     w_s(ijk,m) = w_s(fijk,m)
                  ENDDO

                  CALL set_outflow_fluxes(ijk, fijk)
               ENDIF   ! end if (fluid_at(ip_of(ijk)))


! Fluid cell at South
! --------------------------------------------------------------------//
               IF (fluid_at(jm_of(ijk))) THEN
                  fijk = jm_of(ijk)
                  rvel_g = v_g(fijk)
                  IF(mmax>0) rvel_s(:mmax) = v_s(fijk,:mmax)

                  CALL set_outflow_misc(bcv, ijk, fijk)
                  CALL set_outflow_ep(bcv, ijk, fijk, rvel_g, rvel_s)
                  IF (bc_type_enum(bcv)==p_inflow) &
                     CALL set_pinoutflow(bcv, ijk, fijk, rvel_g, rvel_s)

                  IF (rop_g(ijk) > zero) THEN
                     v_g(ijk) = rop_g(fijk)*v_g(fijk)/rop_g(ijk)
                  ELSE
                     v_g(ijk) = zero
                  ENDIF
                  u_g(ijk) = u_g(fijk)
                  w_g(ijk) = w_g(fijk)

                  IF (mmax > 0) THEN
                      WHERE (rop_s(ijk,:mmax) > zero)
                         v_s(ijk,:mmax) = rop_s(fijk,:mmax)*&
                            v_s(fijk,:mmax)/rop_s(ijk,:mmax)
                      ELSEWHERE
                         v_s(ijk,:mmax) = zero
                      END WHERE
                      u_s(ijk,:mmax) = u_s(fijk,:mmax)
                      w_s(ijk,:mmax) = w_s(fijk,:mmax)
                  ENDIF

                  CALL set_outflow_fluxes(ijk, fijk)
               ENDIF   ! end if (fluid_at(jm_of(ijk)))


! Fluid cell at North
! --------------------------------------------------------------------//
               IF (fluid_at(jp_of(ijk))) THEN
                  fijk = jp_of(ijk)
                  rvel_g = -v_g(ijk)
                  IF (mmax>0) rvel_s(:mmax) = -v_s(ijk,:mmax)

                  CALL set_outflow_misc(bcv, ijk, fijk)
                  CALL set_outflow_ep(bcv, ijk, fijk, rvel_g, rvel_s)
                  IF (bc_type_enum(bcv)==p_inflow) &
                     CALL set_pinoutflow(bcv, ijk, fijk, rvel_g, rvel_s)

                  IF (v_g(ijk) == undefined) THEN
                     IF (rop_g(ijk) > zero) THEN
                        v_g(ijk) = rop_g(fijk)*v_g(fijk)/rop_g(ijk)
                     ELSE
                        v_g(ijk) = zero
                     ENDIF
                  ENDIF
                  u_g(ijk) = u_g(fijk)
                  w_g(ijk) = w_g(fijk)

                  DO m = 1, mmax
                     IF (v_s(ijk,m) == undefined) THEN
                        IF (rop_s(ijk,m) > zero) THEN
                           v_s(ijk,m) = rop_s(fijk,m)*&
                              v_s(fijk,m)/rop_s(ijk,m)
                        ELSE
                           v_s(ijk,m) = zero
                        ENDIF
                     ENDIF
                     u_s(ijk,m) = u_s(fijk,m)
                     w_s(ijk,m) = w_s(fijk,m)
                  ENDDO

                  CALL set_outflow_fluxes(ijk, fijk)
               ENDIF   ! if (fluid_at(jp_of(ijk)))


! Fluid cell at Bottom
! --------------------------------------------------------------------//
               IF (fluid_at(km_of(ijk))) THEN
                  fijk = km_of(ijk)
                  rvel_g = w_g(fijk)
                  IF (mmax>0) rvel_s(:mmax) = w_s(fijk,:mmax)

                  CALL set_outflow_misc(bcv, ijk, fijk)
                  CALL set_outflow_ep(bcv, ijk, fijk, rvel_g, rvel_s)
                  IF (bc_type_enum(bcv)==p_inflow) &
                     CALL set_pinoutflow(bcv, ijk, fijk, rvel_g, rvel_s)

                  IF (rop_g(ijk) > zero) THEN
                     w_g(ijk) = rop_g(fijk)*w_g(fijk)/rop_g(ijk)
                  ELSE
                     w_g(ijk) = zero
                  ENDIF
                  u_g(ijk) = u_g(fijk)
                  v_g(ijk) = v_g(fijk)

                  IF (mmax > 0) THEN
                    WHERE (rop_s(ijk,:mmax) > zero)
                        w_s(ijk,:mmax) = rop_s(fijk,:mmax)*&
                           w_s(fijk,:mmax)/rop_s(ijk,:mmax)
                     ELSEWHERE
                        w_s(ijk,:mmax) = zero
                     END WHERE
                     u_s(ijk,:mmax) = u_s(fijk,:mmax)
                     v_s(ijk,:mmax) = v_s(fijk,:mmax)
                  ENDIF

                  CALL set_outflow_fluxes(ijk, fijk)
               ENDIF   ! if (fluid_at(km_of(ijk)))


! Fluid cell at Top
! --------------------------------------------------------------------//
               IF (fluid_at(kp_of(ijk))) THEN
                  fijk = kp_of(ijk)
                  rvel_g = -w_g(ijk)
                  IF (mmax>0) rvel_s(:mmax) = -w_s(ijk,:mmax)

                  CALL set_outflow_misc(bcv, ijk, fijk)
                  CALL set_outflow_ep(bcv, ijk, fijk, rvel_g, rvel_s)
                  IF (bc_type_enum(bcv)==p_inflow) &
                     CALL set_pinoutflow(bcv, ijk, fijk, rvel_g, rvel_s)

                  IF (w_g(ijk) == undefined) THEN
                     IF (rop_g(ijk) > zero) THEN
                        w_g(ijk) = rop_g(fijk)*w_g(fijk)/rop_g(ijk)
                     ELSE
                        w_g(ijk) = zero
                     ENDIF
                  ENDIF
                  u_g(ijk) = u_g(fijk)
                  v_g(ijk) = v_g(fijk)

                  DO m = 1, mmax
                     IF (w_s(ijk,m) == undefined) THEN
                        IF (rop_s(ijk,m) > zero) THEN
                           w_s(ijk,m) = rop_s(fijk,m)*&
                              w_s(fijk,m)/rop_s(ijk,m)
                        ELSE
                           w_s(ijk,m) = zero
                        ENDIF
                     ENDIF
                     u_s(ijk,m) = u_s(fijk,m)
                     v_s(ijk,m) = v_s(fijk,m)
                  ENDDO

                  CALL set_outflow_fluxes(ijk, fijk)
               ENDIF   ! if (fluid_at(kp_of(ijk)))

            ENDDO   ! end do (i=i1,i2)
         ENDDO   ! end do (j=j1,j2)
      ENDDO   ! end do (k=k1,k2)

      RETURN
      END SUBROUTINE set_outflow


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: SET_OUTFLOW_MISC                                        C
!  Purpose: Set the value of certain variables in the specified        C
!  outflow boundary cell that would not otherwise be set according     C
!  to their value in the adjacent fluid cell.                          C
!                                                                      C
!  READ BEFORE MODIFYING!!                                             C
!  For many of the scalar field variables (e.g., T_g, T_s, X_g, X_s,   C
!  k_turb_g, e_turb_g, theta_m and scalar) their corresponding         C
!  governing equation solver routine sets its own value in the         C
!  boundary cell (see bc_phi). So DO NOT set their value here unless   C
!  a clear explanation to the need is also provided; more likely it    C
!  should simply be dealt with in the governing equation routine.      C
!                                                                      C
!  This routine sets the value of quantites that are derived or with   C
!  unique solver routines (e.g., P_star, P_s, MW_MIX_g, P_g) since no  C
!  other routine will.                                                 C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE set_outflow_misc(BCV, IJK, FIJK)

! Global variables
!---------------------------------------------------------------------//
      use bc
      use run, only: kt_type_enum, ghd_2007
      use fldvar, only: p_g, ro_g, t_g
      use fldvar, only: p_s, p_star
      use physprop, only: smax, mmax
      use physprop, only: ro_g0, mw_mix_g
      use eos, only: eosg

! Global parameters
!---------------------------------------------------------------------//
      use param1, only: undefined
      IMPLICIT NONE

! Dummy arguments
!---------------------------------------------------------------------//
! Boundary condition number
      INTEGER, INTENT(IN) :: BCV
! ijk index for boundary cell
      INTEGER, INTENT(IN) :: IJK
! ijk index for adjacent fluid cell
      INTEGER, INTENT(IN) :: FIJK
!---------------------------------------------------------------------//

      IF (bc_type_enum(bcv) /= p_outflow .AND. &
          bc_type_enum(bcv) /= p_inflow) p_g(ijk) = p_g(fijk)

      mw_mix_g(ijk) = mw_mix_g(fijk)

! T_g (bc_phi), P_g (above depending on bc), and MW_MIX_G (above) have
! now been defined at IJK
      IF (ro_g0 == undefined) ro_g(ijk) = &
         eosg(mw_mix_g(ijk),p_g(ijk),t_g(ijk))

      IF (smax >0) THEN
         p_star(ijk) = p_star(fijk)
         p_s(ijk,:smax) = p_s(fijk,:smax)
      ENDIF

      IF (kt_type_enum == ghd_2007 .AND. mmax>0) &
         p_s(ijk,mmax) = p_s(fijk,mmax)

      RETURN
      END SUBROUTINE set_outflow_misc


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: SET_OUTFLOW_EP                                          C
!  Purpose: Set the volume fraction/bulk density (i.e., ROP_s, EP_g    C
!  ROP_S) in the specified outflow boundary cell that would not        C
!  otherwise be set according to their value in the adjacent fluid     C
!  cell.                                                               C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE set_outflow_ep(BCV, IJK, FIJK, RVEL_G, RVEL_S)

! Global variables
!---------------------------------------------------------------------//
      use bc
      use run, only: kt_type_enum, ghd_2007
      use physprop, only: smax, mmax
      use fldvar, only: rop_g, ro_g, ep_g
      use fldvar, only: rop_s, ep_s
      use discretelement, only: discrete_element, des_mmax
! Global parameters
!---------------------------------------------------------------------//
      use param, only: dimension_m
      use param1, only: undefined, zero, one
      IMPLICIT NONE

! Dummy arguments
!---------------------------------------------------------------------//
! Boundary condition number
      INTEGER, INTENT(IN) :: BCV
! ijk index for boundary cell
      INTEGER, INTENT(IN) :: IJK
! ijk index for adjacent fluid cell
      INTEGER, INTENT(IN) :: FIJK
! the gas or solids velocity in the fluid cell adjacent to the boundary
! cell dot with the outward normal of that bc plane; defines the gas or
! solids velocity component normal to the bc plane as positive when it
! is flowing into the bc cell from the fluid cell. so, for example, for
! an outflow on the eastern boundary this is the u component of velocity
! while for an outflow on the western boundary this is the -u component,
! etc.
      DOUBLE PRECISION, INTENT(IN) :: RVEL_G
      DOUBLE PRECISION, INTENT(IN), DIMENSION(DIMENSION_M) :: RVEL_S

! Local variables
!---------------------------------------------------------------------//
! indices
      INTEGER :: M
! sum of solids phases volume fractions
      DOUBLE PRECISION :: SUM_EPs
! sum of solids phases bulk densities
      DOUBLE PRECISION :: SUM_ROPS
      LOGICAL, SAVE :: FIRST_PASS = .true.
!---------------------------------------------------------------------//

! initializing summation quantities
      sum_rops = zero
      sum_eps = zero

      DO m = 1, smax

         IF(bc_type_enum(bcv) == p_inflow) THEN
            rop_s(ijk,m) = rop_s(fijk,m)
         ELSE
! the outflow type bc do not permit 're-entering' solids, in which
! case the solids are removed
            IF (rvel_s(m) >= zero) THEN
! solids are leaving the domain
               rop_s(ijk,m) = rop_s(fijk,m)
            ELSE
! solids cannot enter the domain at an outflow cell
               rop_s(ijk,m) = zero
            ENDIF
         ENDIF

! if bc_rop_s is defined, set value of rop_s in the ijk boundary cell
! according to user definition
         IF(bc_rop_s(bcv,m)/=undefined) rop_s(ijk,m)=bc_rop_s(bcv,m)

! add to total solids phase bulk density and solids volume fraction
         sum_rops = sum_rops + rop_s(ijk,m)
         sum_eps = sum_eps + ep_s(ijk,m)
      ENDDO   ! end do (m=1,smax)

      IF (kt_type_enum == ghd_2007) rop_s(ijk,mmax) = sum_rops

! this section must be skipped until after the initial setup of the
! discrete element portion of the simulation (set_bc1 is called once
! before the initial dem setup).
      IF (discrete_element .AND. .NOT.first_pass) THEN
         first_pass = .false.
         DO m = mmax+1, des_mmax+mmax
! unlike in the two fluid model, in the discrete element model it is
! possible to actually calculate the bulk density in a flow boundary
! cell. Currently, however, such calculations are not strictly enforced.
! therefore use the bulk density of the adjacent fluid cell
            rop_s(ijk,m) = rop_s(fijk,m)
            sum_rops = sum_rops + rop_s(ijk,m)
            sum_eps = sum_eps + ep_s(ijk,m)
         ENDDO
      ENDIF

! if bc_ep_g undefined, set ep_g accordingly (based on flow condition
! or based on bc_rop_s). if bc_ep_g is defined its set value will be
! maintained (from set_bc0).
      IF (bc_ep_g(bcv) == undefined) ep_g(ijk) = one - sum_eps

! now that ep_g in the boundary cell is known, define the bulk density
! of the gas phase in the boundary cell
      rop_g(ijk) = ro_g(ijk)*ep_g(ijk)

      RETURN
      END SUBROUTINE set_outflow_ep


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Purpose: Update convective fluxes....                               C
!  Set the value of the convective fluxes in the specified boundary    C
!  cell according to their value in the adjacent fluid cell.           C
!                                                                      C
!  Comment/concern:                                                    C
!  Should these be assigned in the same method as the velocity? Note   C
!  if bc_plane is W, S, B then the normal component of velocity may be C
!  assigned a zero value as opposed to value of its neighboring fluid  C
!  cell. This routine would seem to introduce some inconsistency       C
!  between velocity and flux at boundary.                              C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C

      SUBROUTINE set_outflow_fluxes(IJK,FIJK)

! Global variables
!---------------------------------------------------------------------//
      use physprop, only: mmax
      use run, only: kt_type_enum, ghd_2007
      use run, only: added_mass
      use mflux, only: flux_ge, flux_gn, flux_gt
      use mflux, only: flux_se, flux_sn, flux_st
      use mflux, only: flux_ne, flux_nn, flux_nt
      use mflux, only: flux_gse, flux_gsn, flux_gst
      use mflux, only: flux_sse, flux_ssn, flux_sst

! Dummy arguments
!---------------------------------------------------------------------//
! ijk index for boundary cell
      INTEGER, INTENT(IN) :: IJK
! ijk index for adjacent fluid cell
      INTEGER, INTENT(IN) :: FIJK
!---------------------------------------------------------------------//

      flux_ge(ijk) = flux_ge(fijk)
      flux_gn(ijk) = flux_gn(fijk)
      flux_gt(ijk) = flux_gt(fijk)
      IF (mmax >0) THEN
         flux_se(ijk,:mmax) = flux_se(fijk,:mmax)
         flux_sn(ijk,:mmax) = flux_sn(fijk,:mmax)
         flux_st(ijk,:mmax) = flux_st(fijk,:mmax)
      ENDIF

      IF(added_mass) THEN
        flux_gse(ijk) = flux_gse(fijk)
        flux_gsn(ijk) = flux_gsn(fijk)
        flux_gst(ijk) = flux_gst(fijk)
        IF (mmax >0) THEN
           flux_sse(ijk) = flux_sse(fijk)
           flux_ssn(ijk) = flux_ssn(fijk)
           flux_sst(ijk) = flux_sst(fijk)
         ENDIF
      ENDIF

      IF (kt_type_enum == ghd_2007) THEN
         flux_ne(ijk) = flux_ne(fijk)
         flux_nn(ijk) = flux_nn(fijk)
         flux_nt(ijk) = flux_nt(fijk)
      ENDIF

      RETURN
      END SUBROUTINE set_outflow_fluxes


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
! Purpose: Specify the field variable in the bc cell depending on the  C
! flow direction. If flow is into domain apply user specified BC. If   C
! the flow is out of the domain follow outflow type approach wherein   C
! the value of the adjacent fluid cell is applied to the bc cell.      C
!                                                                      C
! WARNING: this routine only serves to specify the field variables     C
! whose governing solver routine invokes bc_phi! do not insert any     C
! other settings here (such as derived quantities) as it is likely     C
! not appropriate!!                                                    C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C

      SUBROUTINE set_pinoutflow(BCV,IJK,FIJK,RVEL_G,RVEL_S)

! Global variables
!---------------------------------------------------------------------//
      use bc, only: bc_t_g, bc_x_g
      use bc, only: bc_scalar, bc_k_turb_g, bc_e_turb_g
      use bc, only: bc_t_s, bc_x_s, bc_theta_m
      use fldvar, only: t_g, t_s, theta_m
      use fldvar, only: x_g, x_s, scalar
      use fldvar, only: k_turb_g, e_turb_g
! needed because of ghd theory
      use fldvar, only: rop_s, ro_s, d_p
      use run, only: kt_type_enum, ghd_2007
      use run, only: k_epsilon
      use physprop, only: nmax, smax, mmax
      use scalars, only: nscalar, phase4scalar

! Global parameters
!---------------------------------------------------------------------//
      use param, only: dimension_m
      use param1, only: zero
      use constant, only: pi
      IMPLICIT NONE

! Dummy arguments
!---------------------------------------------------------------------//
! Boundary condition number
      INTEGER, INTENT(IN) :: BCV
! index for boundary cell
      INTEGER, INTENT(IN) :: IJK
! index for a fluid cell adjacent to the boundary cell
      INTEGER, INTENT(IN) :: FIJK
! gas and solids velocity defined to be positive when flow is leaving
! the domain (outflow!)
      DOUBLE PRECISION, INTENT(IN) :: RVEL_G, RVEL_S(dimension_m)

! Local variables
!---------------------------------------------------------------------//
! number density
      DOUBLE PRECISION :: Nm, NTot
! indices
      INTEGER :: M, N, Ms
!---------------------------------------------------------------------//

! address gas phase
      IF (rvel_g < 0) THEN   ! inflow (apply BC)
         t_g(ijk) = bc_t_g(bcv)
         IF (nmax(0) > 0) &
           x_g(ijk,:nmax(0)) = bc_x_g(bcv,:nmax(0))
         IF (k_epsilon) THEN
           k_turb_g(ijk) = bc_k_turb_g(bcv)
           e_turb_g(ijk) = bc_e_turb_g(bcv)
         ENDIF
         IF (nscalar > 0) THEN
            DO n = 1,nscalar
               ms = phase4scalar(n)
               IF (ms == 0) &
                  scalar(ijk, n) = bc_scalar(bcv,n)
            ENDDO
         ENDIF
      ELSE   ! outflow (apply neighbor fluid)
         t_g(ijk) = t_g(fijk)
         IF (nmax(0) > 0) &
            x_g(ijk,:nmax(0)) = x_g(fijk,:nmax(0))
         IF(k_epsilon) THEN
            k_turb_g(ijk) = k_turb_g(fijk)
            e_turb_g(ijk) = e_turb_g(fijk)
         ENDIF
         IF (nscalar >0) THEN
            DO n = 1, nscalar
               ms = phase4scalar(n)
               IF (ms == 0) &
                  scalar(ijk, n) = scalar(fijk, n)
            ENDDO
         ENDIF
      ENDIF   ! end gas phase

! address solids phases
      DO m = 1, smax
         IF (rvel_s(m) < 0) THEN   ! inflow (apply BC)
            t_s(ijk,m) = bc_t_s(bcv,m)
            theta_m(ijk,m) = bc_theta_m(bcv,m)
            IF (nmax(m) > 0) &
               x_s(ijk,m,:nmax(m)) = bc_x_s(bcv,m,:nmax(m))
            IF (nscalar > 0) THEN
               DO n = 1,nscalar
                  ms = phase4scalar(n)
                  IF (ms == m) &
                     scalar(ijk, n) = bc_scalar(bcv,n)
               ENDDO
            ENDIF
         ELSE   ! outflow (apply neighbor fluid)
            t_s(ijk,m) = t_s(fijk,m)
            theta_m(ijk,m) =  theta_m(fijk,m)
            IF (nmax(m) > 0) &
               x_s(ijk,m,:nmax(m)) = x_s(fijk,m,:nmax(m))
            IF (nscalar > 0) THEN
               DO n = 1,nscalar
                  ms = phase4scalar(n)
                  IF (ms == m) &
                     scalar(ijk, n) = scalar(fijk,n)
               ENDDO
            ENDIF
         ENDIF
      ENDDO


! derived BC field quantity (since it has no explicit user defined BC)
      IF(kt_type_enum == ghd_2007) THEN
         ntot = zero
         DO m = 1, smax
             nm = rop_s(ijk,m)*6d0/ &
                 (pi*d_p(ijk,m)**3*ro_s(ijk,m))
             ntot = ntot + nm
             theta_m(ijk,mmax) = theta_m(ijk,mmax) + &
                nm*theta_m(ijk,m)
         ENDDO
         IF (ntot > zero) &
            theta_m(ijk,mmax) = theta_m(ijk,mmax) / ntot
      ENDIF

      RETURN
      END SUBROUTINE set_pinoutflow