! Copyright (c) 2020 Kannan Masilamani <kannan.masilamani@uni-siegen.de>
! Copyright (c) 2020 Peter Vitt <peter.vitt2@uni-siegen.de>
! Copyright (c) 2020 Harald Klimach <harald.klimach@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: Kannan Masilamani
!! Non-equilibrium extrapolatioon boundary condition treatment routines for
!! fluid simulation. The boundary condition for straight boundary with qVal=0
!! is based on: Guo, Z., & Shi, B. (2002). "Non-equilibrium extrapolation method
!! for velocity and pressure boundary conditions in the lattice Boltzmann
!! method." Chinese Physics.
!! The boundary condition for curved boundary is based on
!! Guo, Z., Shi, B., & Zheng, C. (2002). An extrapolation method for boundary
!! conditions in lattice Boltzmann method. Physics of Fluids, 14(May), 10–14.
!!
!! 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_nonEqExpol_module

  ! include treelm modules
  use env_module,               only: rk
  use tem_param_module,         only: div1_3, div4_3, csInv, cs, cs2, cs2inv, &
    &                                 rho0, rho0Inv
  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
  use mus_derivedQuantities_module2, only: getVelocity, getVelocity_incomp,    &
    &                                      getEqByDensVel,                     &
    &                                      getEquilibrium,                     &
    &                                      getEquilibriumIncomp,               &
    &                                      getEqByDensVelPerDir

  implicit none

  private

  public :: velocity_nonEqExpol, velocity_nonEqExpol_mrt
  public :: velocity_nonEqExpol_incomp, velocity_nonEqExpol_mrt_incomp
  public :: velocity_nonEqExpol_curved, velocity_nonEqExpol_curved_incomp
  public :: pressure_nonEqExpol, pressure_nonEqExpol_mrt
  public :: pressure_nonEqExpol_incomp, pressure_nonEqExpol_mrt_incomp

contains

! **************************************************************************** !
  !> Linkwise Dirichlet velocity non-equilibrium boundary condition for
  !! straight using the subroutine "mus_set_nonEqExpol".
  !! For straight wall, values are extrapolated along
  !! boundary normal instead of along the link and qVal =0.0 for straight wall.
  !!
  !! Notation: b (fictious boundary)   w  (physical bolundary or surface)
  !!           f (local element) ff (first neighbor fluid element)
  !!
  !! Usage
  !! -----
  !!```lua
  !!boundary_condition = {
  !!  { label = 'inlet',
  !!    kind = 'velocity_nonEqExpol',
  !!    velocity = 'inlet_vel,
  !!  }
  !!}
  !!variable = {
  !!  name = 'inlet_vel',
  !!  ncomponents = 3,
  !!  vartype = 'st_fun',
  !!  st_fun = {0.06, 0.0, 0.0}
  !!}
  !!```
  !! This is described in the paper:
  !! Non-equilibrium extrapolation method
  !! for velocity and pressure boundary conditions in the lattice Boltzmann
  !! method." Chinese Physics.
  !!
  !! More informations concerning the nonEqExpol can be found in the
  !! corresponding subroutine.
  !!
  !! 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_nonEqExpol( 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
    ! -------------------------------------------------------------------- !
    ! also determined in mus_bc_header_module
    integer :: iDir
    ! variables for fictious boundary element
    real(kind=rk) :: feq_b, vel_b(3)
    ! variables for surface (density and link-wise velocity)
    real(kind=rk) :: vel_w(globBC%nElems(iLevel)*3)
    ! variables for overnext fluid element
    real(kind=rk) :: feq_ff, rho_ff, rhoInv, vel_ff(3)
    real(kind=rk) :: inv_vel
    ! temporary local pdf values and the ones of overnext fluid
    real(kind=rk) :: pdf_ff(layout%fStencil%QQ)
    real(kind=rk) :: nEqTerm
    integer :: bcVel_pos, iLink, QQ, offset
    ! relaxation parameter
    real(kind=rk) :: omegaKine
    ! --------------------------------------------------------------------------
    QQ = layout%fstencil%QQ
    inv_vel = 1._rk / physics%fac(iLevel)%vel

    ! position of boundary velocity in varSys
    bcVel_pos = me%bc_states%velocity%varPos
    ! get velocity_phy on boundary
    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_w                               )

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

    associate( neighBufferPost => me%neigh(iLevel)%neighBufferPost,        &
      &        posInNeighBuf => me%nonEqExpol(iLevel)%posInNeighBuf,       &
      &        posInBCelems => me%nonEqExpol(iLevel)%posInBCelems,         &
      &        posInState => globBC%elemLvl(iLevel)%elem%val,              &
      &        omegaKineVal => fieldProp%fluid%viscKine%omLvl(iLevel)%val  )

      ! loop over links
      do iLink = 1, me%links(iLevel)%nVals

        ! set link-wise direction
        iDir = me%nonEqExpol(iLevel)%iDir(iLink)

        ! determine needed quantities of the overnext fluid neighbor element x_ff
        pdf_ff = neighBufferPost(1, (posInNeighBuf(iLink)-1)*QQ+1  &
          &                       : (posInNeighBuf(iLink)-1)*QQ+QQ )
        ! calculate density
        rho_ff = sum(pdf_ff)
        rhoInv = 1.0_rk / rho_ff
        ! calulate velocity
?? copy :: velocity(vel_ff, pdf_ff, rhoInv, layout%fStencil)

      ! derive feq from velocity and density per direction
?? copy :: EqDir_comp(fEq_ff, rho_ff, vel_ff, layout, iDir, cs2inv, rk)

        ! need offset for velocity
        offset = (posInBCelems(iLink)-1)*3

        ! determine needed quantities of fictious boundary (_b)
        ! compute velocity on boundary (eq.18)
        ! qVal = 0.0, boundary node is one lattice away from fluid node
        vel_b = vel_w(offset+1:offset+3)
        ! extrapolate density on boundary
        ! derive feq from velocity and density per direction
?? copy :: EqDir_comp(fEq_b, rho_ff, vel_b, layout, iDir, cs2inv, rk)

        ! write into state: feq_b + fneq_ff
        ! fneq is computed from post-collision pdf so we need to divide by
        ! (1-omega)
        omegaKine = omegaKineVal( posInState( posInBCelems(iLink) ) )
        nEqTerm = (pdf_ff(iDir) - feq_ff) / (1.0_rk-omegaKine)

        state( me%links(iLevel)%val(iLink) ) = feq_b + neqTerm

      end do !iLink

    end associate

  end subroutine velocity_nonEqExpol
! **************************************************************************** !

! **************************************************************************** !
  !> Linkwise Dirichlet velocity non-equilibrium boundary condition for
  !! straight wall for mrt collision using the subroutine "mus_set_nonEqExpol".
  !!
  !! 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_nonEqExpol_mrt( 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
    ! -------------------------------------------------------------------- !
    ! also determined in mus_bc_header_module
    integer :: iDir, iDirEq
    ! variables for fictious boundary element
    real(kind=rk) :: feq_b, vel_b(3)
    ! variables for surface (density and link-wise velocity)
    real(kind=rk) :: vel_w(globBC%nElems(iLevel)*3)
    ! variables for overnext fluid element
    real(kind=rk) :: feq_ff(layout%fStencil%QQ), rho_ff, rhoInv, vel_ff(3)
    real(kind=rk) :: inv_vel
    ! temporary local pdf values and the ones of overnext fluid
    real(kind=rk) :: pdf_ff(layout%fStencil%QQ)
    real(kind=rk) :: fneq(layout%fStencil%QQ)
    real(kind=rk) :: mneq(layout%fStencil%QQ)
    real(kind=rk) :: mInvXomega(layout%fStencil%QQ)
    real(kind=rk) :: nonEqScalingFacs(layout%fStencil%QQ)
    real(kind=rk) :: relaxation_mrt(layout%fStencil%QQ)
    real(kind=rk) :: nEqTerm
    integer :: bcVel_pos, iLink, QQ, offset
    ! relaxation parameter
    real(kind=rk) :: omegaKine
    ! --------------------------------------------------------------------------
    QQ = layout%fstencil%QQ
    inv_vel = 1._rk / physics%fac(iLevel)%vel

    ! position of boundary velocity in varSys
    bcVel_pos = me%bc_states%velocity%varPos
    ! get velocity_phy on boundary
    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_w                               )

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

    associate( neighBufferPost => me%neigh(iLevel)%neighBufferPost,        &
      &        posInNeighBuf => me%nonEqExpol(iLevel)%posInNeighBuf,       &
      &        posInBCelems => me%nonEqExpol(iLevel)%posInBCelems,         &
      &        posInState => globBC%elemLvl(iLevel)%elem%val,              &
      &        omegaKineVal => fieldProp%fluid%viscKine%omLvl(iLevel)%val, &
      &        omegaBulk  => fieldProp%fluid%omegaBulkLvl(iLevel),         &
      &        toMoment => layout%moment%toMoments%A,                      &
      &        toPDF => layout%moment%toPDF%A                              )

      ! loop over links
      do iLink = 1, me%links(iLevel)%nVals

        ! set link-wise direction
        iDir = me%nonEqExpol(iLevel)%iDir(iLink)

        ! determine needed quantities of the overnext fluid neighbor element x_ff
        pdf_ff = neighBufferPost(1, (posInNeighBuf(iLink)-1)*QQ+1  &
          &                       : (posInNeighBuf(iLink)-1)*QQ+QQ )
        ! calculate density
        rho_ff = sum(pdf_ff)
        rhoInv = 1.0_rk / rho_ff
        ! calulate velocity
?? copy :: velocity(vel_ff, pdf_ff, rhoInv, layout%fStencil)

        ! derive feq from velocity and density for all direction
        do iDirEq = 1, QQ
?? copy :: EqDir_comp(fEq_ff(iDirEq), rho_ff, vel_ff, layout, iDirEq, cs2inv, rk)
        end do

        ! need offset for velocity
        offset = (posInBCelems(iLink)-1)*3
        ! determine needed quantities of fictious boundary (_b)
        ! compute velocity on boundary (eq.18)
        ! qVal = 0.0, boundary node is one lattice away from fluid node
        vel_b = vel_w(offset+1:offset+3)
        ! extrapolate density on boundary
        ! derive feq from velocity and density per direction
?? copy :: EqDir_comp(fEq_b, rho_ff, vel_b, layout, iDir, cs2inv, rk)

        ! write into state: feq_b + fneq_ff/(1-omega)
        ! fneq is computed from post-collision pdf so we need to divide by
        ! (1-omega)
        omegaKine = omegaKineVal( posInState( posInBCelems(iLink) ) )
        relaxation_mrt = fieldProp%fluid%mrtPtr( omegaKine = omegaKine, &
          &                                      omegaBulk = omegaBulk, &
          &                                      QQ        = QQ         )
        nonEqScalingFacs(:) = 1.0_rk / ( 1.0_rk - relaxation_mrt(:) )

        fneq = pdf_ff - feq_ff
        ! convert fneq to moments neq
        mneq = matmul( toMoment(1:QQ, 1:QQ), fneq(1:QQ) )

        ! compute Minv times nonEqScalingFacs
        mInvXOmega(:) = toPDF(iDir, :) * nonEqScalingFacs(:)

        nEqTerm = dot_product( mInvXomega(1:QQ), mneq(1:QQ) )

        state( me%links(iLevel)%val(iLink) ) = feq_b + nEqTerm

      end do !iLink

    end associate

  end subroutine velocity_nonEqExpol_mrt
! **************************************************************************** !

! **************************************************************************** !
  !> Linkwise Dirichlet velocity non-equilibrium boundary condition for
  !! straight using the subroutine "mus_set_nonEqExpol".
  !! For straight wall, values are extrapolated along
  !! boundary normal instead of along the link
  !!
  !! Notation: b (fictious boundary)   w  (physical bolundary or surface)
  !!           f (local element)      ff  (overnext fluidneighbor)
  !! Usage
  !! -----
  !!```lua
  !!boundary_condition = {
  !!  { label = 'inlet',
  !!    kind = 'velocity_nonEqExpol',
  !!    velocity = 'inlet_vel,
  !!    curved = true
  !!  }
  !!}
  !!variable = {
  !!  name = 'inlet_vel',
  !!  ncomponents = 3,
  !!  vartype = 'st_fun',
  !!  st_fun = {0.06, 0.0, 0.0}
  !!}
  !!```
  !! This is described in the paper:
  !! Guo, Z.; Zheng, C. & Shi B. (2002). An extrapolation method for boundary
  !! conditions in lattice Boltzmann method.
  !! Physics of Fluids 14, 2007 (2002); https://doi.org/10.1063/1.1471914
  !!
  !!
  !! 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_nonEqExpol_curved( 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
    ! -------------------------------------------------------------------- !
    ! coefficients which are calculated in mus_bc_header_module
    real(kind=rk) :: c_w, c_f, c_ff, c_neq_f, c_neq_ff
    ! also determined in mus_bc_header_module
    integer :: iDir, posInBuffer, posInNeighBuf
    ! variables for fictious boundary element
    real(kind=rk) :: feq_b, fneq_b, rho_b, vel_b(3)
    ! variables for surface (density and link-wise velocity)
    real(kind=rk) :: vel_w(me%links(iLevel)%nVals*3)
    ! variables for fluid element
    real(kind=rk) :: feq_f, fneq_f, rho_f, vel_f(3)
    ! variables for overnext fluid element
    real(kind=rk) :: feq_ff, fneq_ff, rho_ff, vel_ff(3)
    real(kind=rk) :: inv_vel
    ! temporary local pdf values and the ones of overnext fluid
    real(kind=rk) :: pdfTmp(layout%fStencil%QQ)
    real(kind=rk) :: pdf_ff(layout%fStencil%QQ)
    integer :: bcVel_pos, iLink, QQ, offset
    ! relaxation parameter
    !real(kind=rk) :: omega
    ! --------------------------------------------------------------------------
    QQ = layout%fstencil%QQ
    inv_vel = 1._rk / physics%fac(iLevel)%vel

    ! position of boundary velocity in varSys
    bcVel_pos = me%bc_states%velocity%varPos
    ! get velocity_phy on boundary (surface)
    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_w                               )

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

    ! loop over links
    do iLink = 1, me%links(iLevel)%nVals

      ! load pre-calculated coefficients from subroutine
      c_w      = me%nonEqExpol(iLevel)%     c_w(iLink)
      c_f      = me%nonEqExpol(iLevel)%     c_f(iLink)
      c_ff     = me%nonEqExpol(iLevel)%    c_ff(iLink)
      c_neq_f  = me%nonEqExpol(iLevel)% c_neq_f(iLink)
      c_neq_ff = me%nonEqExpol(iLevel)%c_neq_ff(iLink)

      ! set link-wise direction
      iDir = me%nonEqExpol(iLevel)%iDir(iLink)

      ! determine needed quantities of the current element (_f)
      ! calculate density
      posInBuffer  = me%nonEqExpol(iLevel)%posInBuffer(iLink)
      pdfTmp(1:QQ) = bcBuffer( (posInBuffer-1)*nScalars+varPos(1) :  &
        &                    ( posInBuffer-1)*nScalars+varPos(QQ) )
      rho_f = sum(pdfTmp)
      ! calulate velocity (facilitate it by using a pure function)
      vel_f = getVelocity( subset  = pdfTmp,          &
                           stencil = layout%fStencil, &
                           varPos  = varPos           )
      ! derive feq from velocity and density per direction (using pure funct.)
      feq_f = getEqByDensVelPerDir(dens   = rho_f,                  &
        &                          vel    = vel_f,                  &
        &                          cx     = layout%fStencil         &
        &                                         %cxDirRK(:,iDir), &
        &                          weight = layout%weight(iDir)     )

      ! determine needed quantities of the overnext fluid neighbor element (_ff)
      ! calculate density
      posInNeighBuf = me%nonEqExpol(iLevel)%posInNeighBuf(iLink)
      pdf_ff = me%neigh(iLevel)%neighBufferPost(1,                    &
        &             (posInNeighBuf-1)*QQ+1: (posInNeighBuf-1)*QQ+QQ )
      rho_ff = sum(pdf_ff)
      ! calulate velocity (facilitate it by using a pure function)
      vel_ff = getVelocity( subset  = pdf_ff,          &
        &                   stencil = layout%fStencil, &
                            varPos  = varPos           )
      ! derive feq from velocity and density per direction (using pure funct.)
      feq_ff = getEqByDensVelPerDir(dens   = rho_ff,                 &
        &                           vel    = vel_ff,                 &
        &                           cx     = layout%fStencil         &
        &                                          %cxDirRK(:,iDir), &
        &                           weight = layout%weight(iDir)     )

      ! need offset for velocity
      offset = (iLink-1)*3

      ! determine needed quantities of fictious boundary (_b)
      ! compute velocity on boundary (eq.18)
      vel_b = c_w*vel_w(offset+1:offset+3) + c_f*vel_f + c_ff*vel_ff
      ! extrapolate density on boundary
      !rho_b = div4_3*rho_f - div1_3*rho_ff
      rho_b = rho_f
      ! derive feq from velocity and density per direction (using pure funct.)
      feq_b = getEqByDensVelPerDir(dens   = rho_b,                  &
        &                          vel    = vel_b,                  &
        &                          cx     = layout%fStencil         &
        &                                         %cxDirRK(:,iDir), &
        &                          weight = layout%weight(iDir)     )

      ! use pdf_b = fEq_b + fneq_b (eq.16) to determine
      ! non-equilibrium components (eq.19)
      fneq_f  = pdfTmp(iDir) - feq_f
      fneq_ff = pdf_ff(iDir) - feq_ff
      fneq_b  = c_nEq_f * fneq_f + c_nEq_ff * fneq_ff

      ! write into state
      ! fneq is computed from post-collision pdf so no need to multiply
      ! (1-omega) as in the literature
      state( me%links(iLevel)%val(iLink) ) = feq_b + fneq_b

    end do !iLink

  end subroutine velocity_nonEqExpol_curved
! **************************************************************************** !

! **************************************************************************** !
  !> Linkwise Dirichlet velocity non-equilibrium boundary condition for
  !! straight using the subroutine "mus_set_nonEqExpol" for incompressible model
  !! For straight wall, values are extrapolated along
  !! boundary normal instead of along the link and qVal=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_nonEqExpol_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
    ! -------------------------------------------------------------------- !
    ! also determined in mus_bc_header_module
    integer :: iDir
    ! variables for fictious boundary element
    real(kind=rk) :: feq_b, vel_b(3)
    ! variables for surface (density and link-wise velocity)
    real(kind=rk) :: vel_w(globBC%nElems(iLevel)*3)
    ! variables for overnext fluid element
    real(kind=rk) :: feq_ff, rho_ff, vel_ff(3)
    real(kind=rk) :: inv_vel
    ! temporary local pdf values
    real(kind=rk) :: pdf_ff(layout%fStencil%QQ)
    real(kind=rk) :: nEqTerm
    integer :: bcVel_pos, iLink, QQ, offset
    ! relaxation parameter
    real(kind=rk) :: omegaKine
    ! --------------------------------------------------------------------------
    QQ = layout%fstencil%QQ
    inv_vel = 1._rk / physics%fac(iLevel)%vel

    ! position of boundary velocity in varSys
    bcVel_pos = me%bc_states%velocity%varPos
    ! get velocity_phy on boundary
    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_w                               )

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

    associate( neighBufferPost => me%neigh(iLevel)%neighBufferPost,        &
      &        posInNeighBuf => me%nonEqExpol(iLevel)%posInNeighBuf,       &
      &        posInBCelems => me%nonEqExpol(iLevel)%posInBCelems,         &
      &        posInState => globBC%elemLvl(iLevel)%elem%val,              &
      &        omegaKineVal => fieldProp%fluid%viscKine%omLvl(iLevel)%val  )

      ! loop over links
      do iLink = 1, me%links(iLevel)%nVals

        ! set link-wise direction
        iDir = me%nonEqExpol(iLevel)%iDir(iLink)

        ! determine needed quantities of the overnext fluid neighbor element x_ff
        pdf_ff = neighBufferPost(1, (posInNeighBuf(iLink)-1)*QQ+1  &
          &                       : (posInNeighBuf(iLink)-1)*QQ+QQ )
        ! calculate density
        rho_ff = sum(pdf_ff)
        ! calulate velocity
?? copy :: velocity(vel_ff, pdf_ff, rho0Inv, layout%fStencil)
        ! derive feq from velocity and density for all direction
?? copy :: EqDir_incomp(fEq_ff, rho_ff, vel_ff, layout, iDir, rho0, cs2inv, rk)

        ! need offset for velocity
        offset = (posInBCelems(iLink)-1)*3
        ! determine needed quantities of fictious boundary (_b)
        ! compute velocity on boundary (eq.18)
        ! qVal = 0.0, boundary node is one lattice away from fluid node
        vel_b = vel_w(offset+1:offset+3)
        ! extrapolate density on boundary
        ! derive feq from velocity and density per direction
?? copy :: EqDir_incomp(fEq_b, rho_ff, vel_b, layout, iDir, rho0, cs2inv, rk)

        ! write into state: feq_b + fneq_f/(1-omega)
        ! fneq is computed from post-collision pdf so we need to divide by
        ! (1-omega)
        omegaKine = omegaKineVal( posInState( posInBCelems(iLink) ) )
        nEqTerm = (pdf_ff(iDir) - feq_ff) / (1.0_rk - omegaKine)

        state( me%links(iLevel)%val(iLink) ) = feq_b + nEqTerm

      end do !iLink

    end associate

  end subroutine velocity_nonEqExpol_incomp
! **************************************************************************** !

! **************************************************************************** !
  !> Linkwise Dirichlet velocity non-equilibrium boundary condition for
  !! straight for mrt collision using the subroutine "mus_set_nonEqExpol" for
  !! incompressible model.
  !! For straight wall, values are extrapolated along
  !! boundary normal instead of along the link and qVal=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_nonEqExpol_mrt_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
    ! -------------------------------------------------------------------- !
    ! also determined in mus_bc_header_module
    integer :: iDir, iDirLoop, posInNeighBuf
    ! variables for fictious boundary element
    real(kind=rk) :: feq_b, vel_b(3)
    ! variables for surface (density and link-wise velocity)
    real(kind=rk) :: vel_w(globBC%nElems(iLevel)*3)
    ! variables for overnext fluid element
    real(kind=rk) :: feq_ff(layout%fStencil%QQ), rho_ff, vel_ff(3)
    real(kind=rk) :: inv_vel
    ! temporary local pdf values
    real(kind=rk) :: pdf_ff(layout%fStencil%QQ)
    real(kind=rk) :: fneq(layout%fStencil%QQ)
    real(kind=rk) :: mneq(layout%fStencil%QQ)
    real(kind=rk) :: mInvXomega(layout%fStencil%QQ)
    real(kind=rk) :: nonEqScalingFacs(layout%fStencil%QQ)
    real(kind=rk) :: relaxation_mrt(layout%fStencil%QQ)
    real(kind=rk) :: nEqTerm
    integer :: bcVel_pos, iLink, QQ, offset
    ! relaxation parameter
    integer :: posInState
    real(kind=rk) :: omegaKine, omegaBulk
    ! --------------------------------------------------------------------------
    QQ = layout%fstencil%QQ
    omegaBulk = fieldProp%fluid%omegaBulkLvl(iLevel)
    inv_vel = 1._rk / physics%fac(iLevel)%vel

    ! position of boundary velocity in varSys
    bcVel_pos = me%bc_states%velocity%varPos
    ! get velocity_phy on boundary
    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_w                               )

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

    associate( neighBufferPost => me%neigh(iLevel)%neighBufferPost,        &
      &        posInNeighBuf => me%nonEqExpol(iLevel)%posInNeighBuf,       &
      &        posInBCelems => me%nonEqExpol(iLevel)%posInBCelems,         &
      &        posInState => globBC%elemLvl(iLevel)%elem%val,              &
      &        omegaKineVal => fieldProp%fluid%viscKine%omLvl(iLevel)%val, &
      &        omegaBulk  => fieldProp%fluid%omegaBulkLvl(iLevel),         &
      &        toMoment => layout%moment%toMoments%A,                      &
      &        toPDF => layout%moment%toPDF%A                              )

      ! loop over links
      do iLink = 1, me%links(iLevel)%nVals

        ! set link-wise direction
        iDir = me%nonEqExpol(iLevel)%iDir(iLink)

        ! determine needed quantities of the overnext fluid neighbor element x_ff
        pdf_ff = neighBufferPost(1, (posInNeighBuf(iLink)-1)*QQ+1  &
          &                       : (posInNeighBuf(iLink)-1)*QQ+QQ )
        ! calculate density
        rho_ff = sum(pdf_ff)
        ! calulate velocity
?? copy :: velocity(vel_ff, pdf_ff, rho0Inv, layout%fStencil)
        ! derive feq from velocity and density for all direction
        do iDirLoop = 1, layout%fStencil%QQ
?? copy :: EqDir_incomp(fEq_ff(iDirLoop), rho_ff, vel_ff, layout, iDirLoop, rho0, cs2inv, rk)
        end do

        ! need offset for velocity
        offset = (posInBCelems(iLink)-1)*3
        ! determine needed quantities of fictious boundary (_b)
        ! compute velocity on boundary (eq.18)
        ! qVal = 0.0, boundary node is one lattice away from fluid node
        vel_b = vel_w(offset+1:offset+3)
        ! extrapolate density on boundary
        ! derive feq from velocity and density per direction
?? copy :: EqDir_incomp(fEq_b, rho_ff, vel_b, layout, iDir, rho0, cs2inv, rk)

        ! write into state: feq_b + fneq_f/(1-omega)
        ! fneq is computed from post-collision pdf so we need to divide by
        ! (1-omega)
        omegaKine = omegaKineVal( posInState( posInBCelems(iLink) ) )
        relaxation_mrt = fieldProp%fluid%mrtPtr( omegaKine = omegaKine, &
          &                                      omegaBulk = omegaBulk, &
          &                                      QQ        = QQ         )
        nonEqScalingFacs(:) = 1.0_rk / ( 1.0_rk - relaxation_mrt(:) )

        fneq = pdf_ff - feq_ff
        ! convert fneq to moments neq
        mneq = matmul( toMoment(1:QQ, 1:QQ), fneq(1:QQ) )

        ! compute Minv times nonEqScalingFacs
        mInvXOmega(:) = toPDF(iDir, :) * nonEqScalingFacs(:)

        nEqTerm = dot_product( mInvXomega(1:QQ), mneq(1:QQ) )

        state( me%links(iLevel)%val(iLink) ) = feq_b + nEqTerm

      end do !iLink

    end associate

  end subroutine velocity_nonEqExpol_mrt_incomp
! **************************************************************************** !

! **************************************************************************** !
  !> Linkwise Dirichlet velocity non-equilibrium boundary condition for
  !! curved wal using the subroutine "mus_set_nonEqExpol" for incompressible
  !! model.
  !!
  !! 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_nonEqExpol_curved_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
    ! -------------------------------------------------------------------- !
    ! coefficients which are calculated in mus_bc_header_module
    real(kind=rk) :: c_w, c_f, c_ff, c_neq_f, c_neq_ff
    ! also determined in mus_bc_header_module
    integer :: iDir, posInBuffer, posInNeighBuf
    ! variables for fictious boundary element
    real(kind=rk) :: feq_b, fneq_b, rho_b, vel_b(3)
    ! variables for surface (density and link-wise velocity)
    real(kind=rk) :: vel_w(me%links(iLevel)%nVals*3)
    ! variables for fluid element
    real(kind=rk) :: feq_f, fneq_f, rho_f, vel_f(3)
    ! variables for overnext fluid element
    real(kind=rk) :: feq_ff, fneq_ff, rho_ff, vel_ff(3)
    real(kind=rk) :: inv_vel
    ! temporary local pdf values and the ones of overnext fluid
    real(kind=rk) :: pdfTmp(layout%fStencil%QQ)
    real(kind=rk) :: pdf_ff(layout%fStencil%QQ)
    integer :: bcVel_pos, iLink, QQ, offset
    ! --------------------------------------------------------------------------
    QQ = layout%fstencil%QQ
    inv_vel = 1._rk / physics%fac(iLevel)%vel

    ! position of boundary velocity in varSys
    bcVel_pos = me%bc_states%velocity%varPos
    ! get velocity_phy on boundary (surface)
    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_w                               )

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

    ! loop over links
    do iLink = 1, me%links(iLevel)%nVals

      ! load pre-calculated coefficients from subroutine
      c_w      = me%nonEqExpol(iLevel)%     c_w(iLink)
      c_f      = me%nonEqExpol(iLevel)%     c_f(iLink)
      c_ff     = me%nonEqExpol(iLevel)%    c_ff(iLink)
      c_neq_f  = me%nonEqExpol(iLevel)% c_neq_f(iLink)
      c_neq_ff = me%nonEqExpol(iLevel)%c_neq_ff(iLink)

      ! set link-wise direction
      iDir = me%nonEqExpol(iLevel)%iDir(iLink)

      ! determine needed quantities of the current element (_f)
      ! calculate density
      posInBuffer  = me%nonEqExpol(iLevel)%posInBuffer(iLink)
      pdfTmp(1:QQ) = bcBuffer( (posInBuffer-1)*nScalars+varPos(1) :  &
        &                      (posInBuffer-1)*nScalars+varPos(QQ) )
      !posInBuffer  = me%nonEqExpol(iLevel)%posInBCelems(iLink)
      !pdfTmp(1:QQ) = me%neigh(iLevel)%computeNeighBuf_now( &
      !  &                    (posInBuffer-1)*nScalars+varPos(1) :  &
      !  &                    (posInBuffer-1)*nScalars+varPos(QQ) )
      rho_f = sum(pdfTmp)
      ! calulate velocity (facilitate it by using a pure function)
      vel_f = getVelocity_incomp( subset  = pdfTmp,          &
                                  stencil = layout%fStencil, &
                                  varPos  = varPos           )
      ! derive feq from velocity and density per direction (using pure funct.)
?? copy :: EqDir_incomp(feq_f, rho_f, vel_f, layout, iDir, rho0, cs2inv, rk)

      ! determine needed quantities of the overnext fluid neighbor element (_ff)
      ! calculate density
      posInNeighBuf = me%nonEqExpol(iLevel)%posInNeighBuf(iLink)
      pdf_ff = me%neigh(iLevel)%neighBufferPost(1,                    &
        &             (posInNeighBuf-1)*QQ+1: (posInNeighBuf-1)*QQ+QQ )
      rho_ff = sum(pdf_ff)
      ! calulate velocity (facilitate it by using a pure function)
      vel_ff = getVelocity_incomp( subset  = pdf_ff,          &
        &                          stencil = layout%fStencil, &
                                   varPos  = varPos           )
      ! derive feq from velocity and density per direction (using pure funct.)
?? copy :: EqDir_incomp(feq_ff, rho_ff, vel_ff, layout, iDir, rho0, cs2inv, rk)

      ! need offset for velocity
      offset = (iLink-1)*3

      ! determine needed quantities of fictious boundary (_b)
      ! compute velocity on boundary (eq.18)
      vel_b = c_w*vel_w(offset+1:offset+3) + c_f*vel_f + c_ff*vel_ff
      ! extrapolate density on boundary
      !rho_b = div4_3*rho_f - div1_3*rho_ff
      rho_b = rho_f
      ! derive feq from velocity and density per direction (using pure funct.)
?? copy :: EqDir_incomp(feq_b, rho_b, vel_b, layout, iDir, rho0, cs2inv, rk)

      ! use pdf_b = fEq_b + fneq_b (eq.16) to determine
      ! non-equilibrium components (eq.19)
      fneq_f  = pdfTmp(iDir) - feq_f
      fneq_ff = pdf_ff(iDir) - feq_ff

      fneq_b  = (c_nEq_f * fneq_f) + (c_nEq_ff * fneq_ff)

      ! write into state
      ! KM: Here dividing nonequilibrium by 1-omega does not work
      state( me%links(iLevel)%val(iLink) ) = feq_b + fneq_b

    end do !iLink

  end subroutine velocity_nonEqExpol_curved_incomp
! **************************************************************************** !


! **************************************************************************** !
  !> Linkwise Dirichlet pressure non-equilibrium boundary condition for
  !! straight using the subroutine "mus_set_nonEqExpol"
  !! For straight wall, values are extrapolated along
  !! boundary normal instead of along the link and qVal=0.
  !!
  !! Notation: b (fictious boundary)   w  (physical bolundary or surface)
  !!           f (local element)      ff  (overnext fluidneighbor)
  !!
  !! Usage
  !! -----
  !!```lua
  !!boundary_condition = {
  !!  { label = 'outlet',
  !!    kind = 'pressure_nonEqExpol',
  !!    pressure = 'press,
  !!    curved = true
  !!  }
  !!}
  !!variable = {
  !!  name = 'press',
  !!  ncomponents = 1,
  !!  vartype = 'st_fun',
  !!  st_fun = 0.0
  !!}
  !!```
  !! It is based on the paper:
  !! Guo, Z.; Zheng, C. & Shi B. (2002). An extrapolation method for boundary
  !! conditions in lattice Boltzmann method.
  !! Physics of Fluids 14, 2007 (2002); https://doi.org/10.1063/1.1471914
  !!
  !! More informations concerning the nonEqExpol can be found in the
  !! corresponding subroutine.
  !!
  !!
  !! 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_nonEqExpol( 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
    ! -------------------------------------------------------------------- !
    ! also determined in mus_bc_header_module
    integer :: iDir
    ! variables for fictious boundary element
    real(kind=rk) :: feq_b, rho_b, vel_b(3)
    ! variables for surface (density and link-wise velocity)
    real(kind=rk) :: rho_w(globBC%nElems(iLevel))
    ! variables for overnext fluid element
    real(kind=rk) :: feq_ff, rho_ff, rhoInv, vel_ff(3)
    ! temporary local pdf values and the ones of overnext fluid
    real(kind=rk) :: pdf_ff(layout%fStencil%QQ)
    real(kind=rk) :: nEqTerm
    integer :: bcPress_pos, iLink, QQ
    ! relaxation parameter
    ! real(kind=rk) :: omega
    real(kind=rk) :: inv_rho_phy
    real(kind=rk) :: omegaKine
    ! --------------------------------------------------------------------------
    QQ = layout%fstencil%QQ
    inv_rho_phy = 1.0_rk / physics%fac(iLevel)%press * cs2inv

    ! position of boundary velocity in varSys
    bcPress_pos = me%bc_states%pressure%varPos
    ! get pressure_phy on boundary (surface) and store it to rho_w
    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_w                                 )

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

    associate( neighBufferPost => me%neigh(iLevel)%neighBufferPost,       &
      &        posInNeighBuf => me%nonEqExpol(iLevel)%posInNeighBuf,      &
      &        posInBCelems => me%nonEqExpol(iLevel)%posInBCelems,        &
      &        posInState => globBC%elemLvl(iLevel)%elem%val,             &
      &        omegaKineVal => fieldProp%fluid%viscKine%omLvl(iLevel)%val )

      ! loop over links
      do iLink = 1, me%links(iLevel)%nVals

        ! set link-wise direction
        iDir = me%nonEqExpol(iLevel)%iDir(iLink)

        ! determine needed quantities of the overnext fluid neighbor element
        ! x_ff
        pdf_ff = neighBufferPost(1, (posInNeighBuf(iLink)-1)*QQ+1  &
          &                       : (posInNeighBuf(iLink)-1)*QQ+QQ )
        ! calculate density
        rho_ff = sum(pdf_ff)
        rhoInv = 1.0_rk / rho_ff
        ! calulate velocity
?? copy :: velocity(vel_ff, pdf_ff, rhoInv, layout%fStencil)

        ! derive feq from velocity and density per direction
?? copy :: EqDir_comp(fEq_ff, rho_ff, vel_ff, layout, iDir, cs2inv, rk)

        ! density on boundary, qVal=0: boundary and fluid node overlaps
        rho_b = rho_w(posInBCelems(iLink))
        ! extrapolate velocity on boundary
        vel_b = vel_ff
        ! derive feq from velocity and density per direction
?? copy :: EqDir_comp(fEq_b, rho_b, vel_b, layout, iDir, cs2inv, rk)

        ! write into state: feq_b + fneq_ff/(1-omega)
        ! fneq is computed from post-collision pdf so we need to divide by
        ! (1-omega)
        omegaKine = omegaKineVal( posInState( posInBCelems(iLink) ) )
        nEqTerm = (pdf_ff(iDir) - feq_ff) / (1.0_rk-omegaKine)

        state( me%links(iLevel)%val(iLink) ) = feq_b + nEqTerm

      end do !iLink

    end associate

  end subroutine pressure_nonEqExpol
! **************************************************************************** !

! **************************************************************************** !
  !> Linkwise Dirichlet pressure non-equilibrium boundary condition for
  !! straight for mrt collision using the subroutine "mus_set_nonEqExpol"
  !!
  !! 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_nonEqExpol_mrt( 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
    ! -------------------------------------------------------------------- !
    ! also determined in mus_bc_header_module
    integer :: iDir, iDirEq
    ! variables for fictious boundary element
    real(kind=rk) :: feq_b, rho_b, vel_b(3)
    ! variables for surface (density and link-wise velocity)
    real(kind=rk) :: rho_w(globBC%nElems(iLevel))
    ! variables for overnext fluid element
    real(kind=rk) :: feq_ff(layout%fStencil%QQ), rho_ff, rhoInv, vel_ff(3)
    ! temporary local pdf values and the ones of overnext fluid
    real(kind=rk) :: pdf_ff(layout%fStencil%QQ)
    real(kind=rk) :: fneq(layout%fStencil%QQ)
    real(kind=rk) :: mneq(layout%fStencil%QQ)
    real(kind=rk) :: mInvXomega(layout%fStencil%QQ)
    real(kind=rk) :: nonEqScalingFacs(layout%fStencil%QQ)
    real(kind=rk) :: relaxation_mrt(layout%fStencil%QQ)
    real(kind=rk) :: nEqTerm
    integer :: bcPress_pos, iLink, QQ
    ! relaxation parameter
    ! real(kind=rk) :: omega
    real(kind=rk) :: inv_rho_phy
    real(kind=rk) :: omegaKine
    ! --------------------------------------------------------------------------
    QQ = layout%fstencil%QQ
    inv_rho_phy = 1.0_rk / physics%fac(iLevel)%press * cs2inv

    ! position of boundary velocity in varSys
    bcPress_pos = me%bc_states%pressure%varPos
    ! get pressure_phy on boundary (surface) and store it to rho_w
    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_w                                 )

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

    associate( neighBufferPost => me%neigh(iLevel)%neighBufferPost,        &
      &        posInNeighBuf => me%nonEqExpol(iLevel)%posInNeighBuf,       &
      &        posInBCelems => me%nonEqExpol(iLevel)%posInBCelems,         &
      &        posInState => globBC%elemLvl(iLevel)%elem%val,              &
      &        omegaKineVal => fieldProp%fluid%viscKine%omLvl(iLevel)%val, &
      &        omegaBulk  => fieldProp%fluid%omegaBulkLvl(iLevel),         &
      &        toMoment => layout%moment%toMoments%A,                      &
      &        toPDF => layout%moment%toPDF%A                              )

      ! loop over links
      do iLink = 1, me%links(iLevel)%nVals

        ! set link-wise direction
        iDir = me%nonEqExpol(iLevel)%iDir(iLink)

        ! determine needed quantities of the overnext fluid neighbor element
        ! x_ff
        pdf_ff = neighBufferPost(1, (posInNeighBuf(iLink)-1)*QQ+1  &
          &                       : (posInNeighBuf(iLink)-1)*QQ+QQ )
        ! calculate density
        rho_ff = sum(pdf_ff)
        rhoInv = 1.0_rk / rho_ff
        ! calulate velocity
?? copy :: velocity(vel_ff, pdf_ff, rhoInv, layout%fStencil)

        ! derive feq from velocity and density for all direction
        do iDirEq = 1, QQ
?? copy :: EqDir_comp(fEq_ff(iDirEq), rho_ff, vel_ff, layout, iDirEq, cs2inv, rk)
        end do

        ! density on boundary, qVal=0: boundary and fluid node overlaps
        rho_b = rho_w(posInBCelems(iLink))
        ! extrapolate velocity on boundary
        vel_b = vel_ff
        ! derive feq from velocity and density per direction
?? copy :: EqDir_comp(fEq_b, rho_b, vel_b, layout, iDir, cs2inv, rk)

        ! write into state: feq_b + fneq_ff/(1-omega)
        ! fneq is computed from post-collision pdf so we need to divide by
        ! (1-omega)
        omegaKine = omegaKineVal( posInState( posInBCelems(iLink) ) )
        relaxation_mrt = fieldProp%fluid%mrtPtr( omegaKine = omegaKine, &
          &                                      omegaBulk = omegaBulk, &
          &                                      QQ        = QQ         )
        nonEqScalingFacs(:) = 1.0_rk / ( 1.0_rk - relaxation_mrt(:) )

        fneq = pdf_ff - feq_ff
        ! convert fneq to moments neq
        mneq = matmul( toMoment(1:QQ, 1:QQ), fneq(1:QQ) )

        ! compute Minv times nonEqScalingFacs
        mInvXOmega(:) = toPDF(iDir, :) * nonEqScalingFacs(:)

        nEqTerm = dot_product( mInvXomega(1:QQ), mneq(1:QQ) )

        state( me%links(iLevel)%val(iLink) ) = feq_b + nEqTerm

      end do !iLink

    end associate

  end subroutine pressure_nonEqExpol_mrt
! **************************************************************************** !

! **************************************************************************** !
  !> Linkwise Dirichlet pressure non-equilibrium boundary condition for
  !! straight using the subroutine "mus_set_nonEqExpol" for incompressible model
  !! For straight wall, values are extrapolated along
  !! boundary normal instead of along the link and qVal=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_nonEqExpol_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
    ! -------------------------------------------------------------------- !
    ! also determined in mus_bc_header_module
    integer :: iDir
    ! variables for fictious boundary element
    real(kind=rk) :: feq_b, rho_b, vel_b(3)
    ! variables for surface (density and link-wise velocity)
    real(kind=rk) :: rho_w(globBC%nElems(iLevel))
    ! variables for overnext fluid element
    real(kind=rk) :: feq_ff, rho_ff, vel_ff(3)
    ! temporary local pdf values and the ones of overnext fluid
    real(kind=rk) :: pdf_ff(layout%fStencil%QQ)
    real(kind=rk) :: nEqTerm
    integer :: bcPress_pos, iLink, QQ
    ! relaxation parameter
    real(kind=rk) :: inv_rho_phy
    real(kind=rk) :: omegaKine
    ! --------------------------------------------------------------------------
    QQ = layout%fstencil%QQ
    inv_rho_phy = 1.0_rk / physics%fac(iLevel)%press * cs2inv

    ! position of boundary velocity in varSys
    bcPress_pos = me%bc_states%pressure%varPos
    ! get pressure_phy on boundary (surface) and store it to rho_w
    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_w                                 )

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

    associate( neighBufferPost => me%neigh(iLevel)%neighBufferPost,       &
      &        posInNeighBuf => me%nonEqExpol(iLevel)%posInNeighBuf,      &
      &        posInBCelems => me%nonEqExpol(iLevel)%posInBCelems,        &
      &        posInState => globBC%elemLvl(iLevel)%elem%val,             &
      &        omegaKineVal => fieldProp%fluid%viscKine%omLvl(iLevel)%val )

      ! loop over links
      do iLink = 1, me%links(iLevel)%nVals

        ! set link-wise direction
        iDir = me%nonEqExpol(iLevel)%iDir(iLink)

        ! determine needed quantities of the overnext fluid neighbor element
        ! x_ff
        pdf_ff = neighBufferPost(1, (posInNeighBuf(iLink)-1)*QQ+1  &
          &                       : (posInNeighBuf(iLink)-1)*QQ+QQ )
        ! fluid density
        rho_ff = sum(pdf_ff)
        ! calulate velocity
?? copy :: velocity(vel_ff, pdf_ff, rho0Inv, layout%fStencil)

        ! derive feq from velocity and density for all direction
?? copy :: EqDir_incomp(fEq_ff, rho_ff, vel_ff, layout, iDir, rho0, cs2inv, rk)

        ! density on boundary, qVal=0: boundary and fluid node overlaps
        rho_b = rho_w(posInBCelems(iLink))
        ! extrapolate velocity on boundary
        vel_b = vel_ff
        ! derive feq from velocity and density per direction
?? copy :: EqDir_incomp(fEq_b, rho_b, vel_b, layout, iDir, rho0, cs2inv, rk)

        ! write into state: feq_b + fneq_ff/(1-omega)
        ! fneq is computed from post-collision pdf so we need to divide by
        ! (1-omega)
        omegaKine = omegaKineVal( posInState( posInBCelems(iLink) ) )
        nEqTerm = (pdf_ff(iDir) - feq_ff) / (1.0_rk - omegaKine)

        state( me%links(iLevel)%val(iLink) ) = feq_b + nEqTerm

      end do !iLink

   end associate

  end subroutine pressure_nonEqExpol_incomp
! **************************************************************************** !


! **************************************************************************** !
  !> Linkwise Dirichlet pressure non-equilibrium boundary condition for
  !! straight for mrt collision using the subroutine "mus_set_nonEqExpol"
  !! for incompressible model
  !!
  !! 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_nonEqExpol_mrt_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
    ! -------------------------------------------------------------------- !
    ! also determined in mus_bc_header_module
    integer :: iDir, iDirLoop
    ! variables for fictious boundary element
    real(kind=rk) :: feq_b, rho_b, vel_b(3)
    ! variables for surface (density and link-wise velocity)
    real(kind=rk) :: rho_w(globBC%nElems(iLevel))
    ! variables for overnext fluid element
    real(kind=rk) :: feq_ff(layout%fStencil%QQ), rho_ff, vel_ff(3)
    ! temporary local pdf values and the ones of overnext fluid
    real(kind=rk) :: pdf_ff(layout%fStencil%QQ)
    real(kind=rk) :: fneq(layout%fStencil%QQ)
    real(kind=rk) :: mneq(layout%fStencil%QQ)
    real(kind=rk) :: mInvXomega(layout%fStencil%QQ)
    real(kind=rk) :: nonEqScalingFacs(layout%fStencil%QQ)
    real(kind=rk) :: relaxation_mrt(layout%fStencil%QQ)
    real(kind=rk) :: nEqTerm
    integer :: bcPress_pos, iLink, QQ
    ! relaxation parameter
    real(kind=rk) :: inv_rho_phy
    real(kind=rk) :: omegaKine
    ! --------------------------------------------------------------------------
    QQ = layout%fstencil%QQ
    inv_rho_phy = 1.0_rk / physics%fac(iLevel)%press * cs2inv

    ! position of boundary velocity in varSys
    bcPress_pos = me%bc_states%pressure%varPos
    ! get pressure_phy on boundary (surface) and store it to rho_w
    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_w                                 )

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

    associate( neighBufferPost => me%neigh(iLevel)%neighBufferPost,        &
      &        posInNeighBuf => me%nonEqExpol(iLevel)%posInNeighBuf,       &
      &        posInBCelems => me%nonEqExpol(iLevel)%posInBCelems,         &
      &        posInState => globBC%elemLvl(iLevel)%elem%val,              &
      &        omegaKineVal => fieldProp%fluid%viscKine%omLvl(iLevel)%val, &
      &        omegaBulk  => fieldProp%fluid%omegaBulkLvl(iLevel),         &
      &        toMoment => layout%moment%toMoments%A,                      &
      &        toPDF => layout%moment%toPDF%A                              )

      ! loop over links
      do iLink = 1, me%links(iLevel)%nVals

        ! set link-wise direction
        iDir = me%nonEqExpol(iLevel)%iDir(iLink)

        ! determine needed quantities of the overnext fluid neighbor element
        ! x_ff
        pdf_ff = neighBufferPost(1, (posInNeighBuf(iLink)-1)*QQ+1  &
          &                       : (posInNeighBuf(iLink)-1)*QQ+QQ )
        ! fluid density
        rho_ff = sum(pdf_ff)
        ! calulate velocity
?? copy :: velocity(vel_ff, pdf_ff, rho0Inv, layout%fStencil)

      ! derive feq from velocity and density for all direction
        do iDirLoop = 1, layout%fStencil%QQ
?? copy :: EqDir_incomp(fEq_ff(iDirLoop), rho_ff, vel_ff, layout, iDirLoop, rho0, cs2inv, rk)
        end do

        ! density on boundary, qVal=0: boundary and fluid node overlaps
        rho_b = rho_w(posInBCelems(iLink))
        ! extrapolate velocity on boundary
        vel_b = vel_ff
        ! derive feq from velocity and density per direction
?? copy :: EqDir_incomp(fEq_b, rho_b, vel_b, layout, iDir, rho0, cs2inv, rk)

        ! write into state: feq_b + fneq_ff/(1-omega)
        ! fneq is computed from post-collision pdf so we need to divide by
        ! (1-omega)
        omegaKine = omegaKineVal( posInState( posInBCelems(iLink) ) )
        relaxation_mrt = fieldProp%fluid%mrtPtr( omegaKine = omegaKine, &
          &                                      omegaBulk = omegaBulk, &
          &                                      QQ        = QQ         )
        nonEqScalingFacs(:) = 1.0_rk / ( 1.0_rk - relaxation_mrt(:) )

        fneq = pdf_ff - feq_ff
        ! convert fneq to moments neq
        mneq = matmul( toMoment(1:QQ, 1:QQ), fneq(1:QQ) )

        ! compute Minv times nonEqScalingFacs
        mInvXOmega(:) = toPDF(iDir, :) * nonEqScalingFacs(:)

        nEqTerm = dot_product( mInvXomega(1:QQ), mneq(1:QQ) )

        state( me%links(iLevel)%val(iLink) ) = feq_b + nEqTerm

      end do !iLink

   end associate

  end subroutine pressure_nonEqExpol_mrt_incomp
! **************************************************************************** !

end module mus_bc_fluid_nonEqExpol_module
