! Copyright (c) 2012-2020 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]]
!!
?? include 'header/lbm_macros.inc'
?? include 'header/lbm_deriveMacros.inc'
?? include 'treelm/source/logMacros.inc'
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 :: 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( ?FETCH?( iDir, iField, elemPos, QQ, nScalars, nSize, neigh )) &
            & = 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( ?FETCH?( iDir, iField, elemPos, QQ, nScalars, nSize, neigh )) &
            & = 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( ?FETCH?( iDir, iField, elemPos, QQ, nScalars, nSize, neigh )) &
            & = 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( ?FETCH?( iDir, iField, elemPos, QQ, nScalars, nSize, neigh )) &
            & = 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( ?FETCH?( iDir, iField, elemPos, QQ, nScalars, nSize, neigh )) = &
          & 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( ?FETCH?( iDir, iField, elemPos, QQ, nScalars,nSize, neigh ))  = &
         ! 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(?FETCH?(invDir, iField, elemPos, QQ, nScalars,nSize,neigh))&
            & + 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( ?FETCH?( iDir, iField, elemPos, QQ, nScalars,nSize, neigh))  = &
        & 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 ------------------------------------------------
?? IF (loglvl >= 10) THEN
write(dbgUnit(1),*) 'iLink: ', iLink, ', iElem: ', iElem
write(dbgUnit(1),*) 'fTmp_1: ', fTmp_1, ', fTmp_2: ', fTmp_2
?? ENDIF
! 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( ?FETCH?( iDir, iField, elemPos, QQ, nScalars, nSize, neigh ) &
          & )  = &
          ! 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 ------------------------------------------------
?? IF (loglvl >= 10) THEN
    do iElem = 1, globBC%nElems(iLevel)
      elemPos = globBC%elemLvl(iLevel)%elem%val( iElem )
write(dbgUnit(1),*) 'iElem ', iElem,' elemID ',levelDesc%total(   &
  &                                            globBC%elemLvl(iLevel)%elem(iElem))
      do iDir = 1, layout%fStencil%QQN
        if( globBC%elemLvl(iLevel)%bitmask%val( iDir, iElem )) then
write(dbgUnit(1),*) iDir, 'state ', &
  & state( ?FETCH?( iDir, iField, elemPos, QQ, nScalars, nSize, neigh )), &
  & ' fEq0 ', fEq0((iElem-1)*QQ+iDir), ' fEq ', fEq( (iElem-1)*QQ+invDir ), &
  & ' fTmp ', fTmp((iElem-1)*QQ+iDir)
        end if
      end do
    end do
?? ENDIF
! 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(?FETCH?(iDir, iField, globBC%elemLvl(iLevel)%elem%val(iElem), QQ, nScalars,nSize, neigh) ) = &
          ! 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(                                                      &
&?FETCH?(iDir,iField,globBC%elemLvl(iLevel)%elem%val(iElem),QQ,nScalars,nSize,neigh) &
&         )  =                                                   &
            &  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(?FETCH?( iDir, iField, elemPos, QQ, nScalars,nSize, neigh )) = &
          ! antibounceback term
! & - state( ?FETCH?( iDir, iField, globBC%elemLvl(iLevel)%elem%val( iElem ), QQ, nScalars, neigh ) )    &
          ! 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(                      &
! & ?FETCH?( iDir, iField, globBC%elemLvl(iLevel)%elem%val( iElem ), QQ, nScalars, nSize,neigh ))
! --------------------------------------------------------------------------

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

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


! ****************************************************************************** !
  !> author: Manuel Hasert
  !! 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( 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) :: fEqPlus, fPlus
    real(kind=rk) :: usq
    real(kind=rk) :: rhoDef(globBC%nElems( iLevel ))
    real(kind=rk) :: Ma, swap, sigma, kappa, K_mod
    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), csMod
    real(kind=rk) :: fTmpPost( nScalars )
    real(kind=rk) :: fTmpPost_2( nScalars )
    real(kind=rk) :: inv_rho_phy
    integer :: iELem, iDir, QQ, elemPos, posInBuffer
    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     = rhoDef                                )

    rhoDef = rhoDef * inv_rho_phy + rho0

    sigma = me%nrbc%sigma
    kappa = me%nrbc%kappa  ! =1 for D2Q9, =1.66666??? for D3Q19
    csMod = me%nrbc%cs_mod

    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
      swap = 1._rk
      ! x-direction
      if( globBC%elemLvl( iLevel )%normal%val( 1, iElem ) /= 0 ) then
        swap = 1._rk*(-real( globBC%elemLvl( iLevel )%normal%val( 1, iElem ),rk))

        ! LODI Treatment
        Ma   = abs(uLodi(1))*csInv ! X direction
        ! Equation (17) for the linear relaxation model constant
        ! lodi_length= distance outlet-inlet
        K_mod = swap*sigma*(1.d0 - Ma*Ma)*cs/(me%nrbc%lodi_length)

        ! 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)
        L1 = K_mod*csMod*csMod*(rhoLodi-rhoDef(iElem))
        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 = (uLodi(1))*csMod*csMod*drhodx*(1.d0-1.d0/kappa)
        ! equation (9) fourth row with EOS equation (3) used
        ! MH: Here rhoLodi should be replaced with rho0
        L5 = ( uLodi(1) + swap * csMod )                                       &
          &* ( csMod * csMod * drhodx / kappa + swap * csMod * dudx )

        ! LODI Equations (18)
        rhoLodi  = rhoLodi  - (L4+0.5d0*(L5+L1))/(csMod*csMod)
        ! Eq. (18b) for x-velocity
        uLodi(1) = uLodi(1) - ((L5-L1))/(swap*csMod)*0.5d0
        uLodi(2) = uLodi(2) - L2  ! equation (18c) for y-direction
        uLodi(3) = uLodi(3) - L3  ! equation (18c) for z-direction
      else
        ! Treat y-direction
        if( globBC%elemLvl( iLevel )%normal%val( 2, iElem ) /= 0 ) then
          swap = 1._rk*(-real( globBC%elemLvl( iLevel )%normal%val(2,iElem),rk))

          ! LODI Treatment
          Ma   = abs(uLodi(2))*csInv ! X direction
          ! lodi_length= distance outlet-inlet
          K_mod = swap*sigma*(1.d0 - Ma*Ma)*cs/(me%nrbc%lodi_length)

          ! calculate derivatives
          ! 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_mod*csMod*csMod*(rhoLodi-rhoDef(iElem))
          L2 = ( uLodi(2) ) * dudx
          L3 = ( uLodi(2) ) * dwdx
          L4 = ( uLodi(2) ) * csMod * csMod * drhodx * ( 1.d0 - 1.d0 / kappa )
          L5 = ( uLodi(2) + swap * csMod )                                     &
            &* ( csMod * csMod * drhodx / kappa                                &
            &+   swap * csMod * dvdx )

          ! LODI Equations
          rhoLodi  = rhoLodi  - (L4+0.5d0*(L5+L1))/(csMod*csMod)
          uLodi(1) = uLodi(1) - L2
          uLodi(2) = uLodi(2) - ((L5-L1))/(swap*csMod)*0.5d0
          uLodi(3) = uLodi(3) - L3
        else
          ! Treat z-direction
          if( globBC%elemLvl( iLevel )%normal%val( 3, iElem ) /= 0 ) then
            swap = 1._rk*(-real( globBC%elemLvl(iLevel)%normal%val(3,iElem),rk ))

            ! LODI Treatment
            Ma   = abs(uLodi(3))*csInv ! X direction
            !lodi_length= distance outlet-inlet
            K_mod = swap*sigma*(1.d0 - Ma*Ma)*cs/(me%nrbc%lodi_length)

            ! calculate derivatives
           ! 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_mod*csMod*csMod*(rhoLodi-rhoDef(iElem))
            L2 = ( uLodi(3) ) * dudx
            L3 = ( uLodi(3) ) * dvdx
            L4 = ( uLodi(3) ) * csMod * csMod * drhodx * ( 1.d0 - 1.d0 / kappa )
            L5 = ( uLodi(3) +swap * csMod )                                    &
              &* ( csMod * csMod * drhodx / kappa                              &
              &+   swap * csMod * dwdx )

            ! LODI Equations
            rhoLodi  = rhoLodi  - (L4+0.5d0*(L5+L1))/(csMod*csMod)
            uLodi(1) = uLodi(1) - L2
            uLodi(2) = uLodi(2) - L3
            uLodi(3) = uLodi(3) - ((L5-L1))/(swap*csMod)*0.5d0
          end if
        end if
      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

      ! density and velocity of my first neighbor
      ! MH: Velocity and Density should be calculated with access as SAVE
      ! For pull, this is fine however, as IDX and SAVE coincide
      posInBuffer = globBC%elemLvl( iLevel )%posInBcElemBuf%val( iElem )
      fTmpPost(1:QQ) = bcBuffer( (posInBuffer-1)*nScalars+varPos(1) : &
        &                        (posInBuffer-1)*nScalars+varPos(QQ)  )
      ! @todo: fTmpPost_2 is probably wrong
      do iDir = 1,QQ
        fTmpPost_2(iDir) = me%neigh( iLevel )%neighBufferPost( 1,iDir+(iElem-1)*QQ )
      end do

      ! density
      me%elemLvl( iLevel )%lodi( 1, 2, iElem ) = sum( fTmpPost(1:QQ) )
      ! velocity
      call derVarPos%velFromState( state  = fTmpPost,   &
        &                          iField = iField,     &
        &                          nElems = 1,      &
        &                          varSys = varSys, &
        &                          layout = layout, &
        &                          res    = me%elemLvl(iLevel)%lodi(2:4,2,iElem) )

      ! density and velocity of my secnd neighbor
      me%elemLvl( iLevel )%lodi( 1, 3, iElem )  = sum( fTmpPost_2(1:QQ) )
      call derVarPos%velFromState( state  = fTmpPost_2, &
        &                          iField = iField,     &
        &                          nElems = 1,          &
        &                          varSys = varSys,     &
        &                          layout = layout,     &
        &                          res    = me%elemLvl(iLevel)%lodi(2:4,3,iElem) )

      usq = uLodi(1)*uLodi(1) + uLodi(2)*uLodi(2) + uLodi(3)*uLodi(3)

      elemPos = globBC%elemLvl(iLevel)%elem%val( iElem )
      ! Now do the anti-bounce back
      do iDir = 1,layout%fStencil%QQN
        if( globBC%elemLvl(iLevel)%bitmask%val(iDir, iElem )) then
          ! MH: This calculation must be the same as for the pressure_antiBounceBack,
          !     except that we are using the Lodi values
          ! Calc the correction term and the Dirichlet (Eq) term
          ! MH: It makes no difference to use either iDir or inv(iDir)
          fEqPlus = layout%weight(iDir)*rhoLodi                                &
            &     + 4.5_rk*layout%weight(iDir)*((         &
            &       layout%fStencil%cxDirRK( 1, layout%fStencil%             &
            &                                        cxDirInv(iDir ))*uLodi(1) &
            &     + layout%fStencil%cxDirRK( 2, layout%fStencil%             &
            &                                        cxDirInv(iDir ))*uLodi(2) &
            &     + layout%fStencil%cxDirRK( 3, layout%fStencil%             &
            &                                        cxDirInv(iDir ))*uLodi(3) &
            &     )**2 - div1_3*usq)

          fPlus = 0.5_rk*(                                                     &
&  state( ?FETCH?( iDir, iField, elemPos, QQ, nScalars, nSize,neigh ) ) &
&+ state( ?FETCH?( layout%fStencil%cxDirInv(iDir), iField, elemPos, QQ, nScalars, nSize,neigh ) &
&            ))

          ! Now assign the values
state( ?FETCH?( iDir, iField, elemPos, QQ, nScalars, nSize,neigh ) )  =     &
          ! antibounceback term
  & - state( ?FETCH?( iDir, iField, elemPos, QQ, nScalars, nSize,neigh ) )  &
  & + 2._rk*fEqPlus ! Dirichlet pressure term
       ! MH: Correciton term increases error!
!       + (2._rk - omega )*( fPlus - fEqPlus ) ! Correction term
        end if ! bitMask
      end do ! iDir
    end do

  end subroutine outlet_nrbc
! ****************************************************************************** !


! ****************************************************************************** !
  !> author: Manuel Hasert
  !! Characteristic-based non-reflective open boundary conditions
  !!
  !! 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_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) :: rhoDef(globBC%nElems( iLevel ))
    real(kind=rk) :: Ma, swap, sigma, kappa, K_mod
    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), csMod
    real(kind=rk) :: fEq( layout%fStencil%QQ )
    real(kind=rk) :: fTmpPost( nScalars )
    real(kind=rk) :: fTmpPost_2( nScalars )
    real(kind=rk) :: inv_rho_phy
    integer :: iELem, iDir, QQ
    integer :: posInBuffer, 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     = rhoDef                                )

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

    sigma = me%nrbc%sigma
    kappa = me%nrbc%kappa  ! =1 for D2Q9, D3Q19
    csMod = me%nrbc%cs_mod

    do iElem = 1, globBC%nElems( iLevel )
      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
      swap = 1._rk
      ! x-direction
      if( globBC%elemLvl( iLevel )%normal%val( 1, iElem ) /= 0 ) then
        ! Use a swapping variable to treat +x and -x directions
        swap = 1._rk*(-real( globBC%elemLvl( iLevel )%normal%val( 1, iElem ),kind=rk))

        ! LODI Treatment
        Ma   = abs(uLodi(1))*csInv ! X direction
        ! Equation (17) for the linear relaxation model constant
        ! lodi_length= distance outlet-inlet
        K_mod = swap*sigma*(1.d0 - Ma*Ma)*cs/(me%nrbc%lodi_length)

        ! 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)
        L1 = K_mod*csMod*csMod*(rhoLodi-rhoDef(iElem))
        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 = (uLodi(1))*csMod*csMod*drhodx*(1.d0-1.d0/kappa)
        ! equation (9) fourth row with EOS equation (3) used
        ! MH: Here rhoLodi should be replaced with rho0
        L5 = ( uLodi(1) + swap * csMod )                                       &
          &* ( csMod * csMod * drhodx / kappa                                  &
          &+   swap * csMod * dudx )

        ! LODI Equations (18)
        rhoLodi  = rhoLodi  - (L4+0.5d0*(L5+L1))/(csMod*csMod)
        ! Eq. (18b) for x-velocity
        uLodi(1) = uLodi(1) - ((L5-L1))/(swap*csMod)*0.5d0
        uLodi(2) = uLodi(2) - L2  ! equation (18c) for y-direction
        uLodi(3) = uLodi(3) - L3  ! equation (18c) for z-direction
      else
        ! Treat y-direction
        if( globBC%elemLvl( iLevel )%normal%val( 2, iElem ) /= 0 ) then
          swap = 1._rk*(-real( globBC%elemLvl( iLevel )%normal%val(2,iElem),kind=rk ))

          ! LODI Treatment
          Ma   = abs(uLodi(2))*csInv ! X direction
          ! lodi_length= distance outlet-inlet
          K_mod = swap*sigma*(1.d0 - Ma*Ma)*cs/(me%nrbc%lodi_length)

          ! calculate derivatives
          ! 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_mod*csMod*csMod*(rhoLodi-rhoDef(iElem))
          L2 = (uLodi(2))*dudx
          L3 = (uLodi(2))*dwdx
          L4 = (uLodi(2))*csMod*csMod*drhodx*(1.d0-1.d0/kappa)
          L5 = ( uLodi(2) + swap * csMod )                                     &
            &* ( csMod * csMod * drhodx / kappa                                &
            &+   swap * csMod * dvdx )

          ! LODI Equations
          rhoLodi  = rhoLodi  - (L4+0.5d0*(L5+L1))/(csMod*csMod)
          uLodi(1) = uLodi(1) - L2
          uLodi(2) = uLodi(2) - ((L5-L1))/(swap*csMod)*0.5d0
          uLodi(3) = uLodi(3) - L3
        else
          ! Treat z-direction
          if( globBC%elemLvl( iLevel )%normal%val( 3, iElem ) /= 0 ) then
            swap = 1._rk*(-real(globBC%elemLvl( iLevel )%normal%val(3,iElem),kind=rk))

            ! LODI Treatment
            Ma   = abs(uLodi(3))*csInv ! X direction
            !lodi_length= distance outlet-inlet
            K_mod = swap*sigma*(1.d0 - Ma*Ma)*cs/(me%nrbc%lodi_length)

            ! calculate derivatives
           ! 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_mod*csMod*csMod*(rhoLodi-rhoDef(iElem))
            L2 = (uLodi(3))*dudx
            L3 = (uLodi(3))*dvdx
            L4 = (uLodi(3))*csMod*csMod*drhodx*(1.d0-1.d0/kappa)
            L5 = ( uLodi(3) + swap * csMod )                                   &
              &* ( csMod * csMod * drhodx / kappa                              &
              &+   swap * csMod * dwdx )

            ! LODI Equations
            rhoLodi  = rhoLodi  - (L4+0.5d0*(L5+L1))/(csMod*csMod)
            uLodi(1) = uLodi(1) - L2
            uLodi(2) = uLodi(2) - L3
            uLodi(3) = uLodi(3) - ((L5-L1))/(swap*csMod)*0.5d0
          end if
        end if
      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
      ! density and velocity of my first neighbor
      ! MH: Velocity and Density should be calculated with access as SAVE
      ! For pull, this is fine however, as IDX and SAVE coincide
      posInBuffer = globBC%elemLvl( iLevel )%posInBcElemBuf%val( iElem )
      fTmpPost(1:QQ) = bcBuffer( (posInBuffer-1)*nScalars+varPos(1) : &
        &                        (posInBuffer-1)*nScalars+varPos(QQ)  )
      do iDir = 1, QQ
        fTmpPost_2(iDir) = me%neigh( iLevel )%neighBufferPost( 1,iDir+(iElem-1)*QQ)
      end do
      !density
      me%elemLvl( iLevel )%lodi( 1, 2, iElem ) = sum( fTmpPost(1:QQ) )
      !velocity
      call derVarPos%velFromState( state  = fTmpPost,   &
        &                          iField = iField,     &
        &                          nElems = 1,      &
        &                          varSys = varSys, &
        &                          layout = layout, &
        &                          res    = me%elemLvl(iLevel)%lodi(2:4,2,iElem) )

      ! density and velocity of my secnd neighbor
      me%elemLvl( iLevel )%lodi( 1, 3, iElem )  = sum( fTmpPost_2(1:QQ) )
      call derVarPos%velFromState( state  = fTmpPost_2, &
        &                          iField = iField,     &
        &                          nElems = 1,          &
        &                          varSys = varSys,     &
        &                          layout = layout,     &
        &                          res    = me%elemLvl(iLevel)%lodi(2:4,3,iElem) )

      ! Calculate the equilibrium distribution
      call derVarPos%equilFromState( state  = fTmpPost_2, &
        &                            iField = iField,     &
        &                            nElems = 1,          &
        &                            varSys = varSys,     &
        &                            layout = layout,     &
        &                            res     = fEq )

      ! Now do the anti-bounce back
      do iDir = 1, layout%fStencil%QQN
        if( globBC%elemLvl(iLevel)%bitmask%val(iDir, iElem )) then
          ! Now assign the values
          state(                                                                &
& ?FETCH?( iDir, iField, globBC%elemLvl( iLevel )%elem%val( iElem ), QQ, nScalars, nSize,neigh ) &
&         ) = fEq( iDir )
        end if
      end do ! iDir
    end do

  end subroutine outlet_nrbc_eq
! ****************************************************************************** !


! ****************************************************************************** !
  !> author: Kannan Masilamani
  !! Moment based open boundary condition from Sam Bennent PhD thesis
  !! "A Lattice Boltzmann Model for Diffusion of Binary Gas Mixtures"
  !! for weakly compressible LBM model
  !!
  !! Usage
  !! -----
  !!```lua
  !!boundary_condition = {
  !! { label = 'inlet',
  !!    kind = 'pressure_momentsbased',
  !!    pressure = 'p0' }
  !!}
  !! variable = {
  !!   name = 'p0',
  !!   ncomponents = 1,
  !!   vartype = 'st_fun',
  !!   st_fun = 1.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 pressure_momentsbased( 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) :: pressB( globBC%nElems(iLevel) ), rho, vel(3), inv_press
    integer :: iELem, iDir
    integer :: nLinks, iLink, elemPos, QQ
    integer, allocatable :: missing_links(:)
    real(kind=rk), allocatable :: rhs(:)
    real(kind=rk), allocatable :: unKnown_fTmp(:)
    real(kind=rk) :: fTmp( layout%fStencil%QQ )
    real(kind=rk) :: moments( layout%fStencil%QQ )
    real(kind=rk) :: first_moments(3), second_moments(6), third_moments(6), &
      &              fourth_moments(3)
    integer :: nMom(4)
    integer :: posInBuffer, bcPress_pos
    ! ---------------------------------------------------------------------------
    QQ = layout%fStencil%QQ
    inv_press = 1.0_rk / physics%fac( iLevel )%press

!    write(*,*) 'Boundary label ', trim(me%label)
    ! 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     = pressB                                )

    ! If physical quantities are given, transform to lattice units by division
    ! with the conversion factor
    pressB = pressB * inv_press

    nMom(1) = 1 + size(layout%moment%first_moments)
    nMom(2) = nMom(1) + size(layout%moment%second_moments)
    nMom(3) = nMom(2) + size(layout%moment%third_moments)
    nMom(4) = nMom(3) + size(layout%moment%fourth_moments)

    do iElem = 1, globBC%nElems( iLevel )
      fTmp = 0.0_rk
      nLinks = me%elemLvl(iLevel)%moments(iElem)%nUnKnownPdfs
      allocate(missing_links(nLinks))
      allocate(unKnown_fTmp(nLinks))
      allocate(rhs(nLinks))
      elemPos = globBC%elemLvl(iLevel)%elem%val( iElem )

      iLink = 0
      do iDir = 1, layout%fStencil%QQN
        if( globBC%elemLvl(iLevel)%bitmask%val( iDir, iElem )) then
          iLink = iLink + 1
          missing_links(iLink) = iDir
        else
          fTmp( iDir ) = &
            &state( ?FETCH?( iDir, iField, elemPos, QQ, nScalars, nSize, neigh ) )
        end if
      end do
      fTmp( layout%fStencil%restPosition ) = &
&state( ?FETCH?( layout%fStencil%restPosition, iField, elemPos, QQ, nScalars, nSize, neigh ) )

      ! density from defined pressure
      rho = cs2inv*pressB(iElem) + rho0

      !local velocity
      posInBuffer = globBC%elemLvl( iLevel )%posInBcElemBuf%val( iElem )
      vel(1) = sum ( bcBuffer( (posInBuffer-1)*nScalars+varPos(1) :  &
        &                      (posInBuffer-1)*nScalars+varPos(QQ) ) &
        &            * layout%fStencil%cxDir(1,:) )
      vel(2) = sum ( bcBuffer( (posInBuffer-1)*nScalars+varPos(1) :  &
        &                      (posInBuffer-1)*nScalars+varPos(QQ) ) &
        &            * layout%fStencil%cxDir(2,:) )
      vel(3) = sum ( bcBuffer( (posInBuffer-1)*nScalars+varPos(1) :  &
        &                      (posInBuffer-1)*nScalars+varPos(QQ) ) &
        &            * layout%fStencil%cxDir(3,:) )
      vel = vel/rho

      ! moments for in-compressible model
      first_moments = rho*vel !momX, momY, momZ
      second_moments = (/ pressB(iElem) + rho*vel(1)*vel(1), & !momXX
        &                 pressB(iElem) + rho*vel(2)*vel(2), & !momYY
        &                 pressB(iElem) + rho*vel(3)*vel(3), & !momZZ
        &                 rho*vel(1)*vel(2),                 & !momXY
        &                 rho*vel(2)*vel(3),                 & !momYZ
        &                 rho*vel(3)*vel(1) /)                 !momXZ
      third_moments = (/ cs2*rho*vel(2), & !momXXY
        &                cs2*rho*vel(3), & !momXXZ
        &                cs2*rho*vel(1), & !momYYX
        &                cs2*rho*vel(3), & !momYYZ
        &                cs2*rho*vel(1), & !momZZX
        &                cs2*rho*vel(2) /) !momZZY
      !fourth_moments = cs2*pressB(iElem) !momXXYY, momYYZZ, momZZXX
      fourth_moments = (/ &
        & cs2*(pressB(iElem)+rho*( vel(1)**2 + vel(2)**2 - 0.5_rk*vel(3)**2 )), &!momXXYY
        & cs2*(pressB(iElem)+rho*(-0.5_rk*vel(1)**2 + vel(2)**2 + vel(3)**2 )), &!momYYZZ
        & cs2*(pressB(iElem)+rho*( vel(1)**2 - 0.5_rk*vel(2)**2 + vel(3)**2 ))  &!momZZXX
        & /)

      ! moments for pressure-wall corner
      ! KM: for pressure-velocity, velocity at boundary node must be defined
      moments = 0.0_rk
      moments(1) = rho
      moments(2:nMom(1)) = &
        & first_moments(layout%moment%first_moments(:))
      moments(nMom(1)+1:nMom(2)) = &
        & second_moments(layout%moment%second_moments(:))
      moments(nMom(2)+1:nMom(3)) = &
        & third_moments(layout%moment%third_moments(:))
      moments(nMom(3)+1:nMom(4)) = &
        & fourth_moments(layout%moment%fourth_moments(:))

      ! compute rhs of lse to compute unknown pdfs
      do iLink = 1, nLinks
        rhs(iLink) = &
          & moments( me%elemLvl(iLevel)%moments(iElem)%knownMom_pos(iLink) ) &
          ! convert pdf to moments
          & - dot_product(layout%moment%toMoments%A(                        &
          &    me%elemLvl(iLevel)%moments(iElem)%knownMom_pos(iLink), :), fTmp)
      end do

      ! unknown pdfs
      unKnown_fTmp = matmul( &
        & me%elemLvl(iLevel)%moments(iElem)%unKnownPdfs_MatInv, rhs )

      do iLink = 1, nLinks
        state(                                                                 &
          & ?FETCH?( missing_links(iLink), iField, elemPos, QQ, nScalars, nSize,neigh ) &
          & )  = unKnown_fTmp(iLink)
      end do

     deallocate(missing_links)
     deallocate(unKnown_fTmp)
     deallocate(rhs)
   end do

  end subroutine pressure_momentsbased
! ****************************************************************************** !



! ****************************************************************************** !
  !> author: Kannan Masilamani
  !! Moment based open boundary condition from Sam Bennent PhD thesis
  !! "A Lattice Boltzmann Model for Diffusion of Binary Gas Mixtures"
  !! for incompressible LBM model
  !!
  !! Usage
  !! -----
  !!```lua
  !!boundary_condition = {
  !! { label = 'inlet',
  !!    kind = 'pressure_momentsbased',
  !!    pressure = 'p0'}
  !!}
  !! variable = {
  !!   name = 'p0',
  !!   ncomponents = 1,
  !!   vartype = 'st_fun',
  !!   st_fun = 1.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 pressure_momentsbased_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) :: pressB( globBC%nElems(iLevel) ), rho, vel(3)
    real(kind=rk) :: inv_press
    integer :: iELem, iDir
    integer :: nLinks, iLink, QQ, elemPos
    integer, allocatable :: missing_links(:)
    real(kind=rk), allocatable :: rhs(:)
    real(kind=rk), allocatable :: unKnown_fTmp(:)
    real(kind=rk) :: fTmp( layout%fStencil%QQ )
    real(kind=rk) :: moments( layout%fStencil%QQ )
    real(kind=rk) :: first_moments(3), second_moments(6), third_moments(6), &
      &              fourth_moments(3)
    integer :: nMom(4)
    integer :: posInBuffer, bcPress_pos
    ! ---------------------------------------------------------------------------
    QQ = layout%fStencil%QQ
    inv_press = 1.0_rk / physics%fac( iLevel )%press

    ! 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     = pressB                                )

    ! If physical quantities are given, transform to lattice units by division
    ! with the conversion factor
    pressB = pressB * inv_press

    nMom(1) = 1 + size(layout%moment%first_moments)
    nMom(2) = nMom(1) + size(layout%moment%second_moments)
    nMom(3) = nMom(2) + size(layout%moment%third_moments)
    nMom(4) = nMom(3) + size(layout%moment%fourth_moments)

    do iElem = 1, globBC%nElems( iLevel )
      fTmp = 0.0_rk
      nLinks = me%elemLvl(iLevel)%moments(iElem)%nUnKnownPdfs
      allocate(missing_links(nLinks))
      allocate(unKnown_fTmp(nLinks))
      allocate(rhs(nLinks))
      elemPos = globBC%elemLvl(iLevel)%elem%val( iElem )

      iLink = 0
      do iDir = 1, layout%fStencil%QQN
        if( globBC%elemLvl(iLevel)%bitmask%val( iDir, iElem )) then
          iLink = iLink + 1
          missing_links(iLink) = iDir
        else
          fTmp( iDir ) = &
            & state( ?FETCH?( iDir, iField, elemPos, QQ, nScalars, nSize,neigh ) )
        end if
      end do
      fTmp( layout%fStencil%restPosition ) = &
& state( ?FETCH?( layout%fStencil%restPosition, iField, elemPos, QQ, nScalars, nSize,neigh ) )

      !local velocity
      posInBuffer = globBC%elemLvl( iLevel )%posInBcElemBuf%val( iElem )
      vel(1) = sum ( bcBuffer( (posInBuffer-1)*nScalars+varPos(1) :  &
        &                      (posInBuffer-1)*nScalars+varPos(QQ) ) &
        &            * layout%fStencil%cxDir(1,:) )
      vel(2) = sum ( bcBuffer( (posInBuffer-1)*nScalars+varPos(1) :  &
        &                      (posInBuffer-1)*nScalars+varPos(QQ) ) &
        &            * layout%fStencil%cxDir(2,:) )
      vel(3) = sum ( bcBuffer( (posInBuffer-1)*nScalars+varPos(1) :  &
        &                      (posInBuffer-1)*nScalars+varPos(QQ) ) &
        &            * layout%fStencil%cxDir(3,:) )
      vel = vel/rho0

      ! density from defined pressure
      rho = cs2inv*pressB(iElem) + rho0

      ! moments for in-compressible model
      !first_moments = (/ 0.0_rk, 0.0_rk, 0.0_rk /) !momX, momY, momZ
      first_moments = rho0*vel !momX, momY, momZ
      second_moments = (/ pressB(iElem) + rho0*vel(1)*vel(1), & !momXX
        &                 pressB(iElem) + rho0*vel(2)*vel(2), & !momYY
        &                 pressB(iElem) + rho0*vel(3)*vel(3), & !momZZ
        &                 rho0*vel(1)*vel(2),                 & !momXY
        &                 rho0*vel(2)*vel(3),                 & !momYZ
        &                 rho0*vel(3)*vel(1) /)                 !momXZ
      third_moments = (/ cs2*rho0*vel(2), & !momXXY
        &                cs2*rho0*vel(3), & !momXXZ
        &                cs2*rho0*vel(1), & !momYYX
        &                cs2*rho0*vel(3), & !momYYZ
        &                cs2*rho0*vel(1), & !momZZX
        &                cs2*rho0*vel(2) /) !momZZY
      !fourth_moments = cs2*pressB(iElem) !momXXYY, momYYZZ, momZZXX
      fourth_moments = (/ &
        & cs2*(pressB(iElem)+rho0*( vel(1)**2 + vel(2)**2 - 0.5_rk*vel(3)**2 )), &!momXXYY
        & cs2*(pressB(iElem)+rho0*(-0.5_rk*vel(1)**2 + vel(2)**2 + vel(3)**2 )), &!momYYZZ
        & cs2*(pressB(iElem)+rho0*( vel(1)**2 - 0.5_rk*vel(2)**2 + vel(3)**2 ))  &!momZZXX
        & /)

      ! moments for pressure-wall corner
      ! KM: for pressure-velocity, velocity at boundary node must be defined
      moments = 0.0_rk
      moments(1) = rho
      moments(2:nMom(1)) = &
        & first_moments(layout%moment%first_moments(:))
      moments(nMom(1)+1:nMom(2)) = &
        & second_moments(layout%moment%second_moments(:))
      moments(nMom(2)+1:nMom(3)) = &
        & third_moments(layout%moment%third_moments(:))
      moments(nMom(3)+1:nMom(4)) = &
        & fourth_moments(layout%moment%fourth_moments(:))

      ! compute rhs of lse to compute unknown pdfs
      do iLink = 1, nLinks
        rhs(iLink) = &
          & moments( me%elemLvl(iLevel)%moments(iElem)%knownMom_pos(iLink) ) &
          ! convert pdf to moments
          & - dot_product(layout%moment%toMoments%A(                        &
          &    me%elemLvl(iLevel)%moments(iElem)%knownMom_pos(iLink), :), fTmp)
      end do

      ! unknown pdfs
      unKnown_fTmp = matmul( &
        & me%elemLvl(iLevel)%moments(iElem)%unKnownPdfs_MatInv, rhs )

      do iLink = 1, nLinks
        state(                                                                 &
          & ?FETCH?( missing_links(iLink), iField, elemPos, QQ, nScalars, nSize,neigh ) &
          & )  = unKnown_fTmp(iLink)
      end do

     deallocate(missing_links)
     deallocate(unKnown_fTmp)
     deallocate(rhs)
   end do

  end subroutine pressure_momentsbased_incomp
! ****************************************************************************** !


! ****************************************************************************** !
  !> author: Kannan Masilamani
  !! Moment based wall boundary condition from Sam Bennent PhD thesis
  !! "A Lattice Boltzmann Model for Diffusion of Binary Gas Mixtures"
  !! Usage
  !! -----
  !!```lua
  !!boundary_condition = {
  !! { label = 'inlet',
  !!    kind = 'moments_wall'}
  !!```
  !!
  !! 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 moments_wall( 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) :: rhoB, press
    integer :: iELem, iDir, tmpDir
    integer :: nLinks, iLink, QQ, elemPos
    integer, allocatable :: missing_links(:)
    real(kind=rk), allocatable :: rhs(:)
    real(kind=rk), allocatable :: unKnown_fTmp( : )
    real(kind=rk) :: fTmp( layout%fStencil%QQ )
    real(kind=rk) :: moments( layout%fStencil%QQ )
    real(kind=rk) :: first_moments(3), second_moments(6), third_moments(6), &
      &              fourth_moments(3)
    integer :: nMom(4)
    ! ---------------------------------------------------------------------------
    QQ = layout%fStencil%QQ

    nMom(1) = 1 + size(layout%moment%first_moments)
    nMom(2) = nMom(1) + size(layout%moment%second_moments)
    nMom(3) = nMom(2) + size(layout%moment%third_moments)
    nMom(4) = nMom(3) + size(layout%moment%fourth_moments)

    do iElem = 1, globBC%nElems( iLevel )
      fTmp = 0.0_rk
      nLinks = me%elemLvl(iLevel)%moments(iElem)%nUnKnownPdfs
      allocate(missing_links(nLinks))
      allocate(rhs(nLinks))
      allocate(unKnown_fTmp(nLinks))
      elemPos = globBC%elemLvl(iLevel)%elem%val( iElem )

      iLink = 0
      do iDir = 1, layout%fStencil%QQN
        if( globBC%elemLvl(iLevel)%bitmask%val( iDir, iElem )) then
          iLink = iLink + 1
          missing_links(iLink) = iDir
        else
          fTmp( iDir ) = &
            &state(?FETCH?( iDir, iField, elemPos, QQ, nScalars, nSize,neigh ) )
        end if
      end do
      fTmp( layout%fStencil%restPosition ) = &
& state( ?FETCH?( layout%fStencil%restPosition, iField, elemPos, QQ, nScalars, nSize,neigh ) )

      rhoB = fTmp(layout%fStencil%restPosition)
      do iDir = 1, layout%fStencil%QQN
        if( globBC%elemLvl(iLevel)%bitmask%val( iDir, iElem )) then
          tmpDir = layout%fStencil%cxDirInv(iDir)
          if( globBC%elemLvl(iLevel)%bitmask%val( tmpDir, iElem )) then
            tmpDir = layout%fStencil%cxDirInv(globBC%elemLvl(iLevel)%          &
                &                                         normalInd%val(iElem))
          end if
        else
          tmpDir = iDir
        end if
        rhoB  = rhoB + fTmp(tmpDir)
      end do
      press = cs2 * rhoB

      ! moments for in-compressible model
      first_moments = (/ 0.0_rk, 0.0_rk, 0.0_rk /) !momX, momY, momZ
      second_moments = (/ press, & !momXX
        &                 press, & !momYY
        &                 press, & !momZZ
        &                 0.0_rk,& !momXY
        &                 0.0_rk,& !momYZ
        &                 0.0_rk /)!momXZ
      third_moments = 0.0_rk
      fourth_moments = cs2*press !momXXYY, momYYZZ, momZZXX

      ! moments for pressure-wall corner
      moments = 0.0_rk
      moments(1) = rhoB
      moments(2:nMom(1)) = &
        & first_moments(layout%moment%first_moments(:))
      moments(nMom(1)+1:nMom(2)) = &
        & second_moments(layout%moment%second_moments(:))
      moments(nMom(2)+1:nMom(3)) = &
        & third_moments(layout%moment%third_moments(:))
      moments(nMom(3)+1:nMom(4)) = &
        & fourth_moments(layout%moment%fourth_moments(:))

!     write(*,*) 'Moments ', moments

      do iLink = 1, nLinks
        rhs(iLink) = &
          & moments( me%elemLvl(iLevel)%moments(iElem)%knownMom_pos(iLink) ) &
          ! convert pdf to moments
          & - dot_product(layout%moment%toMoments%A(                        &
          &    me%elemLvl(iLevel)%moments(iElem)%knownMom_pos(iLink), :), fTmp)
      end do
     !write(*,*) 'rhs ', rhs

     unKnown_fTmp = matmul( &
       & me%elemLvl(iLevel)%moments(iElem)%unKnownPdfs_MatInv, rhs )

     !write(*,*) 'unknown fTmp ', unKnown_fTmp

      do iLink = 1, nLinks
        state(                                                         &
&?FETCH?( missing_links(iLink), iField, elemPos, QQ, nScalars, nSize,neigh ) &
         & )  = unKnown_fTmp(iLink)
      end do

      deallocate(missing_links)
      deallocate(rhs)
      deallocate(unKnown_fTmp)
    end do

  end subroutine moments_wall
! ****************************************************************************** !


! ****************************************************************************** !
  !> author: Kannan Masilamani
  !! Moment based velocity boundary condition for weakly compressible LBM model
  !! Based on Sam Bennent PhD thesis
  !! "A Lattice Boltzmann Model for Diffusion of Binary Gas Mixtures"
  !!  Usage
  !! -----
  !!```lua
  !!boundary_condition = {
  !! { label = 'inlet',
  !!    kind = 'velocity_momentsbased'
  !!    velocityX = 0.06,
  !!    velocityY = 0.0,
  !!    velocityZ = 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_momentsbased( 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) :: rho, press, dens_correction, vel(3)
    real(kind=rk) :: vel_b(globBC%nElems(iLevel)*3), inv_vel
    integer :: iELem, iDir, tmpDir
    integer :: nLinks, iLink, QQ, elemPos
    integer, allocatable :: missing_links(:)
    real(kind=rk), allocatable :: rhs(:)
    real(kind=rk), allocatable :: unKnown_fTmp( : )
    real(kind=rk) :: fTmp( layout%fStencil%QQ )
    real(kind=rk) :: moments( layout%fStencil%QQ )
    real(kind=rk) :: first_moments(3), second_moments(6), third_moments(6), &
      &              fourth_moments(3)
    integer :: nMom(4)
    integer :: bcVel_pos
    ! ---------------------------------------------------------------------------
    QQ = layout%fStencil%QQ
    inv_vel = 1.0_rk / physics%fac( iLevel )%vel
    !write(*,*) 'Boundary label ', trim(me%label)

    ! 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                               )

    ! If physical quantities are given, transform to lattice units by division
    ! with the conversion factor
    !write(*,*) 'phy Vel', ux, uy, uz
    vel_b = vel_b * inv_vel

    nMom(1) = 1 + size(layout%moment%first_moments)
    nMom(2) = nMom(1) + size(layout%moment%second_moments)
    nMom(3) = nMom(2) + size(layout%moment%third_moments)
    nMom(4) = nMom(3) + size(layout%moment%fourth_moments)

    do iElem = 1, globBC%nElems( iLevel )

      vel(:) = vel_b( (iElem-1)*3+1 : iElem*3 )

      fTmp = 0.0_rk
      dens_correction = 0.0_rk
      nLinks = me%elemLvl(iLevel)%moments(iElem)%nUnKnownPdfs
      allocate(missing_links(nLinks))
      allocate(rhs(nLinks))
      allocate(unKnown_fTmp(nLinks))
      elemPos = globBC%elemLvl(iLevel)%elem%val( iElem )

      iLink = 0
      do iDir = 1,layout%fStencil%QQN
        if( globBC%elemLvl(iLevel)%bitmask%val( iDir, iElem )) then
          iLink = iLink +1
          missing_links(iLink) = iDir
        else
          fTmp( iDir ) = &
            &state(?FETCH?( iDir, iField, elemPos, QQ, nScalars, nSize, neigh ) )
        end if
      end do
      fTmp( layout%fStencil%restPosition ) = &
& state( ?FETCH?( layout%fStencil%restPosition, iField, elemPos, QQ, nScalars, nSize,neigh ) )

      rho = fTmp(layout%fStencil%restPosition)
      do iDir = 1,layout%fStencil%QQN
        if( globBC%elemLvl(iLevel)%bitmask%val( iDir, iElem )) then
          tmpDir = layout%fStencil%cxDirInv(iDir)
          if( globBC%elemLvl(iLevel)%bitmask%val( tmpDir, iElem )) then
            tmpDir = layout%fStencil%cxDirInv(globBC%elemLvl(iLevel)%          &
                &                                         normalInd%val(iElem))
          end if
        else
          tmpDir = iDir
        end if
        rho  = rho + fTmp(tmpDir)
      end do

      ! add small correction term to density which comes from momentum in the
      ! flow direction while substituing unknowns pdfs terms obtained
      ! from moments BC derivation to compute density.
      ! With this correction term density is recovered correctly at this node
      dens_correction = dot_product(vel, globBC%elemLvl(iLevel)%normal%val(:, iElem))
      ! for compressible model, rho_new = rho_old + rho_new*dens_correction
      ! thus, rho_new = rho_old/(1-dens_correction)
      rho = rho/(1.0_rk-dens_correction)

      press = cs2 * rho

      ! moments for in-compressible model
      first_moments = vel !momX, momY, momZ
      second_moments = (/ press + rho*vel(1)*vel(1), & !momXX
        &                 press + rho*vel(2)*vel(2), & !momYY
        &                 press + rho*vel(3)*vel(3), & !momZZ
        &                 rho*vel(1)*vel(2),                 & !momXY
        &                 rho*vel(2)*vel(3),                 & !momYZ
        &                 rho*vel(3)*vel(1) /)                 !momXZ
      third_moments = (/ cs2*rho*vel(2), & !momXXY
        &                cs2*rho*vel(3), & !momXXZ
        &                cs2*rho*vel(1), & !momYYX
        &                cs2*rho*vel(3), & !momYYZ
        &                cs2*rho*vel(1), & !momZZX
        &                cs2*rho*vel(2) /) !momZZY
      !fourth_moments = cs2*pressB(iElem) !momXXYY, momYYZZ, momZZXX
      fourth_moments = (/ &
        & cs2*(press+rho*( vel(1)**2 + vel(2)**2 - 0.5_rk*vel(3)**2 )), &!momXXYY
        & cs2*(press+rho*(-0.5_rk*vel(1)**2 + vel(2)**2 + vel(3)**2 )), &!momYYZZ
        & cs2*(press+rho*( vel(1)**2 - 0.5_rk*vel(2)**2 + vel(3)**2 ))  &!momZZXX
        & /)

      ! moments for pressure-wall corner
      ! KM: for pressure-velocity, velocity at boundary node must be defined
      moments = 0.0_rk
      moments(1) = rho
      moments(2:nMom(1)) = &
        & first_moments(layout%moment%first_moments(:))
      moments(nMom(1)+1:nMom(2)) = &
        & second_moments(layout%moment%second_moments(:))
      moments(nMom(2)+1:nMom(3)) = &
        & third_moments(layout%moment%third_moments(:))
      moments(nMom(3)+1:nMom(4)) = &
        & fourth_moments(layout%moment%fourth_moments(:))

      do iLink = 1, nLinks
        rhs(iLink) = &
          & moments( me%elemLvl(iLevel)%moments(iElem)%knownMom_pos(iLink) ) &
          ! convert pdf to moments
          & - dot_product(layout%moment%toMoments%A(                        &
          &    me%elemLvl(iLevel)%moments(iElem)%knownMom_pos(iLink), :), fTmp)
      end do

      unKnown_fTmp = matmul( &
        & me%elemLvl(iLevel)%moments(iElem)%unKnownPdfs_MatInv, rhs )

      do iLink = 1, nLinks
        state(                                                                 &
&?FETCH?( missing_links(iLink), iField, elemPos, QQ, nScalars, nSize,neigh ) &
         & )  = unKnown_fTmp(iLink)
      end do

      deallocate(missing_links)
      deallocate(rhs)
      deallocate(unKnown_fTmp)

    end do

  end subroutine velocity_momentsbased
! ****************************************************************************** !


! ****************************************************************************** !
  !> author: Kannan Masilamani
  !! Moment based velocity boundary condition for incompressible LBM model
  !! Based on Sam Bennent PhD thesis
  !! "A Lattice Boltzmann Model for Diffusion of Binary Gas Mixtures"
  !!  Usage
  !! -----
  !!```lua
  !!boundary_condition = {
  !! { label = 'inlet',
  !!    kind = 'velocity_momentsbased'
  !!    velocityX = 0.06,
  !!    velocityY = 0.0,
  !!    velocityZ = 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_momentsbased_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) :: rho, press, dens_correction, vel(3)
    real(kind=rk) :: vel_b(globBC%nElems(iLevel)*3), inv_vel
    integer :: iELem, iDir, tmpDir
    integer :: nLinks, iLink, QQ, elemPos
    integer, allocatable :: missing_links(:)
    real(kind=rk), allocatable :: rhs(:)
    real(kind=rk), allocatable :: unKnown_fTmp( : )
    real(kind=rk) :: fTmp( layout%fStencil%QQ )
    real(kind=rk) :: moments( layout%fStencil%QQ )
    real(kind=rk) :: first_moments(3), second_moments(6), third_moments(6), &
      &              fourth_moments(3)
    integer :: nMom(4)
    integer :: bcVel_pos
    ! ---------------------------------------------------------------------------
    QQ = layout%fStencil%QQ
    inv_vel = 1.0_rk / physics%fac( iLevel )%vel
    !write(*,*) 'Boundary label ', trim(me%label)

    ! 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                               )

    ! If physical quantities are given, transform to lattice units by division
    ! with the conversion factor
    vel_b = vel_b * inv_vel

    nMom(1) = 1 + size(layout%moment%first_moments)
    nMom(2) = nMom(1) + size(layout%moment%second_moments)
    nMom(3) = nMom(2) + size(layout%moment%third_moments)
    nMom(4) = nMom(3) + size(layout%moment%fourth_moments)

    do iElem = 1, globBC%nElems( iLevel )

      vel(:) = vel_b( (iElem-1)*3+1 : iElem*3 )

      fTmp = 0.0_rk
      dens_correction = 0.0_rk
      nLinks = me%elemLvl(iLevel)%moments(iElem)%nUnKnownPdfs
      allocate(missing_links(nLinks))
      allocate(rhs(nLinks))
      allocate(unKnown_fTmp(nLinks))
      elemPos = globBC%elemLvl(iLevel)%elem%val( iElem )

      iLink = 0
      do iDir = 1,layout%fStencil%QQN
        if( globBC%elemLvl(iLevel)%bitmask%val( iDir, iElem )) then
          iLink = iLink +1
          missing_links(iLink) = iDir
        else
          fTmp( iDir ) = &
            &  state(?FETCH?( iDir, iField, elemPos, QQ, nScalars, nSize,neigh ) )
        end if
      end do
      fTmp( layout%fStencil%restPosition ) = &
& state( ?FETCH?( layout%fStencil%restPosition, iField, elemPos, QQ, nScalars, nSize,neigh ) )

      rho = fTmp( layout%fStencil%restPosition )
      do iDir = 1,layout%fStencil%QQN
        if( globBC%elemLvl(iLevel)%bitmask%val( iDir, iElem )) then
          tmpDir = layout%fStencil%cxDirInv(iDir)
          if( globBC%elemLvl(iLevel)%bitmask%val( tmpDir, iElem )) then
            tmpDir = layout%fStencil%cxDirInv(globBC%elemLvl(iLevel)%          &
                &                                         normalInd%val(iElem))
          end if
        else
          tmpDir = iDir
        end if
        rho  = rho + fTmp(tmpDir)
      end do

      ! add small correction term to density which comes from momentum in the
      ! flow direction while substituing unknowns pdfs terms obtained
      ! from moments BC derivation to compute density.
      ! With this correction term density is recovered correctly at this node
      dens_correction = dot_product(vel, globBC%elemLvl(iLevel)%normal%val(:, iElem))
      ! for compressible model, rho_new = rho_old + rho_new*dens_correction
      ! thus, rho_new = rho_old/(1-dens_correction)
      ! correction term for in-compressible mode
      rho = rho + rho0*dens_correction

      press = cs2 * rho

      ! moments for in-compressible model
      first_moments = vel !momX, momY, momZ
      second_moments = (/ press + rho0*vel(1)*vel(1), & !momXX
        &                 press + rho0*vel(2)*vel(2), & !momYY
        &                 press + rho0*vel(3)*vel(3), & !momZZ
        &                 rho0*vel(1)*vel(2),                 & !momXY
        &                 rho0*vel(2)*vel(3),                 & !momYZ
        &                 rho0*vel(3)*vel(1) /)                 !momXZ
      third_moments = (/ cs2*rho0*vel(2), & !momXXY
        &                cs2*rho0*vel(3), & !momXXZ
        &                cs2*rho0*vel(1), & !momYYX
        &                cs2*rho0*vel(3), & !momYYZ
        &                cs2*rho0*vel(1), & !momZZX
        &                cs2*rho0*vel(2) /) !momZZY
      !fourth_moments = cs2*pressB(iElem) !momXXYY, momYYZZ, momZZXX
      fourth_moments = (/ &
        & cs2*(press+rho0*( vel(1)**2 + vel(2)**2 - 0.5_rk*vel(3)**2 )), &!momXXYY
        & cs2*(press+rho0*(-0.5_rk*vel(1)**2 + vel(2)**2 + vel(3)**2 )), &!momYYZZ
        & cs2*(press+rho0*( vel(1)**2 - 0.5_rk*vel(2)**2 + vel(3)**2 ))  &!momZZXX
        & /)

      ! moments for pressure-wall corner
      ! KM: for pressure-velocity, velocity at boundary node must be defined
      moments = 0.0_rk
      moments(1) = rho
      moments(2:nMom(1)) = &
        & first_moments(layout%moment%first_moments(:))
      moments(nMom(1)+1:nMom(2)) = &
        & second_moments(layout%moment%second_moments(:))
      moments(nMom(2)+1:nMom(3)) = &
        & third_moments(layout%moment%third_moments(:))
      moments(nMom(3)+1:nMom(4)) = &
        & fourth_moments(layout%moment%fourth_moments(:))

      do iLink = 1, nLinks
        rhs(iLink) = &
          & moments( me%elemLvl(iLevel)%moments(iElem)%knownMom_pos(iLink) ) &
          ! convert pdf to moments
          & - dot_product(layout%moment%toMoments%A(                        &
          &    me%elemLvl(iLevel)%moments(iElem)%knownMom_pos(iLink), :), fTmp)
      end do

      unKnown_fTmp = matmul( &
        & me%elemLvl(iLevel)%moments(iElem)%unKnownPdfs_MatInv, rhs )

      do iLink = 1, nLinks
        state(                                                         &
&?FETCH?( missing_links(iLink), iField, elemPos, QQ, nScalars, nSize,neigh ) &
         & )  = unKnown_fTmp(iLink)
      end do

      deallocate(missing_links)
      deallocate(rhs)
      deallocate(unKnown_fTmp)
    end do

  end subroutine velocity_momentsbased_incomp
! ****************************************************************************** !


! ****************************************************************************** !
  !> Boundary condition for APESMATE that supports pdfs as input variable
  !! Makes it possible to use pdf from, for example, the left to the right domain
  !!
  !! Usage
  !! -----
  !!```lua
  !!boundary_condition = {
  !! { label = 'some_label',
  !!   kind  = 'bc_pdf',
  !!   pdf   = 'some_value'
  !! },
  !! ...
  !!}
  !!```
  !!
  !! 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 bc_pdf( 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) :: in_pdf( layout%fStencil%QQ * globBC%nElems(iLevel) )
    !> get number of directions from stencil
    integer :: QQ
    !> position of pdf in variable system
    integer :: bcPdf_pos
    integer :: iELem, iDir, elemPos
    ! ---------------------------------------------------------------------------

    QQ = layout%fStencil%QQ

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

    do iElem = 1, globBC%nElems( iLevel )
      elemPos = globBC%elemLvl(iLevel)%elem%val( iElem )
      do iDir = 1,layout%fStencil%QQN
        if( globBC%elemLvl(iLevel)%bitmask%val( iDir, iElem )) then
          state( ?FETCH?( iDir, iField, elemPos, QQ, nScalars, nSize, neigh ) &
            & )  = in_pdf((iElem-1)*QQ+iDir)
        end if
      end do
    end do

  end subroutine bc_pdf
! ****************************************************************************** !


end module mus_bc_fluid_module
! ****************************************************************************** !
