mus_bc_fluid_module.f90 Source File


This file depends on

sourcefile~~mus_bc_fluid_module.f90~~EfferentGraph sourcefile~mus_bc_fluid_module.f90 mus_bc_fluid_module.f90 sourcefile~mus_scheme_layout_module.f90 mus_scheme_layout_module.f90 sourcefile~mus_bc_fluid_module.f90->sourcefile~mus_scheme_layout_module.f90 sourcefile~mus_bc_header_module.f90 mus_bc_header_module.f90 sourcefile~mus_bc_fluid_module.f90->sourcefile~mus_bc_header_module.f90 sourcefile~mus_mixture_module.f90 mus_mixture_module.f90 sourcefile~mus_bc_fluid_module.f90->sourcefile~mus_mixture_module.f90 sourcefile~mus_dervarpos_module.f90 mus_derVarPos_module.f90 sourcefile~mus_bc_fluid_module.f90->sourcefile~mus_dervarpos_module.f90 sourcefile~mus_param_module.f90 mus_param_module.f90 sourcefile~mus_bc_fluid_module.f90->sourcefile~mus_param_module.f90 sourcefile~mus_physics_module.f90 mus_physics_module.f90 sourcefile~mus_bc_fluid_module.f90->sourcefile~mus_physics_module.f90 sourcefile~mus_field_prop_module.f90 mus_field_prop_module.f90 sourcefile~mus_bc_fluid_module.f90->sourcefile~mus_field_prop_module.f90 sourcefile~mus_moments_type_module.f90 mus_moments_type_module.f90 sourcefile~mus_scheme_layout_module.f90->sourcefile~mus_moments_type_module.f90 sourcefile~mus_bc_header_module.f90->sourcefile~mus_scheme_layout_module.f90 sourcefile~mus_bc_header_module.f90->sourcefile~mus_mixture_module.f90 sourcefile~mus_bc_header_module.f90->sourcefile~mus_dervarpos_module.f90 sourcefile~mus_bc_header_module.f90->sourcefile~mus_physics_module.f90 sourcefile~mus_bc_header_module.f90->sourcefile~mus_field_prop_module.f90 sourcefile~mus_pdf_module.f90 mus_pdf_module.f90 sourcefile~mus_bc_header_module.f90->sourcefile~mus_pdf_module.f90 sourcefile~mus_mixture_module.f90->sourcefile~mus_physics_module.f90 sourcefile~mus_enrtl_module.f90 mus_eNRTL_module.f90 sourcefile~mus_mixture_module.f90->sourcefile~mus_enrtl_module.f90 sourcefile~mus_scheme_header_module.f90 mus_scheme_header_module.f90 sourcefile~mus_mixture_module.f90->sourcefile~mus_scheme_header_module.f90 sourcefile~mus_dervarpos_module.f90->sourcefile~mus_scheme_layout_module.f90 sourcefile~mus_param_module.f90->sourcefile~mus_physics_module.f90 sourcefile~mus_abortcriteria_module.f90 mus_abortCriteria_module.f90 sourcefile~mus_param_module.f90->sourcefile~mus_abortcriteria_module.f90 sourcefile~mus_field_prop_module.f90->sourcefile~mus_physics_module.f90 sourcefile~mus_fluid_module.f90 mus_fluid_module.f90 sourcefile~mus_field_prop_module.f90->sourcefile~mus_fluid_module.f90 sourcefile~mus_poisson_module.f90 mus_poisson_module.f90 sourcefile~mus_field_prop_module.f90->sourcefile~mus_poisson_module.f90 sourcefile~mus_field_prop_module.f90->sourcefile~mus_scheme_header_module.f90 sourcefile~mus_species_module.f90 mus_species_module.f90 sourcefile~mus_field_prop_module.f90->sourcefile~mus_species_module.f90 sourcefile~mus_fluid_module.f90->sourcefile~mus_physics_module.f90 sourcefile~mus_fluid_module.f90->sourcefile~mus_pdf_module.f90 sourcefile~mus_fluid_module.f90->sourcefile~mus_scheme_header_module.f90 sourcefile~mus_turb_viscosity_module.f90 mus_turb_viscosity_module.f90 sourcefile~mus_fluid_module.f90->sourcefile~mus_turb_viscosity_module.f90 sourcefile~mus_relaxationparam_module.f90 mus_relaxationParam_module.f90 sourcefile~mus_fluid_module.f90->sourcefile~mus_relaxationparam_module.f90 sourcefile~mus_turbulence_module.f90 mus_turbulence_module.f90 sourcefile~mus_fluid_module.f90->sourcefile~mus_turbulence_module.f90 sourcefile~mus_nonnewtonian_module.f90 mus_nonNewtonian_module.f90 sourcefile~mus_fluid_module.f90->sourcefile~mus_nonnewtonian_module.f90 sourcefile~mus_mrtrelaxation_module.f90 mus_mrtRelaxation_module.f90 sourcefile~mus_fluid_module.f90->sourcefile~mus_mrtrelaxation_module.f90 sourcefile~mus_poisson_module.f90->sourcefile~mus_physics_module.f90 sourcefile~mus_species_module.f90->sourcefile~mus_physics_module.f90 sourcefile~mus_turb_viscosity_module.f90->sourcefile~mus_scheme_layout_module.f90 sourcefile~mus_turb_viscosity_module.f90->sourcefile~mus_scheme_header_module.f90 sourcefile~mus_turb_viscosity_module.f90->sourcefile~mus_turbulence_module.f90 sourcefile~mus_vreman_module.f90 mus_Vreman_module.f90 sourcefile~mus_turb_viscosity_module.f90->sourcefile~mus_vreman_module.f90 sourcefile~mus_smagorinsky_module.f90 mus_Smagorinsky_module.f90 sourcefile~mus_turb_viscosity_module.f90->sourcefile~mus_smagorinsky_module.f90 sourcefile~mus_graddata_module.f90 mus_gradData_module.f90 sourcefile~mus_turb_viscosity_module.f90->sourcefile~mus_graddata_module.f90 sourcefile~mus_wale_module.f90 mus_WALE_module.f90 sourcefile~mus_turb_viscosity_module.f90->sourcefile~mus_wale_module.f90 sourcefile~mus_relaxationparam_module.f90->sourcefile~mus_scheme_layout_module.f90 sourcefile~mus_relaxationparam_module.f90->sourcefile~mus_dervarpos_module.f90 sourcefile~mus_relaxationparam_module.f90->sourcefile~mus_physics_module.f90 sourcefile~mus_relaxationparam_module.f90->sourcefile~mus_moments_type_module.f90 sourcefile~mus_relaxationparam_module.f90->sourcefile~mus_scheme_header_module.f90 sourcefile~mus_relaxationparam_module.f90->sourcefile~mus_turbulence_module.f90 sourcefile~mus_relaxationparam_module.f90->sourcefile~mus_nonnewtonian_module.f90 sourcefile~mus_relaxationparam_module.f90->sourcefile~mus_graddata_module.f90 sourcefile~mus_turbulence_module.f90->sourcefile~mus_scheme_layout_module.f90 sourcefile~mus_turbulence_module.f90->sourcefile~mus_graddata_module.f90 sourcefile~mus_nonnewtonian_module.f90->sourcefile~mus_scheme_layout_module.f90 sourcefile~mus_nonnewtonian_module.f90->sourcefile~mus_physics_module.f90 sourcefile~mus_nonnewtonian_module.f90->sourcefile~mus_scheme_header_module.f90 sourcefile~mus_mrtrelaxation_module.f90->sourcefile~mus_moments_type_module.f90 sourcefile~mus_mrtrelaxation_module.f90->sourcefile~mus_scheme_header_module.f90

Files dependent on this one

sourcefile~~mus_bc_fluid_module.f90~~AfferentGraph sourcefile~mus_bc_fluid_module.f90 mus_bc_fluid_module.f90 sourcefile~mus_bc_general_module.f90 mus_bc_general_module.f90 sourcefile~mus_bc_general_module.f90->sourcefile~mus_bc_fluid_module.f90 sourcefile~mus_debug_module.f90 mus_debug_module.f90 sourcefile~mus_debug_module.f90->sourcefile~mus_bc_general_module.f90 sourcefile~mus_dynloadbal_module.f90 mus_dynLoadBal_module.f90 sourcefile~mus_dynloadbal_module.f90->sourcefile~mus_bc_general_module.f90 sourcefile~mus_construction_module.f90 mus_construction_module.f90 sourcefile~mus_dynloadbal_module.f90->sourcefile~mus_construction_module.f90 sourcefile~mus_control_module.f90 mus_control_module.f90 sourcefile~mus_control_module.f90->sourcefile~mus_bc_general_module.f90 sourcefile~mus_control_module.f90->sourcefile~mus_debug_module.f90 sourcefile~mus_program_module.f90 mus_program_module.f90 sourcefile~mus_program_module.f90->sourcefile~mus_bc_general_module.f90 sourcefile~mus_program_module.f90->sourcefile~mus_dynloadbal_module.f90 sourcefile~mus_program_module.f90->sourcefile~mus_control_module.f90 sourcefile~mus_program_module.f90->sourcefile~mus_construction_module.f90 sourcefile~mus_construction_module.f90->sourcefile~mus_debug_module.f90 sourcefile~musubi.f90 musubi.f90 sourcefile~musubi.f90->sourcefile~mus_control_module.f90 sourcefile~musubi.f90->sourcefile~mus_program_module.f90 sourcefile~mus_harvesting.f90 mus_harvesting.f90 sourcefile~mus_harvesting.f90->sourcefile~mus_construction_module.f90 sourcefile~mus_hvs_construction_module.f90 mus_hvs_construction_module.f90 sourcefile~mus_hvs_construction_module.f90->sourcefile~mus_construction_module.f90

Contents


Source Code

! Copyright (c) 2012-2021 Kannan Masilamani <kannan.masilamani@uni-siegen.de>
! Copyright (c) 2012, 2014-2016 Jiaxing Qi <jiaxing.qi@uni-siegen.de>
! Copyright (c) 2012-2013 Manuel Hasert <m.hasert@grs-sim.de>
! Copyright (c) 2012-2014 Simon Zimny <s.zimny@grs-sim.de>
! Copyright (c) 2012, 2014 Kartik Jain <kartik.jain@uni-siegen.de>
! Copyright (c) 2013, 2019 Harald Klimach <harald.klimach@uni-siegen.de>
! Copyright (c) 2015-2016 Tobias Schneider <tobias1.schneider@student.uni-siegen.de>
! Copyright (c) 2016 Verena Krupp <verena.krupp@uni-siegen.de>
! Copyright (c) 2017-2018 Raphael Haupt <raphael.haupt@uni-siegen.de>
! Copyright (c) 2018 Jana Gericke <jana.gericke@uni-siegen.de>
! Copyright (c) 2019-2020 Peter Vitt <peter.vitt2@uni-siegen.de>
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions are met:
!
! 1. Redistributions of source code must retain the above copyright notice,
! this list of conditions and the following disclaimer.
!
! 2. Redistributions in binary form must reproduce the above copyright notice,
! this list of conditions and the following disclaimer in the documentation
! and/or other materials provided with the distribution.
!
! THIS SOFTWARE IS PROVIDED BY THE UNIVERSITY OF SIEGEN “AS IS” AND ANY EXPRESS
! OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
! OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
! IN NO EVENT SHALL UNIVERSITY OF SIEGEN OR CONTRIBUTORS BE LIABLE FOR ANY
! DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
! (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
! ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
! ****************************************************************************** !
!> author: Manuel Hasert
!! author: Kannan Masilamani
!! Boundary condition treatment routines for fluid simulation
!! Details on implementation can be found in [[tem_bc_module]]
!!
! Copyright (c) 2011-2013 Manuel Hasert <m.hasert@grs-sim.de>
! Copyright (c) 2011 Harald Klimach <harald.klimach@uni-siegen.de>
! Copyright (c) 2011 Konstantin Kleinheinz <k.kleinheinz@grs-sim.de>
! Copyright (c) 2011-2012 Simon Zimny <s.zimny@grs-sim.de>
! Copyright (c) 2012, 2014-2016 Jiaxing Qi <jiaxing.qi@uni-siegen.de>
! Copyright (c) 2012 Kartik Jain <kartik.jain@uni-siegen.de>
! Copyright (c) 2013-2015, 2019 Kannan Masilamani <kannan.masilamani@uni-siegen.de>
! Copyright (c) 2016 Tobias Schneider <tobias1.schneider@student.uni-siegen.de>
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions are met:
!
! 1. Redistributions of source code must retain the above copyright notice,
! this list of conditions and the following disclaimer.
!
! 2. Redistributions in binary form must reproduce the above copyright notice,
! this list of conditions and the following disclaimer in the documentation
! and/or other materials provided with the distribution.
!
! THIS SOFTWARE IS PROVIDED BY THE UNIVERSITY OF SIEGEN “AS IS” AND ANY EXPRESS
! OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
! OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
! IN NO EVENT SHALL UNIVERSITY OF SIEGEN OR CONTRIBUTORS BE LIABLE FOR ANY
! DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
! (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
! ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
! Copyright (c) 2013 Harald Klimach <harald.klimach@uni-siegen.de>
! Copyright (c) 2016 Tobias Schneider <tobias1.schneider@student.uni-siegen.de>
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions are met:
!
! 1. Redistributions of source code must retain the above copyright notice, this
! list of conditions and the following disclaimer.
!
! 2. Redistributions in binary form must reproduce the above copyright notice,
! this list of conditions and the following disclaimer in the documentation
! and/or other materials provided with the distribution.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
! AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
! IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
! DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
! FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
! SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
! OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
! OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
! Make sure loglvl is defined to an integer value.
! Usually this should be defined on the command line with -Dloglvl=
! to some value.
module mus_bc_fluid_module

  ! include treelm modules
  use env_module,               only: rk
  use tem_param_module,         only: div1_3, div4_3, csInv, cs, cs2, cs2inv, &
    &                                 rho0
  use tem_time_module,          only: tem_time_type
  use treelmesh_module,         only: treelmesh_type
  use tem_varSys_module,        only: tem_varSys_type
  use tem_property_module,      only: prp_hasQVal
  ! use tem_subTree_type_module,  only: tem_subTree_type
  use tem_geometry_module,      only: tem_BaryOfId,tem_ElemSizeLevel
  use tem_logging_module,       only: logUnit
  use tem_debug_module,         only: dbgUnit
  use tem_construction_module,  only: tem_levelDesc_type

  ! include musubi modules
  use mus_bc_header_module,          only: boundary_type, glob_boundary_type
  use mus_scheme_layout_module,      only: mus_scheme_layout_type
  use mus_field_prop_module,         only: mus_field_prop_type
  use mus_derVarPos_module,          only: mus_derVarPos_type
  use mus_param_module,              only: mus_param_type
  use mus_physics_module,            only: mus_physics_type
  use mus_mixture_module,            only: mus_mixture_type

  implicit none

  private

  public :: velocity_eq
  public :: mfr_bounceback
  public :: mfr_eq
  public :: velocity_bounceback_incomp, velocity_bounceback
  public :: velocity_bfl, velocity_bfl_incomp
  public :: vel_neq

  public :: bc_pdf

  public :: press_neq
  public :: pressure_antiBounceBack
  public :: pressure_eq
  public :: outlet_dnt
  public :: pressure_expol
  public :: outlet_zero_prsgrd

  public :: inlet_nrbc_incomp, inlet_nrbc
  public :: outlet_nrbc_incomp, outlet_nrbc
  public :: outlet_nrbc_eq

  ! moments based BC
  public :: velocity_momentsbased
  public :: pressure_momentsbased
  public :: moments_wall

  public :: velocity_momentsbased_incomp
  public :: pressure_momentsbased_incomp

contains

! ****************************************************************************** !
  !> Inlet Velocity Equilibrium type boundary conditions for weakly
  !! compressible lbm scheme
  !!
  !! The incoming densities are set to the equilibrium distribution with
  !! macroscopic values:
  !! - density is extrapolated (0 or 1st order) from the flow domain
  !!   + 0-order: so we simply use the density of the boundary element itself
  !! - velocity is used as defined in the configuration file
  !!
  !! Usage
  !! -----
  !!```lua
  !!boundary_condition = {
  !!  { label = 'inlet',
  !!    kind = 'velocity_eq',
  !!    velocity = 'inlet_vel',
  !!  }
  !!}
  !!variable = {
  !!  name = 'inlet_vel',
  !!  ncomponents = 3,
  !!  vartype = 'st_fun',
  !!  st_fun = {0.06, 0.0, 0.0}
  !!}
  !!```
  !!
  !! This subroutine's interface must match the abstract interface definition
  !! [[boundaryRoutine]] in bc/[[mus_bc_header_module]].f90 in order to be
  !! callable via [[boundary_type:fnct]] function pointer.
  subroutine velocity_eq( me, state, bcBuffer, globBC, levelDesc, tree, nSize, &
    &                     iLevel, sim_time, neigh, layout, fieldProp, varPos,  &
    &                     nScalars, varSys, derVarPos, physics, iField,        &
    &                     mixture                                              )
    ! -------------------------------------------------------------------- !
    !> global boundary type
    class(boundary_type) :: me
    !> Current state vector of iLevel
    real(kind=rk), intent(inout) :: state(:)
    !> size of state array ( in terms of elements )
    integer, intent(in) :: nSize
    !> state values of boundary elements of all fields of iLevel
    real(kind=rk), intent(in) :: bcBuffer(:)
    !> iLevel descriptor
    type(tem_levelDesc_type), intent(in) :: levelDesc
    !> Treelm Mesh
    type(treelmesh_type), intent(in) :: tree
    !> fluid parameters and properties
    type(mus_field_prop_type), intent(in) :: fieldProp
    !> stencil layout information
    type(mus_scheme_layout_type), intent(in) :: layout
    !> the level On which this boundary was invoked
    integer, intent(in) :: iLevel
    !> connectivity array corresponding to state vector
    integer, intent(in) :: neigh(:)
    !> global time information
    type(tem_time_type), intent(in)  :: sim_time
    !> pointer to field variable in the state vector
    integer, intent(in) :: varPos(:)
    !> number of Scalars in the scheme var system
    integer, intent(in) :: nScalars
    !> scheme variable system
    type(tem_varSys_type), intent(in) :: varSys
    !> position of derived quantities in varsys
    type(mus_derVarPos_type), intent(in) :: derVarPos
    !> scheme global boundary type
    type(glob_boundary_type), intent(in) :: globBC
    !> scheme global boundary type
    type(mus_physics_type), intent(in) :: physics
    !> current field
    integer, intent(in) :: iField
    !> mixture info
    type(mus_mixture_type), intent(in) :: mixture
    ! -------------------------------------------------------------------- !
    real(kind=rk) :: fEq( layout%fStencil%QQ )
    real(kind=rk) :: rho(1), vel(3,1)
    real(kind=rk) :: vel_b(3*globBC%nElems(iLevel)), inv_vel
    integer :: iELem, iDir, QQ, elemPos, posInBuffer
    integer :: bcVel_pos
    ! ---------------------------------------------------------------------------

    QQ = layout%fStencil%QQ
    inv_vel = 1.0_rk / physics%fac( iLevel )%vel

    ! position of boundary velocity in varSys
    bcVel_pos = me%bc_states%velocity%varPos
    ! Get velocity
    call varSys%method%val(bcVel_pos)%get_valOfIndex( &
      & varSys  = varSys,                             &
      & time    = sim_time,                           &
      & iLevel  = iLevel,                             &
      & idx     = me%bc_states%velocity               &
      &           %pntIndex%indexLvl(iLevel)          &
      &           %val(1:globBC%nElems(iLevel)),      &
      & nVals   = globBC%nElems(iLevel),              &
      & res     = vel_b                               )

    ! convert physical velocity into LB velocity
    vel_b = vel_b * inv_vel

    ! @todo: these if condition can only be done once.
    do iElem = 1, globBC%nElems(iLevel)
      ! Calculate the density of current element
      posInBuffer = globBC%elemLvl( iLevel )%posInBcElemBuf%val( iElem )
      rho(1) = sum( bcBuffer( (posInBuffer-1)*nScalars+varPos(1) :  &
        &                  (posInBuffer-1)*nScalars+varPos(QQ) ) )
      vel(:,1) = vel_b((iElem-1)*3+1:iElem*3)

      ! compute equilibrium
      call derVarPos%equilFromMacro( density  = rho,    &
        &                            velocity = vel,    &
        &                            iField   = iField, &
        &                            nElems   = 1,      &
        &                            varSys   = varSys, &
        &                            layout   = layout, &
        &                            res      = fEq     )

      elemPos = globBC%elemLvl(iLevel)%elem%val( iElem )
      do iDir = 1, layout%fStencil%QQN
        ! Write the values
        ! The bitmask points into the incoming direction into the flow domain,
        !  which actually we want to update
        ! * For PUSH, we write to the incoming position,
        !   as the kernel reads it from there without propagation.
        ! * For PULL, we need to write to the inverse direction, as the kernel
        !   performs a bounce back before reading it.
        !   However, this bounced back direction actually comes from the
        !   non-existent boundary element and would point into the incoming
        !   direction, so the value has to be treated and set as if it points
        !   really into the incoming direction.
        if( globBC%elemLvl(iLevel)%bitmask%val( iDir, iElem )) then
          ! Depending on PUSH or pull, use + or - for cxDir, because directions
          ! are inverted
          state(  neigh (( idir-1)* nsize+ elempos)+( ifield-1)* qq+ nscalars*0) &
            & = fEq( iDir )
        end if ! bitMask
      end do ! iDir
    end do ! iElem

  end subroutine velocity_eq
! ****************************************************************************** !

! ****************************************************************************** !
  !> Velocity Non-Equilibrium type boundary conditions from
  !! Guo, Z., & Shi, B. (2002). "Non-equilibrium extrapolation method for
  !! velocity and pressure boundary conditions in the lattice Boltzmann method."
  !! Chinese Physics, (November 2016).
  !!
  !! The incoming densities are updated as
  !! $$f_i(\vec{x}_b,t) = f^{eq}_i(\vec{x}_b,t)
  !!          + (1-\omega)(f_i(\vec{x}_f,t) - f^{eq}_i(\vec{x}_f,t)) $$
  !! macroscopic values for $$f^{eq}_i(\vec{x}_b,t)$$:
  !! - density is extrapolated (0 or 1st order) from the flow domain
  !!   + 0-order: so we simply use the density of the boundary element itself
  !! - velocity is used as defined in the configuration file
  !!
  !! Usage
  !! -----
  !!```lua
  !!boundary_condition = {
  !!  { label = 'inlet',
  !!    kind = 'vel_neq',
  !!    velocity = 'vel',
  !!  }
  !!}
  !!variable = {
  !!  name = 'inlet_vel',
  !!  ncomponents = 3,
  !!  vartype = 'st_fun',
  !!  st_fun = {0.06, 0.0, 0.0}
  !!}
  !!```
  !!
  !! This subroutine's interface must match the abstract interface definition
  !! [[boundaryRoutine]] in bc/[[mus_bc_header_module]].f90 in order to be
  !! callable via [[boundary_type:fnct]] function pointer.
  subroutine vel_neq( me, state, bcBuffer, globBC, levelDesc, tree, nSize,  &
    &                 iLevel, sim_time, neigh, layout, fieldProp, varPos,   &
    &                 nScalars, varSys, derVarPos, physics, iField, mixture )
    ! -------------------------------------------------------------------- !
    !> global boundary type
    class(boundary_type) :: me
    !> Current state vector of iLevel
    real(kind=rk), intent(inout) :: state(:)
    !> size of state array ( in terms of elements )
    integer, intent(in) :: nSize
    !> state values of boundary elements of all fields of iLevel
    real(kind=rk), intent(in) :: bcBuffer(:)
    !> iLevel descriptor
    type(tem_levelDesc_type), intent(in) :: levelDesc
    !> Treelm Mesh
    type(treelmesh_type), intent(in) :: tree
    !> fluid parameters and properties
    type(mus_field_prop_type), intent(in) :: fieldProp
    !> stencil layout information
    type(mus_scheme_layout_type), intent(in) :: layout
    !> the level On which this boundary was invoked
    integer, intent(in) :: iLevel
    !> connectivity array corresponding to state vector
    integer, intent(in) :: neigh(:)
    !> global time information
    type(tem_time_type), intent(in)  :: sim_time
    !> pointer to field variable in the state vector
    integer, intent(in) :: varPos(:)
    !> number of Scalars in the scheme var system
    integer, intent(in) :: nScalars
    !> scheme variable system
    type(tem_varSys_type), intent(in) :: varSys
    !> position of derived quantities in varsys
    type(mus_derVarPos_type), intent(in) :: derVarPos
    !> scheme global boundary type
    type(glob_boundary_type), intent(in) :: globBC
    !> scheme global boundary type
    type(mus_physics_type), intent(in) :: physics
    !> current field
    integer, intent(in) :: iField
    !> mixture info
    type(mus_mixture_type), intent(in) :: mixture
    ! -------------------------------------------------------------------- !
    real(kind=rk) :: fEq( layout%fStencil%QQ )
    real(kind=rk) :: fEq_Bnd( layout%fStencil%QQ )
    real(kind=rk) :: fTmp( layout%fStencil%QQ )
    real(kind=rk) :: rho(1), vel(3,1), omega
    real(kind=rk) :: vel_b(3*globBC%nElems(iLevel)), inv_vel
    integer :: iELem, iDir, QQ, elemPos, posInBuffer
    integer :: bcVel_pos
    ! ---------------------------------------------------------------------------

    QQ = layout%fStencil%QQ
    inv_vel = 1.0_rk / physics%fac( iLevel )%vel

    ! position of boundary velocity in varSys
    bcVel_pos = me%bc_states%velocity%varPos
    ! Get velocity
    call varSys%method%val(bcVel_pos)%get_valOfIndex( &
      & varSys  = varSys,                             &
      & time    = sim_time,                           &
      & iLevel  = iLevel,                             &
      & idx     = me%bc_states%velocity               &
      &           %pntIndex%indexLvl(iLevel)          &
      &           %val(1:globBC%nElems(iLevel)),      &
      & nVals   = globBC%nElems(iLevel),              &
      & res     = vel_b                               )

    ! convert physical velocity into LB velocity
    vel_b = vel_b * inv_vel

    ! @todo: these if condition can only be done once.
    do iElem = 1, globBC%nElems(iLevel)
      ! Calculate the density of current element
      posInBuffer = globBC%elemLvl( iLevel )%posInBcElemBuf%val( iElem )
      fTmp(1:QQ) = bcBuffer( (posInBuffer-1)*nScalars+varPos(1) : &
        &                    (posInBuffer-1)*nScalars+varPos(QQ)  )
      rho(1) = sum( fTmp )
      vel(:,1) = vel_b((iElem-1)*3+1:iElem*3)

      ! compute equilibrium on Boundary element
      call derVarPos%equilFromMacro( density  = rho,    &
        &                            velocity = vel,    &
        &                            iField   = iField, &
        &                            nElems   = 1,      &
        &                            varSys   = varSys, &
        &                            layout   = layout, &
        &                            res      = fEq_Bnd )

      ! compute equilibrium on Fluid element
      call derVarPos%equilFromState( state  = fTmp,   &
        &                            iField = iField, &
        &                            nElems = 1,      &
        &                            varSys = varSys, &
        &                            layout = layout, &
        &                            res    = fEq     )

      elemPos = globBC%elemLvl(iLevel)%elem%val( iElem )
      ! omega is STfun so get omega for current element
      omega = fieldProp%fluid%viscKine%omLvl(iLevel)%val(elemPos)

      do iDir = 1, layout%fStencil%QQN
        ! Write the values
        ! The bitmask points into the incoming direction into the flow domain,
        !  which actually we want to update
        ! * For PUSH, we write to the incoming position,
        !   as the kernel reads it from there without propagation.
        ! * For PULL, we need to write to the inverse direction, as the kernel
        !   performs a bounce back before reading it.
        !   However, this bounced back direction actually comes from the
        !   non-existent boundary element and would point into the incoming
        !   direction, so the value has to be treated and set as if it points
        !   really into the incoming direction.
        if( globBC%elemLvl(iLevel)%bitmask%val( iDir, iElem )) then
          ! Depending on PUSH or pull, use + or - for cxDir, because directions
          ! are inverted
          state(  neigh (( idir-1)* nsize+ elempos)+( ifield-1)* qq+ nscalars*0) &
            & = fEq_Bnd( iDir ) + (1.0_rk-omega)*(fTmp(iDir) - fEq(iDir))
        end if ! bitMask
      end do ! iDir
    end do ! iElem

  end subroutine vel_neq
! ****************************************************************************** !

! ****************************************************************************** !
  !> Pressure Non-Equilibrium type boundary conditions from
  !! Guo, Z., & Shi, B. (2002). "Non-equilibrium extrapolation method for
  !! velocity and pressure boundary conditions in the lattice Boltzmann method."
  !! Chinese Physics, (November 2016).
  !!
  !! The incoming densities are updated as
  !! $$f_i(\vec{x}_b,t) = f^{eq}_i(\vec{x}_b,t)
  !!          + (1-\omega)(f_i(\vec{x}_f,t) - f^{eq}_i(\vec{x}_f,t)) $$
  !! macroscopic values for $$f^{eq}_i(\vec{x}_b,t)$$:
  !! - density is defined in config file as pressure
  !! - velocity is extrapolated from fluid node
  !!
  !! Usage
  !! -----
  !!```lua
  !!boundary_condition = {
  !!  { label = 'outlet',
  !!    kind = 'pressure_neq',
  !!    pressure = 'press',
  !!  }
  !!}
  !!variable = {
  !!  name = 'press',
  !!  ncomponents = 1,
  !!  vartype = 'st_fun',
  !!  st_fun = 0.0
  !!}
  !!```
  !!
  !! This subroutine's interface must match the abstract interface definition
  !! [[boundaryRoutine]] in bc/[[mus_bc_header_module]].f90 in order to be
  !! callable via [[boundary_type:fnct]] function pointer.
  subroutine press_neq( me, state, bcBuffer, globBC, levelDesc, tree, nSize,  &
    &                   iLevel, sim_time, neigh, layout, fieldProp, varPos,   &
    &                   nScalars, varSys, derVarPos, physics, iField, mixture )
    ! -------------------------------------------------------------------- !
    !> global boundary type
    class(boundary_type) :: me
    !> Current state vector of iLevel
    real(kind=rk), intent(inout) :: state(:)
    !> size of state array ( in terms of elements )
    integer, intent(in) :: nSize
    !> state values of boundary elements of all fields of iLevel
    real(kind=rk), intent(in) :: bcBuffer(:)
    !> iLevel descriptor
    type(tem_levelDesc_type), intent(in) :: levelDesc
    !> Treelm Mesh
    type(treelmesh_type), intent(in) :: tree
    !> fluid parameters and properties
    type(mus_field_prop_type), intent(in) :: fieldProp
    !> stencil layout information
    type(mus_scheme_layout_type), intent(in) :: layout
    !> the level On which this boundary was invoked
    integer, intent(in) :: iLevel
    !> connectivity array corresponding to state vector
    integer, intent(in) :: neigh(:)
    !> global time information
    type(tem_time_type), intent(in)  :: sim_time
    !> pointer to field variable in the state vector
    integer, intent(in) :: varPos(:)
    !> number of Scalars in the scheme var system
    integer, intent(in) :: nScalars
    !> scheme variable system
    type(tem_varSys_type), intent(in) :: varSys
    !> position of derived quantities in varsys
    type(mus_derVarPos_type), intent(in) :: derVarPos
    !> scheme global boundary type
    type(glob_boundary_type), intent(in) :: globBC
    !> scheme global boundary type
    type(mus_physics_type), intent(in) :: physics
    !> current field
    integer, intent(in) :: iField
    !> mixture info
    type(mus_mixture_type), intent(in) :: mixture
    ! -------------------------------------------------------------------- !
    real(kind=rk) :: fEq( layout%fStencil%QQ * globBC%nElems(iLevel) )
    real(kind=rk) :: fEq_Bnd( layout%fStencil%QQ * globBC%nElems(iLevel) )
    real(kind=rk) :: fTmp( layout%fStencil%QQ * globBC%nElems(iLevel) )
    real(kind=rk) :: fTmp_loc( layout%fStencil%QQ )
    real(kind=rk) :: vel(3,globBC%nElems(iLevel)), omega
    real(kind=rk) :: rhoB(globBC%nElems(iLevel)), inv_rho_phy
    integer :: iELem, iDir, QQ, elemPos, posInBuffer, offset
    integer :: bcPress_pos
    ! ---------------------------------------------------------------------------

    QQ = layout%fStencil%QQ
    inv_rho_phy = 1.0_rk / physics%fac(iLevel)%press * cs2inv

    ! position of boundary pressure in varSys
    bcPress_pos = me%bc_states%pressure%varPos
    ! get pressure variable from spacetime function
    call varSys%method%val(bcPress_pos)%get_valOfIndex( &
      & varSys  = varSys,                               &
      & time    = sim_time,                             &
      & iLevel  = iLevel,                               &
      & idx     = me%bc_states%pressure                 &
      &           %pntIndex%indexLvl(iLevel)            &
      &           %val(1:globBC%nElems(iLevel)),        &
      & nVals   = globBC%nElems(iLevel),                &
      & res     = rhoB                                  )

    ! convert physical pressure to LB density
    rhoB =  rhoB* inv_rho_phy + rho0

    do iElem = 1, globBC%nElems(iLevel)
      ! Calculate the density of current element
      posInBuffer = globBC%elemLvl( iLevel )%posInBcElemBuf%val( iElem )
      fTmp_loc = bcBuffer( (posInBuffer-1)*nScalars+varPos(1) : &
        &                  (posInBuffer-1)*nScalars+varPos(QQ)  )

      fTmp((iElem-1)*QQ+varPos(1):(iElem-1)*QQ+varPos(QQ)) = fTmp_loc

      ! compute velocity on Fluid element
      call derVarPos%velFromState( state  = fTmp_loc,     &
        &                          iField = iField,       &
        &                          nElems = 1,            &
        &                          varSys = varSys,       &
        &                          layout = layout,       &
        &                          res    = vel(:, iElem) )
    end do

    ! compute equilibrium on Boundary element
    call derVarPos%equilFromMacro( density  = rhoB,                  &
      &                            velocity = vel,                   &
      &                            iField   = iField,                &
      &                            nElems   = globBC%nElems(iLevel), &
      &                            varSys   = varSys,                &
      &                            layout   = layout,                &
      &                            res      = fEq_Bnd                )

    ! compute equilibrium on Fluid element
    call derVarPos%equilFromState( state  = fTmp,                  &
      &                            iField = iField,                &
      &                            nElems = globBC%nElems(iLevel), &
      &                            varSys = varSys,                &
      &                            layout = layout,                &
      &                            res    = fEq                    )

    do iElem = 1, globBC%nElems(iLevel)
      elemPos = globBC%elemLvl(iLevel)%elem%val( iElem )
      ! omega is STfun so get omega for current element
      omega = fieldProp%fluid%viscKine%omLvl(iLevel)%val(elemPos)

      do iDir = 1, layout%fStencil%QQN
        ! Write the values
        ! The bitmask points into the incoming direction into the flow domain,
        !  which actually we want to update
        ! * For PUSH, we write to the incoming position,
        !   as the kernel reads it from there without propagation.
        ! * For PULL, we need to write to the inverse direction, as the kernel
        !   performs a bounce back before reading it.
        !   However, this bounced back direction actually comes from the
        !   non-existent boundary element and would point into the incoming
        !   direction, so the value has to be treated and set as if it points
        !   really into the incoming direction.
        if( globBC%elemLvl(iLevel)%bitmask%val( iDir, iElem )) then
          offset = (iElem-1)*QQ
          ! Depending on PUSH or pull, use + or - for cxDir, because directions
          ! are inverted
          state(  neigh (( idir-1)* nsize+ elempos)+( ifield-1)* qq+ nscalars*0) &
            & = fEq_Bnd( offset+iDir )                                         &
            & + (1.0_rk-omega)*(fTmp(offset+iDir) - fEq(offset+iDir))
        end if ! bitMask
      end do ! iDir
    end do ! iElem

  end subroutine press_neq
! ****************************************************************************** !


! ****************************************************************************** !
  !> Inlet Velocity Equilibrium type boundary conditions with mass flow rate
  !! as input
  !!
  !! The incoming densities are set to the equilibrium distribution with
  !! macroscopic values:
  !! - density is extrapolated (0 or 1st order) from the flow domain
  !!   + 0-order: so we simply use the density of the boundary element itself
  !! - velocity is used as defined in the configuration file
  !!
  !! Usage
  !! -----
  !!```lua
  !!boundary_condition = {
  !! { label = 'inlet',
  !!    kind = 'mfr_eq',
  !!    massflowrate = 'mfr' }
  !!}
  !!variable = {
  !!  name = 'mfr',
  !!  ncomponents = 1,
  !!  vartype = 'st_fun',
  !!  st_fun = 0.1
  !!}
  !!```
  !!
  !! This subroutine's interface must match the abstract interface definition
  !! [[boundaryRoutine]] in bc/[[mus_bc_header_module]].f90 in order to be
  !! callable via [[boundary_type:fnct]] function pointer.
  subroutine mfr_eq( me, state, bcBuffer, globBC, levelDesc, tree, nSize,  &
    &                iLevel, sim_time, neigh, layout, fieldProp, varPos,   &
    &                nScalars, varSys, derVarPos, physics, iField, mixture )
    ! -------------------------------------------------------------------- !
    !> global boundary type
    class(boundary_type) :: me
    !> Current state vector of iLevel
    real(kind=rk), intent(inout) :: state(:)
    !> size of state array ( in terms of elements )
    integer, intent(in) :: nSize
    !> state values of boundary elements of all fields of iLevel
    real(kind=rk), intent(in) :: bcBuffer(:)
    !> iLevel descriptor
    type(tem_levelDesc_type), intent(in) :: levelDesc
    !> Treelm Mesh
    type(treelmesh_type), intent(in) :: tree
    !> fluid parameters and properties
    type(mus_field_prop_type), intent(in) :: fieldProp
    !> stencil layout information
    type(mus_scheme_layout_type), intent(in) :: layout
    !> the level On which this boundary was invoked
    integer, intent(in) :: iLevel
    !> connectivity array corresponding to state vector
    integer, intent(in) :: neigh(:)
    !> global time information
    type(tem_time_type), intent(in)  :: sim_time
    !> pointer to field variable in the state vector
    integer, intent(in) :: varPos(:)
    !> number of Scalars in the scheme var system
    integer, intent(in) :: nScalars
    !> scheme variable system
    type(tem_varSys_type), intent(in) :: varSys
    !> position of derived quantities in varsys
    type(mus_derVarPos_type), intent(in) :: derVarPos
    !> scheme global boundary type
    type(glob_boundary_type), intent(in) :: globBC
    !> scheme global boundary type
    type(mus_physics_type), intent(in) :: physics
    !> current field
    integer, intent(in) :: iField
    !> mixture info
    type(mus_mixture_type), intent(in) :: mixture
    ! -------------------------------------------------------------------- !
    real(kind=rk) :: massFlowRate(globBC%nElems(iLevel))
    real(kind=rk) :: massFlowRateToVel
    real(kind=rk) :: rho(globBC%nElems(iLevel))
    real(kind=rk) :: velocity(3,globBC%nElems(iLevel))
    integer :: iELem, iDir, iLvl, QQ, elemPos
    real(kind=rk) :: fEq( layout%fStencil%QQ*globBC%nElems(iLevel) )
    reaL(kind=rk) :: area
    integer :: posInBuffer
    integer :: bcMfr_pos
    ! ---------------------------------------------------------------------------
    QQ = layout%fStencil%QQ

    ! position of mass flow rate spacetime function variable in varSys
    bcMfr_pos = me%bc_states%massFlowRate%varPos
    ! get mass flow rate
    call varSys%method%val(bcMfr_pos)%get_valOfIndex( &
      & varSys  = varSys,                             &
      & time    = sim_time,                           &
      & iLevel  = iLevel,                             &
      & idx     = me%bc_states%massFlowRate           &
      &           %pntIndex%indexLvl(iLevel)          &
      &           %val(1:globBC%nElems(iLevel)),      &
      & nVals   = globBC%nElems(iLevel),              &
      & res     = massFlowRate                        )

    ! convertion factor to convert mass flow rate to velocity
    ! vel_lattice = massFlowRate/(density*areaOfCrossSection)
    ! area of cross section = sum(nElems(iLevel)*dx(iLevel)**2)
    area = 0.0_rk
    do iLvl = tree%global%minLevel, tree%global%maxLevel
      area = area + globBC%nElems_totalLevel(iLvl) * physics%dxLvl(iLvl)**2
    end do
    massFlowRateToVel = 1.0_rk / ( physics%rho0 * area ) &
      &               / physics%fac( iLevel )%vel

    ! If physical quantities are given, transform to lattice units by division
    ! with the conversion factor
!write(*,*) 'velocity_P ', massFlowRateToVel, 1.0/massFlowRateToVel
!write(*,*) 'velocity_L ', massFlowRateToVel

    ! Calculate the density of current element
    do iElem = 1, globBC%nElems(iLevel)
      posInBuffer = globBC%elemLvl( iLevel )%posInBcElemBuf%val( iElem )
      rho(iElem) = sum(bcBuffer( (posInBuffer-1)*nScalars+varPos(1) :  &
        &                        (posInBuffer-1)*nScalars+varPos(QQ) ) )
    end do

    ! compute velocity from massFlowRate
    do iElem = 1, globBC%nElems(iLevel)
      velocity(:,iElem) = massFlowRate(iElem) * massFlowRateToVel              &
        & * layout%fStencil%cxDirRK(:, globBC%elemLvl(iLevel)%                 &
        &                                                  normalInd%val(iElem))
    end do

    ! compute equilibrium
    call derVarPos%equilFromMacro( density  = rho,                   &
      &                            velocity = velocity,              &
      &                            iField   = iField,                &
      &                            nElems   = globBC%nElems(iLevel), &
      &                            varSys   = varSys,                &
      &                            layout   = layout,                &
      &                            res      = fEq                    )

    do iElem = 1, globBC%nElems(iLevel)
      ! bc element position in total
      elemPos = globBC%elemLvl(iLevel)%elem%val( iElem )
      do iDir = 1, layout%fStencil%QQN
        ! Write the values
        ! The bitmask points into the incoming direction into the flow domain,
        !  which actually we want to update
        ! * For PUSH, we write to the incoming position,
        !   as the kernel reads it from there without propagation.
        ! * For PULL, we need to write to the inverse direction, as the kernel
        !   performs a bounce back before reading it.
        !   However, this bounced back direction actually comes from the
        !   non-existent boundary element and would point into the incoming
        !   direction, so the value has to be treated and set as if it points
        !   really into the incoming direction.
        if( globBC%elemLvl(iLevel)%bitmask%val( iDir, iElem )) then
          ! Depending on PUSH or pull, use + or - for cxDir, because directions
          ! are inverted
          state(  neigh (( idir-1)* nsize+ elempos)+( ifield-1)* qq+ nscalars*0) &
            & = fEq( (iElem-1)*QQ + iDir )
        end if ! bitMask
      end do ! iDir
    end do ! iElem

  end subroutine mfr_eq
! ****************************************************************************** !


! ****************************************************************************** !
  !> Outlet Pressure Equilibrium type boundary conditions
  !!
  !! The incoming densities are set to the equilibrium distribution with
  !! macroscopic values used from:
  !! * Velocity: Linear extrapolation of velocity of precollision of current
  !!   time step
  !! * Pressure as defined in the configuration file
  !!
  !! Usage:\n
  !!```lua
  !!  boundary_condition = {
  !!    { label = 'outlet',
  !!       kind = 'pressure_eq',
  !!       pressure = 'pressure_out',
  !!    }
  !!  }
  !! variable = {
  !!   name = 'pressure_out',
  !!   ncomponents = 1,
  !!   vartype = 'st_fun',
  !!   st_fun = 1.0
  !! }
  !!```
  !!
  !! This is described in the paper:\n
  !! S. Izquierdo, P. Martinez-Lera, and N. Fueyo,\n
  !! Analysis of open boundary effects in unsteady lattice Boltzmann simulations
  !! Comput. Math. Applications, vol. 58, no. 5, pp. 914–921, Sep. 2009.\n
  !!
  !! The velocity is extrapolated (1st order) from the flow domain
  !! See equation (6)
  !!
  !! DISCLAIMER: This BC requires pre-collision PDF from 1st and 2nd neighbor
  !! which are not available if those neighbors are boundary halo on different
  !! process. So, if you encounter a problem with this BC, please check the
  !! domain distribution to verify whether 2 neighbors of boundary elements
  !! are in one process. The domain distribution can be changed by shifting the
  !! boundary by one element or by changing the origin of bounding cube or
  !! by using different number of MPI processes.
  !!
  !! \todo 20200410, KM: Check if this BC works if neighbufferPre is replaced
  !! by neighbufferPost
  !!
  !! This subroutine's interface must match the abstract interface definition
  !! [[boundaryRoutine]] in bc/[[mus_bc_header_module]].f90 in order to be
  !! callable via [boundary_type:fnct]] function pointer.
  subroutine pressure_eq( me, state, bcBuffer, globBC, levelDesc, tree, nSize, &
    &                     iLevel, sim_time, neigh, layout, fieldProp, varPos,  &
    &                     nScalars, varSys, derVarPos, physics, iField,        &
    &                     mixture                                              )
    ! -------------------------------------------------------------------- !
    !> global boundary type
    class(boundary_type) :: me
    !> Current state vector of iLevel
    real(kind=rk), intent(inout) :: state(:)
    !> size of state array ( in terms of elements )
    integer, intent(in) :: nSize
    !> state values of boundary elements of all fields of iLevel
    real(kind=rk), intent(in) :: bcBuffer(:)
    !> iLevel descriptor
    type(tem_levelDesc_type), intent(in) :: levelDesc
    !> Treelm Mesh
    type(treelmesh_type), intent(in) :: tree
    !> fluid parameters and properties
    type(mus_field_prop_type), intent(in) :: fieldProp
    !> stencil layout information
    type(mus_scheme_layout_type), intent(in) :: layout
    !> the level On which this boundary was invoked
    integer, intent(in) :: iLevel
    !> connectivity array corresponding to state vector
    integer, intent(in) :: neigh(:)
    !> global time information
    type(tem_time_type), intent(in)  :: sim_time
    !> pointer to field variable in the state vector
    integer, intent(in) :: varPos(:)
    !> number of Scalars in the scheme var system
    integer, intent(in) :: nScalars
    !> scheme variable system
    type(tem_varSys_type), intent(in) :: varSys
    !> position of derived quantities in varsys
    type(mus_derVarPos_type), intent(in) :: derVarPos
    !> scheme global boundary type
    type(glob_boundary_type), intent(in) :: globBC
    !> scheme global boundary type
    type(mus_physics_type), intent(in) :: physics
    !> current field
    integer, intent(in) :: iField
    !> mixture info
    type(mus_mixture_type), intent(in) :: mixture
    ! -------------------------------------------------------------------- !
    real(kind=rk) :: fEq( layout%fStencil%QQ*globBC%nElems(iLevel) )
    ! need to compute precollision velocity of neighbor element
    ! Density on boundary element
    real(kind=rk) :: rho( globBC%nElems(iLevel))
    ! Velocity on boundary element
    real(kind=rk) :: velocity(3,globBC%nElems(iLevel))
    ! Vel on frst neighbor element
    real(kind=rk) :: uxB_1(3*globBC%nElems(iLevel))
    ! Vel on scnd neighbor element
    real(kind=rk) :: uxB_2(3*globBC%nElems(iLevel))
    integer :: iDir, iElem, QQ, elemPos
    integer :: bcPress_pos
    ! ---------------------------------------------------------------------------
    QQ = layout%fStencil%QQ

    ! position of boundary pressure in varSys
    bcPress_pos = me%bc_states%pressure%varPos
    ! get pressure variable from spacetime function
    call varSys%method%val(bcPress_pos)%get_valOfIndex( &
      & varSys  = varSys,                               &
      & time    = sim_time,                             &
      & iLevel  = iLevel,                               &
      & idx     = me%bc_states%pressure                 &
      &           %pntIndex%indexLvl(iLevel)            &
      &           %val(1:globBC%nElems(iLevel)),        &
      & nVals   = globBC%nElems(iLevel),                &
      & res     = rho                                   )

    ! convert physical pressure to LB density
    rho = rho / physics%fac( iLevel )%press * cs2inv + rho0

    ! KM: velocity is taken from precollision of current
    ! Calculate velocity of two neighbors
    call derVarPos%velFromState(                           &
      &     state  = me%neigh(iLevel)%neighBufferPre(1,:), &
      &     iField = iField,                               &
      &     nElems = globBC%nElems(iLevel),                &
      &     varSys = varSys,                               &
      &     layout = layout,                               &
      &     res    = uxB_1                                 )
    call derVarPos%velFromState(                           &
      &     state  = me%neigh(iLevel)%neighBufferPre(2,:), &
      &     iField = iField,                               &
      &     nElems = globBC%nElems(iLevel),                &
      &     varSys = varSys,                               &
      &     layout = layout,                               &
      &     res    = uxB_2                                 )

    ! caluclate my velocity by interpolation
    do iElem = 1, globBC%nElems(iLevel)
      velocity(:,iElem) =   1.5_rk * uxB_1((iElem-1)*3+1:iElem*3) &
        &                  - 0.5_rk * uxB_2((iElem-1)*3+1:iElem*3)
    end do

    ! compute equilibrium
    call derVarPos%equilFromMacro( density  = rho,                   &
      &                            velocity = velocity,              &
      &                            iField   = iField,                &
      &                            nElems   = globBC%nElems(iLevel), &
      &                            varSys   = varSys,                &
      &                            layout   = layout,                &
      &                            res      = fEq                    )

    ! Run over all elements in list for this boundary
    do iElem = 1, globBC%nElems( iLevel )
      elemPos = globBC%elemLvl(iLevel)%elem%val( iElem )
      do iDir = 1, layout%fStencil%QQN
        ! the bitmask points into the incoming direction, which must
        ! be updated. As the kernel reads this value with FETCH,
        ! we store it to the FETCH position.
        ! This results for
        ! - PUSH in storing the density to where the bitmask
        !   points, i.e. incoming
        ! - PULL storing the density pointing outwards, but as
        !   this density is being bounced back before collision in
        !   the kernel, this is the correct position
        if( globBC%elemLvl(iLevel)%bitmask%val( iDir, iElem )) then
          ! Now assign the values
state(  neigh (( idir-1)* nsize+ elempos)+( ifield-1)* qq+ nscalars*0) = &
          & fEq( (iElem-1)*QQ + iDir )
        end if ! bitMask
      end do ! iDir
    end do ! iElem

  end subroutine pressure_eq
! ****************************************************************************** !


! ****************************************************************************** !
  !> Outlet Pressure do-nothing boundary is the open boundary condition for
  !! incompressible model.
  !! This BC sets reference density at boundary so the background pressure
  !! should be zero.
  !! Here, the normal velocity is extrapolated from 1st fluid node and
  !! tangential velocity is extrapolated from 2nd fluid node in normal
  !! direction.
  !! Algorithm used in this boundary condition:
  !! fEq(1,u) and fEq(rho,u) are computed using macroscopic values from current
  !! element.
  !! In fNeq, post-collision of current time step is used for normal direction.
  !!
  !! This is taken from the paper:
  !! M. Junk and Z. Yang,
  !! Asymptotic Analysis of Lattice Boltzmann Outflow Treatments
  !! Commun. Comput. Phys., pp. 1–11, 2011.
  !!
  !! This subroutine's interface must match the abstract interface definition
  !! [[boundaryRoutine]] in bc/[[mus_bc_header_module]].f90 in order to be
  !! callable via [[boundary_type:fnct]] function pointer.
  subroutine outlet_dnt( me, state, bcBuffer, globBC, levelDesc, tree, nSize,  &
    &                    iLevel, sim_time, neigh, layout, fieldProp, varPos,   &
    &                    nScalars, varSys, derVarPos, physics, iField, mixture )
    ! -------------------------------------------------------------------- !
    !> global boundary type
    class(boundary_type) :: me
    !> Current state vector of iLevel
    real(kind=rk), intent(inout) :: state(:)
    !> size of state array ( in terms of elements )
    integer, intent(in) :: nSize
    !> state values of boundary elements of all fields of iLevel
    real(kind=rk), intent(in) :: bcBuffer(:)
    !> iLevel descriptor
    type(tem_levelDesc_type), intent(in) :: levelDesc
    !> Treelm Mesh
    type(treelmesh_type), intent(in) :: tree
    !> fluid parameters and properties
    type(mus_field_prop_type), intent(in) :: fieldProp
    !> stencil layout information
    type(mus_scheme_layout_type), intent(in) :: layout
    !> the level On which this boundary was invoked
    integer, intent(in) :: iLevel
    !> connectivity array corresponding to state vector
    integer, intent(in) :: neigh(:)
    !> global time information
    type(tem_time_type), intent(in)  :: sim_time
    !> pointer to field variable in the state vector
    integer, intent(in) :: varPos(:)
    !> number of Scalars in the scheme var system
    integer, intent(in) :: nScalars
    !> scheme variable system
    type(tem_varSys_type), intent(in) :: varSys
    !> position of derived quantities in varsys
    type(mus_derVarPos_type), intent(in) :: derVarPos
    !> scheme global boundary type
    type(glob_boundary_type), intent(in) :: globBC
    !> scheme global boundary type
    type(mus_physics_type), intent(in) :: physics
    !> current field
    integer, intent(in) :: iField
    !> mixture info
    type(mus_mixture_type), intent(in) :: mixture
    ! -------------------------------------------------------------------- !
    ! fTmp uses AOS layout
    real(kind=rk) :: fTmp( nScalars*globBC%nElems(iLevel) )
    ! fEq, fEq0 use AOS layout
    real(kind=rk) :: fEq ( layout%fStencil%QQ*globBC%nElems(iLevel) )
    real(kind=rk) :: fEq0( layout%fStencil%QQ*globBC%nElems(iLevel) )
    real(kind=rk) :: fNEq
    real(kind=rk) :: rho( globBC%nElems(iLevel) ) ! my density
    real(kind=rk) :: uxB(3)  ! Velocity on boundary element
    real(kind=rk) :: uxB_n(3), uxB_n2(3), uxB_t(3), normal(3)
    integer :: normalDir
    real(kind=rk) :: velocity(3,globBC%nElems(iLevel)) ! my velocity
    real(kind=rk) :: uxF_2(globBC%nElems(iLevel)* 3) ! first  neigh velocity
    real(kind=rk) :: omega, visc, facFNeq
    integer :: iDir, invDir, iElem, QQ, elemPos, posInBuffer
    ! ---------------------------------------------------------------------------

    QQ    = layout%fStencil%QQ

    do iElem = 1, globBC%nElems(iLevel)

      posInBuffer = globBC%elemLvl( iLevel )%posInBcElemBuf%val( iElem )
      fTmp( (iElem-1)*nScalars+1: (iElem-1)*nScalars+QQ ) &
        &       = bcBuffer( (posInBuffer-1)*nScalars+varPos(1) : &
        &                   (posInBuffer-1)*nScalars+varPos(QQ)  )

      rho( iElem ) = sum( fTmp( (iElem-1)*QQ+1:iElem*QQ ) )

      ! Get my velocity
      call derVarPos%velFromState( state  = fTmp( (iElem-1)*nScalars+1:iElem*nScalars ), &
        &                          iField = iField,            &
        &                          nElems = 1,                 &
        &                          varSys = varSys,            &
        &                          layout = layout,            &
        &                          res    = velocity(:,iElem) )
    end do ! iElem


    ! compute my equilibrium using rho
    call derVarPos%equilFromMacro( density  = rho,                   &
      &                            velocity = velocity,              &
      &                            iField   = iField,                &
      &                            nElems   = globBC%nElems(iLevel), &
      &                            varSys   = varSys,                &
      &                            layout   = layout,                &
      &                            res      = fEq                    )

    ! Calculate my equilibrium distribution using rho0
    rho = 1.0_rk

    ! compute my equilibrium using rho0
    call derVarPos%equilFromMacro( density  = rho,                   &
      &                            velocity = velocity,              &
      &                            iField   = iField,                &
      &                            nElems   = globBC%nElems(iLevel), &
      &                            varSys   = varSys,                &
      &                            layout   = layout,                &
      &                            res      = fEq0                   )

    ! Get my neighbors velocity
    !KM: velocity is taken from precollision of current
    call derVarPos%velFromState(                            &
      &     state  = me%neigh(iLevel)%neighBufferPost(1,:), &
      &     iField = iField,                               &
      &     nElems = globBC%nElems(iLevel),                &
      &     varSys = varSys,                               &
      &     layout = layout,                               &
      &     res    = uxF_2                                 )

    ! Run all elements in list for this boundary
    do iElem = 1, globBC%nElems(iLevel)
      elemPos = globBC%elemLvl(iLevel)%elem%val( iElem )

      ! Extrapolate normal velocity from first neighbor
      ! and tangential velocity from second neighbor as per paper
      normalDir = globBC%elemLvl( iLevel )%normalInd%val( iElem )
      normal = layout%fStencil%cxDirRK(:, normalDir)
      ! Extract normal velocity
      uxB_n = dot_product( velocity(:, iElem), abs(normal) ) &
        &   * abs(normal)
      uxB_n2 = dot_product( uxF_2((iElem-1)*3+1 : iElem*3),   &
        &                   abs(normal)                     ) &
        &    * abs(normal)
      ! Extract tangential velocity
      uxB_t = uxF_2((iElem-1)*3+1 : iElem*3) - uxB_n2
      ! velocity on boundary
      uxB = uxB_n + uxB_t


      ! Treat all directions, but actually we only want to treat
      ! non-normal directions
      ! Equation (3.2b)
      do iDir = 1, layout%fStencil%QQN
        ! the bitmask points into the incoming direction, which must
        ! be updated. As the kernel reads this value with FETCH,
        ! we store it to the FETCH position.
        ! This results for
        ! - PUSH in storing the density to where the bitmask
        !   points, i.e. incoming
        ! - PULL storing the density pointing outwards, but as
        !   this density is being bounced back before collision in
        !   the kernel, this is the correct position
        if ( globBC%elemLvl(iLevel)%bitmask%val( iDir, iElem )) then
          invDir = layout%fStencil%cxDirInv(iDir)
          ! update the incoming velocity direction
          state(  neigh (( idir-1)*nsize+ elempos)+( ifield-1)* qq+ nscalars*0)  = &
         ! This bc is stable only using FETCH with inverse direction of bitmask
         ! For PULL this means, get the outgoing one from upstream neighbor
         ! For PUSH this means, get the outgoing one from local element
         ! KM: Using FETCH gives stable result
         ! @todo test for different testcase for stability
            & state(neigh((invdir-1)*nsize+ elempos)+( ifield-1)* qq+ nscalars*0)&
            & + 6._rk*layout%weight( iDir )                  &
            !extrapolated velocity is used here to remove fluctuations
            & * ( layout%fStencil%cxDirRK( 1, iDir )*uxB(1)  &
            & +   layout%fStencil%cxDirRK( 2, iDir )*uxB(2)  &
            & +   layout%fStencil%cxDirRK( 3, iDir )*uxB(3)  )
        end if ! bitMask
      end do ! iDir

      ! omega is STfun so get omega for current element
      omega = fieldProp%fluid%viscKine%omLvl(iLevel)%val(elemPos)
      visc  = fieldProp%fluid%viscKine%dataOnLvl(iLevel)%val(elemPos)
      facFNeq = visc * omega - 1._rk

      ! then overwrite the normal direction with special treatment
      ! Equation (3.2a)
      iDir = globBC%elemLvl( iLevel )%normalInd%val( iElem )
      invDir = layout%fStencil%cxDirInv(iDir)
      fNEq =  fTmp( (iElem-1)*QQ + invDir ) - fEq( (iElem-1)*QQ + invDir )
      state(  neigh(( idir-1)*nsize+ elempos)+( ifield-1)* qq+ nscalars*0)  = &
        & fEq0( iDir+(iElem-1)*QQ ) - facFNeq * fNeq
    end do

  end subroutine outlet_dnt
! ****************************************************************************** !


  ! ****************************************************************************** !
  !> Outlet Pressure extrapolation boundary.
  !! qVal for this boundary must be 0.0
  !!
  !! This is taken from the paper:
  !! M. Junk and Z. Yang,
  !! Asymptotic Analysis of Lattice Boltzmann Outflow Treatments
  !! Commun. Comput. Phys., pp. 1-11, 2011.
  !!
  !! This boundary condition prescribes density = 1.0 at outlet by forcing the
  !! equlibirium part of \( f_i \) in normal direction (equation 3.13a):
  !! \[
  !!  f_i(n+1,j_0) = \left[ F^{eq}_{i}(1,u) + f^{neq}_i \right] (n,j_0)
  !! \]
  !! For non-normal direction, \( f_i \) is extrapolated from neighrbors by
  !! (equation 3.13b):
  !! \[
  !!  f_i(n+1,j_0) = 2f_i(n+1,j_0-n) - f_i(n+1,j_0-2n)
  !! \]
  !!
  !! Usage:
  !!----
  !!```lua
  !!  boundary_condition = {
  !!    { label = 'outlet',
  !!       kind = 'pressure_expol',
  !!       pressure = 'pressure_out',
  !!    }
  !!  }
  !!  variable = {
  !!    name = 'pressure_out',
  !!    ncomponents = 1,
  !!    vartype = 'st_fun',
  !!    st_fun = 1.0
  !!  }
  !!```
  !! DISCLAIMER: This BC requires pre-collision PDF from 1st and 2nd neighbor
  !! which are not available if those neighbors are boundary halo on different
  !! process. So, if you encounter a problem with this BC, please check the
  !! domain distribution to verify whether 2 neighbors of boundary elements
  !! are in one process. The domain distribution can be changed by shifting the
  !! boundary by one element or by changing the origin of bounding cube or
  !! by using different number of MPI processes.
  !!
  !! This subroutine's interface must match the abstract interface definition
  !! [[boundaryRoutine]] in bc/[[mus_bc_header_module]].f90 in order to be
  !! callable via [[boundary_type:fnct]] function pointer.
  subroutine pressure_expol( me, state, bcBuffer, globBC, levelDesc, tree,   &
    &                        nSize, iLevel, sim_time, neigh, layout,         &
    &                        fieldProp, varPos, nScalars, varSys, derVarPos, &
    &                        physics, iField, mixture                        )
    ! -------------------------------------------------------------------- !
    !> global boundary type
    class(boundary_type) :: me
    !> Current state vector of iLevel
    real(kind=rk), intent(inout) :: state(:)
    !> size of state array ( in terms of elements )
    integer, intent(in) :: nSize
    !> state values of boundary elements of all fields of iLevel
    real(kind=rk), intent(in) :: bcBuffer(:)
    !> iLevel descriptor
    type(tem_levelDesc_type), intent(in) :: levelDesc
    !> Treelm Mesh
    type(treelmesh_type), intent(in) :: tree
    !> fluid parameters and properties
    type(mus_field_prop_type), intent(in) :: fieldProp
    !> stencil layout information
    type(mus_scheme_layout_type), intent(in) :: layout
    !> the level On which this boundary was invoked
    integer, intent(in) :: iLevel
    !> connectivity array corresponding to state vector
    integer, intent(in) :: neigh(:)
    !> global time information
    type(tem_time_type), intent(in)  :: sim_time
    !> pointer to field variable in the state vector
    integer, intent(in) :: varPos(:)
    !> number of Scalars in the scheme var system
    integer, intent(in) :: nScalars
    !> scheme variable system
    type(tem_varSys_type), intent(in) :: varSys
    !> position of derived quantities in varsys
    type(mus_derVarPos_type), intent(in) :: derVarPos
    !> scheme global boundary type
    type(glob_boundary_type), intent(in) :: globBC
    !> scheme global boundary type
    type(mus_physics_type), intent(in) :: physics
    !> current field
    integer, intent(in) :: iField
    !> mixture info
    type(mus_mixture_type), intent(in) :: mixture
    ! -------------------------------------------------------------------- !
    real(kind=rk) :: fTmp( nScalars * globBC%nElems(iLevel) )
    real(kind=rk) :: fTmp_1 ! first  neighbor pdf
    real(kind=rk) :: fTmp_2 ! second neighbor pdf
    real(kind=rk) :: fEq(  layout%fStencil%QQ * globBC%nElems(iLevel) )
    real(kind=rk) :: fEq0( layout%fStencil%QQ * globBC%nElems(iLevel) )
    real(kind=rk) :: rho( globBC%nElems(iLevel) ), inv_rho_phy
    real(kind=rk) :: velocity( 3, globBC%nElems(iLevel) )
    real(kind=rk) :: vel_res( globBC%nElems(iLevel)*3 )
    integer :: iDir, iElem, QQ, elemPos, invDir, posInBuffer
    integer :: iLink, bcPress_pos
    logical :: axisNormal
    ! ---------------------------------------------------------------------------
!write(dbgUnit(1),*) 'bclabel ', trim(me%label)
!write(dbgUnit(1),*) 'sim time ', sim_time%iter
    ! Here starts the copied part of pressure_expol
    QQ = layout%fStencil%QQ
    inv_rho_phy = 1.0_rk / physics%fac(iLevel)%press * cs2inv

    ! calculate rho and velocity from PDF
    do iElem = 1, globBC%nElems( iLevel )
      posInBuffer = globBC%elemLvl( iLevel )%posInBcElemBuf%val( iElem )
      fTmp( (iElem-1)*nScalars+1: (iElem-1)*nScalars+QQ ) &
        &       = bcBuffer( (posInBuffer-1)*nScalars+varPos(1) : &
        &                   (posInBuffer-1)*nScalars+varPos(QQ) )
      rho( iElem ) = sum( fTmp( (iElem-1)*nScalars+1 : (iElem-1)*nScalars+QQ ) )
    end do

    ! Get my velocity
    ! All state variables have to be put inside fTmp
    ! (e.g. pdf and omegas for nonNewtonian fluid)
    call derVarPos%velFromState( state  = fTmp(:),               &
      &                          iField = iField,                &
      &                          nElems = globBC%nElems(iLevel), &
      &                          varSys = varSys,                &
      &                          layout = layout,                &
      &                          res    = vel_res(:)             )

    ! convert  velocity arrangement
    do iElem = 1, globBC%nElems( iLevel )
      velocity( 1, iElem ) = vel_res( (iElem-1)*3 + 1 )
      velocity( 2, iElem ) = vel_res( (iElem-1)*3 + 2 )
      velocity( 3, iElem ) = vel_res( (iElem-1)*3 + 3 )
    end do

    ! compute my equilibrium (fEq) using calculated rho
    call derVarPos%equilFromMacro( density  = rho,                   &
      &                            velocity = velocity,              &
      &                            iField   = iField,                &
      &                            nElems   = globBC%nElems(iLevel), &
      &                            varSys   = varSys,                &
      &                            layout   = layout,                &
      &                            res      = fEq                    )

    ! position of boundary pressure in varSys
    bcPress_pos = me%bc_states%pressure%varPos
    ! get pressure variable from spacetime function
    call varSys%method%val(bcPress_pos)%get_valOfIndex( &
      & varSys  = varSys,                               &
      & time    = sim_time,                             &
      & iLevel  = iLevel,                               &
      & idx     = me%bc_states%pressure                 &
      &           %pntIndex%indexLvl(iLevel)            &
      &           %val(1:globBC%nElems(iLevel)),        &
      & nVals   = globBC%nElems(iLevel),                &
      & res     = rho                                   )

    ! convert physical pressure to LB density
    rho = rho * inv_rho_phy + rho0

    ! compute my equilibrium (fEq0)using user defined rho
    call derVarPos%equilFromMacro( density  = rho,                   &
      &                            velocity = velocity,              &
      &                            iField   = iField,                &
      &                            nElems   = globBC%nElems(iLevel), &
      &                            varSys   = varSys,                &
      &                            layout   = layout,                &
      &                            res      = fEq0                   )
    ! Here ends the copied part of pressure_expol

    ! linkwise treatment for links that have boundary
    do iLink = 1, me%links(iLevel)%nVals
      iElem = me%outletExpol(iLevel)%iElem(iLink)

      ! if any of the neighbor has property qVal then fall back to simple
      ! equilibrium boundary condition
      ! elements that has property qVal
      ! @todo: create a list for these special links.
      if( btest( levelDesc%property(                        &
        &        me%neigh(iLevel)%posInState(1, iElem)), prp_hasQVal )   &
        & .or. btest( levelDesc%property(                   &
        &        me%neigh(iLevel)%posInState(2, iElem)), prp_hasQVal ) ) &
        & then
        fTmp_1 = 0.0_rk
        fTmp_2 = 0.0_rk
        state( me%links(iLevel)%val(iLink) ) &
          &  = fEq0( me%outletExpol(iLevel)%statePos(iLink) )
!write(dbgUnit(1),*) 'fEq0 ', fEq0( me%outletExpol(iLevel)%statePos(iLink) )
      else
        fTmp_1 = me%neigh(iLevel)%neighBufferPre(1, me%outletExpol(iLevel) &
          &                                               %statePos(iLink) )
        fTmp_2 = me%neigh(iLevel)%neighBufferPre(2, me%outletExpol(iLevel) &
          &                                               %statePos(iLink) )
        state( me%links(iLevel)%val(iLink) ) &
          & = 1.5_rk*fTmp_1 - 0.5_rk*fTmp_2

! debug output ------------------------------------------------
! debug output ------------------------------------------------

      end if
    end do ! iLink, linkwise treatment ends here

    ! then overwrite the normal direction with special treatment
    ! Equation (3.13a)
    do iElem = 1, globBC%nElems(iLevel)
      elemPos = globBC%elemLvl(iLevel)%elem%val( iElem )
      iDir   = globBC%elemLvl( iLevel )%normalInd%val( iElem )

      axisNormal = ( ( abs(layout%fStencil%cxDir( 1, iDir ))        &
        &            + abs(layout%fStencil%cxDir( 2, iDir ))        &
        &            + abs(layout%fStencil%cxDir( 3, iDir )) ) == 1 )

      if (axisNormal) then
        invDir = layout%fStencil%cxDirInv( iDir )
        state(  neigh (( idir-1)* nsize+ elempos)+( ifield-1)* qq+ nscalars*0 &
          & )  = &
          ! KM:we use fEq0 since this equilibrium has been computed from
          ! defined density at boundary
          & fEq0((iElem-1)*QQ+iDir) + (fTmp( (iElem-1)*nScalars+invDir ) &
          &                       -  fEq( (iElem-1)*QQ+invDir ))
      end if
    end do ! iElem

! debug output ------------------------------------------------
! debug output ------------------------------------------------

  end subroutine pressure_expol
  ! ****************************************************************************** !


  ! ****************************************************************************** !
  !> Inlet Velocity Bounce Back boundary condition with qvalues
  !!
  !! This is taken from the paper:
  !! [1] S. Izquierdo and N. Fueyo, "Characteristic  non-reflecting
  !! boundary  conditions  for  open  boundaries  in  lattice
  !!  Boltzmann  methods," Physical Review E, vol. 78, no. 46707, 2008
  !! The incoming densities \( \bar{f_{alpha}} \) are reconstructed from the
  !! already reflected densities, enforcing the given velocity me\%ubb\%velocity.
  !! [2] M. Junk and Z. Yang. One-point boundary condition for the
  !! lattice Boltzmann method, Physical Review E, vol 72, issue 8, year 2005.
  !!
  !! This boundary condition has error of 1st order pressure and 2nd order velocity
  !! only if wall is located exactly at q=1/2 else the accuracy of both pressure
  !! and velocity reduce by order 1.
  !!
  !! Usage
  !! -----
  !!```lua
  !!boundary_condition = {
  !!  { label = 'inlet',
  !!    kind = 'velocity_bounceback',
  !!    velocity = 'inlet_vel',
  !!  }
  !!}
  !!variable = {
  !!  name = 'inlet_vel',
  !!  ncomponents = 3,
  !!  vartype = 'st_fun',
  !!  st_fun = {0.06, 0.0, 0.0}
  !!}
  !!```
  !!
  !! This subroutine's interface must match the abstract interface definition
  !! [[boundaryRoutine]] in bc/[[mus_bc_header_module]].f90 in order to be
  !! callable via [[boundary_type:fnct]] function pointer.
  subroutine velocity_bounceback_incomp( me, state, bcBuffer, globBC,         &
    &                                    levelDesc, tree, nSize, iLevel,      &
    &                                    sim_time, neigh, layout, fieldProp,  &
    &                                    varPos, nScalars, varSys, derVarPos, &
    &                                    physics, iField, mixture             )
    ! -------------------------------------------------------------------- !
    !> global boundary type
    class(boundary_type) :: me
    !> Current state vector of iLevel
    real(kind=rk), intent(inout) :: state(:)
    !> size of state array ( in terms of elements )
    integer, intent(in) :: nSize
    !> state values of boundary elements of all fields of iLevel
    real(kind=rk), intent(in) :: bcBuffer(:)
    !> iLevel descriptor
    type(tem_levelDesc_type), intent(in) :: levelDesc
    !> Treelm Mesh
    type(treelmesh_type), intent(in) :: tree
    !> fluid parameters and properties
    type(mus_field_prop_type), intent(in) :: fieldProp
    !> stencil layout information
    type(mus_scheme_layout_type), intent(in) :: layout
    !> the level On which this boundary was invoked
    integer, intent(in) :: iLevel
    !> connectivity array corresponding to state vector
    integer, intent(in) :: neigh(:)
    !> global time information
    type(tem_time_type), intent(in)  :: sim_time
    !> pointer to field variable in the state vector
    integer, intent(in) :: varPos(:)
    !> number of Scalars in the scheme var system
    integer, intent(in) :: nScalars
    !> scheme variable system
    type(tem_varSys_type), intent(in) :: varSys
    !> position of derived quantities in varsys
    type(mus_derVarPos_type), intent(in) :: derVarPos
    !> scheme global boundary type
    type(glob_boundary_type), intent(in) :: globBC
    !> scheme global boundary type
    type(mus_physics_type), intent(in) :: physics
    !> current field
    integer, intent(in) :: iField
    !> mixture info
    type(mus_mixture_type), intent(in) :: mixture
    ! -------------------------------------------------------------------- !
    integer       :: iLink, iDir
    real(kind=rk) :: fOut, eqPlus, inv_vel
    real(kind=rk) :: vel_b(me%links(iLevel)%nVals*3)
    integer :: bcVel_pos, offset
    ! ---------------------------------------------------------------------------

    inv_vel = 1.0_rk / physics%fac( iLevel )%vel

    ! position of boundary velocity in varSys
    bcVel_pos = me%bc_states%velocity%varPos
    ! Multiply transient velocity with spatial velocity. Spatial term includes
    ! maxvalue
    call varSys%method%val(bcVel_pos)%get_valOfIndex( &
      & varSys  = varSys,                             &
      & time    = sim_time,                           &
      & iLevel  = iLevel,                             &
      & idx     = me%bc_states%velocity               &
      &           %pntIndex%indexLvl(iLevel)          &
      &           %val(1:me%links(iLevel)%nVals),     &
      & nVals   = me%links(iLevel)%nVals,             &
      & res     = vel_b                               )

    !NEC$ ivdep
    !DIR$ ivdep
    do iLink = 1, me%links(iLevel)%nVals
      offset = (iLink-1)*3

      fOut = bcBuffer( me%inletUbbQVal(iLevel)%outPos(iLink) )

      ! According to the bounce back rule,
      ! Subtract the equlibrium if its computed from alpha-(inverse of
      ! bitmask) or add the equlibrium if its computed in the direction of
      ! bitmask.
      ! In the Izquiredo paper, this is rho0, the reference density
      iDir   = me%inletUbbQVal(iLevel)%iDir(iLink)
      eqPlus = layout%weight(iDir) * 6.0_rk * rho0  &
        &    * (   layout%fStencil%cxDirRK(1, iDir) * vel_b(offset+1) * inv_vel &
        &        + layout%fStencil%cxDirRK(2, iDir) * vel_b(offset+2) * inv_vel &
        &        + layout%fStencil%cxDirRK(3, iDir) * vel_b(offset+3) * inv_vel )

      state( me%links(iLevel)%val(iLink) ) = fOut + eqPlus

    end do ! iLink

  end subroutine velocity_bounceback_incomp
! ****************************************************************************** !


  ! ****************************************************************************** !
  !> Inlet Velocity Bounce Back boundary condition with qvalues for compressible
  !! flows. It is similar to velocity_bounceback except the density is extrapolated
  !! from fluid element.
  !!
  !! This subroutine's interface must match the abstract interface definition
  !! [[boundaryRoutine]] in bc/[[mus_bc_header_module]].f90 in order to be
  !! callable via [[boundary_type:fnct]] function pointer.
  subroutine velocity_bounceback( me, state, bcBuffer, globBC, levelDesc,      &
    &                             tree, nSize, iLevel, sim_time, neigh,        &
    &                             layout, fieldProp, varPos, nScalars, varSys, &
    &                             derVarPos, physics, iField, mixture          )
    ! -------------------------------------------------------------------- !
    !> global boundary type
    class(boundary_type) :: me
    !> Current state vector of iLevel
    real(kind=rk), intent(inout) :: state(:)
    !> size of state array ( in terms of elements )
    integer, intent(in) :: nSize
    !> state values of boundary elements of all fields of iLevel
    real(kind=rk), intent(in) :: bcBuffer(:)
    !> iLevel descriptor
    type(tem_levelDesc_type), intent(in) :: levelDesc
    !> Treelm Mesh
    type(treelmesh_type), intent(in) :: tree
    !> fluid parameters and properties
    type(mus_field_prop_type), intent(in) :: fieldProp
    !> stencil layout information
    type(mus_scheme_layout_type), intent(in) :: layout
    !> the level On which this boundary was invoked
    integer, intent(in) :: iLevel
    !> connectivity array corresponding to state vector
    integer, intent(in) :: neigh(:)
    !> global time information
    type(tem_time_type), intent(in)  :: sim_time
    !> pointer to field variable in the state vector
    integer, intent(in) :: varPos(:)
    !> number of Scalars in the scheme var system
    integer, intent(in) :: nScalars
    !> scheme variable system
    type(tem_varSys_type), intent(in) :: varSys
    !> position of derived quantities in varsys
    type(mus_derVarPos_type), intent(in) :: derVarPos
    !> scheme global boundary type
    type(glob_boundary_type), intent(in) :: globBC
    !> scheme global boundary type
    type(mus_physics_type), intent(in) :: physics
    !> current field
    integer, intent(in) :: iField
    !> mixture info
    type(mus_mixture_type), intent(in) :: mixture
    ! -------------------------------------------------------------------- !
    integer       :: iLink, iDir, QQ
    real(kind=rk) :: fOut, eqPlus, inv_vel, rho
    real(kind=rk) :: vel_b(me%links(iLevel)%nVals*3)
    integer :: bcVel_pos, offset, posInBuffer
    ! ---------------------------------------------------------------------------
    QQ = layout%fstencil%QQ

    inv_vel = 1.0_rk / physics%fac( iLevel )%vel

    ! position of boundary velocity in varSys
    bcVel_pos = me%bc_states%velocity%varPos
    ! Multiply transient velocity with spatial velocity. Spatial term includes
    ! maxvalue
    call varSys%method%val(bcVel_pos)%get_valOfIndex( &
      & varSys  = varSys,                             &
      & time    = sim_time,                           &
      & iLevel  = iLevel,                             &
      & idx     = me%bc_states%velocity               &
      &           %pntIndex%indexLvl(iLevel)          &
      &           %val(1:me%links(iLevel)%nVals),     &
      & nVals   = me%links(iLevel)%nVals,             &
      & res     = vel_b                               )

    !NEC$ ivdep
    !DIR$ ivdep
    do iLink = 1, me%links(iLevel)%nVals
      offset = (iLink-1)*3

      fOut = bcBuffer( me%inletUbbQVal(iLevel)%outPos(iLink) )

      ! calculate density
      posInBuffer  = me%inletUbbQVal(iLevel)%posInBuffer(iLink)
      rho = sum(bcBuffer( (posInBuffer-1)*nScalars+varPos(1) :  &
        &                 (posInBuffer-1)*nScalars+varPos(QQ) ))

      ! According to the bounce back rule,
      ! Subtract the equlibrium if its computed from alpha-(inverse of
      ! bitmask) or add the equlibrium if its computed in the direction of
      ! bitmask.
      ! In the Izquiredo paper, this is rho0, the reference density
      iDir   = me%inletUbbQVal(iLevel)%iDir(iLink)
      eqPlus = layout%weight(iDir) * 6.0_rk * rho  &
        &    * (   layout%fStencil%cxDirRK(1, iDir) * vel_b(offset+1) * inv_vel &
        &        + layout%fStencil%cxDirRK(2, iDir) * vel_b(offset+2) * inv_vel &
        &        + layout%fStencil%cxDirRK(3, iDir) * vel_b(offset+3) * inv_vel )

      state( me%links(iLevel)%val(iLink) ) = fOut + eqPlus

    end do ! iLink

  end subroutine velocity_bounceback
! ****************************************************************************** !


  ! ************************************************************************ !
  !> Inlet Velocity BFL rule  boundary condition
  !!
  !! This is taken from the paper:
  !! "One-point boundary condition for the lattice Boltzmann method", by
  !! M. Junk and Z. Yang published in Physical Review E 72 (2005) and leads to
  !! the chapter 'B. BFL rule'.
  !! This boundary condition shall be an improvement of the already existing
  !! velocity bounce back condition. It introduces a term called
  !! \( f_{i^{\star}}^{b} \).
  !! This one represents the particle distribution:
  !! \[ \hat{f}_{i}(n+1,\mathbf{j}) = \hat{f}_{i^{\star}}^{c}(n,\mathbf{j})
  !! + 6hf_{i}^{\star}\mathbf{c}_{i}
  !! \cdot \mathbf{\Phi}(t_{n},\mathbf{x}_{\mathbf{j}i}) \]
  !!
  !! This boundary condition has error of 1st order velocity and less 1st order
  !! pressure tested for wall is located exactly at q=1/2.
  !!
  !! Like for velocity_bounceback_qval we compute the values at the boundary
  !! points.
  !!
  !! ---------------------------------------------------------------------------
  !! Main differences to velocity_bounceback:
  !!  - implemented "link-wise", only incoming links at the boundary are updated
  !!    for each time step
  !!  - uses q-Values in pre-processing (see mus_bc_header_module)
  !!  - computes state values at the boundary nodes instead of barycentre
  !!    position
  !!  - can be an improvement of the bounce back rule for q-values ~= 0.5
  !! ---------------------------------------------------------------------------
  !!
  !! This subroutine's interface must match the abstract interface definition
  !! [[boundaryRoutine]] in bc/[[mus_bc_header_module]].f90 in order to be
  !! callable via [[boundary_type:fnct]] function pointer.
  subroutine velocity_bfl( me, state, bcBuffer, globBC, levelDesc, tree,      &
    &                      nSize, iLevel, sim_time, neigh, layout, fieldProp, &
    &                      varPos, nScalars, varSys, derVarPos, physics,      &
    &                      iField, mixture                                    )
    ! -------------------------------------------------------------------- !
    !> global boundary type
    class(boundary_type) :: me
    !> Current state vector of iLevel
    real(kind=rk), intent(inout) :: state(:)
    !> size of state array ( in terms of elements )
    integer, intent(in) :: nSize
    !> state values of boundary elements of all fields of iLevel
    real(kind=rk), intent(in) :: bcBuffer(:)
    !> iLevel descriptor
    type(tem_levelDesc_type), intent(in) :: levelDesc
    !> Treelm Mesh
    type(treelmesh_type), intent(in) :: tree
    !> fluid parameters and properties
    type(mus_field_prop_type), intent(in) :: fieldProp
    !> stencil layout information
    type(mus_scheme_layout_type), intent(in) :: layout
    !> the level On which this boundary was invoked
    integer, intent(in) :: iLevel
    !> connectivity array corresponding to state vector
    integer, intent(in) :: neigh(:)
    !> global time information
    type(tem_time_type), intent(in)  :: sim_time
    !> pointer to field variable in the state vector
    integer, intent(in) :: varPos(:)
    !> number of Scalars in the scheme var system
    integer, intent(in) :: nScalars
    !> scheme variable system
    type(tem_varSys_type), intent(in) :: varSys
    !> position of derived quantities in varsys
    type(mus_derVarPos_type), intent(in) :: derVarPos
    !> scheme global boundary type
    type(glob_boundary_type), intent(in) :: globBC
    !> scheme global boundary type
    type(mus_physics_type), intent(in) :: physics
    !> current field
    integer, intent(in) :: iField
    !> mixture info
    type(mus_mixture_type), intent(in) :: mixture
    ! -------------------------------------------------------------------- !
    integer       :: iLink, iDir, QQ, posInBuffer
    real(kind=rk) :: weight
    real(kind=rk) :: fIn, fOut, fNgh, eqPlus, rho
    real(kind=rk) :: cIn, cOut, cNgh, cVel
    real(kind=rk) :: vel_b(me%links(iLevel)%nVals*3), inv_vel
    integer :: bcVel_pos, offset
    ! ---------------------------------------------------------------------------
    QQ = layout%fstencil%QQ

    inv_vel = 1.0_rk / physics%fac( iLevel )%vel

    ! position of boundary velocity in varSys
    bcVel_pos = me%bc_states%velocity%varPos
    ! Multiply transient velocity with spatial velocity. Spatial term includes
    ! maxvalue
    call varSys%method%val(bcVel_pos)%get_valOfIndex( &
      & varSys  = varSys,                             &
      & time    = sim_time,                           &
      & iLevel  = iLevel,                             &
      & idx     = me%bc_states%velocity               &
      &           %pntIndex%indexLvl(iLevel)          &
      &           %val(1:me%links(iLevel)%nVals),     &
      & nVals   = me%links(iLevel)%nVals,             &
      & res     = vel_b                               )

    do iLink = 1, me%links(iLevel)%nVals
      offset = (iLink-1)*3

      ! QValues and coefficients
      cIn  = me%inletBfl(iLevel)% cIn(iLink)
      cOut = me%inletBfl(iLevel)%cOut(iLink)
      cNgh = me%inletBfl(iLevel)%cNgh(iLink)
      cVel = me%inletBfl(iLevel)%cVel(iLink)

      ! positions of element
      fIn  = bcBuffer(me%inletBfl(iLevel)% inPos(iLink))
      fOut = bcBuffer(me%inletBfl(iLevel)%outPos(iLink))
      fNgh = me%neigh(iLevel)%computeNeighBuf(me%inletbfl(iLevel)%nghPos(iLink))

      ! local density
      posInBuffer  = me%inletBfl(iLevel)%posInBuffer(iLink)
      rho = sum(bcBuffer( (posInBuffer-1)*nScalars+varPos(1) :  &
        &                 (posInBuffer-1)*nScalars+varPos(QQ) ))

      iDir   = me%inletBfl(iLevel)%iDir(iLink)
      weight = layout%weight(iDir)
      eqPlus = weight * 6._rk * rho  &
        &    * (  layout%fStencil%cxDirRK(1, iDir) * vel_b(offset+1) * inv_vel &
        &       + layout%fStencil%cxDirRK(2, iDir) * vel_b(offset+2) * inv_vel &
        &       + layout%fStencil%cxDirRK(3, iDir) * vel_b(offset+3) * inv_vel )

      ! We need to get post-collision pdf in direction
      ! alpha-, which is the inverse direction of bitmask
      ! For PULL this means, get the outgoing one, as this is the one which will be bounced back
      ! For PUSH this means, get the already bounced back pdf back, so take the incoming

      state( me%links(iLevel)%val(iLink) ) &
        &  = cIn*fIn + cOut*fOut + cNgh*fNgh + cVel*eqPlus

    end do ! iLink

  end subroutine velocity_bfl
! ****************************************************************************** !

! ****************************************************************************** !
  !> No comment yet!
  !!
  !! @TODO add comment
  !!
  !! This subroutine's interface must match the abstract interface definition
  !! [[boundaryRoutine]] in bc/[[mus_bc_header_module]].f90 in order to be
  !! callable via [[boundary_type:fnct]] function pointer.
  subroutine velocity_bfl_incomp( me, state, bcBuffer, globBC, levelDesc,      &
    &                             tree, nSize, iLevel, sim_time, neigh,        &
    &                             layout, fieldProp, varPos, nScalars, varSys, &
    &                             derVarPos, physics, iField, mixture          )
    ! -------------------------------------------------------------------- !
    !> global boundary type
    class(boundary_type) :: me
    !> Current state vector of iLevel
    real(kind=rk), intent(inout) :: state(:)
    !> size of state array ( in terms of elements )
    integer, intent(in) :: nSize
    !> state values of boundary elements of all fields of iLevel
    real(kind=rk), intent(in) :: bcBuffer(:)
    !> iLevel descriptor
    type(tem_levelDesc_type), intent(in) :: levelDesc
    !> Treelm Mesh
    type(treelmesh_type), intent(in) :: tree
    !> fluid parameters and properties
    type(mus_field_prop_type), intent(in) :: fieldProp
    !> stencil layout information
    type(mus_scheme_layout_type), intent(in) :: layout
    !> the level On which this boundary was invoked
    integer, intent(in) :: iLevel
    !> connectivity array corresponding to state vector
    integer, intent(in) :: neigh(:)
    !> global time information
    type(tem_time_type), intent(in)  :: sim_time
    !> pointer to field variable in the state vector
    integer, intent(in) :: varPos(:)
    !> number of Scalars in the scheme var system
    integer, intent(in) :: nScalars
    !> scheme variable system
    type(tem_varSys_type), intent(in) :: varSys
    !> position of derived quantities in varsys
    type(mus_derVarPos_type), intent(in) :: derVarPos
    !> scheme global boundary type
    type(glob_boundary_type), intent(in) :: globBC
    !> scheme global boundary type
    type(mus_physics_type), intent(in) :: physics
    !> current field
    integer, intent(in) :: iField
    !> mixture info
    type(mus_mixture_type), intent(in) :: mixture
    ! -------------------------------------------------------------------- !
    integer       :: iLink, iDir
    real(kind=rk) :: weight
    real(kind=rk) :: fIn, fOut, fNgh, eqPlus
    real(kind=rk) :: cIn, cOut, cNgh, cVel
    real(kind=rk) :: vel_b(me%links(iLevel)%nVals*3), inv_vel
    integer :: bcVel_pos, offset
    ! ---------------------------------------------------------------------------

    inv_vel = 1.0_rk / physics%fac( iLevel )%vel

    ! position of boundary velocity in varSys
    bcVel_pos = me%bc_states%velocity%varPos
    ! Multiply transient velocity with spatial velocity. Spatial term includes
    ! maxvalue
    call varSys%method%val(bcVel_pos)%get_valOfIndex( &
      & varSys  = varSys,                             &
      & time    = sim_time,                           &
      & iLevel  = iLevel,                             &
      & idx     = me%bc_states%velocity               &
      &           %pntIndex%indexLvl(iLevel)          &
      &           %val(1:me%links(iLevel)%nVals),     &
      & nVals   = me%links(iLevel)%nVals,             &
      & res     = vel_b                               )

    do iLink = 1, me%links(iLevel)%nVals
      offset = (iLink-1)*3

      ! QValues and coefficients
      cIn  = me%inletBfl(iLevel)% cIn(iLink)
      cOut = me%inletBfl(iLevel)%cOut(iLink)
      cNgh = me%inletBfl(iLevel)%cNgh(iLink)
      cVel = me%inletBfl(iLevel)%cVel(iLink)

      ! positions of element
      fIn  = bcBuffer(me%inletBfl(iLevel)% inPos(iLink))
      fOut = bcBuffer(me%inletBfl(iLevel)%outPos(iLink))
      fNgh = me%neigh(iLevel)%computeNeighBuf(me%inletbfl(iLevel)%nghPos(iLink))

      iDir   = me%inletBfl(iLevel)%iDir(iLink)
      weight = layout%weight(iDir)
      eqPlus = weight * 6._rk * rho0  &
        &    * (  layout%fStencil%cxDirRK(1, iDir) * vel_b(offset+1) * inv_vel &
        &       + layout%fStencil%cxDirRK(2, iDir) * vel_b(offset+2) * inv_vel &
        &       + layout%fStencil%cxDirRK(3, iDir) * vel_b(offset+3) * inv_vel )

      ! We need to get post-collision pdf in direction
      ! alpha-, which is the inverse direction of bitmask
      ! For PULL this means, get the outgoing one, as this is the one which will be bounced back
      ! For PUSH this means, get the already bounced back pdf back, so take the incoming

      state( me%links(iLevel)%val(iLink) ) &
        &  = cIn*fIn + cOut*fOut + cNgh*fNgh + cVel*eqPlus

    end do ! iLink

  end subroutine velocity_bfl_incomp
! ****************************************************************************** !


  ! **************************************************************************** !
  !> Inlet Velocity Bounce Back boundary condition with mass flow rate as input
  !!
  !! This is taken from the paper:
  !! [1] S. Izquierdo and N. Fueyo, “Characteristic  non-reflecting
  !! boundary  conditions  for  open  boundaries  in  lattice
  !!  Boltzmann  methods,” Physical Review E, vol. 78, no. 46707, 2008
  !! The incoming densities \( \bar{f_{alpha}} \) are reconstructed from the
  !! already reflected densities, enforcing the given velocity
  !! me\%ubb\%velocity.
  !! [2] M. Junk and Z. Yang. One-point boundary condition for the
  !! lattice Boltzmann method,
  !! Physical Review E, vol 72, issue 8, year 2005.
  !!
  !! This boundary condition has error of 1st order pressure and 2nd order
  !! velocity only if wall is located exactly at q=1/2 else the accuracy of
  !! both pressure and velocity reduce by order 1.
  !!
  !! Usage
  !! -----
  !!```lua
  !!boundary_condition = {
  !! { label = 'inlet',
  !!    kind = 'mfr_bounceback',
  !!    massflowrate = 'mfr }
  !!}
  !! variable = {
  !!   name = 'mfr',
  !!   ncomponents = 1,
  !!   vartype = 'st_fun',
  !!   st_fun = 0.1
  !! }
  !!```
  !!
  !! This subroutine's interface must match the abstract interface definition
  !! [[boundaryRoutine]] in bc/[[mus_bc_header_module]].f90 in order to be
  !! callable via [[boundary_type:fnct]] function pointer.
  subroutine mfr_bounceback( me, state, bcBuffer, globBC, levelDesc, tree,   &
    &                        nSize, iLevel, sim_time, neigh, layout,         &
    &                        fieldProp, varPos, nScalars, varSys, derVarPos, &
    &                        physics, iField, mixture                        )
    ! -------------------------------------------------------------------- !
    !> global boundary type
    class(boundary_type) :: me
    !> Current state vector of iLevel
    real(kind=rk), intent(inout) :: state(:)
    !> size of state array ( in terms of elements )
    integer, intent(in) :: nSize
    !> state values of boundary elements of all fields of iLevel
    real(kind=rk), intent(in) :: bcBuffer(:)
    !> iLevel descriptor
    type(tem_levelDesc_type), intent(in) :: levelDesc
    !> Treelm Mesh
    type(treelmesh_type), intent(in) :: tree
    !> fluid parameters and properties
    type(mus_field_prop_type), intent(in) :: fieldProp
    !> stencil layout information
    type(mus_scheme_layout_type), intent(in) :: layout
    !> the level On which this boundary was invoked
    integer, intent(in) :: iLevel
    !> connectivity array corresponding to state vector
    integer, intent(in) :: neigh(:)
    !> global time information
    type(tem_time_type), intent(in)  :: sim_time
    !> pointer to field variable in the state vector
    integer, intent(in) :: varPos(:)
    !> number of Scalars in the scheme var system
    integer, intent(in) :: nScalars
    !> scheme variable system
    type(tem_varSys_type), intent(in) :: varSys
    !> position of derived quantities in varsys
    type(mus_derVarPos_type), intent(in) :: derVarPos
    !> scheme global boundary type
    type(glob_boundary_type), intent(in) :: globBC
    !> scheme global boundary type
    type(mus_physics_type), intent(in) :: physics
    !> current field
    integer, intent(in) :: iField
    !> mixture info
    type(mus_mixture_type), intent(in) :: mixture
    ! -------------------------------------------------------------------- !
    real(kind=rk) :: eqPlus
    real(kind=rk) :: fTmp( layout%fStencil%QQ )
    real(kind=rk) :: massFlowRate(globBC%nElems(iLevel))
    real(kind=rk) :: velocity(3), massFlowRateToVel
    integer :: iELem, iDir, iLvl, QQ, posInBuffer
    reaL(kind=rk) :: area
    integer :: bcMfr_pos
    ! ---------------------------------------------------------------------------

    QQ = layout%fStencil%QQ

    ! position of mass flow rate spacetime function variable in varSys
    bcMfr_pos = me%bc_states%massFlowRate%varPos
    ! get mass flow rate
    call varSys%method%val(bcMfr_pos)%get_valOfIndex( &
      & varSys  = varSys,                             &
      & time    = sim_time,                           &
      & iLevel  = iLevel,                             &
      & idx     = me%bc_states%massFlowRate           &
      &           %pntIndex%indexLvl(iLevel)          &
      &           %val(1:globBC%nElems(iLevel)),      &
      & nVals   = globBC%nElems(iLevel),              &
      & res     = massFlowRate                        )

    ! convertion factor to convert mass flow rate to velocity
    ! vel_lattice = massFlowRate/(density*nElems*dx^2)
    ! area of cross section = sum(nElems(iLevel)*dx(iLevel)**2)
    ! @todo: element area can be computed once and saved
    area = 0.0_rk
    do iLvl = tree%global%minLevel, tree%global%maxLevel
      area = area + globBC%nElems_totalLevel(iLvl) * physics%dxLvl(iLvl)**2
    end do

    massFlowRateToVel = 1.0_rk / ( physics%rho0 * area )

    ! If physical quantities are given, transform to lattice units by division
    ! with the conversion factor
!write(*,*) 'velocity_P ', massFlowRateToVel, 1.0/massFlowRateToVel
    massFlowRateToVel = massFlowRateToVel / physics%fac( iLevel )%vel
!write(*,*) 'velocity_L ', massFlowRateToVel

    do iElem = 1, globBC%nElems( iLevel )
      ! compute velocity from massFlowRate
      velocity = massFlowRate(iElem) * massFlowRateToVel                     &
        & * layout%fStencil%cxDir(:, globBC%elemLvl(iLevel)%normalInd%val(iElem))

      !store state in temp array
      !KM:produce no checker board effect if we use pre-collision values of
      !current time step
      posInBuffer = globBC%elemLvl( iLevel )%posInBcElemBuf%val( iElem )
      fTmp(1:QQ) = bcBuffer( (posInBuffer-1)*nScalars+varPos(1) : &
        &                    (posInBuffer-1)*nScalars+varPos(QQ)  )

      do iDir = 1, layout%fStencil%QQN
        ! Write the values
        if( globBC%elemLvl(iLevel)%bitmask%val(iDir, iElem )) then
          ! Accourding to the bounce back rule,
          ! Subtract the equlibrium if its computed from alpha-(inverse of
          ! bitmask) or add the equlibrium if its computed in the direction of
          ! bitmask.
          ! In the Izquiredo paper, this is rho0, the reference density
          eqPlus =    layout%weight( iDir ) * 6._rk * rho0             &
            &       * ( layout%fStencil%cxDir(1,iDir) * velocity(1)   &
            &       +   layout%fStencil%cxDir(2,iDir) * velocity(2)   &
            &       +   layout%fStencil%cxDir(3,iDir) * velocity(3)   )
 state( neigh((idir-1)*nsize+ globbc%elemlvl(ilevel)%elem%val(ielem))+( ifield-1)* qq+ nscalars*0 ) = &
          ! We need to get post-collision pdf in direction
          ! alpha-, which is the inverse direction of bitmask
          ! For PULL this means, get the outgoing one, as this is the one which will be bounced back
          ! For PUSH this means, get the already bounced back pdf back, so take the incoming
            & fTmp( layout%fStencil%cxDirInv( iDir ) ) + eqPlus
        end if ! bitMask
      end do ! iDir
    end do

  end subroutine mfr_bounceback
! ****************************************************************************** !


  ! ****************************************************************************** !
  !> author: Kannan Masilamani
  !! Outlet boundary conditions with zero pressure gradient.
  !!
  !! These boundary conditions use the neighbor cells density and velocity in the
  !! equilibrium function. It is not necessary to specify density at boundary in
  !! the lua configuration file \( f = f^{eq}(\rho_{b-1},u_{b-1}) \)
  !!
  !! MH: Kannan, please update them
  !!
  !! This subroutine's interface must match the abstract interface definition
  !! [[boundaryRoutine]] in bc/[[mus_bc_header_module]].f90 in order to be
  !! callable via [[boundary_type:fnct]] function pointer.
  subroutine outlet_zero_prsgrd( me, state, bcBuffer, globBC, levelDesc, tree, &
    &                            nSize, iLevel, sim_time, neigh, layout,       &
    &                            fieldProp, varPos, nScalars, varSys,          &
    &                            derVarPos, physics, iField, mixture           )
    ! -------------------------------------------------------------------- !
    !> global boundary type
    class(boundary_type) :: me
    !> Current state vector of iLevel
    real(kind=rk), intent(inout) :: state(:)
    !> size of state array ( in terms of elements )
    integer, intent(in) :: nSize
    !> state values of boundary elements of all fields of iLevel
    real(kind=rk), intent(in) :: bcBuffer(:)
    !> iLevel descriptor
    type(tem_levelDesc_type), intent(in) :: levelDesc
    !> Treelm Mesh
    type(treelmesh_type), intent(in) :: tree
    !> fluid parameters and properties
    type(mus_field_prop_type), intent(in) :: fieldProp
    !> stencil layout information
    type(mus_scheme_layout_type), intent(in) :: layout
    !> the level On which this boundary was invoked
    integer, intent(in) :: iLevel
    !> connectivity array corresponding to state vector
    integer, intent(in) :: neigh(:)
    !> global time information
    type(tem_time_type), intent(in)  :: sim_time
    !> pointer to field variable in the state vector
    integer, intent(in) :: varPos(:)
    !> number of Scalars in the scheme var system
    integer, intent(in) :: nScalars
    !> scheme variable system
    type(tem_varSys_type), intent(in) :: varSys
    !> position of derived quantities in varsys
    type(mus_derVarPos_type), intent(in) :: derVarPos
    !> scheme global boundary type
    type(glob_boundary_type), intent(in) :: globBC
    !> scheme global boundary type
    type(mus_physics_type), intent(in) :: physics
    !> current field
    integer, intent(in) :: iField
    !> mixture info
    type(mus_mixture_type), intent(in) :: mixture
    ! -------------------------------------------------------------------- !
    ! defining local variables
    real(kind=rk) :: rho, ux(globBC%nElems(iLevel),3), usq1_3
    integer :: iElem, iDir, QQ
    real(kind=rk) :: fTmp( nScalars * globBC%nElems(iLevel) )
    integer :: posInBuffer
    ! ---------------------------------------------------------------------------

    QQ   = layout%fStencil%QQ

    do iElem = 1, globBC%nElems(iLevel)
      posInBuffer = globBC%elemLvl( iLevel )%posInBcElemBuf%val( iElem )
      fTmp( (iElem-1)*nScalars+1: (iElem-1)*nScalars+QQ ) &
        &       = bcBuffer( (posInBuffer-1)*nScalars+varPos(1) : &
        &                   (posInBuffer-1)*nScalars+varPos(QQ)  )

      ! Get my velocity
      call derVarPos%velFromState( state  = fTmp( (iElem-1)*nScalars+1 &
        &                                        :iElem*nScalars ),    &
        &                          iField = iField,                    &
        &                          nElems = 1,                         &
        &                          varSys = varSys,                    &
        &                          layout = layout,                    &
        &                          res    = ux(iElem, :)               )
    end do

    do iElem = 1, globBC%nElems(iLevel)
      ! Calculate local density and velocity moments
      rho = sum( fTmp( (iElem-1)*QQ+1:iElem*QQ ) )

      usq1_3 =  sum( ux(iElem,:) ) * div1_3

      do iDir = 1, layout%fStencil%QQN
        ! Write the values
        if( globBC%elemLvl(iLevel)%bitmask%val(iDir, iElem )) then
          state(                                                      &
&neigh((idir-1)*nsize+globbc%elemlvl(ilevel)%elem%val(ielem))+(ifield-1)*qq+nscalars*0 &
&         )  =                                                   &
            &  layout%weight( iDir )*( 2._rk*rho + 9._rk*rho0* ( (    &
            &            layout%fStencil%cxDirRK(1,iDir)*ux(iElem,1)  &
            &         +  layout%fStencil%cxDirRK(2,iDir)*ux(iElem,2)  &
            &         +  layout%fStencil%cxDirRK(3,iDir)*ux(iElem,3)  &
            &         )**2                                                 &
            &         -  usq1_3 ) )                                        &
          ! We need to get post-collision pdf in direction
          ! alpha-, which is the inverse direction of bitmask
          ! For PULL this means, get the outgoing one, as this is the one which will be bounced back
          ! For PUSH this means, get the already bounced back pdf back, so take the incoming
            &         - fTmp( (iElem-1)*QQ + layout%fStencil%cxDirInv(iDir) )
        end if
      end do
      ! end if
    end do

  end subroutine outlet_zero_prsgrd
! ****************************************************************************** !


  ! ************************************************************************ !
  !> author: Manuel Hasert
  !! Outlet Pressure Bounce Back boundary condition
  !!
  !! An anti-bounce back is performed for the outgoing links, so that the
  !! incoming densities are set with the defined reference density while
  !! maintaining the correct non-equilibrium part.
  !! The velocity is extrapolated from two neighboring fluid elements (boundary
  !! element itself + one fluid neighbor). The position of the wall is hence
  !! between the boundary element and the (non-existing) solid element
  !! Based on the Bouzidi Boundary conditions, the current implementation is
  !! taken from the paper:
  !! S. Izquierdo and N. Fueyo, "Characteristic  non-reflecing  boundary
  !!  conditions  for  open  boundaries  in  lattice  Boltzmann
  !!  methods," Physical Review E, vol. 78, no. 46707, 2008.
  !!
  !! MH: I deactivated the correction term as it involves information from the
  !! old time step. In parallel, this can not be provided.
  !!
  !! DISCLAIMER: This BC requires pre-collision PDF from 1st and 2nd neighbor
  !! which are not available if those neighbors are boundary halo on different
  !! process. So, if you encounter a problem with this BC, please check the
  !! domain distribution to verify whether 2 neighbors of boundary elements
  !! are in one process. The domain distribution can be changed by shifting the
  !! boundary by one element or by changing the origin of bounding cube or
  !! by using different number of MPI processes.
  !!
  !! \todo 20200410, KM: Check if this BC works if neighbufferPre is replaced
  !! by neighbufferPost
  !!
  !! This subroutine's interface must match the abstract interface definition
  !! [[boundaryRoutine]] in bc/[[mus_bc_header_module]].f90 in order to be
  !! callable via [[boundary_type:fnct] function pointer.
  subroutine pressure_antiBounceBack( me, state, bcBuffer, globBC, levelDesc, &
    &                                 tree, nSize, iLevel, sim_time, neigh,   &
    &                                 layout, fieldProp, varPos, nScalars,    &
    &                                 varSys, derVarPos, physics, iField,     &
    &                                 mixture                                 )
    ! -------------------------------------------------------------------- !
    !> global boundary type
    class(boundary_type) :: me
    !> Current state vector of iLevel
    real(kind=rk), intent(inout) :: state(:)
    !> size of state array ( in terms of elements )
    integer, intent(in) :: nSize
    !> state values of boundary elements of all fields of iLevel
    real(kind=rk), intent(in) :: bcBuffer(:)
    !> iLevel descriptor
    type(tem_levelDesc_type), intent(in) :: levelDesc
    !> Treelm Mesh
    type(treelmesh_type), intent(in) :: tree
    !> fluid parameters and properties
    type(mus_field_prop_type), intent(in) :: fieldProp
    !> stencil layout information
    type(mus_scheme_layout_type), intent(in) :: layout
    !> the level On which this boundary was invoked
    integer, intent(in) :: iLevel
    !> connectivity array corresponding to state vector
    integer, intent(in) :: neigh(:)
    !> global time information
    type(tem_time_type), intent(in)  :: sim_time
    !> pointer to field variable in the state vector
    integer, intent(in) :: varPos(:)
    !> number of Scalars in the scheme var system
    integer, intent(in) :: nScalars
    !> scheme variable system
    type(tem_varSys_type), intent(in) :: varSys
    !> position of derived quantities in varsys
    type(mus_derVarPos_type), intent(in) :: derVarPos
    !> scheme global boundary type
    type(glob_boundary_type), intent(in) :: globBC
    !> scheme global boundary type
    type(mus_physics_type), intent(in) :: physics
    !> current field
    integer, intent(in) :: iField
    !> mixture info
    type(mus_mixture_type), intent(in) :: mixture
    ! -------------------------------------------------------------------- !
    real(kind=rk) :: rhoDef(globBC%nElems(iLevel)) ! Density on boundary element
    real(kind=rk) :: rhoF    ! density on fluid element
    real(kind=rk) :: uxF(globBC%nElems(iLevel)*3)  ! Velocity on fluid element
    real(kind=rk) :: uxN(globBC%nElems(iLevel)*3)  ! Velocity on neighbor
    real(kind=rk) :: uxB(3)  ! Velocity on boundary surface
    real(kind=rk) :: usqB, usqF ! Squared velocity magnitude
    real(kind=rk) :: fPlusFluid ! f Plus for anti-bounceback correction
    ! Equilibrium Plus for anti-bounceback correction
    real(kind=rk) :: fEqPlus
    ! Equlibrium Plus from the old time step
    real(kind=rk) :: fEqPlusFluid
    real(kind=rk) :: inv_rho_phy
    real(kind=rk) :: fTmp( nScalars * globBC%nElems(iLevel) )
    real(kind=rk) :: omega       ! Relaxation parameter
    integer :: iDir, iElem, QQ, invDir, elemPos
    integer :: posInBuffer, bcPress_pos
    ! ---------------------------------------------------------------------------

    ! write(dbgUnit(6), *) ''
    ! write(dbgUnit(6), *) ' do pressure_antiBounceBack for bc label: '//trim(globBC%label)

    QQ = layout%fStencil%QQ
    inv_rho_phy = 1.0_rk / physics%fac(iLevel)%press * cs2inv

    ! position of boundary pressure in varSys
    bcPress_pos = me%bc_states%pressure%varPos
    ! get pressure variable from spacetime function
    call varSys%method%val(bcPress_pos)%get_valOfIndex( &
      & varSys  = varSys,                               &
      & time    = sim_time,                             &
      & iLevel  = iLevel,                               &
      & idx     = me%bc_states%pressure                 &
      &           %pntIndex%indexLvl(iLevel)            &
      &           %val(1:globBC%nElems(iLevel)),        &
      & nVals   = globBC%nElems(iLevel),                &
      & res     = rhoDef                                )

    ! convert physical pressure into LB density
    rhoDef = rhoDef * inv_rho_phy + rho0

    do iElem = 1, globBC%nElems( iLevel )
      posInBuffer = globBC%elemLvl( iLevel )%posInBcElemBuf%val( iElem )
      fTmp( (iElem-1)*nScalars+1: (iElem-1)*nScalars+QQ ) &
        &       = bcBuffer( (posInBuffer-1)*nScalars+varPos(1) : &
        &                   (posInBuffer-1)*nScalars+varPos(QQ)  )
    end do

    ! Get my velocity
    call derVarPos%velFromState( state  = fTmp,                  &
      &                          iField = iField,                &
      &                          nElems = globBC%nElems(iLevel), &
      &                          varSys = varSys,                &
      &                          layout = layout,                &
      &                          res    = uxF                    )

    ! Calculate my neighbor velocity
    call derVarPos%velFromState( state  = me%neigh(iLevel)%neighBufferPre(1,:), &
      &                          iField = iField,                &
      &                          nElems = globBC%nElems(iLevel), &
      &                          varSys = varSys,                &
      &                          layout = layout,                &
      &                          res    = uxN                    )

    do iElem = 1, globBC%nElems( iLevel )

      rhoF = sum(fTmp( (iElem-1)*nScalars+1: (iElem-1)*nScalars+QQ ))

      uxB =   1.5_rk * uxF( (iElem-1)*3+1 : iElem*3 ) &
        &   - 0.5_rk * uxN( (iElem-1)*3+1 : iElem*3 )

      usqB = uxB(1)*uxB(1) + uxB(2)*uxB(2) + uxB(3)*uxB(3)
      usqF =   uxF((iElem-1)*3+1) * uxF((iElem-1)*3+1) &
        &    + uxF((iElem-1)*3+2) * uxF((iElem-1)*3+2) &
        &    + uxF((iElem-1)*3+3) * uxF((iElem-1)*3+3)

! -----------  debug output    --------------------------------------------
! write(dbgUnit(6), *) iElem, ' treeID: ',&
!   &     levelDesc%total(globBC%elemLvl(iLevel)%elem%val(iElem)), &
!   &     'elemBuffer: ', globBC%elemLvl( iLevel )%posInBcElemBuf%val(iElem)
! write(dbgUnit(6), *) 'rhoDef: ', rhoDef( iElem )
! write(dbgUnit(6), *) '  rhoF: ', rhoF
! write(dbgUnit(6), *) '   uxF: ', uxF( (iElem-1)*3+1:(iElem-1)*3+3)
! write(dbgUnit(6), *) '   uxN: ', uxN( (iElem-1)*3+1:(iElem-1)*3+3)
! write(dbgUnit(6), *) '   uxB: ', uxB(1:3)
! write(dbgUnit(6), *) '  usqB: ', usqB
! write(dbgUnit(6), *) '  usqF: ', usqF
! write(dbgUnit(6), *) ' origin pdf:'
! do iDir = 1, QQ
!   write( dbgUnit(6), *) iDir, fTmp( (iElem-1)*nScalars + iDir )
! end do
! -------------------------------------------------------------------------
      elemPos = globBC%elemLvl(iLevel)%elem%val( iElem )
      ! omega is STfun so get omega for current element
      omega = fieldProp%fluid%viscKine%omLvl(iLevel)%val(elemPos)

      do iDir = 1, layout%fStencil%QQN
        if( globBC%elemLvl(iLevel)%bitmask%val(iDir, iElem )) then
          ! Calc the correction term and the Dirichlet (Eq) term

          invDir = layout%fStencil%cxDirInv(iDir)

          fEqPlusFluid = layout%weight( iDir )*rhoF                  &
            &          + 4.5_rk*layout%weight( iDir )*rho0*((        &
            &            layout%fStencil%cxDirRK( 1, invDir )*uxF((iElem-1)*3+1) &
            &          + layout%fStencil%cxDirRK( 2, invDir )*uxF((iElem-1)*3+2) &
            &          + layout%fStencil%cxDirRK( 3, invDir )*uxF((iElem-1)*3+3) &
            &          )**2                                       &
            &          - div1_3*usqF)

          fEqPlus = layout%weight( iDir )*rhoDef(iElem)        &
            &     + 4.5_rk*layout%weight( iDir )*rho0*((       &
            &       layout%fStencil%cxDirRK( 1, InvDir)*uxB(1) &
            &     + layout%fStencil%cxDirRK( 2, InvDir)*uxB(2) &
            &     + layout%fStencil%cxDirRK( 3, InvDir)*uxB(3) &
            &     )**2                                         &
            &     - div1_3*usqB)

          fPlusFluid = 0.5_rk*(   fTmp( (iElem-1)*nScalars + iDir )  &
            &                   + fTmp( (iElem-1)*nScalars + invDir ) )

          ! Now assign the values
          ! Actually, we should take fPlusLAst and fEqPlusLast.
          ! But this introduces errors in parallel
          state( neigh (( idir-1)*nsize+ elempos)+( ifield-1)* qq+ nscalars*0) = &
          ! antibounceback term
! & - state(  )    &
          ! We need to get post-collision pdf in direction
          ! alpha-, which is the inverse direction of bitmask
          ! For PULL this means, get the outgoing one, as this is the one which will be bounced back
          ! For PUSH this means, get the already bounced back pdf back, so take the incoming
    &   - fTmp( (iElem-1)*nScalars + invDir ) &
    &   + 2._rk*fEqPlus     &                          ! Dirichlet pressure term
    &   + (2._rk - omega )*( fPlusFluid - fEqPlusFluid ) ! Correction term

! -----------   DEBUG output     -------------------------------------------
! write(dbgUnit(6), *) '      fEqPlus: ', fEqPlus
! write(dbgUnit(6), *) ' fEqPlusFluid: ', fEqPlusFluid
! write(dbgUnit(6), *) '   fPlusFluid: ', fPlusFluid
! write(dbgUnit(6), *) '  updated pdf:'
! write( dbgUnit(6), *) 'iDir', iDir, 'invDir', invDir, state(                      &
! & neigh (( idir-1)* nsize+ globbc%elemlvl(ilevel)%elem%val( ielem ))+( ifield-1)* qq+ nscalars*0)
! --------------------------------------------------------------------------

        end if ! bitMask
      end do ! iDir
    end do ! iElem

  end subroutine pressure_antiBounceBack
! ****************************************************************************** !

! ****************************************************************************** !
  !> author: Kannan Masilamani
  !! Characteristic-based non-reflective inlet boundary conditions for
  !! incompressible flows
  !! @todo: add explaination of these steps
  !!
  !! These boundary conditions are taken from the paper:
  !!   S. Izquierdo and N. Fueyo,
  !!   "Characteristic nonreflecting boundary
  !!   conditions for open boundaries in lattice Boltzmann methods,"
  !!   Physical Review E, vol. 78, no. 46707, 2008.
  !!
  !! > Note: unstable behavior for high omega relaxation parameters.
  !! >       It is recommended to use a sponge layer with a
  !! >       \ref mus_aux_module::mus_setspatialomega "spatial omega"
  !! >       where the omega is decreased towards the boundaries to increase
  !! >       stability
  !!
  !! This subroutine's interface must match the abstract interface definition
  !! [[boundaryRoutine]] in bc/[[mus_bc_header_module]].f90 in order to be
  !! callable via [[boundary_type:fnct]] function pointer.
  subroutine inlet_nrbc_incomp( me, state, bcBuffer, globBC, levelDesc, &
    &                           tree, nSize, iLevel, sim_time, neigh,   &
    &                           layout, fieldProp, varPos, nScalars,    &
    &                           varSys, derVarPos, physics, iField,     &
    &                           mixture                                 )
    ! -------------------------------------------------------------------- !
    !> global boundary type
    class(boundary_type) :: me
    !> Current state vector of iLevel
    real(kind=rk), intent(inout) :: state(:)
    !> size of state array ( in terms of elements )
    integer, intent(in) :: nSize
    !> state values of boundary elements of all fields of iLevel
    real(kind=rk), intent(in) :: bcBuffer(:)
    !> iLevel descriptor
    type(tem_levelDesc_type), intent(in) :: levelDesc
    !> Treelm Mesh
    type(treelmesh_type), intent(in) :: tree
    !> fluid parameters and properties
    type(mus_field_prop_type), intent(in) :: fieldProp
    !> stencil layout information
    type(mus_scheme_layout_type), intent(in) :: layout
    !> the level On which this boundary was invoked
    integer, intent(in) :: iLevel
    !> connectivity array corresponding to state vector
    integer, intent(in) :: neigh(:)
    !> global time information
    type(tem_time_type), intent(in)  :: sim_time
    !> pointer to field variable in the state vector
    integer, intent(in) :: varPos(:)
    !> number of Scalars in the scheme var system
    integer, intent(in) :: nScalars
    !> scheme variable system
    type(tem_varSys_type), intent(in) :: varSys
    !> position of derived quantities in varsys
    type(mus_derVarPos_type), intent(in) :: derVarPos
    !> scheme global boundary type
    type(glob_boundary_type), intent(in) :: globBC
    !> scheme global boundary type
    type(mus_physics_type), intent(in) :: physics
    !> current field
    integer, intent(in) :: iField
    !> mixture info
    type(mus_mixture_type), intent(in) :: mixture
    ! -------------------------------------------------------------------- !
    real(kind=rk) :: rhoDef(globBC%nElems(iLevel)) ! Density on boundary element
    real(kind=rk) :: rhoF    ! density on fluid element
    real(kind=rk) :: uxF(globBC%nElems(iLevel)*3)  ! Velocity on fluid element
    real(kind=rk) :: uxN(globBC%nElems(iLevel)*3)  ! Velocity on neighbor
    ! Equilibrium Plus for bounceback
    real(kind=rk) :: fEqPlus
    real(kind=rk) :: inv_rho_phy
    real(kind=rk) :: fTmp( nScalars )
    real(kind=rk) :: fTmpPre( nScalars )
    real(kind=rk) :: fTmpPre_2( nScalars )
    real(kind=rk) :: omega       ! Relaxation parameter
    integer :: iDir, iElem, QQ, invDir, elemPos
    integer :: posInBuffer, bcPress_pos
    real(kind=rk) :: swap
    real(kind=rk) :: L1, L2, L3, L4, L5
    real(kind=rk) :: dudx, dvdx, dwdx, drhodx
    real(kind=rk) :: rhoN1Prev, uN1Prev(3), rhoN2Prev, uN2Prev(3)
    real(kind=rk) :: rhoLodi, uLodi(3)
    integer :: normal(3)
    real(kind=rk) :: vel_b(globBC%nElems(iLevel)*3), inv_vel, vel_bE(3)
    integer :: bcVel_pos
    ! ---------------------------------------------------------------------------

    ! write(dbgUnit(6), *) ''
    ! write(dbgUnit(6), *) ' do characteristic velocity BounceBack for ' &
    !   &                //'bc label: '//trim(globBC%label)

    QQ = layout%fStencil%QQ
    inv_vel = 1.0_rk / physics%fac( iLevel )%vel

    ! position of boundary velocity in varSys
    bcVel_pos = me%bc_states%velocity%varPos
    ! Get velocity
    call varSys%method%val(bcVel_pos)%get_valOfIndex( &
      & varSys  = varSys,                             &
      & time    = sim_time,                           &
      & iLevel  = iLevel,                             &
      & idx     = me%bc_states%velocity               &
      &           %pntIndex%indexLvl(iLevel)          &
      &           %val(1:globBC%nElems(iLevel)),      &
      & nVals   = globBC%nElems(iLevel),              &
      & res     = vel_b                               )

    ! convert physical velocity into LB velocity
    vel_b = vel_b * inv_vel

    do iElem = 1, globBC%nElems( iLevel )

      ! Dirichlet velocity of current element
      vel_bE(:) = vel_b((iElem-1)*3+1:iElem*3)

      ! my density and velocity
      rhoLodi   = me%elemLvl( iLevel )%lodi( 1,   1, iElem )
      uLodi     = me%elemLvl( iLevel )%lodi( 2:4, 1, iElem )
      ! density and velocity of my first neighbor from last iteration
      rhoN1Prev = me%elemLvl( iLevel )%lodi( 1,   2, iElem )
      uN1Prev   = me%elemLvl( iLevel )%lodi( 2:4, 2, iElem )
      ! density and velocity of my secnd neighbor from last iteration
      rhoN2Prev = me%elemLvl( iLevel )%lodi( 1,   3, iElem )
      uN2Prev   = me%elemLvl( iLevel )%lodi( 2:4, 3, iElem )
      ! Use a swapping variable to treat +x and -x directions
      ! Incoming wave amplitudes L2, L3 and L4 in Eq (9) point inwards
      ! and outgoing wave amplitude L1 points outwards.
      ! The derivatives at +x and +y direction are computed by second-order
      ! backward stencil i.e. Eq (15) and at -x and -y, the derivatives
      ! are computed by second-order forward stencil i.e. Eq (21).
      !
      ! Swap is used to switch between forward and backward stencil.
      ! It is positive for backward stencil and negative for forward stencil.
      ! NormalInd refers to incoming direction so use cxDirInv to get normal
      ! pointing away from boundary.
      normal = layout%fStencil%cxDir(:,                                   &
        &                      layout%fStencil%cxDirInv(globBC%normalInd) )
!write(dbgUnit(6),*) 'normal ', normal
      ! x-direction
      if (abs(normal(1)) == 1) then
        swap = -1.0_rk*(real(normal(1),rk))

        ! LODI Treatment
        ! calculate spatial derivatives according to (15) which will be used to
        ! compute the characteristics
        ! This is a quadratic extrapolation
        drhodx =  swap*div1_3*(-8.d0*rhoLodi  + 9.d0*rhoN1Prev  - rhoN2Prev)
        dudx   =  swap*div1_3*(-8.d0*uLodi(1) + 9.d0*uN1Prev(1) - uN2Prev(1))
        dvdx   =  swap*div1_3*(-8.d0*uLodi(2) + 9.d0*uN1Prev(2) - uN2Prev(2))
        dwdx   =  swap*div1_3*(-8.d0*uLodi(3) + 9.d0*uN1Prev(3) - uN2Prev(3))

        ! wave amplitudes
        ! Outgoioming wave L1 Eq (9)
        L1 = (uLodi(1) - swap * cs) &
          & * (cs2 * drhodx - swap * cs * rho0 * dudx)
        ! incoming wave amplitudes L2, L3, L4 and L5 is approached with a
        ! linear relaxation model Eq (22)
        L2 = - (vel_bE(2)-uLodi(2))
        L3 = - (vel_bE(3)-uLodi(3))
        ! L5 is valid only if cs = sqrt(kappa*p/rho)
        L4 = 0.0_rk
        L5 = L1 - 2.0_rk*rho0*swap*cs*(vel_bE(1)-uLodi(1))

        ! LODI Equations (18)
        rhoLodi  = rhoLodi  - ( L4 + 0.5d0*(L5+L1) ) * cs2inv
        ! Eq. (18b) for x-velocity
        uLodi(1) = uLodi(1) - (L5-L1)/(rho0*swap*cs)*0.5d0
        uLodi(2) = uLodi(2) - L2  ! equation (18c) for y-direction
        uLodi(3) = uLodi(3) - L3  ! equation (18c) for z-direction
      else if( abs(normal(2)) == 1 ) then
        ! Treat y-direction
        swap = -1.0_rk*(real(normal(2),rk))

        ! LODI Treatment
        ! calculate derivativies along y-direction
        ! This is a quadratic extrapolation
        drhodx =  swap*div1_3*(-8.d0*rhoLodi  + 9.d0*rhoN1Prev  - rhoN2Prev)
        dudx   =  swap*div1_3*(-8.d0*uLodi(1) + 9.d0*uN1Prev(1) - uN2Prev(1))
        dvdx   =  swap*div1_3*(-8.d0*uLodi(2) + 9.d0*uN1Prev(2) - uN2Prev(2))
        dwdx   =  swap*div1_3*(-8.d0*uLodi(3) + 9.d0*uN1Prev(3) - uN2Prev(3))

        ! wave amplitudes
        ! Outgoioming wave L1 Eq (9)
        L1 = (uLodi(2) - swap * cs) &
          & * (cs2 * drhodx - swap * cs * rho0 * dvdx)
        ! Outgoing wave amplitudes L2, L3, L4 and L5 are calculated from Eq (6)
        L2 = - ( vel_bE(1) - uLodi(1) )
        L3 = - ( vel_bE(3) - uLodi(3) )
        L4 = 0.0_rk
        L5 = L1 - 2.0_rk*rho0*swap*cs*(vel_bE(2)-uLodi(2))


        ! LODI Equations
        rhoLodi  = rhoLodi  - ( L4 + 0.5d0*(L5+L1) ) * cs2inv
        uLodi(1) = uLodi(1) - L2
        uLodi(2) = uLodi(2) - (L5-L1)/(rho0*swap*cs)*0.5d0
        uLodi(3) = uLodi(3) - L3
      else if( abs(normal(3)) == 1 ) then
        ! Treat z-direction
        swap = -1.0_rk*(real( normal(3),rk))

        ! LODI Treatment
        ! calculate derivatives in z-direction
        ! This is a quadratic extrapolation
        drhodx =  swap*div1_3*(-8.d0*rhoLodi  + 9.d0*rhoN1Prev  - rhoN2Prev)
        dudx   =  swap*div1_3*(-8.d0*uLodi(1) + 9.d0*uN1Prev(1) - uN2Prev(1))
        dvdx   =  swap*div1_3*(-8.d0*uLodi(2) + 9.d0*uN1Prev(2) - uN2Prev(2))
        dwdx   =  swap*div1_3*(-8.d0*uLodi(3) + 9.d0*uN1Prev(3) - uN2Prev(3))

        ! wave amplitudes
        ! Outgoioming wave L1 Eq (9)
        L1 = (uLodi(3) - swap * cs) &
          & * (cs2 * drhodx - swap * cs * rho0 * dwdx)
        ! Outgoing wave amplitudes L2, L3, L4 and L5 is approached with a
        ! linear relaxation model Eq (22)
        L2 = - ( vel_bE(1) - uLodi(1) )
        L3 = - ( vel_bE(2) - uLodi(2) )
        L4 = 0.0_rk
        L5 = L1 - 2.0_rk*rho0*swap*cs*(vel_bE(3)-uLodi(3))

        ! LODI Equations
        rhoLodi  = rhoLodi  - ( L4 + 0.5d0*(L5+L1) ) * cs2inv
        uLodi(1) = uLodi(1) - L2
        uLodi(2) = uLodi(2) - L3
        uLodi(3) = uLodi(3) - (L5-L1)/(rho0*swap*cs)*0.5d0
      end if


      ! Now store the previous variables for next iteration
      me%elemLvl(iLevel)%lodi(1,   1, iElem) = rhoLodi
      me%elemLvl(iLevel)%lodi(2:4, 1, iElem) = uLodi

      do iDir = 1,QQ
        fTmpPre(iDir) = me%neigh(iLevel)%neighBufferPre(1, iDir+(iElem-1)*QQ)
        fTmpPre_2(iDir) = me%neigh(iLevel)%neighBufferPre(2, iDir+(iElem-1)*QQ)
      end do

      ! 1st neighbor
      ! density
      me%elemLvl(iLevel)%lodi(1, 2, iElem) = sum( fTmpPre(1:QQ) )
      ! Calculate my neighbor velocity for next iteration
      call derVarPos%velFromState( state  = fTmpPre,             &
        &                          iField = iField,              &
        &                          nElems = 1,                   &
        &                          varSys = varSys,              &
        &                          layout = layout,              &
        &                          res    = me%elemLvl(iLevel)   &
        &                                     %lodi(2:4,2,iElem) )

      ! 2nd neighbor
      ! density
      me%elemLvl(iLevel)%lodi(1, 3, iElem) = sum( fTmpPre_2(1:QQ) )
      ! Calculate my neighbor velocity for next iteration
      call derVarPos%velFromState( state  = fTmpPre_2,           &
        &                          iField = iField,              &
        &                          nElems = 1,                   &
        &                          varSys = varSys,              &
        &                          layout = layout,              &
        &                          res    = me%elemLvl(iLevel)   &
        &                                     %lodi(2:4,3,iElem) )

      ! Position in state array
      elemPos = globBC%elemLvl(iLevel)%elem%val( iElem )
      ! Position in bcBuffer
      posInBuffer = globBC%elemLvl(iLevel)%posInBcElemBuf%val( iElem )
      fTmp(1:QQ) = bcBuffer( (posInBuffer-1)*nScalars+varPos(1) : &
        &                    (posInBuffer-1)*nScalars+varPos(QQ)  )

      do iDir = 1, layout%fStencil%QQN
        if( globBC%elemLvl(iLevel)%bitmask%val(iDir, iElem )) then

          ! According to the bounce back rule,
          ! Subtract the equlibrium if its computed from alpha-(inverse of
          ! bitmask) or add the equlibrium if its computed in the direction of
          ! bitmask.
          ! In the Izquiredo paper, this is rho0, the reference density
          fEqPlus = layout%weight(iDir) * 6.0_rk * rho0             &
            &     * (   layout%fStencil%cxDirRK(1, iDir) * uLodi(1) &
            &         + layout%fStencil%cxDirRK(2, iDir) * uLodi(2) &
            &         + layout%fStencil%cxDirRK(3, iDir) * uLodi(3) )

          ! Now assign the values according to bounceback rule
          state( neigh (( idir-1)*nsize+ elempos)+( ifield-1)* qq+ nscalars*0) &
            & = fTmp(layout%fStencil%cxDirInv(iDir)) + fEqPlus
        end if ! bitMask
      end do ! iDir
    end do ! iElem


  end subroutine inlet_nrbc_incomp
! **************************************************************************** !

! ****************************************************************************** !
  !> author: Kannan Masilamani
  !! Characteristic-based non-reflective inlet boundary conditions
  !! @todo: add explaination of these steps
  !!
  !! These boundary conditions are taken from the paper:
  !!   S. Izquierdo and N. Fueyo,
  !!   "Characteristic nonreflecting boundary
  !!   conditions for open boundaries in lattice Boltzmann methods,"
  !!   Physical Review E, vol. 78, no. 46707, 2008.
  !!
  !! > Note: unstable behavior for high omega relaxation parameters.
  !! >       It is recommended to use a sponge layer with a
  !! >       \ref mus_aux_module::mus_setspatialomega "spatial omega"
  !! >       where the omega is decreased towards the boundaries to increase
  !! >       stability
  !!
  !! This subroutine's interface must match the abstract interface definition
  !! [[boundaryRoutine]] in bc/[[mus_bc_header_module]].f90 in order to be
  !! callable via [[boundary_type:fnct]] function pointer.
  subroutine inlet_nrbc( me, state, bcBuffer, globBC, levelDesc, &
    &                    tree, nSize, iLevel, sim_time, neigh,   &
    &                    layout, fieldProp, varPos, nScalars,    &
    &                    varSys, derVarPos, physics, iField,     &
    &                    mixture                                 )
    ! -------------------------------------------------------------------- !
    !> global boundary type
    class(boundary_type) :: me
    !> Current state vector of iLevel
    real(kind=rk), intent(inout) :: state(:)
    !> size of state array ( in terms of elements )
    integer, intent(in) :: nSize
    !> state values of boundary elements of all fields of iLevel
    real(kind=rk), intent(in) :: bcBuffer(:)
    !> iLevel descriptor
    type(tem_levelDesc_type), intent(in) :: levelDesc
    !> Treelm Mesh
    type(treelmesh_type), intent(in) :: tree
    !> fluid parameters and properties
    type(mus_field_prop_type), intent(in) :: fieldProp
    !> stencil layout information
    type(mus_scheme_layout_type), intent(in) :: layout
    !> the level On which this boundary was invoked
    integer, intent(in) :: iLevel
    !> connectivity array corresponding to state vector
    integer, intent(in) :: neigh(:)
    !> global time information
    type(tem_time_type), intent(in)  :: sim_time
    !> pointer to field variable in the state vector
    integer, intent(in) :: varPos(:)
    !> number of Scalars in the scheme var system
    integer, intent(in) :: nScalars
    !> scheme variable system
    type(tem_varSys_type), intent(in) :: varSys
    !> position of derived quantities in varsys
    type(mus_derVarPos_type), intent(in) :: derVarPos
    !> scheme global boundary type
    type(glob_boundary_type), intent(in) :: globBC
    !> scheme global boundary type
    type(mus_physics_type), intent(in) :: physics
    !> current field
    integer, intent(in) :: iField
    !> mixture info
    type(mus_mixture_type), intent(in) :: mixture
    ! -------------------------------------------------------------------- !
    real(kind=rk) :: rhoDef(globBC%nElems(iLevel)) ! Density on boundary element
    real(kind=rk) :: rhoF    ! density on fluid element
    real(kind=rk) :: uxF(globBC%nElems(iLevel)*3)  ! Velocity on fluid element
    real(kind=rk) :: uxN(globBC%nElems(iLevel)*3)  ! Velocity on neighbor
    ! Equilibrium Plus for bounceback
    real(kind=rk) :: fEqPlus
    real(kind=rk) :: inv_rho_phy
    real(kind=rk) :: fTmp( nScalars )
    real(kind=rk) :: fTmpPre( nScalars )
    real(kind=rk) :: fTmpPre_2( nScalars )
    real(kind=rk) :: omega       ! Relaxation parameter
    integer :: iDir, iElem, QQ, invDir, elemPos
    integer :: posInBuffer, bcPress_pos
    real(kind=rk) :: swap
    real(kind=rk) :: L1, L2, L3, L4, L5
    real(kind=rk) :: dudx, dvdx, dwdx, drhodx
    real(kind=rk) :: rhoN1Prev, uN1Prev(3), rhoN2Prev, uN2Prev(3)
    real(kind=rk) :: rhoLodi, uLodi(3)
    integer :: normal(3)
    real(kind=rk) :: vel_b(globBC%nElems(iLevel)*3), inv_vel, vel_bE(3)
    integer :: bcVel_pos
    ! ---------------------------------------------------------------------------

    ! write(dbgUnit(6), *) ''
    ! write(dbgUnit(6), *) ' do characteristic velocity BounceBack for ' &
    !   &                //'bc label: '//trim(globBC%label)

    QQ = layout%fStencil%QQ
    inv_vel = 1.0_rk / physics%fac( iLevel )%vel

    ! position of boundary velocity in varSys
    bcVel_pos = me%bc_states%velocity%varPos
    ! Get velocity
    call varSys%method%val(bcVel_pos)%get_valOfIndex( &
      & varSys  = varSys,                             &
      & time    = sim_time,                           &
      & iLevel  = iLevel,                             &
      & idx     = me%bc_states%velocity               &
      &           %pntIndex%indexLvl(iLevel)          &
      &           %val(1:globBC%nElems(iLevel)),      &
      & nVals   = globBC%nElems(iLevel),              &
      & res     = vel_b                               )

    ! convert physical velocity into LB velocity
    vel_b = vel_b * inv_vel

    do iElem = 1, globBC%nElems( iLevel )

      ! Dirichlet velocity of current element
      vel_bE(:) = vel_b((iElem-1)*3+1:iElem*3)

      ! my density and velocity
      rhoLodi   = me%elemLvl( iLevel )%lodi( 1,   1, iElem )
      uLodi     = me%elemLvl( iLevel )%lodi( 2:4, 1, iElem )
      ! density and velocity of my first neighbor from last iteration
      rhoN1Prev = me%elemLvl( iLevel )%lodi( 1,   2, iElem )
      uN1Prev   = me%elemLvl( iLevel )%lodi( 2:4, 2, iElem )
      ! density and velocity of my secnd neighbor from last iteration
      rhoN2Prev = me%elemLvl( iLevel )%lodi( 1,   3, iElem )
      uN2Prev   = me%elemLvl( iLevel )%lodi( 2:4, 3, iElem )
      ! Use a swapping variable to treat +x and -x directions
      ! Incoming wave amplitudes L2, L3 and L4 in Eq (9) point inwards
      ! and outgoing wave amplitude L1 points outwards.
      ! The derivatives at +x and +y direction are computed by second-order
      ! backward stencil i.e. Eq (15) and at -x and -y, the derivatives
      ! are computed by second-order forward stencil i.e. Eq (21).
      !
      ! Swap is used to switch between forward and backward stencil.
      ! It is positive for backward stencil and negative for forward stencil.
      ! NormalInd refers to incoming direction so use cxDirInv to get normal
      ! pointing away from boundary.
      normal = layout%fStencil%cxDir(:,                                   &
        &                      layout%fStencil%cxDirInv(globBC%normalInd) )
!write(dbgUnit(6),*) 'normal ', normal
      ! x-direction
      if (abs(normal(1)) == 1) then
        swap = -1.0_rk*(real(normal(1),rk))

        ! LODI Treatment
        ! calculate spatial derivatives according to (15) which will be used to
        ! compute the characteristics
        ! This is a quadratic extrapolation
        drhodx =  swap*div1_3*(-8.d0*rhoLodi  + 9.d0*rhoN1Prev  - rhoN2Prev)
        dudx   =  swap*div1_3*(-8.d0*uLodi(1) + 9.d0*uN1Prev(1) - uN2Prev(1))
        dvdx   =  swap*div1_3*(-8.d0*uLodi(2) + 9.d0*uN1Prev(2) - uN2Prev(2))
        dwdx   =  swap*div1_3*(-8.d0*uLodi(3) + 9.d0*uN1Prev(3) - uN2Prev(3))

        ! wave amplitudes
        ! Outgoioming wave L1 Eq (9)
        L1 = (uLodi(1) - swap * cs) &
          & * (cs2 * drhodx - swap * cs * rhoLodi * dudx)
        ! incoming wave amplitudes L2, L3, L4 and L5 is approached with a
        ! linear relaxation model Eq (22)
        L2 = - (vel_bE(2)-uLodi(2))
        L3 = - (vel_bE(3)-uLodi(3))
        ! L5 is valid only if cs = sqrt(kappa*p/rho)
        L4 = 0.0_rk
        L5 = L1 - 2.0_rk*rhoLodi*swap*cs*(vel_bE(1)-uLodi(1))

        ! LODI Equations (18)
        rhoLodi  = rhoLodi  - ( L4 + 0.5d0*(L5+L1) ) * cs2inv
        ! Eq. (18b) for x-velocity
        uLodi(1) = uLodi(1) - (L5-L1)/(rhoLodi*swap*cs)*0.5d0
        uLodi(2) = uLodi(2) - L2  ! equation (18c) for y-direction
        uLodi(3) = uLodi(3) - L3  ! equation (18c) for z-direction
      else if( abs(normal(2)) == 1 ) then
        ! Treat y-direction
        swap = -1.0_rk*(real(normal(2),rk))

        ! LODI Treatment
        ! calculate derivativies along y-direction
        ! This is a quadratic extrapolation
        drhodx =  swap*div1_3*(-8.d0*rhoLodi  + 9.d0*rhoN1Prev  - rhoN2Prev)
        dudx   =  swap*div1_3*(-8.d0*uLodi(1) + 9.d0*uN1Prev(1) - uN2Prev(1))
        dvdx   =  swap*div1_3*(-8.d0*uLodi(2) + 9.d0*uN1Prev(2) - uN2Prev(2))
        dwdx   =  swap*div1_3*(-8.d0*uLodi(3) + 9.d0*uN1Prev(3) - uN2Prev(3))

        ! wave amplitudes
        ! Outgoioming wave L1 Eq (9)
        L1 = (uLodi(2) - swap * cs) &
          & * (cs2 * drhodx - swap * cs * rhoLodi * dvdx)
        ! Outgoing wave amplitudes L2, L3, L4 and L5 are calculated from Eq (6)
        L2 = - ( vel_bE(1) - uLodi(1) )
        L3 = - ( vel_bE(3) - uLodi(3) )
        L4 = 0.0_rk
        L5 = L1 - 2.0_rk*rhoLodi*swap*cs*(vel_bE(2)-uLodi(2))


        ! LODI Equations
        rhoLodi  = rhoLodi  - ( L4 + 0.5d0*(L5+L1) ) * cs2inv
        uLodi(1) = uLodi(1) - L2
        uLodi(2) = uLodi(2) - (L5-L1)/(rhoLodi*swap*cs)*0.5d0
        uLodi(3) = uLodi(3) - L3
      else if( abs(normal(3)) == 1 ) then
        ! Treat z-direction
        swap = -1.0_rk*(real( normal(3),rk))

        ! LODI Treatment
        ! calculate derivatives in z-direction
        ! This is a quadratic extrapolation
        drhodx =  swap*div1_3*(-8.d0*rhoLodi  + 9.d0*rhoN1Prev  - rhoN2Prev)
        dudx   =  swap*div1_3*(-8.d0*uLodi(1) + 9.d0*uN1Prev(1) - uN2Prev(1))
        dvdx   =  swap*div1_3*(-8.d0*uLodi(2) + 9.d0*uN1Prev(2) - uN2Prev(2))
        dwdx   =  swap*div1_3*(-8.d0*uLodi(3) + 9.d0*uN1Prev(3) - uN2Prev(3))

        ! wave amplitudes
        ! Outgoioming wave L1 Eq (9)
        L1 = (uLodi(3) - swap * cs) &
          & * (cs2 * drhodx - swap * cs * rhoLodi * dwdx)
        ! Outgoing wave amplitudes L2, L3, L4 and L5 is approached with a
        ! linear relaxation model Eq (22)
        L2 = - ( vel_bE(1) - uLodi(1) )
        L3 = - ( vel_bE(2) - uLodi(2) )
        L4 = 0.0_rk
        L5 = L1 - 2.0_rk*rhoLodi*swap*cs*(vel_bE(3)-uLodi(3))

        ! LODI Equations
        rhoLodi  = rhoLodi  - ( L4 + 0.5d0*(L5+L1) ) * cs2inv
        uLodi(1) = uLodi(1) - L2
        uLodi(2) = uLodi(2) - L3
        uLodi(3) = uLodi(3) - (L5-L1)/(rhoLodi*swap*cs)*0.5d0
      end if

      ! Now store the previous variables for next iteration
      me%elemLvl(iLevel)%lodi(1,   1, iElem) = rhoLodi
      me%elemLvl(iLevel)%lodi(2:4, 1, iElem) = uLodi

      do iDir = 1,QQ
        fTmpPre(iDir) = me%neigh(iLevel)%neighBufferPre(1, iDir+(iElem-1)*QQ)
        fTmpPre_2(iDir) = me%neigh(iLevel)%neighBufferPre(2, iDir+(iElem-1)*QQ)
      end do

      ! 1st neighbor
      ! density
      me%elemLvl(iLevel)%lodi(1, 2, iElem) = sum( fTmpPre(1:QQ) )
      ! Calculate my neighbor velocity for next iteration
      call derVarPos%velFromState( state  = fTmpPre,             &
        &                          iField = iField,              &
        &                          nElems = 1,                   &
        &                          varSys = varSys,              &
        &                          layout = layout,              &
        &                          res    = me%elemLvl(iLevel)   &
        &                                     %lodi(2:4,2,iElem) )

      ! 2nd neighbor
      ! density
      me%elemLvl(iLevel)%lodi(1, 3, iElem) = sum( fTmpPre_2(1:QQ) )
      ! Calculate my neighbor velocity for next iteration
      call derVarPos%velFromState( state  = fTmpPre_2,           &
        &                          iField = iField,              &
        &                          nElems = 1,                   &
        &                          varSys = varSys,              &
        &                          layout = layout,              &
        &                          res    = me%elemLvl(iLevel)   &
        &                                     %lodi(2:4,3,iElem) )

      ! Position in state array
      elemPos = globBC%elemLvl(iLevel)%elem%val( iElem )
      ! Position in bcBuffer
      posInBuffer = globBC%elemLvl(iLevel)%posInBcElemBuf%val( iElem )
      fTmp(1:QQ) = bcBuffer( (posInBuffer-1)*nScalars+varPos(1) : &
        &                    (posInBuffer-1)*nScalars+varPos(QQ)  )

      do iDir = 1, layout%fStencil%QQN
        if( globBC%elemLvl(iLevel)%bitmask%val(iDir, iElem )) then

          ! According to the bounce back rule,
          ! Subtract the equlibrium if its computed from alpha-(inverse of
          ! bitmask) or add the equlibrium if its computed in the direction of
          ! bitmask.
          ! In the Izquiredo paper, this is rho0, the reference density
          fEqPlus = layout%weight(iDir) * 6.0_rk * rhoLodi  &
            &     * (   layout%fStencil%cxDirRK(1, iDir) * uLodi(1) &
            &         + layout%fStencil%cxDirRK(2, iDir) * uLodi(2) &
            &         + layout%fStencil%cxDirRK(3, iDir) * uLodi(3) )

          ! Now assign the values according to bounceback rule
          state( neigh (( idir-1)*nsize+ elempos)+( ifield-1)* qq+ nscalars*0) &
            & = fTmp(layout%fStencil%cxDirInv(iDir)) + fEqPlus
        end if ! bitMask
      end do ! iDir
    end do ! iElem


  end subroutine inlet_nrbc
! **************************************************************************** !


! ****************************************************************************** !
  !> author: Kannan Masilamani for incompressible flows
  !! Characteristic-based non-reflective open boundary conditions
  !! @todo: add explaination of these steps
  !!
  !! These boundary conditions are taken from the paper:
  !!   S. Izquierdo and N. Fueyo,
  !!   "Characteristic nonreflecting boundary
  !!   conditions for open boundaries in lattice Boltzmann methods,"
  !!   Physical Review E, vol. 78, no. 46707, 2008.
  !!
  !! > Note: unstable behavior for high omega relaxation parameters.
  !! >       It is recommended to use a sponge layer with a
  !! >       \ref mus_aux_module::mus_setspatialomega "spatial omega"
  !! >       where the omega is decreased towards the boundaries to increase
  !! >       stability
  !!
  !! This subroutine's interface must match the abstract interface definition
  !! [[boundaryRoutine]] in bc/[[mus_bc_header_module]].f90 in order to be
  !! callable via [[boundary_type:fnct]] function pointer.
  subroutine outlet_nrbc_incomp( me, state, bcBuffer, globBC, levelDesc, &
    &                            tree, nSize, iLevel, sim_time, neigh,   &
    &                            layout, fieldProp, varPos, nScalars,    &
    &                            varSys, derVarPos, physics, iField,     &
    &                            mixture                                 )
    ! -------------------------------------------------------------------- !
    !> global boundary type
    class(boundary_type) :: me
    !> Current state vector of iLevel
    real(kind=rk), intent(inout) :: state(:)
    !> size of state array ( in terms of elements )
    integer, intent(in) :: nSize
    !> state values of boundary elements of all fields of iLevel
    real(kind=rk), intent(in) :: bcBuffer(:)
    !> iLevel descriptor
    type(tem_levelDesc_type), intent(in) :: levelDesc
    !> Treelm Mesh
    type(treelmesh_type), intent(in) :: tree
    !> fluid parameters and properties
    type(mus_field_prop_type), intent(in) :: fieldProp
    !> stencil layout information
    type(mus_scheme_layout_type), intent(in) :: layout
    !> the level On which this boundary was invoked
    integer, intent(in) :: iLevel
    !> connectivity array corresponding to state vector
    integer, intent(in) :: neigh(:)
    !> global time information
    type(tem_time_type), intent(in)  :: sim_time
    !> pointer to field variable in the state vector
    integer, intent(in) :: varPos(:)
    !> number of Scalars in the scheme var system
    integer, intent(in) :: nScalars
    !> scheme variable system
    type(tem_varSys_type), intent(in) :: varSys
    !> position of derived quantities in varsys
    type(mus_derVarPos_type), intent(in) :: derVarPos
    !> scheme global boundary type
    type(glob_boundary_type), intent(in) :: globBC
    !> scheme global boundary type
    type(mus_physics_type), intent(in) :: physics
    !> current field
    integer, intent(in) :: iField
    !> mixture info
    type(mus_mixture_type), intent(in) :: mixture
    ! -------------------------------------------------------------------- !
    real(kind=rk) :: rhoDef(globBC%nElems(iLevel)) ! Density on boundary element
    real(kind=rk) :: rhoF    ! density on fluid element
    real(kind=rk) :: uxF(globBC%nElems(iLevel)*3)  ! Velocity on fluid element
    real(kind=rk) :: usqB, usqF ! Squared velocity magnitude
    real(kind=rk) :: fPlusFluid ! f Plus for anti-bounceback correction
    ! Equilibrium Plus for anti-bounceback correction
    real(kind=rk) :: fEqPlus
    ! Equlibrium Plus from the old time step
    real(kind=rk) :: fEqPlusFluid
    real(kind=rk) :: inv_rho_phy
    real(kind=rk) :: fTmp( nScalars * globBC%nElems(iLevel) )
    real(kind=rk) :: fTmpPre( nScalars )
    real(kind=rk) :: fTmpPre_2( nScalars )
    real(kind=rk) :: omega       ! Relaxation parameter
    integer :: iDir, iElem, QQ, invDir, elemPos
    integer :: posInBuffer, bcPress_pos
    real(kind=rk) :: sigma, kappa, csMod, Ma_L, lodi_length, k_i, swap
    real(kind=rk) :: L1, L2, L3, L4, L5
    real(kind=rk) :: dudx, dvdx, dwdx, drhodx
    real(kind=rk) :: rhoN1Prev, uN1Prev(3), rhoN2Prev, uN2Prev(3)
    real(kind=rk) :: rhoLodi, uLodi(3)
    integer :: normal(3)
    ! ---------------------------------------------------------------------------

    ! write(dbgUnit(6), *) ''
    ! write(dbgUnit(6), *) ' do non-reflecting pressure_antiBounceBack for ' &
    !   &                //'bc label: '//trim(globBC%label)

    sigma = me%nrbc%sigma
    kappa = me%nrbc%kappa  ! =1 for D2Q9, =1.66666??? for D3Q19
    csMod = me%nrbc%cs_mod
    Ma_L  = me%nrbc%Ma_L
    lodi_length = me%nrbc%lodi_length
    ! Eq (17)
    k_i = sigma*(1.0_rk - Ma_L**2)*cs/lodi_length
!write(dngUnit(6),*) 'k_i ', k_i

    QQ = layout%fStencil%QQ
    inv_rho_phy = 1.0_rk / physics%fac(iLevel)%press * cs2inv

    ! position of boundary pressure in varSys
    bcPress_pos = me%bc_states%pressure%varPos
    ! get pressure variable from spacetime function
    call varSys%method%val(bcPress_pos)%get_valOfIndex( &
      & varSys  = varSys,                               &
      & time    = sim_time,                             &
      & iLevel  = iLevel,                               &
      & idx     = me%bc_states%pressure                 &
      &           %pntIndex%indexLvl(iLevel)            &
      &           %val(1:globBC%nElems(iLevel)),        &
      & nVals   = globBC%nElems(iLevel),                &
      & res     = rhoDef                                )

    ! convert physical pressure into LB density
    rhoDef = rhoDef * inv_rho_phy + rho0

    do iElem = 1, globBC%nElems( iLevel )
      posInBuffer = globBC%elemLvl( iLevel )%posInBcElemBuf%val( iElem )
      fTmp( (iElem-1)*nScalars+1: (iElem-1)*nScalars+QQ ) &
        &       = bcBuffer( (posInBuffer-1)*nScalars+varPos(1) : &
        &                   (posInBuffer-1)*nScalars+varPos(QQ)  )
    end do

    ! Get my velocity
    call derVarPos%velFromState( state  = fTmp,                  &
      &                          iField = iField,                &
      &                          nElems = globBC%nElems(iLevel), &
      &                          varSys = varSys,                &
      &                          layout = layout,                &
      &                          res    = uxF                    )


    do iElem = 1, globBC%nElems( iLevel )

      ! my density and velocity
      rhoLodi   = me%elemLvl( iLevel )%lodi( 1,   1, iElem )
      uLodi     = me%elemLvl( iLevel )%lodi( 2:4, 1, iElem )
      ! density and velocity of my first neighbor from last iteration
      rhoN1Prev = me%elemLvl( iLevel )%lodi( 1,   2, iElem )
      uN1Prev   = me%elemLvl( iLevel )%lodi( 2:4, 2, iElem )
      ! density and velocity of my secnd neighbor from last iteration
      rhoN2Prev = me%elemLvl( iLevel )%lodi( 1,   3, iElem )
      uN2Prev   = me%elemLvl( iLevel )%lodi( 2:4, 3, iElem )
      ! Use a swapping variable to treat +x and -x directions
      ! Outgoing wave amplitudes L2, L3 and L4 in Eq (9) point outwards
      ! and incoming wave amplitude L1 points inwards.
      ! The derivatives at +x and +y direction are computed by second-order
      ! backward stencil i.e. Eq (15) and at -x and -y, the derivatives
      ! are computed by second-order forward stencil i.e. Eq (21).
      !
      ! Swap is used to switch between forward and backward stencil.
      ! It is positive for backward stencil and negative for forward stencil.
      ! NormalInd refers to incoming direction so use cxDirInv to get normal
      ! pointing away from boundary.
      normal = layout%fStencil%cxDir(:,                                   &
        &                      layout%fStencil%cxDirInv(globBC%normalInd) )
!write(dbgUnit(6),*) 'normal ', normal
      ! x-direction
      if (abs(normal(1)) == 1) then
        swap = real(normal(1),rk)

        ! LODI Treatment
        ! calculate spatial derivatives according to (15) which will be used to
        ! compute the characteristics
        ! This is a quadratic extrapolation
        drhodx =  swap*div1_3*(8.d0*rhoLodi  - 9.d0*rhoN1Prev  + rhoN2Prev)
        dudx   =  swap*div1_3*(8.d0*uLodi(1) - 9.d0*uN1Prev(1) + uN2Prev(1))
        dvdx   =  swap*div1_3*(8.d0*uLodi(2) - 9.d0*uN1Prev(2) + uN2Prev(2))
        dwdx   =  swap*div1_3*(8.d0*uLodi(3) - 9.d0*uN1Prev(3) + uN2Prev(3))

        ! wave amplitudes
        ! Incoming wave L1 is approached with a linear relaxation model (16)
        ! \todo 21012021 KM:check whether swap needs to be multiplied in L1
        L1 = k_i * cs2 * ( rhoLodi-rhoDef(iElem) )
        ! Outgoing wave amplitudes L2, L3, L4 and L5
        L2 = uLodi(1) * dvdx  ! equation (9) second row
        L3 = uLodi(1) * dwdx  ! equation (9) second row adapted for z-direction
        ! equation (9) third row with EOS equation (3) used
        L4 = 0.0_rk
        ! equation (9) fourth row with EOS equation (3) used
        ! MH: Here rhoLodi should be replaced with rho0
        L5 = ( uLodi(1) + swap * cs )                     &
          & * ( cs2 * drhodx + swap * cs * rho0 * dudx )

        ! LODI Equations (18)
        rhoLodi  = rhoLodi  - ( L4 + 0.5d0*(L5+L1) ) * cs2inv
        ! Eq. (18b) for x-velocity
        uLodi(1) = uLodi(1) - (L5-L1)/(rho0*swap*cs)*0.5d0
        uLodi(2) = uLodi(2) - L2  ! equation (18c) for y-direction
        uLodi(3) = uLodi(3) - L3  ! equation (18c) for z-direction
      else if( abs(normal(2)) == 1 ) then
        ! Treat y-direction
        swap = real(normal(2),rk)

        ! LODI Treatment
        ! calculate derivativies along y-direction
        ! This is a quadratic extrapolation
        drhodx =  swap*div1_3*(8.d0*rhoLodi  - 9.d0*rhoN1Prev  + rhoN2Prev)
        dudx   =  swap*div1_3*(8.d0*uLodi(1) - 9.d0*uN1Prev(1) + uN2Prev(1))
        dvdx   =  swap*div1_3*(8.d0*uLodi(2) - 9.d0*uN1Prev(2) + uN2Prev(2))
        dwdx   =  swap*div1_3*(8.d0*uLodi(3) - 9.d0*uN1Prev(3) + uN2Prev(3))

        ! wave amplitudes
        L1 = k_i * cs2 * ( rhoLodi-rhoDef(iElem) )
        L2 = ( uLodi(2) ) * dudx
        L3 = ( uLodi(2) ) * dwdx
        L4 = 0.0_rk
        L5 = ( uLodi(2) + swap * cs )                         &
          &* ( cs2 * drhodx + swap * cs * rho0 * dvdx )

        ! LODI Equations
        rhoLodi  = rhoLodi  - ( L4 + 0.5d0*(L5+L1) ) * cs2inv
        uLodi(1) = uLodi(1) - L2
        uLodi(2) = uLodi(2) - (L5-L1)/(rho0*swap*cs)*0.5d0
        uLodi(3) = uLodi(3) - L3
      else if( abs(normal(3)) == 1 ) then
        ! Treat z-direction
        swap = real( normal(3),rk )

        ! LODI Treatment
        ! calculate derivatives in z-direction
        ! This is a quadratic extrapolation
        drhodx =  swap*div1_3*(8.d0*rhoLodi  - 9.d0*rhoN1Prev  + rhoN2Prev)
        dudx   =  swap*div1_3*(8.d0*uLodi(1) - 9.d0*uN1Prev(1) + uN2Prev(1))
        dvdx   =  swap*div1_3*(8.d0*uLodi(2) - 9.d0*uN1Prev(2) + uN2Prev(2))
        dwdx   =  swap*div1_3*(8.d0*uLodi(3) - 9.d0*uN1Prev(3) + uN2Prev(3))

        ! wave amplitudes
        L1 = k_i * cs2 * ( rhoLodi-rhoDef(iElem) )
        L2 = ( uLodi(3) ) * dudx
        L3 = ( uLodi(3) ) * dvdx
        L4 = 0.0_rk
        L5 = ( uLodi(3) + swap * cs )                    &
          &* ( cs2 * drhodx + swap * cs * rho0 * dwdx )

        ! LODI Equations
        rhoLodi  = rhoLodi  - ( L4 + 0.5d0*(L5+L1) ) * cs2inv
        uLodi(1) = uLodi(1) - L2
        uLodi(2) = uLodi(2) - L3
        uLodi(3) = uLodi(3) - (L5-L1)/(rho0*swap*cs)*0.5d0
      end if

      ! Now store the previous variables for next iteration
      me%elemLvl(iLevel)%lodi(1,   1, iElem) = rhoLodi
      me%elemLvl(iLevel)%lodi(2:4, 1, iElem) = uLodi

      do iDir = 1,QQ
        fTmpPre(iDir) = me%neigh(iLevel)%neighBufferPre(1, iDir+(iElem-1)*QQ)
        fTmpPre_2(iDir) = me%neigh(iLevel)%neighBufferPre(2, iDir+(iElem-1)*QQ)
      end do

      ! 1st neighbor
      ! density
      me%elemLvl(iLevel)%lodi(1, 2, iElem) = sum( fTmpPre(1:QQ) )
      ! Calculate my neighbor velocity for next iteration
      call derVarPos%velFromState( state  = fTmpPre,             &
        &                          iField = iField,              &
        &                          nElems = 1,                   &
        &                          varSys = varSys,              &
        &                          layout = layout,              &
        &                          res    = me%elemLvl(iLevel)