set_bc0.f

Go to the documentation of this file.

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: set_bc0                                                 C
!  Purpose: This subroutine does the initial setting of all boundary   C
!  conditions. The user specifications of the boundary conditions are  C
!  checked for veracity in various check_data/ routines:               C
!  (e.g., check_boundary_conditions).                                  C
!                                                                      C
!  Author: M. Syamlal                                 Date: 29-JAN-92  C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE set_bc0

! Modules
!--------------------------------------------------------------------//
      use bc, only: ijk_p_g
      use bc
      use fldvar, only: x_g, t_g, p_g
      use mms, only: calculate_mms, calculate_mms_source, use_mms

      use param, only: dimension_bc
      use param1, only: undefined_i

      use mpi_utility
      use sendrecv
      IMPLICIT NONE

! Local variables
!--------------------------------------------------------------------//
! Local index for boundary condition
      INTEGER ::  L
!--------------------------------------------------------------------//

! Incompressible cases require that Ppg specified for one cell.
! The following attempts to pick an appropriate cell.
      CALL set_ijk_p_g

      IF(use_mms) THEN
! IJK_P_G is set as UNDEFINED for MMS since current IJK_P_G setting
! is not a second order operation.
         ijk_p_g = undefined_i
! Calculate MMS variables. Better place might be inside interate before
! every time-step. (?)
         CALL calculate_mms
         CALL calculate_mms_source
      END IF

      DO l = 1, dimension_bc
         IF (bc_defined(l)) THEN

            SELECT CASE (bc_type_enum(l))

            CASE (free_slip_wall)
               CALL set_bc0_walls(l)
            CASE (no_slip_wall)
               CALL set_bc0_walls(l)
            CASE (par_slip_wall)
               CALL set_bc0_walls(l)
            CASE (p_outflow)
               CALL set_bc0_init_bcdt_calcs(l)
               CALL set_bc0_outflow(l)
               CALL set_bc0_init_bcscalars(l)
            CASE (mass_outflow)
               CALL set_bc0_init_bcdt_calcs(l)
               CALL set_bc0_inflow(l)
               CALL set_bc0_init_bcscalars(l)
            CASE (outflow)
               CALL set_bc0_init_bcdt_calcs(l)
               CALL set_bc0_outflow(l)
               CALL set_bc0_init_bcscalars(l)
            CASE (mass_inflow)
               CALL set_bc0_init_jet(l)
               CALL set_bc0_inflow(l)
            CASE (p_inflow)
               CALL set_bc0_inflow(l)
            END SELECT
         ENDIF
      ENDDO

! Make T_g nonzero in k=0,1 ghost layers when k-decomposition employed
      call send_recv(t_g,2)
      call send_recv(p_g,2)
      call send_recv(x_g,2)

      RETURN
      END SUBROUTINE set_bc0


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: set_bc0_walls                                           C
!  Purpose: Define certain field variables at the wall boundaries      C
!  according to the user specifications. These are not the real        C
!  values in the wall cells, only initial guesses.                     C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE set_bc0_walls(BCV)

! Modules
!--------------------------------------------------------------------//
      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 bc, only: bc_tw_g, bc_tw_s, bc_thetaw_m
      use bc, only: bc_xw_g, bc_xw_s, bc_scalarw

      use fldvar, only: x_g, x_s, scalar
      use fldvar, only: t_g, t_s, theta_m

      use param1, only: undefined
      use physprop, only: smax, mmax, nmax
      use scalars, only: nscalar

      use functions, only: is_on_mype_plus2layers, bound_funijk
      use compar, only: dead_cell_at
      IMPLICIT NONE

! Dummy arguments
!--------------------------------------------------------------------//
! index of boundary
      INTEGER, INTENT(IN) :: BCV

! Local variables
!--------------------------------------------------------------------//
! indices
      INTEGER :: I, J, K, IJK, M
!--------------------------------------------------------------------//

      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)

         IF (.NOT.is_on_mype_plus2layers(i,j,k)) cycle
         IF (dead_cell_at(i,j,k)) cycle  ! skip dead cells
         ijk = bound_funijk(i,j,k)

         IF(bc_tw_g(bcv) /= undefined) t_g(ijk) = bc_tw_g(bcv)

         IF(nmax(0) > 0) &
            WHERE (bc_xw_g(bcv,:nmax(0)) /= undefined) &
            x_g(ijk,:nmax(0)) = bc_xw_g(bcv,:nmax(0))

         IF(smax > 0) &
            WHERE(bc_tw_s(bcv,:smax) /= undefined) &
            t_s(ijk,:smax) = bc_tw_s(bcv,:smax)

         IF(mmax > 0) &
            WHERE(bc_thetaw_m(bcv,:mmax) /= undefined) &
            theta_m(ijk,:mmax) = bc_thetaw_m(bcv,:mmax)

         DO m = 1, smax
            IF(nmax(m) > 0) &
               WHERE (bc_xw_s(bcv,m,:nmax(m)) /= undefined) &
               x_s(ijk,m,:nmax(m)) = bc_xw_s(bcv,m,:nmax(m))
         ENDDO

         IF(nscalar > 0) &
            WHERE (bc_scalarw(bcv,:nscalar) /= undefined) &
            scalar(ijk,:nscalar) = bc_scalarw(bcv,:nscalar)

      ENDDO   ! end do (i=bc_i_w(l),bc_i_e(l))
      ENDDO   ! end do (j=bc_j_s(l),bc_j_n(l))
      ENDDO   ! end do (k=bc_k_b(l),bc_k_t(l))

      RETURN
      END SUBROUTINE set_bc0_walls


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: set_bc0_outflow                                         C
!  Purpose: Set the initial settings of the boundary conditions (if    C
!  they are defined) for pressure outflow (PO) or outflow (O)          C
!  boundary types.                                                     C
!                                                                      C
!  Comments: For a new run the field variables are undefined in the    C
!  boundary cell locations, while for a restart run the field variable C
!  may have an existing value based on the preceding simulation.       C
!  Regardless, a user defined BC value will supercede any existing     C
!  value.                                                              C
!                                                                      C
!  Note: Scalar field quantities (i.e., T_G, T_s, x_g, x_s, Theta_m,   C
!  scalar, k_turb_g, e_turb_g) do not (should not) need to be defined  C
!  in any type of outflow boundary as the boundary values are unused   C
!  by the corresponding matrix equation (bc_phi). The only reason      C
!  this is potentially necessary is during a calculation another cell  C
!  references the value of this boundary cell before the corresponding C
!  governing equation solver is called.                                C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE set_bc0_outflow(BCV)

! Modules
!--------------------------------------------------------------------//
      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 bc, only: bc_p_g
      use bc, only: bc_t_g, bc_t_s, bc_theta_m
      use bc, only: bc_x_g, bc_x_s, bc_scalar
      use bc, only: bc_ep_g, bc_rop_s
      use bc, only: bc_k_turb_g, bc_e_turb_g

      use constant, only: pi

      use fldvar, only: x_g, x_s, scalar
      use fldvar, only: t_g, t_s, theta_m
      use fldvar, only: p_g, p_star
      use fldvar, only: k_turb_g, e_turb_g
      use fldvar, only: ep_g, rop_s
      use fldvar, only: d_p, ro_s

      use mms, only: use_mms
      use mms, only: mms_p_g, mms_t_g, mms_ep_g
      use mms, only: mms_theta_m, mms_t_s, mms_rop_s

      use param1, only: undefined, zero
      use physprop, only: smax, mmax, nmax
      use run, only: k_epsilon, kt_type_enum, ghd_2007
      use scalars, only: nscalar
      use scales, only: scale_pressure
      use toleranc, only: tmin

      use functions, only: is_on_mype_plus2layers, bound_funijk
      use compar, only: dead_cell_at
      IMPLICIT NONE

! Dummy arguments
!--------------------------------------------------------------------//
! index of boundary condition
      INTEGER, INTENT(IN) :: BCV

! Local variables
!--------------------------------------------------------------------//
! indices
      INTEGER :: I, J, K, IJK, M
! number densities for use in GHD theory only
      DOUBLE PRECISION :: nM, nTOT
!--------------------------------------------------------------------//


      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)
         IF (.NOT.is_on_mype_plus2layers(i,j,k)) cycle
         IF (dead_cell_at(i,j,k)) cycle  ! skip dead cells
         ijk = bound_funijk(i,j,k)

         p_star(ijk) = zero
         p_g(ijk) = scale_pressure(bc_p_g(bcv))
         IF (bc_ep_g(bcv) /= undefined) ep_g(ijk) = bc_ep_g(bcv)

         t_g(ijk)= merge(bc_t_g(bcv), tmin,&
            bc_t_g(bcv) /= undefined)

         IF (nmax(0) > 0) &
            WHERE (bc_x_g(bcv,:nmax(0)) /= undefined) &
            x_g(ijk,:nmax(0)) = bc_x_g(bcv,:nmax(0))

         IF (nscalar > 0) &
            WHERE (bc_scalar(bcv,:nscalar) /= undefined) &
            scalar(ijk,:nscalar) = bc_scalar(bcv,:nscalar)

         IF (k_epsilon) THEN
            IF (bc_k_turb_g(bcv) /= undefined) &
               k_turb_g(ijk) = bc_k_turb_g(bcv)
            IF (bc_e_turb_g(bcv) /= undefined) &
                e_turb_g(ijk) = bc_e_turb_g(bcv)
         ENDIF

         DO m = 1, smax
            IF (bc_rop_s(bcv,m) /= undefined) &
               rop_s(ijk,m) = bc_rop_s(bcv,m)
            IF(bc_t_s(bcv,m) /= undefined) &
               t_s(ijk,m)=bc_t_s(bcv,m)
            IF (bc_theta_m(bcv,m) /= undefined) &
               theta_m(ijk,m) = bc_theta_m(bcv,m)

            IF (nmax(m) > 0) &
               WHERE (bc_x_s(bcv,m,:nmax(m)) /= undefined) &
               x_s(ijk,m,:nmax(m)) = bc_x_s(bcv,m,:nmax(m))
         ENDDO

! for GHD theory to compute mixture BC of velocity and density
         IF(kt_type_enum == ghd_2007) THEN
            rop_s(ijk,mmax) = zero
            ntot = zero
            theta_m(ijk,mmax) = zero
            DO m = 1, smax
               IF (bc_rop_s(bcv,m) /= undefined) THEN
                  rop_s(ijk,mmax) = rop_s(ijk,mmax) + &
                     bc_rop_s(bcv,m)
                  nm = bc_rop_s(bcv,m)*6d0 / &
                     (pi*d_p(ijk,m)**3*ro_s(ijk,m))
                  ntot = ntot + nm
                  IF (bc_theta_m(bcv,m) /= undefined) &
                     theta_m(ijk,mmax) = theta_m(ijk,mmax) + &
                     nm*bc_theta_m(bcv,m)
               ENDIF
            ENDDO
            IF(rop_s(ijk,mmax) > zero) &
               theta_m(ijk,mmax) = theta_m(ijk,mmax) / ntot
         ENDIF   ! end if(kt_type_enum== ghd_2007)

! Set MMS BCs when PO boundary condition is used.
         IF (use_mms) THEN
            p_g(ijk) = scale_pressure(mms_p_g(ijk))
            ep_g(ijk) = mms_ep_g(ijk)
            t_g(ijk) = mms_t_g(ijk)

            m = 1 ! Single solid phase for MMS cases
            rop_s(ijk,m) = mms_rop_s(ijk)
            t_s(ijk,m) = mms_t_s(ijk)
               theta_m(ijk,m) = mms_theta_m(ijk)
         ENDIF ! end if(USE_MMS)

      ENDDO   ! do i
      ENDDO   ! do j
      ENDDO   ! do k

      RETURN
      END SUBROUTINE set_bc0_outflow


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: set_bc0_inflow                                          C
!  Purpose: Set the initial settings of the boundary conditions        C
!  pressure inflow (PI) and. mass inflow (MI) boundary types. Also do  C
!  mass outflow (MO) boundary types... due to velocity...              C
!                                                                      C
!  Comments: Unlike the treament of PO or O boundary types no checks   C
!  are made for these boundary types to determine whether the given    C
!  BC value is defined before it is assigned to the field variable.    C
!  However, the corresponding check routines generally ensure such BC  C
!  quantities are defined for MI or PI boundaries if they are needed   C
!  for the simulation.                                                 C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE set_bc0_inflow(BCV)

! Modules
!--------------------------------------------------------------------//
      use bc, only: bc_plane
      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 bc, only: bc_u_g, bc_v_g, bc_w_g
      use bc, only: bc_u_s, bc_v_s, bc_w_s
      use bc, only: bc_p_g
      use bc, only: bc_t_g, bc_t_s, bc_theta_m
      use bc, only: bc_x_g, bc_x_s, bc_scalar
      use bc, only: bc_ep_g, bc_rop_s
      use bc, only: bc_k_turb_g, bc_e_turb_g

      use constant, only: pi

      use fldvar, only: u_g, v_g, w_g
      use fldvar, only: u_s, v_s, w_s
      use fldvar, only: x_g, x_s, scalar
      use fldvar, only: t_g, t_s, theta_m
      use fldvar, only: p_g, p_star
      use fldvar, only: k_turb_g, e_turb_g
      use fldvar, only: ep_g, rop_s
      use fldvar, only: d_p, ro_s

      use mms, only: use_mms
      use mms, only: mms_u_g, mms_v_g, mms_w_g
      use mms, only: mms_u_s, mms_v_s, mms_w_s
      use mms, only: mms_p_g, mms_t_g, mms_ep_g
      use mms, only: mms_theta_m, mms_t_s, mms_rop_s

      use param1, only: zero
      use physprop, only: smax, mmax, nmax
      use run, only: k_epsilon, kt_type_enum, ghd_2007
      use scalars, only: nscalar
      use scales, only: scale_pressure

      use indices, only: im1, jm1, km1
      use functions, only: is_on_mype_plus2layers, bound_funijk
      use compar, only: dead_cell_at
      IMPLICIT NONE

! Dummy arguments
!--------------------------------------------------------------------//
! index for boundary condition
      INTEGER, INTENT(IN) :: BCV

! Local variables
!--------------------------------------------------------------------//
! indices
      INTEGER :: I, J, K, IJK, M
! ijk index for setting normal component of velocity
      INTEGER :: FIJK
! number densities for use in GHD theory only
      DOUBLE PRECISION :: nM, nTOT
! calculation for normal component of mixture velocity
      DOUBLE PRECISION :: lvel_s
!--------------------------------------------------------------------//

      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)
         IF (.NOT.is_on_mype_plus2layers(i,j,k)) cycle
         IF (dead_cell_at(i,j,k)) cycle  ! skip dead cells
         ijk = bound_funijk(i,j,k)

         p_star(ijk) = zero
         p_g(ijk) = scale_pressure(bc_p_g(bcv))
         ep_g(ijk) = bc_ep_g(bcv)
         t_g(ijk) = bc_t_g(bcv)

         IF (nmax(0) > 0) &
           x_g(ijk,:nmax(0)) = bc_x_g(bcv,:nmax(0))

         IF (nscalar > 0) &
           scalar(ijk,:nscalar) = bc_scalar(bcv,:nscalar)

         IF (k_epsilon) THEN
           k_turb_g(ijk) = bc_k_turb_g(bcv)
           e_turb_g(ijk) = bc_e_turb_g(bcv)
         ENDIF

         DO m = 1, smax
           rop_s(ijk,m) = bc_rop_s(bcv,m)
           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))
         ENDDO

         u_g(ijk) = bc_u_g(bcv)
         v_g(ijk) = bc_v_g(bcv)
         w_g(ijk) = bc_w_g(bcv)
         IF (smax > 0) THEN
            u_s(ijk,:smax) = bc_u_s(bcv,:smax)
            v_s(ijk,:smax) = bc_v_s(bcv,:smax)
            w_s(ijk,:smax) = bc_w_s(bcv,:smax)
         ENDIF

! When the boundary plane is located on the E, N, T side of the domain
! (fluid cell is located w, s, b), set the component of velocity normal
! to the boundary plane of the adjacent fluid cell
         SELECT CASE (trim(bc_plane(bcv)))
            CASE ('W')
               fijk = bound_funijk(im1(i),j,k)
               u_g(fijk) = bc_u_g(bcv)
               IF (smax >0) u_s(fijk,:smax) = bc_u_s(bcv,:smax)
            CASE ('S')
               fijk = bound_funijk(i,jm1(j),k)
               v_g(fijk) = bc_v_g(bcv)
               IF(smax>0) v_s(fijk,:smax) = bc_v_s(bcv,:smax)
            CASE ('B')
               fijk = bound_funijk(i,j,km1(k))
               w_g(fijk) = bc_w_g(bcv)
               IF (smax>0) w_s(fijk,:smax) = bc_w_s(bcv,:smax)
         END SELECT

! for GHD theory to compute mixture BC of velocity and density
         IF(kt_type_enum == ghd_2007) THEN
            rop_s(ijk,mmax) = zero
            ntot = zero
            theta_m(ijk,mmax) = zero
            u_s(ijk,mmax) =  zero
            v_s(ijk,mmax) =  zero
            w_s(ijk,mmax) =  zero
            lvel_s = zero
            DO m = 1, smax
               rop_s(ijk,mmax) = rop_s(ijk,mmax) + &
                  bc_rop_s(bcv,m)
               nm = bc_rop_s(bcv,m)*6d0/ &
                  (pi*d_p(ijk,m)**3*ro_s(ijk,m))
               ntot = ntot + nm
               theta_m(ijk,mmax) = theta_m(ijk,mmax) + &
                  nm*bc_theta_m(bcv,m)
               u_s(ijk,mmax) = u_s(ijk,mmax) + &
                  bc_rop_s(bcv,m)*bc_u_s(bcv,m)
               v_s(ijk,mmax) = v_s(ijk,mmax) + &
                  bc_rop_s(bcv,m)*bc_v_s(bcv,m)
               w_s(ijk,mmax) = w_s(ijk,mmax) + &
                  bc_rop_s(bcv,m)*bc_w_s(bcv,m)
! set velocity component normal to plane in adjacent fluid cell
               SELECT CASE (trim(bc_plane(bcv)))
                  CASE ('W')
                     lvel_s = lvel_s + bc_rop_s(bcv,m)*bc_u_s(bcv,m)
                  CASE ('S')
                     lvel_s = lvel_s + bc_rop_s(bcv,m)*bc_v_s(bcv,m)
                  CASE ('B')
                     lvel_s = lvel_s + bc_rop_s(bcv,m)*bc_w_s(bcv,m)
               END SELECT
            ENDDO
            IF(rop_s(ijk,mmax) > zero) THEN
               theta_m(ijk,mmax) = theta_m(ijk,mmax) / ntot
               u_s(ijk,mmax) = u_s(ijk,mmax) / rop_s(ijk,mmax)
               v_s(ijk,mmax) = v_s(ijk,mmax) / rop_s(ijk,mmax)
               w_s(ijk,mmax) = w_s(ijk,mmax) / rop_s(ijk,mmax)
               lvel_s = lvel_s/rop_s(ijk,mmax)
            ENDIF
            SELECT CASE (trim(bc_plane(bcv)))
               CASE ('W')
                 fijk = bound_funijk(im1(i),j,k)
                 u_s(fijk,mmax) = lvel_s
               CASE ('S')
                 fijk = bound_funijk(i,jm1(j),k)
                 v_s(fijk,mmax) = lvel_s
               CASE ('B')
                 fijk = bound_funijk(i,j,km1(k))
                 w_s(fijk,mmax) = lvel_s
            END SELECT

         ENDIF  ! end if (trim(kt_type) ='ghd')

! Set MMS BCs when MI boundary condition is used.
         IF (use_mms) THEN
            p_g(ijk) = scale_pressure(mms_p_g(ijk))
            ep_g(ijk) = mms_ep_g(ijk)
            t_g(ijk) = mms_t_g(ijk)

            DO m = 1, smax
               rop_s(ijk,m) = mms_rop_s(ijk)
               t_s(ijk,m) = mms_t_s(ijk)
               theta_m(ijk,m) = mms_theta_m(ijk)
            ENDDO

            u_g(ijk) = mms_u_g(ijk)
            v_g(ijk) = mms_v_g(ijk)
            w_g(ijk) = mms_w_g(ijk)
            IF (smax > 0) THEN
               u_s(ijk,:smax) = mms_u_s(ijk)
               v_s(ijk,:smax) = mms_v_s(ijk)
               w_s(ijk,:smax) = mms_w_s(ijk)
            ENDIF
! When the boundary plane is W, S, or B the normal component of velocity
! needs to be set for both sides of the boundary cell.
            SELECT CASE (trim(bc_plane(bcv)))
               CASE ('W')
                 fijk = bound_funijk(im1(i),j,k)
                 u_g(fijk) = mms_u_g(fijk)
                 IF(smax>0) u_s(fijk,:smax) = mms_u_s(fijk)
               CASE ('S')
                 fijk = bound_funijk(i,jm1(j),k)
                 v_g(fijk) = mms_v_g(fijk)
                 IF(smax>0) v_s(fijk,:smax) = mms_v_s(fijk)
               CASE ('B')
                 fijk = bound_funijk(i,j,km1(k))
                 w_g(fijk) = mms_w_g(fijk)
                 IF(smax>0) w_s(fijk,:smax) = mms_w_s(fijk)
            END SELECT
         ENDIF ! end if(USE_MMS)

      ENDDO   ! do i
      ENDDO   ! do j
      ENDDO   ! do k

      RETURN
      END SUBROUTINE set_bc0_inflow


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: set_bc0_init_jet                                        C
!  Purpose: initializing time dependent jet conditions                 C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE set_bc0_init_jet(BCV)

! Modules
!--------------------------------------------------------------------//
      use bc, only: bc_plane
      use bc, only: bc_jet_g, bc_jet_g0
      use bc, only: bc_dt_0, bc_time
      use bc, only: bc_u_g, bc_v_g, bc_w_g
      use param1, only: undefined
      use run, only: time
      IMPLICIT NONE

! Dummy arguments
!--------------------------------------------------------------------//
! index of boundary
      INTEGER, INTENT(IN) :: BCV
!--------------------------------------------------------------------//

      bc_jet_g(bcv) = undefined
      IF (bc_dt_0(bcv) /= undefined) THEN
         bc_time(bcv) = time + bc_dt_0(bcv)
         bc_jet_g(bcv) = bc_jet_g0(bcv)
         IF (bc_jet_g(bcv) /= undefined) THEN
            SELECT CASE (trim(bc_plane(bcv)))
            CASE ('W', 'E')
               bc_u_g(bcv) = bc_jet_g(bcv)
            CASE ('S', 'N')
               bc_v_g(bcv) = bc_jet_g(bcv)
            CASE ('B', 'T')
               bc_w_g(bcv) = bc_jet_g(bcv)
            END SELECT
         ENDIF
      ELSE
         bc_time(bcv) = undefined
      ENDIF
      RETURN
      END SUBROUTINE set_bc0_init_jet


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvC
!                                                                      C
!  Subroutine: set_bc0_init_bcdt_calcs                                 C
!  Purpose: initializing time dependent outflow calculations for       C
!  modifying outflow conditions (MO type) or simple reporting          C
!  outflow conditions (PO or O types)                                  C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^C
      SUBROUTINE set_bc0_init_bcdt_calcs(BCV)

! Modules
!--------------------------------------------------------------------//
      use bc, only: bc_dt_0, bc_time
      use bc, only: bc_out_n
      use bc, only: bc_mout_g, bc_mout_s
      use bc, only: bc_vout_g, bc_vout_s
      use physprop, only: smax
      use param1, only: undefined, zero
      use run, only: time
      IMPLICIT NONE
! Dummy arguments
!--------------------------------------------------------------------//
! index of boundary
      INTEGER, INTENT(IN) :: BCV
!--------------------------------------------------------------------//

! initializing for time dependent outflow reporting calculation
      IF (bc_dt_0(bcv) /= undefined) THEN
         bc_time(bcv) = time + bc_dt_0(bcv)
         bc_out_n(bcv) = 0
         bc_mout_g(bcv) = zero
         bc_vout_g(bcv) = zero
         IF (smax > 0) THEN
            bc_mout_s(bcv,:smax) = zero
            bc_vout_s(bcv,:smax) = zero
         ENDIF
      ELSE
         bc_time(bcv) = undefined
      ENDIF
      RETURN
      END SUBROUTINE set_bc0_init_bcdt_calcs


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      C
!  Subroutine: set_bc0_init_scalars                                    C
!  Purpose: Provide a mechanism to initialize certain field variables  C
!  in the boundary cell to the value of the neighboring fluid cell.    C
!  It is not clear to what extent initialization is necessary as the   C
!  governing matrix equation for each of the field variables is set    C
!  such that the bc cells are not actually solved but are set to the   C
!  value of the adjacent fluid cell(see bc_phi). However, there may be C
!  some instance in the code where the bc cells are referenced before  C
!  the code has called the governing equation.                         C
!  This subroutine illustrates how these particular variables should   C
!  not need (or use) bc values defined by the user!                    C
!                                                                      C
!  Comments: this call replaces some initializations that were         C
!  'inadvertently' occurring in set_bc1                                C
!                                                                      C
!  WARNING: this routine only serves to initialize field variables     C
!  whose governing solver routine invokes bc_phi! do not insert any    C
!  other initializations here as it is likely not appropriate!!        C
!                                                                      C
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE set_bc0_init_bcscalars(BCV)

! Modules
!--------------------------------------------------------------------//
      use bc, only: bc_plane
      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: t_g, t_s, theta_m
      use fldvar, only: x_g, x_s, scalar
      use fldvar, only: k_turb_g, e_turb_g
      use physprop, only: smax, mmax
      use physprop, only: nmax
      use run, only: k_epsilon, kt_type_enum, ghd_2007
      use scalars, only: nscalar

      use indices, only: im1, ip1, jm1, jp1, km1, kp1
      use functions, only: is_on_mype_plus2layers, bound_funijk
      use compar, only: dead_cell_at
      IMPLICIT NONE

! Dummy arguments
!--------------------------------------------------------------------//
! index of boundary
      INTEGER, INTENT(IN) :: BCV

! Local variables
!--------------------------------------------------------------------//
! indices of current boundary cell
      INTEGER :: IJK, I, J, K
! index of neighboring fluid cell
      INTEGER :: FIJK
! local indices
      INTEGER :: M
!--------------------------------------------------------------------//

      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)
         IF (.NOT.is_on_mype_plus2layers(i,j,k)) cycle
         IF (dead_cell_at(i,j,k)) cycle  ! skip dead cells
         ijk = bound_funijk(i,j,k)

         SELECT CASE (trim(bc_plane(bcv)))
            CASE ('W')   ! fluid cell at west
               fijk = bound_funijk(im1(i),j,k)
            CASE ('E')
               fijk = bound_funijk(ip1(i),j,k)
            CASE ('S')   ! fluid cell at south
               fijk = bound_funijk(i,jm1(j),k)
            CASE ('N')
               fijk = bound_funijk(i,jp1(j),k)
            CASE ('B')   ! fluid cell at bottom
               fijk = bound_funijk(i,j,km1(k))
            CASE ('T')
               fijk = bound_funijk(i,j,kp1(k))
         END SELECT

         t_g(ijk) = t_g(fijk)
         IF (nmax(0) > 0) &
            x_g(ijk,:nmax(0)) = x_g(fijk,:nmax(0))

! setting scalar quantities
         IF (nscalar >0) &
            scalar(ijk, :nscalar) = scalar(fijk, :nscalar)

! setting turbulence quantities
         IF(k_epsilon) THEN
            k_turb_g(ijk) = k_turb_g(fijk)
            e_turb_g(ijk) = e_turb_g(fijk)
         ENDIF

! this should only loop of tfm phases
         DO m = 1, smax
            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))
         ENDDO

         IF(kt_type_enum == ghd_2007) &
             theta_m(ijk,mmax) =  theta_m(fijk,mmax)

      ENDDO
      ENDDO
      ENDDO
      RETURN
      END SUBROUTINE set_bc0_init_bcscalars


!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Subroutine: SET_IJK_P_G                                             !
!  Purpose: Pick an appropriate control volume to specify Ppg.         !
!                                                                      !
!  Author: J. Musser                                  Date: 07-Nov-13  !
!  Reviewer:                                          Date:            !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE set_ijk_p_g

! IJK location where Ppg is fixed.
      use bc, only: ijk_p_g
! Specified constant gas density.
      use physprop, only: ro_g0

      use bc
      use exit, only: mfix_exit
      use funits, only: dmp_log
      use geometry, only: cyclic_x, cyclic_x_pd, cyclic_x_mf
      use geometry, only: cyclic_y, cyclic_y_pd, cyclic_y_mf
      use geometry, only: cyclic_z, cyclic_z_pd, cyclic_z_mf
      use geometry, only: do_k
      use geometry, only: imax1, imin1
      use geometry, only: jmax1, jmin1
      use geometry, only: kmax1, kmin1

! MFIX Runtime parameters:
      use param, only: dimension_bc
      use param1, only: undefined
      use param1, only: undefined_i

      use mpi_utility

      implicit none
!--------------------------------------------------------------------//
      INTEGER :: BCV

      CHARACTER(len=7) :: Map
      CHARACTER(len=128) :: lMsg

      INTEGER :: l3
      INTEGER :: l2, u2
      INTEGER :: l1, u1

      LOGICAL, parameter :: setDBG = .false.
      LOGICAL :: dFlag
      INTEGER :: iErr

!--------------------------------------------------------------------//

      dflag = (dmp_log .AND. setdbg)

! Initialize.
      ierr = 0
      ijk_p_g = undefined_i

! This is not needed for compressible cases.
      IF(ro_g0 == undefined) THEN
         IF(dflag) write(*,"(3x,A)")                                   &
            'Compressible: IJK_P_g remaining undefined.'
         return
      ELSEIF(ro_g0 == 0.0d0) THEN
         IF(dflag) write(*,"(3x,A)")                                   &
            'No gas phase: IJK_P_g remaining undefined.'
         return
      ENDIF

! If there are no cyclic boundaries, look for a pressure outflow.
      lpbcv: DO bcv = 1, dimension_bc
         IF (.NOT.bc_defined(bcv)) cycle lpbcv
         IF (bc_type_enum(bcv) == p_outflow .OR. &
             bc_type_enum(bcv) == p_inflow) THEN
            IF(dflag) write(*,"(3x,A)")                                &
               'Outflow PC defined: IJK_P_g remaining undefined.'
            RETURN
         ENDIF
      ENDDO lpbcv

! Initialize.
         l3 = undefined_i
         l2 = undefined_i; u2=l2
         l1 = undefined_i; u1=l1

! If there are cyclic boundaries, flag a cell along the positive
! domain extreme in the cyclic direction (e.g., JMAX1).
      IF(cyclic_y .OR. cyclic_y_pd .OR. cyclic_y_mf) THEN

         map = 'JKI_MAP'
         l3 = jmax1
         l2 = kmin1;  u2 = kmax1
         l1 = imin1;  u1 = imax1
         lmsg='Cyclic in Y'

      ELSEIF(cyclic_x .OR. cyclic_x_pd .OR. cyclic_x_mf) THEN

         map = 'IKJ_MAP'
         l3 = imax1
         l2 = kmin1;  u2 = kmax1
         l1 = jmin1;  u1 = jmax1
         lmsg='Cyclic in X'

      ELSEIF(cyclic_z .OR. cyclic_z_pd .OR. cyclic_z_mf) THEN

         map = 'KIJ_MAP'
         l3 = kmax1
         l2 = imin1;  u2 = imax1
         l1 = jmin1;  u1 = jmax1
         lmsg='Cyclic in Z'

      ENDIF

! No cyclic boundaries or pressure outflows. The IJ plane is used in
! this case to maximize search region for 2D problems.
      IF(l3 == undefined_i) THEN
         map = 'KIJ_MAP'
         l3 = merge(max((kmax1-kmin1)/2+1,2), kmin1, do_k)
         l2 = imin1;  u2 = imax1
         l1 = jmin1;  u1 = jmax1
         lmsg='Center of domain'
      ENDIF

! Debugging messages.
      IF(dflag) THEN
         write(*,"(3/,3x,'Map: ',A)") map
         write(*,"(/5x,'l3:',2x,I4)") l3
         write(*,"( 5x,'l2:',2(2x,I4))") l2, u2
         write(*,"( 5x,'l1:',2(2x,I4))") l1, u1
         write(*,"( 5x,'Msg: ',A)") trim(lmsg)
      ENDIF

! Invoke the search routine.
      CALL ijk_pg_search(l3, l2, u2, l1, u1, map, dflag, ierr)

      IF(ierr == 0) RETURN

! Error management.
      IF(dmp_log) THEN
         SELECT CASE (ierr)
         CASE ( 1001);  WRITE(unit_log, 1001); WRITE(*,1001)
         CASE ( 2000);  WRITE(unit_log, 2000); WRITE(*,2000)
         CASE ( 2001);  WRITE(unit_log, 2001); WRITE(*,2001)
         CASE ( 2002);  WRITE(unit_log, 2002); WRITE(*,2002)
         CASE DEFAULT
            WRITE(unit_log, 1000) ierr
            WRITE(*,1000) ierr
         END SELECT

         WRITE(unit_log, 9000) map(1:1), l3, map(2:2),                 &
            l2, u2, map(3:3), l1, u1
         WRITE(*, 9000) map(1:1), l3, map(2:2),                        &
            l2, u2, map(3:3), l1, u1

         WRITE(*, 9999)
         WRITE(unit_log, 9999)

      ENDIF


      CALL mfix_exit(mype)


 1000 FORMAT(//1x,70('*')/' From: SET_IJK_Pg',/,                       &
         ' Error 1000: Unknown error reported. x', i4.4)

 1001 FORMAT(//1x,70('*')/' From: SET_IJK_Pg',/,                       &
         ' Error 1001: Invalid mapping function.')

 2000 FORMAT(//1x,70('*')/' From: SET_IJK_Pg > IJK_Pg_SEARCH',/,       &
         ' Error 2000: Unknown error reported from IJK_Pg_SEARCH.')

 2001 FORMAT(//1x,70('*')/' From: SET_IJK_Pg > IJK_Pg_SEARCH',/,       &
         ' Error 2001: Unable to locate fluid cell in search region.')

 2002 FORMAT(//1x,70('*')/' From: SET_IJK_Pg > IJK_Pg_SEARCH',/,       &
         ' Error 2002: Unable to locate fluid cell owner.')

 9000 FORMAT(/' Search plane information:',/,3x,a1,': ',i8,            &
          2(/3x,a1,': ',i8,' x ',i8))

 9999 FORMAT(/' Fatal Error --> Invoking MFIX_EXIT',/1x,70('*'),2/)

      END SUBROUTINE set_ijk_p_g

!vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv!
!                                                                      !
!  Author: J. Musser                                  Date: 07-Nov-13  !
!  Reviewer:                                          Date:            !
!                                                                      !
!^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^!
      SUBROUTINE ijk_pg_search(ll3, ll2, lu2, ll1, lu1, lMAP,          &
         ldflag, ierr)

! Modules
!--------------------------------------------------------------------//
! IJK location where Ppg is fixed.
      USE exit, only: mfix_exit
      use bc, only: ijk_p_g
      use functions
      use indices
      use mpi_utility
      use param1, only: undefined_i
      implicit none

! Dummy arguments
!--------------------------------------------------------------------//
      INTEGER, INTENT(IN)  :: ll3
      INTEGER, INTENT(IN)  :: ll2, lu2
      INTEGER, INTENT(IN)  :: ll1, lu1
      LOGICAL, INTENT(IN) :: ldFlag
      INTEGER, INTENT(OUT)  :: iErr
      CHARACTER(len=*), INTENT(IN) :: lMAP

! Local variables
!--------------------------------------------------------------------//
      INTEGER :: lc2, lS2, lE2
      INTEGER :: lc1, lS1, lE1
      INTEGER :: I, J, K, IJK
      LOGICAL :: recheck
      INTEGER :: IJK_Pg_Owner, proc
      INTEGER :: gIJK(0:numpes-1)
      INTEGER :: I_J_K_Pg(3)
      INTEGER :: lpCnt

      CHARACTER(len=32) :: cInt

!--------------------------------------------------------------------//

! Initialize Error Flag
      ierr = 2000

! Initialize the Owner ID
      ijk_pg_owner = undefined_i

! Set the initial search region, a single cell.
      ls1 = ll1 + (lu1-ll1)/2 + 1; le1 = ls1
      ls2 = ll2 + (lu2-ll2)/2 + 1; le2 = ls2

      lpcnt = 1
      recheck = .true.
      do while(recheck)

! Initialize the global IJK array to zero. Resetting this array inside
! this do-loop is most likely overkill. This loop should only cycle
! if gIJK is zero.
         gijk = 0

! Report debugging information for the search region.
         if(ldflag) then
            write(*,"(/5x,'Pass: ',I4)") lpcnt
            write(*,"( 5x,'lp2 bounds:',2(2x,I4))")ls2, le2
            write(*,"( 5x,'lp1 bounds:',2(2x,I4))")ls1, le1
         endif

         lp2: do lc2 = ls2, le2
         lp1: do lc1 = ls1, le1
! Map the loop counters to I/J/K indices.
            SELECT CASE (lmap)
            CASE ('JKI_MAP')
               i=lc1; j=ll3; k=lc2
            CASE ('IKJ_MAP')
               i=ll3; j=lc1; k=lc2
            CASE ('KIJ_MAP')
               i=lc2; j=lc1; k=ll3
            CASE DEFAULT
               ierr = 1001
            END SELECT

! Only the rank that owns this I/J/K proceeds.
            if(.NOT.is_on_mype_owns(i,j,k)) cycle
! Calculate the triple loop index.
            ijk = funijk(i,j,k)
! If there is fluid at this location, store the IJK and exit loops.
            if(fluid_at(ijk)) then
               gijk(mype) = ijk
               exit lp2
            endif
         enddo lp1
         enddo lp2

! Sync gIJK across all processes. Select the lowest ranked process that
! has gIJK defined. This choice is arbitray and doesn't really matter.
! It just needs to be consistent.
         CALL global_all_sum(gijk)
         proc_lp: do proc=0, numpes-1
            if(gijk(proc) /= 0) then
               ijk_p_g = gijk(proc)
               ijk_pg_owner = proc
               recheck = .false.
               exit proc_lp
            endif
         enddo proc_lp

! If the proceeding section did not find a fluid cell, expand the search
! area and try again.
         if(recheck) then
            if(ls1 > ll1 .OR. le1 < lu1 .OR.                           &
               ls2 > ll2 .OR. le2 < lu2) then
! Expand the 1-axis
               ls1 = max((ls1-1), ll1)
               le1 = min((le1+1), lu1)
! Expand the 2-axis
               ls2 = max((ls2-1), ll2)
               le2 = min((le2+1), lu2)
! The entire seach plane was checked with no fluid cell identified.
! Force IJK_P_g to undefined for later error checking.
            else
               recheck = .false.
               ijk_p_g = undefined_i
            endif
         endif
      enddo

! Verify that one fluid cell was detected. Otherwise flag the possible
! errors and return.
      if(ijk_p_g == undefined_i) then
         ierr = 2001
         return
      elseif(ijk_pg_owner == undefined_i) then
         ierr = 2002
         return
      endif


! The owner if the IJK_Pg gets the global I/J/K values and sends
! them to all ranks.
      if(mype == ijk_pg_owner) then
         i_j_k_pg(1) = i_of(ijk_p_g)
         i_j_k_pg(2) = j_of(ijk_p_g)
         i_j_k_pg(3) = k_of(ijk_p_g)
      endif
      CALL bcast(i_j_k_pg, ijk_pg_owner)

      i = i_j_k_pg(1)
      j = i_j_k_pg(2)
      k = i_j_k_pg(3)

! If debugging, have PE_IO report some information before the
! data is overwritten.
      if(ldflag) then
         write(*,"(/3x,'IJK_P_g successfully identified!')")
         cint=''; write(cint,*) ijk_pg_owner
         write(*,"( 5x,'Owner Rank: ',A)")trim(adjustl(cint))
         cint=''; write(cint,*) ijk_p_g
         write(*,"(5x, 'IJK: ',A)", advance='no') trim(adjustl(cint))
         write(*,"(' :: ')", advance='no')
         cint=''; write(cint,*) i
         write(*,"('(',A)",advance='no') trim(adjustl(cint))
         cint=''; write(cint,*) j
         write(*,"(',',A)",advance='no') trim(adjustl(cint))
         cint=''; write(cint,*) k
         write(*,"(',',A,')',2/)") trim(adjustl(cint))
      endif

! Ranks that 'see' IJK_P_g store their local IJK value. Everyone else
! resets IJK_P_g to UNDEFINED_I. This removes the need for getting
! I/J/K values later on in source_PPg.
!      IJK_P_g = merge(funijk(I,J,K), UNDEFINED_I,                      &
!         IS_ON_myPE_plus2layers(I,J,K))

      IF(is_on_mype_plus2layers(i,j,k)) THEN
         ijk_p_g = funijk(i,j,k)
      ELSE
         ijk_p_g = undefined_i
      ENDIF

      ierr = 0
      RETURN
      END SUBROUTINE ijk_pg_search