! 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-2017, 2019-2020 Kannan Masilamani <kannan.masilamani@uni-siegen.de>
! Copyright (c) 2013-2014 Kartik Jain <kartik.jain@uni-siegen.de>
! Copyright (c) 2013 Harald Klimach <harald.klimach@uni-siegen.de>
! Copyright (c) 2014-2016 Jiaxing Qi <jiaxing.qi@uni-siegen.de>
! Copyright (c) 2016 Verena Krupp <verena.krupp@uni-siegen.de>
! Copyright (c) 2016 Tobias Schneider <tobias1.schneider@student.uni-siegen.de>
! Copyright (c) 2017-2018 Raphael Haupt <raphael.haupt@uni-siegen.de>
! Copyright (c) 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.
?? include 'header/lbm_macros.inc'
! ****************************************************************************** !
!> Boundary condition treatment routines for fluid simulation
!!
!! A detailed description on the implementation details are given in
!! \ref boundary_implementation
!!
module mus_bc_fluid_experimental_module
  use iso_c_binding, only: c_ptr, c_f_pointer

  ! include treelm modules
  use env_module,               only: rk
  use tem_param_module,         only: 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_debug_module,         only: dbgUnit
  use tem_logging_module,       only: logUnit
  use tem_construction_module,  only: tem_levelDesc_type
  use tem_property_module,      only: prp_hasQVal
  use tem_aux_module,           only: tem_abort

  ! 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_varSys_module,         only: mus_varSys_data_type


  implicit none

  private

  public :: pressure_expol_slow

?? IF( .not. SOA) then
  public :: moments_inflow, moments_outflow
  public :: spc_moments_outflow
  public :: spc_bb_wall
  public :: spc_bb_vel_test
?? ENDIF


contains

?? IF( .not. SOA) then
! ****************************************************************************** !
  !> species bounce back velocity boundary
  !! Usage
  !! -----
  !!```lua
  !!boundary_condition = {
  !! { label = 'outlet',
  !!   kind = 'spc_bb_vel_test',
  !!   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 spc_bb_vel_test( 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( layout%fStencil%QQ )
    real(kind=rk) :: fEq( layout%fStencil%QQ, varSys%nStateVars )
    real(kind=rk) :: fTmp_all( layout%fStencil%QQ*varSys%nStateVars )
    integer :: iELem, iDir, iFieldLoc, nFields, pos, invPos
    integer :: iNeigh, elemPos, fluidNeighbor, QQ
    real(kind=rk) :: mass_dens(varSys%nStateVars), vel(3, varSys%nStateVars)
    real(kind=rk) :: fTmp_allNeigh( layout%fStencil%QQ*varSys%nStateVars )
    real(kind=rk) :: fTmp_Next( layout%fStencil%QQ*varSys%nStateVars )
    !real(kind=rk) :: fEqNeigh( layout%fStencil%QQ )
    real(kind=rk) :: fEqNeigh( layout%fStencil%QQ, varSys%nStateVars )
    real(kind=rk) :: fEqNext( layout%fStencil%QQ )
    real(kind=rk) :: ucx
    real(kind=rk) :: vel_b(globBC%nElems(iLevel)*3), inv_vel
    integer :: posInBuffer, bcVel_pos, offset
    ! ------------------------------------------------------------------------
    write(logUnit(1),*) 'WARNING: Experimental BC spc_bb_vel_test'
!!write(dbgUnit(1),*) 'Boundary label ', trim(me%label)
!!write(dbgUnit(1),*) 'iField ', iField

    nFields = varSys%nStateVars
    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                               )

    ! 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

    do iElem = 1, globBC%nElems(iLevel)
      posInBuffer = globBC%elemLvl( iLevel )%posInBcElemBuf%val( iElem )
!!write(dbgUnit(1),*) 'iElem ', iElem
      fTmp_all = 0.0_rk
      ! Calculate the density of current element
      do iFieldLoc = 1, nFields
        do iDir = 1, QQ
          pos = varSys%method%val(iFieldLoc)%state_varPos(iDir)
          fTmp_all(pos) = bcBuffer( pos + ( posInBuffer-1)*nScalars )
        end do
      end do !iField

      fTmp(1:QQ) = bcBuffer( (posInBuffer-1)*nScalars+varPos(1) : &
        &                    (posInBuffer-1)*nScalars+varPos(1)+QQ-1 )

!!write(dbgUnit(1),*) 'fTmp_all ', fTmp_all
      fTmp_Next = 0.0_rk
      do iFieldLoc = 1, nFields
        call derVarPos%equilFromState( state  = fTmp_All,         &
          &                            iField = iFieldLoc,        &
          &                            nElems = 1,                &
          &                            varSys = varSys,           &
          &                            layout = layout,           &
          &                            res    = fEq(:, iFieldLoc) )

        !write(*,*) 'iField ', iFieldLoc, 'local fEq ', fEq
        call tem_abort('Spc_vel_bb_test: bitmast is allcoated with QQN')
        do iDir = 1,layout%fStencil%QQ
          if( globBC%elemLvl(iLevel)%bitmask%val( iDir, iElem ) .or. &
            & iDir == layout%fStencil%restPosition ) then
            pos = varSys%method%val(iFieldLoc)%state_varPos(iDir)
            invPos = varSys%method%val(iFieldLoc)%state_varPos( &
              & layout%fStencil%cxDirInv(iDir) )
            !write(*,*) 'iDir ', iDir, pos, invPos
            fTmp_Next(pos) = ( fTmp_all( invPos )                      &
              &            + mixture%omega_diff * 0.5_rk *             &
              &              fEq( layout%fStencil%cxDirInv(iDir), iFieldLoc ) ) &
              &            / (1.0_rk+mixture%omega_diff*0.5_rk)
          end if
        end do
      end do

!!write(dbgUnit(1),*) 'local fTmp_Next ', fTmp_Next

      ! current element position in state array
      elemPos = globBC%elemLvl(iLevel)%elem%val(iElem)
      offset = (iElem-1)*3
!!write(dbgUnit(1),*) 'elemPos', elemPos, 'treeID ', levelDesc%total(elemPos)

      ! Transform g to f in each neighbor node
      do iNeigh = 1,layout%fStencil%QQN
        ! fluid neighbor in inverse direction of iNeigh (x-c_i)
        ! @todo: fluid neighbor should be found before
        fluidNeighbor = &
 & ?ElemIDX?(neigh( ?NGPOS?( iNeigh, elemPos, nSize )), nScalars, nSize)

        !if (elemPos /= fluidNeighbor) then
        if( .not. globBC%elemLvl(iLevel)%bitmask%val( iNeigh, iElem ) ) then
          if (elemPos == fluidNeighbor) then
            fTmp_allNeigh =  fTmp_all
          else
            ! get pdf of neighbor in iNeigh
            do iFieldLoc = 1, nFields
              do iDir = 1, QQ
                pos = varSys%method%val(iFieldLoc)%state_varPos(iDir)
                fTmp_allNeigh(pos) = &
& state( ?SAVE?( iDir, iFieldLoc, fluidNeighbor, QQ, nScalars,nSize, neigh ) )
              end do
            end do !iFieldLoc
          end if

!!write(dbgUnit(1),*) 'iDir ', iNeigh, 'fluidNeighbor ', fluidNeighbor
!!write(dbgUnit(1),*) 'neighID ', levelDesc%total(fluidNeighbor)
!!write(dbgUnit(1),*) 'fTmp_allNeigh ', fTmp_allNeigh

          ! compute equilibrium of neighbor
          fEqNeigh = 0.0_rk
          do iFieldLoc = 1, nFields
            call derVarPos%equilFromState( state  = fTmp_AllNeigh,         &
              &                            iField = iFieldLoc,             &
              &                            nElems = 1,                     &
              &                            varSys = varSys,                &
              &                            layout = layout,                &
              &                            res    = fEqNeigh(:, iFieldLoc) )
          !write(*,*) 'iFieldLoc ', iFieldLoc, 'fEqNeigh ', fEqNeigh(:, iFieldLoc)
          !write(*,*) 'iFieldLoc ', iFieldLoc, 'sum(fEqNeigh) ', sum(fEqNeigh(:, iFieldLoc))
          end do

          do iFieldLoc = 1, nFields
            pos = varSys%method%val(iFieldLoc)%state_varPos(iNeigh)
            fTmp_Next(pos) = ( fTmp_allNeigh(pos) &
              &               + mixture%omega_diff*0.5_rk*fEqNeigh(iNeigh, iFieldLoc) ) &
              &               / (1.0_rk+mixture%omega_diff*0.5_rk)
          !fTmp_Next(varSys%method%val(iFieldLoc)%state_varPos(tmpDir)) = &
          !  & ( fTmp_allNeigh(varSys%method%val(iFieldLoc)%state_varPos(iNeigh)) &
          !  &               + mixture%omega_diff*0.5_rk*fEqNeigh(iNeigh, iFieldLoc) ) &
          !  &               / (1.0_rk+mixture%omega_diff*0.5_rk)
          end do
        end if
      end do !iNeigh
!write(dbgUnit(1),*) 'fTmp_Next ', fTmp_Next

      mass_dens = 0.0_rk
      do iFieldLoc = 1, nFields
        do iDir = 1,QQ
          pos = varSys%method%val(iFieldLoc)%state_varPos(iDir)
          mass_dens(iFieldLoc) = mass_dens(iFieldLoc) + fTmp_Next(pos)
        end do
      end do
!write(dbgUnit(1),*) 'mass_dens ', mass_dens
!write(dbgUnit(1),*) 'vel ', ux(iElem), uy(iElem), uz(iElem)

      ! update untransformed f with velocity bounce back
      do iDir = 1,QQ
        if( globBC%elemLvl(iLevel)%bitmask%val( iDir, iElem ) ) then
!write(dbgUnit(1),*) 'iDir', iDir
          ucx = layout%fStencil%cxDir(1,iDir) * vel_b(offset+1)   &
            & + layout%fStencil%cxDir(2,iDir) * vel_b(offset+2)   &
            & + layout%fStencil%cxDir(3,iDir) * vel_b(offset+3)
          do iFieldLoc = 1, nFields
            pos = varSys%method%val(iFieldLoc)%state_varPos(iDir)
            fTmp_Next(pos) = fTmp_Next(pos) &
              & + layout%weight(iDir) * 2._rk * cs2inv * mass_Dens(iFieldLoc) &
              & * ucx
          end do
        end if
      end do
!write(dbgUnit(1),*) 'after ubb fTmp_Next ', fTmp_Next

      mass_dens = 0.0_rk
      vel = 0.0_rk
      do iFieldLoc = 1, nFields
        do iDir = 1,QQ
          pos = varSys%method%val(iFieldLoc)%state_varPos(iDir)
          mass_dens(iFieldLoc) = mass_dens(iFieldLoc) + fTmp_Next(pos)
          vel(:, iFieldLoc)  = vel(:, iFieldLoc) &
            & + real(layout%fStencil%cxDir( :, iDir), kind=rk) * fTmp_Next(pos)
        end do
        !vel(:, iFieldLoc)  = (/ ux(iElem), uy(iElem), uz(iElem) /)
      end do
      do iFieldLoc = 1, nFields
        vel(:, iFieldLoc) = vel(:,iFieldLoc)/mass_dens(iFieldLoc)
      end do

!write(dbgUnit(1),*) 'after ubb mass_dens ', mass_dens
!write(dbgUnit(1),*) 'after ubb vel ', vel

!write(dbgUnit(1),*) 'fEqNext ', fEqNext
      ! Calculate the equilibrium distribution
      call derVarPos%equilFromMacro( density  = mass_dens, &
        &                            velocity = vel,       &
        &                            iField   = iField,    &
        &                            nElems   = 1,         &
        &                            varSys   = varSys,    &
        &                            layout   = layout,    &
        &                            res      = fEqNext    )

!KM!      ! density
!KM!      rho_spc(1) = sum(fTmp_all(varSys%method%val(iField)%state_varPos(:)))
!KM!
!KM!      ! equilibrium velocity
!KM!      call derVarPos%velocitiesFromState( state  = fTmp_Next, &
!KM!        &                                 iField = iField,    &
!KM!        &                                 nElems = 1,         &
!KM!        &                                 varSys = varSys,    &
!KM!        &                                 layout = layout,    &
!KM!        &                                 res    = velocities )
!KM!
!KM!      eqVel = equilVelFromMacro( iField, moleFraction, velocities, nFields, &
!KM!        &                        paramBInv, phi, resi_coeff )
!KM!
      ! set boundary
      do iDir = 1, QQ
        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, globBC%elemLvl(iLevel)%elem%val(iElem), QQ, nScalars,nSize, neigh ) ) &
            !! simple Bounce Back
            !& = fTmp(layout%fStencil%cxDirInv( iDir ))

            ! using f in next time step
            & = fTmp_Next(varPos(iDir)) &
            & + mixture%omega_diff*0.5_rk*( fTmp_Next(varPos(iDir)) &
            & - fEqNext(iDir) )

            ! Bounce back with correction term for g with fEq in next time step
            !& = fTmp(layout%fStencil%cxDirInv( iDir ))  &
            !& + mixture%omega_diff*0.5_rk*(fEq(layout%fStencil%cxDirInv( iDir ), iField)&
            !& - fEqNext(iDir) )

            !& = fTmp(layout%fStencil%cxDirInv( iDir ))  &
            !& - layout%weight( iDir )&
            !& * mixture%omega_diff*cs2inv*rho_spc(1)                        &
            !&    * ( layout%fStencil%cxDir( 1, iDir )*eqVel(1)            &
            !&    +   layout%fStencil%cxDir( 2, iDir )*eqVel(2)            &
            !&    +   layout%fStencil%cxDir( 3, iDir )*eqVel(3))

        end if
      end do
    end do

!KM!    ! debug output
!KM!    if (iField == nFields) then
!KM!    do iElem = 1, globBC%nElems( iLevel )
!KM!      !ux = 0.0_rk
!KM!      do iFieldLoc = 1, nFields
!KM!        do iDir = 1,layout%fStencil%QQ
!KM!          pos = varSys%method%val(iFieldLoc)%state_varPos(iDir)
!KM!          fTmp_all(pos) = state(                                                  &
!KM!&?FETCH?( iDir, iFieldLoc, globBC%elemLvl(iLevel)%elem%val( iElem ), QQ, nScalars, neigh ) &
!KM!        & )
!KM!        end do
!KM!      end do
!KM!      bound(1)=iElem
!KM!      bound(2)=iElem
!KM!      write(dbgUnit(1),*) 'iElem ', iElem
!KM!      write(dbgUnit(1),*) 'fTmp_all ', fTmp_all
!KM!      ! density
!KM!      do iFieldLoc = 1, nFields
!KM!        rhoTmp = sum(fTmp_all(varSys%method%val(iFieldLoc)%state_varPos(:)))
!KM!
!KM!        ! mass flux
!KM!        call derVarPos%momFromState( state  = fTmp_all,  &
!KM!          &                          iField = iFieldLoc, &
!KM!          &                          nElems = 1,         &
!KM!          &                          varSys = varSys,    &
!KM!          &                          layout = layout,    &
!KM!          &                          res    = uxTmp      )
!KM!        write(dbgUnit(1),*) 'iField ', iFieldLoc, 'density ', rhoTmp
!KM!        write(dbgUnit(1),*) 'iField ', iFieldLoc, 'velocity ', uxTmp
!KM!      end do
!KM!    end do
!KM!    end if
  end subroutine spc_bb_vel_test
! ****************************************************************************** !


! ****************************************************************************** !
  !> species bounce back wall boundary
  !! Usage
  !! -----
  !!```lua
  !!boundary_condition = {
  !! { label = 'outlet',
  !!    kind = 'spc_bb_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 spc_bb_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) :: fEq( layout%fStencil%QQ*globBC%nElems(iLevel) )
!    real(kind=rk) :: fEq
    real(kind=rk) :: fTmp( layout%fStencil%QQ )
    real(kind=rk) :: fEq( layout%fStencil%QQ, varSys%nStateVars )
    real(kind=rk) :: fTmp_all( layout%fStencil%QQ*varSys%nStateVars )
!    real(kind=rk) :: eqVel(3)
    integer :: iELem, iDir, iFieldLoc, nFields, pos, invPos, QQ
    integer :: iNeigh, elemPos, fluidNeighbor
    real(kind=rk) :: mass_dens(varSys%nStateVars), vel(3, varSys%nStateVars)
    real(kind=rk) :: fTmp_allNeigh( layout%fStencil%QQ*varSys%nStateVars )
    real(kind=rk) :: fTmp_Next( layout%fStencil%QQ*varSys%nStateVars )
    !real(kind=rk) :: fEqNeigh( layout%fStencil%QQ )
    real(kind=rk) :: fEqNeigh( layout%fStencil%QQ, varSys%nStateVars )
    real(kind=rk) :: fEqNext( layout%fStencil%QQ )
    integer :: posInBuffer
    ! ------------------------------------------------------------------------
    write(logUnit(1),*) 'WARNING: Experimental BC spc_bb_wall'
!write(dbgUnit(1),*) 'Boundary label ', trim(me%label)
!write(dbgUnit(1),*) 'iField ', iField
    nFields = varSys%nStateVars
    QQ = layout%fStencil%QQ

    do iElem = 1, globBC%nElems(iLevel)

      posInBuffer = globBC%elemLvl( iLevel )%posInBcElemBuf%val( iElem )
!write(dbgUnit(1),*) 'iElem ', iElem
      fTmp_all = 0.0_rk
      ! Calculate the density of current element
      do iFieldLoc = 1, nFields
        do iDir = 1,QQ
          pos = varSys%method%val(iFieldLoc)%state_varPos(iDir)
          fTmp_all(pos) = bcBuffer( pos+(posInBuffer-1)*nScalars )
        end do
      end do !iField

      fTmp(1:QQ) = bcBuffer( (posInBuffer-1)*nScalars+varPos(1) : &
        &                    (posInBuffer-1)*nScalars+varPos(1)+QQ-1 )

!write(dbgUnit(1),*) 'fTmp_all ', fTmp_all
      fTmp_Next = 0.0_rk
      do iFieldLoc = 1, nFields
        call derVarPos%equilFromState( state  = fTmp_All,         &
          &                            iField = iFieldLoc,        &
          &                            nElems = 1,                &
          &                            varSys = varSys,           &
          &                            layout = layout,           &
          &                            res    = fEq(:, iFieldLoc) )

!write(dbgUnit(1),*) 'iFieldLoc ', iFieldLoc, 'fEqLocal ', fEq(:, iFieldLoc)

        !write(*,*) 'iField ', iFieldLoc, 'local fEq ', fEq
        call tem_abort('Spc_vel_wall: bitmast is allcoated with QQN')
        do iDir = 1,QQ
          if( globBC%elemLvl(iLevel)%bitmask%val( iDir, iElem ) .or. &
            & iDir == layout%fStencil%restPosition ) then
            pos = varSys%method%val(iFieldLoc)%state_varPos(iDir)
            invPos = varSys%method%val(iFieldLoc)%state_varPos( &
              & layout%fStencil%cxDirInv(iDir) )
            !write(*,*) 'iDir ', iDir, pos, invPos
            fTmp_Next(pos) = ( fTmp_all( invPos )                              &
              &            + mixture%omega_diff * 0.5_rk *                     &
              &              fEq( layout%fStencil%cxDirInv(iDir), iFieldLoc ) )&
              &            / (1.0_rk+mixture%omega_diff*0.5_rk)
          end if
        end do
      end do

!write(dbgUnit(1),*) 'local fTmp_Next ', fTmp_Next

      ! current element position in state array
      elemPos = globBC%elemLvl(iLevel)%elem%val(iElem)
!write(dbgUnit(1),*) 'elemPos', elemPos, 'treeID ', levelDesc%total(elemPos)

      ! Transform g to f in each neighbor node
      do iNeigh = 1,layout%fStencil%QQN
        ! fluid neighbor in inverse direction of iNeigh (x-c_i)
        fluidNeighbor = &
 & ?ElemIDX?(neigh( ?NGPOS?( iNeigh, elemPos, nSize )), nScalars, nSize)

        !if (elemPos /= fluidNeighbor) then
        if( .not. globBC%elemLvl(iLevel)%bitmask%val( iNeigh, iElem ) ) then
          if (elemPos == fluidNeighbor) then
            fTmp_allNeigh =  fTmp_all
          else
            ! get pdf of neighbor in iNeigh
            do iFieldLoc = 1, nFields
              do iDir = 1,QQ
                pos = varSys%method%val(iFieldLoc) &
                  &         %state_varPos(iDir)
                fTmp_allNeigh(pos) =  &
& state( ?SAVE?( iDir, iFieldLoc, fluidNeighbor, QQ, nScalars, nSize, neigh ) )
              end do
            end do !iFieldLoc
          end if

!write(dbgUnit(1),*) 'iDir ', iNeigh, 'fluidNeighbor ', fluidNeighbor
!write(dbgUnit(1),*) 'neighID ', levelDesc%total(fluidNeighbor)
!write(dbgUnit(1),*) 'fTmp_allNeigh ', fTmp_allNeigh

          ! compute equilibrium of neighbor
          fEqNeigh = 0.0_rk
          do iFieldLoc = 1, nFields
            call derVarPos%equilFromState( state  = fTmp_AllNeigh,         &
              &                            iField = iFieldLoc,             &
              &                            nElems = 1,                     &
              &                            varSys = varSys,                &
              &                            layout = layout,                &
              &                            res    = fEqNeigh(:, iFieldLoc) )

!write(dbgUnit(1),*) 'iFieldLoc ', iFieldLoc, 'fEqNeigh ', fEqNeigh(:, iFieldLoc)
!write(dbgUnit(1),*) 'iFieldLoc ', iFieldLoc, 'sum(fEqNeigh) ', sum(fEqNeigh(:, iFieldLoc))
          end do

          do iFieldLoc = 1, nFields
            pos = varSys%method%val(iFieldLoc)%state_varPos(iNeigh)
            !fTmp_Next(pos) = fTmp_allNeigh(pos)
            fTmp_Next(pos) = ( fTmp_allNeigh(pos) &
              &               + mixture%omega_diff*0.5_rk*fEqNeigh(iNeigh, iFieldLoc) ) &
              &               / (1.0_rk+mixture%omega_diff*0.5_rk)
          !fTmp_Next(varSys%method%val(iFieldLoc)%state_varPos(tmpDir)) = &
          !  & ( fTmp_allNeigh(varSys%method%val(iFieldLoc)%state_varPos(iNeigh)) &
          !  &               + mixture%omega_diff*0.5_rk*fEqNeigh(iNeigh, iFieldLoc) ) &
          !  &               / (1.0_rk+mixture%omega_diff*0.5_rk)
          end do
        end if
      end do !iNeigh
!write(dbgUnit(1),*) 'fTmp_Next ', fTmp_Next

      mass_dens = 0.0_rk
      vel = 0.0_rk
      do iFieldLoc = 1, nFields
        do iDir = 1,QQ
          pos = varSys%method%val(iFieldLoc)%state_varPos(iDir)
          mass_dens(iFieldLoc) = mass_dens(iFieldLoc) + fTmp_Next(pos)
          vel(:, iFieldLoc)  = vel(:, iFieldLoc) &
            & + real(layout%fStencil%cxDir( :, iDir), kind=rk) * fTmp_Next(pos)
        end do
      end do
      do iFieldLoc = 1, nFields
        vel(:, iFieldLoc) = vel(:,iFieldLoc)/mass_dens(iFieldLoc)
      end do

!write(dbgUnit(1),*) 'mass_dens ', mass_dens
!write(dbgUnit(1),*) 'vel ', vel

      ! Calculate the equilibrium distribution
      call derVarPos%equilFromMacro( density  = mass_dens, &
        &                            velocity = vel,       &
        &                            iField   = iField,    &
        &                            nElems   = 1,         &
        &                            varSys   = varSys,    &
        &                            layout   = layout,    &
        &                            res      = fEqNext    )

!write(dbgUnit(1),*) 'fEqNext ', fEqNext

!KM!      ! density
!KM!      rho_spc(1) = sum(fTmp_all(varSys%method%val(iField)%state_varPos(:)))
!KM!
!KM!      ! equilibrium velocity
!KM!      call derVarPos%velocitiesFromState( state  = fTmp_Next, &
!KM!        &                                 iField = iField,    &
!KM!        &                                 nElems = 1,         &
!KM!        &                                 varSys = varSys,    &
!KM!        &                                 layout = layout,    &
!KM!        &                                 res    = velocities )
!KM!
!KM!      eqVel = equilVelFromMacro( iField, moleFraction, velocities, nFields, &
!KM!        &                        paramBInv, phi, resi_coeff )
!KM!
      ! set boundary
      do iDir = 1,QQ
        if( globBC%elemLvl(iLevel)%bitmask%val( iDir, iElem )) then
          state(                                                               &
&?FETCH?( iDir, iField, globBC%elemLvl( iLevel )%elem%val( iElem ), QQ, nScalars,nSize, neigh ) ) &
              !! simple Bounce Back
              !& = fTmp(layout%fStencil%cxDirInv( iDir ))

              ! using f in next time step
              & = fTmp_Next(varPos(iDir)) &
              & + mixture%omega_diff*0.5_rk*( fTmp_Next(varPos(iDir)) &
              & - fEqNext(iDir) )

              ! Bounce back with correction term for g with fEq in next time step
              !& = fTmp(layout%fStencil%cxDirInv( iDir ))  &
              !& + mixture%omega_diff*0.5_rk*(fEq(layout%fStencil%cxDirInv( iDir ), iField)&
              !& - fEqNext(iDir) )

              !& = fTmp(layout%fStencil%cxDirInv( iDir ))  &
              !& + mixture%omega_diff*0.5_rk*(fEq(layout%fStencil%cxDirInv( iDir ), iField)&
              !& - fEq(iDir, iField) )

              !& = fTmp(layout%fStencil%cxDirInv( iDir ))  &
              !& - layout%weight( iDir )&
              !& * mixture%omega_diff*cs2inv*rho_spc(1)                        &
              !&    * ( layout%fStencil%cxDir( 1, iDir )*eqVel(1)            &
              !&    +   layout%fStencil%cxDir( 2, iDir )*eqVel(2)            &
              !&    +   layout%fStencil%cxDir( 3, iDir )*eqVel(3))

        end if
      end do
    end do

    ! debug output
!KM!    if (iField == nFields) then
!KM!    do iElem = 1, globBC%nElems( iLevel )
!KM!      !ux = 0.0_rk
!KM!      do iFieldLoc = 1, nFields
!KM!        do iDir = 1,layout%fStencil%QQ
!KM!          pos = varSys%method%val(iFieldLoc)%state_varPos(iDir)
!KM!          fTmp_all(pos) = state(                                                  &
!KM!&?FETCH?( iDir, iFieldLoc, globBC%elemLvl(iLevel)%elem%val( iElem ), nScalars, neigh ) &
!KM!        & )
!KM!        end do
!KM!      end do
!KM!      bound(1)=iElem
!KM!      bound(2)=iElem
!KM!!write(dbgUnit(1),*) 'iElem ', iElem
!KM!!write(dbgUnit(1),*) 'fTmp_all ', fTmp_all
!KM!      ! density
!KM!      do iFieldLoc = 1, nFields
!KM!        rhoTmp = sum(fTmp_all(varSys%method%val(iFieldLoc)%state_varPos(:)))
!KM!
!KM!        ! mass flux
!KM!        call derVarPos%momFromState( state  = fTmp_all,  &
!KM!          &                          iField = iFieldLoc, &
!KM!          &                          nElems = 1,         &
!KM!          &                          varSys = varSys,    &
!KM!          &                          layout = layout,    &
!KM!          &                          res    = uxTmp      )
!KM!!write(dbgUnit(1),*) 'iField ', iFieldLoc, 'density ', rhoTmp
!KM!!write(dbgUnit(1),*) 'iField ', iFieldLoc, 'velocity ', uxTmp
!KM!      end do
!KM!    end do
!KM!    end if
  end subroutine spc_bb_wall
! ****************************************************************************** !


! ****************************************************************************** !
  !> molar flux boundary condition like moments velocity bc type
  !! Usage
  !! -----
  !!```lua
  !!boundary_condition = {
  !! { label = 'outlet',
  !!    kind = 'spc_moments_outflow',
  !!    moleflux = 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 spc_moments_outflow( 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 :: iELem, iDir, iFieldLoc, nFields, pos
    real(kind=rk) :: rho, press
    integer :: nLinks, iLink, elemPos, neighPos, 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) :: fTmp_all( layout%fStencil%QQ*varSys%nStateVars )
    real(kind=rk) :: toMoments( layout%fStencil%QQ )
    real(kind=rk) :: moments( layout%fStencil%QQ )
    real(kind=rk) :: vel_inG(3)
    real(kind=rk) :: mom(3,varSys%nStateVars), mom_inG(3)
    real(kind=rk) :: mass_dens(varSys%nStateVars), num_dens(varSys%nStateVars)
    real(kind=rk) :: totMassDens, velAvg(3)
    real(kind=rk) :: moleFrac(varSys%nStateVars)
    real(kind=rk) :: resi_coeff( varSys%nStateVars, varSys%nStateVars )
    real(kind=rk) :: molWeight(varSys%nStateVars), phi(varSys%nStateVars)
    type(mus_varSys_data_type), pointer :: fPtr
    ! ------------------------------------------------------------------------
    write(logUnit(1),*) 'WARNING: Experimental BC spc_moments_outflow'
    !write(*,*) 'Boundary label ', trim(me%label)

    nFields = varSys%nStateVars
    QQ = layout%fStencil%QQ

    call C_F_POINTER( varSys%method%val(iField)%method_Data, fPtr )

    do iFieldLoc = 1, nFields
      ! species properties
      ! molecular weight inverse
      molWeight(iFieldLoc) =                                      &
        &  fPtr%solverData%scheme%field(iFieldLoc)%fieldProp%species%molWeight
      ! molecular weight ratio
      phi(iFieldLoc) = &
        & fPtr%solverData%scheme%field(iFieldLoc)%fieldProp%species%molWeigRatio
      ! resistivity coefficients
      resi_coeff(iFieldLoc, :) =                                     &
        & fPtr%solverData%scheme%field(iFieldLoc)%fieldProp%species%resi_coeff(:)
    end do

!    write(*,*) 'iField ', iField
    ! Calculate the density of current element
    do iElem = 1, globBC%nElems(iLevel)
!      write(*,*) 'iElem ', iElem

      fTmp = 0.0_rk
      rho = 0
      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 )
      neighPos = me%neigh(iLevel)%posInState(1,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 state vector to compute density and velocity of other species
      do iFieldLoc = 1, varSys%nStateVars
        do iDir = 1,QQ
          pos = varSys%method%val(iFieldLoc)%state_varPos(iDir)
          fTmp_all( pos ) = &
!& state(?FETCH?( iDir, iFieldLoc, elemPos, QQ, nScalars, neigh ) )
& state(?SAVE?( iDir, iFieldLoc, neighPos, QQ, nScalars,nSize, neigh ) )
        end do
      end do

      do iFieldLoc = 1, nFields
        mass_dens(iFieldLoc) =                                          &
          & sum(fTmp_all( varSys%method%val(iFieldLoc)%state_varPos(:) ))
        num_dens(iFieldLoc) = mass_dens(iFieldLoc) &
          &                 / molWeight(iFieldLoc)

        ! mass flux
        call derVarPos%momFromState( state  = fTmp_all,         &
          &                          iField = iFieldLoc,        &
          &                          nElems = 1,                &
          &                          varSys = varSys,           &
          &                          layout = layout,           &
          &                          res    = mom(:, iFieldLoc) )
      end do

      ! mass density
      mass_dens = num_dens * molWeight
      totMassDens = sum(mass_dens)

      rho = mass_dens(iField)
!      write(*,*) 'density ', rho

      ! mole fraction
      moleFrac = num_dens/sum(num_dens)

!      do iFieldLoc = 1, nFields
!        write(*,*) iFieldLoc, 'vel ', vel(:, iFieldLoc), &
!          & 'dens ', mass_dens(iFieldLoc), 'moleFrac ', moleFrac(iFieldLoc)
!      end do

      ! compute velocity of transformed variable g. Eq. which is normally
      ! used to compute actual velocity by solving LSE
      mom_inG = mom(:, iField)
      do iFieldLoc = 1, nFields
        mom_inG = mom_inG + ( mixture%omega_diff * 0.5_rk                    &
            &               * resi_coeff(iField, iFieldLoc) * phi(iField)    &
            &               * moleFrac(iFieldLoc) / mixture%paramB )         &
            &             * mom(:, iField)
      end do
      !set non-diagonal part
      do iFieldLoc = 1, nFields
        mom_inG = mom_inG - ( mixture%omega_diff * 0.5_rk                    &
          &                 * resi_coeff(iField, iFieldLoc) * phi(iFieldLoc) &
          &                 * moleFrac(iField) / mixture%paramB )            &
          &               * mom(:,iFieldLoc)
      end do

      vel_inG = mom_inG/rho
!      write(*,*) 'transformed momentum ', mom_inG

      ! mass averaged mixture velocity
      velAvg(1) = sum(mom(1,:)) / totMassDens
      velAvg(2) = sum(mom(2,:)) / totMassDens
      velAvg(3) = sum(mom(3,:)) / totMassDens

!      write(*,*) 'velAvg ', velAvg

      press = cs2*rho*fieldProp%species%molWeigRatio
      !write(*,*) 'rho ', rho, 'press', press

      toMoments = matmul(layout%moment%toMoments%A, fTmp)
      !write(*,*) 'fTmp ', fTmp
      !write(*,*) 'toMoments ', toMoments

      ! moments for liquid mixture
      moments = (/ rho, mom_inG(1), mom_inG(2), press+rho*velAvg(1)**2,    &
        & press+rho*velAvg(2)**2, rho*velAvg(1)*velAvg(2), cs2*mom_inG(2), &
        & cs2*mom_inG(1), cs2*(press+rho*velAvg(1)**2+rho*velAvg(2)**2) /)

      !write(*,*) 'Moments ', moments
      !write(*,*) 'knownMomPos ', me%elemLvl(iLevel)%moments(iElem)%knownMom_pos

      do iLink = 1, nLinks
        rhs(iLink) = &
          & moments( me%elemLvl(iLevel)%moments(iElem)%knownMom_pos(iLink) ) &
          & - toMoments( me%elemLvl(iLevel)%moments(iElem)%knownMom_pos(iLink) )
      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 spc_moments_outflow
! ****************************************************************************** !


! ****************************************************************************** !
  !> author: Kannan Masilamani
  !! Moment based velocity boundary condition from Sam Bennent PhD thesis
  !! "A Lattice Boltzmann Model for Diffusion of Binary Gas Mixtures"
  !!
  !! 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_inflow( 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, vel(3), dens_correction
    real(kind=rk) :: vel_b(globBC%nElems(iLevel)*3), &
      &              pressB(globBC%nElems(iLevel))
    integer :: iELem, iDir, tmpDir
    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) :: fTmp_2( layout%fStencil%QQ )
    real(kind=rk) :: rhoTmp, uxTmp(3), inv_vel, inv_press
    integer :: bcVel_pos, bcPress_pos
    ! ---------------------------------------------------------------------------
    write(logUnit(1),*) 'WARNING: Experimental BC moments_inflow'
    !write(*,*) 'Boundary label ', trim(me%label)
    QQ = layout%fStencil%QQ
    inv_vel = 1.0_rk / physics%fac( iLevel )%vel
    inv_press = 1.0_rk / physics%fac( iLevel )%press

    ! 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

    ! position of boundary pressure in varSys
    bcPress_pos = me%bc_states%pressure%varPos
    ! Get pressure
    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

    do iElem = 1, globBC%nElems( iLevel )

      write(*,*) 'iElem ', iElem
      vel(:) = vel_b((iElem-1)*3+1:iElem*3)
      write(*,*) 'velocity ', vel

      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 ) = &
!           & bcBuffer( ?IDX?( varPos(iDir), globBC%elemLvl( iLevel )%posInBcElemBuf%val(iElem), nScalars ))
            &state(?FETCH?( iDir, iField, elemPos, QQ, nScalars, nSize, neigh ))
        end if
!        write(*,*) 'fTmp0 ',  fTmp(iDir)
      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

write(*,*) 'dens before correction ', rho

      ! 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
write(*,*) 'dens after correction ', rho, rho*cs2
write(*,*) 'dens from press', pressB(iElem)*cs2inv, pressB(iElem)

      rho = pressB(iElem) * cs2inv + rho0

      !!! 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))
      !!! correction term for Multispecies
      !!! rho_new = rho_old/(1.0 - dens_correction)
      !!rho = rho + rho0*dens_correction

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

      ! moments for compressible model
      !moments = (/ rho, rho*ux(iElem), rho*uy(iElem), press+rho*ux(iElem)**2, &
      !  & press+rho*uy(iElem)**2, rho*ux(iElem)*uy(iElem), cs2*rho*uy(iElem), &
      !  & cs2*rho*ux(iElem), cs2*(press+rho*ux(iElem)**2+rho*uy(iElem)**2) /)
      ! moments for in-compressible model
      moments = (/ rho, rho0*vel(1), rho0*vel(2), pressB(iElem)+rho0*vel(1)**2,&
        & pressB(iElem)+rho0*vel(2)**2, rho0*vel(1)*vel(2), cs2*rho0*vel(2),   &
        & cs2*rho0*vel(1), cs2*(pressB(iElem)+rho0*vel(1)**2+rho0*vel(2)**2) /)
      !write(*,*) 'Moments ', moments
      !write(*,*) 'knownMomPos ', me%elemLvl(iLevel)%moments(iElem)%knownMom_pos

      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

    do iElem = 1, globBC%nElems( iLevel )
      !ux = 0.0_rk
      do iDir = 1,QQ
        fTmp_2(iDir) = state(                                                  &
&?FETCH?( iDir, iField, globBC%elemLvl(iLevel)%elem%val( iElem ), QQ, nScalars,nSize, neigh ) &
        & )
      end do
      rhoTmp = sum(fTmp_2)
      call derVarPos%velFromState( state  = fTmp_2, &
        &                          iField = iField, &
        &                          nElems = 1,      &
        &                          varSys = varSys, &
        &                          layout = layout, &
        &                          res    = uxTmp   )
      write(*,*) 'iElem ', iElem
      write(*,*) 'density ', rhoTmp, 'press ', cs2*rhoTmp, cs2*rhoTmp*physics%fac( iLevel )%press
      write(*,*) 'velocity ', uxTmp
!      write(*,*) 'velocity ', uxTmp*rhoTmp
    end do
!
!  stop

  end subroutine moments_inflow
! ****************************************************************************** !


! ****************************************************************************** !
  !> author: Kannan Masilamani
  !! Moment based velocity boundary condition from Sam Bennent PhD thesis
  !! "A Lattice Boltzmann Model for Diffusion of Binary Gas Mixtures"
  !!
  !! 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_outflow( 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, vel(3)
    integer :: iELem, iDir, QQ
    integer :: nLinks, iLink, 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 )
    ! Vel on frst neighbor element
    real(kind=rk) :: uxB_1(3*globBC%nElems(iLevel))
    ! density on first neighbor element
    real(kind=rk) :: rho_1(globBC%nElems(iLevel))
    ! ---------------------------------------------------------------------------
    write(logUnit(1),*) 'WARNING: Experimental BC moments_outflow'
    !write(*,*) 'Boundary label ', trim(me%label)

    QQ = layout%fStencil%QQ

    ! velocity is taken from precollision of current
    ! @todo: what layout is used inside this routine
    call derVarPos%velFromState(                            &
      &     state  = me%neigh(iLevel)%neighBufferPre(1, :), &
      &     iField = iField,                                &
      &     nElems = globBC%nElems(iLevel),                 &
      &     varSys = varSys,                                &
      &     layout = layout,                                &
      &     res    = uxB_1                                  )

    rho_1 = 0.0_rk
    ! density is taken from precollision of current
    do iElem = 1, globBC%nElems( iLevel )
      do iDir = 1,QQ
        rho_1 = rho_1 + me%neigh(iLevel)%neighBufferPre( 1, (iElem-1)*QQ+iDir )
      end do
    end do

    do iElem = 1, globBC%nElems( iLevel )

      !write(*,*) 'iElem ', iElem
      rho = rho_1(iElem)
      vel = (/ uxB_1((iElem-1)*3+1), uxB_1((iElem-1)*3+2), uxB_1((iElem-1)*3+3) /)
      !write(*,*) 'dens', rho, 'velocity ', vel

      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
!        write(*,*) 'fTmp0 ',  fTmp(iDir)
      end do
      fTmp( layout%fStencil%restPosition ) = &
& state(?FETCH?( layout%fstencil%restPosition, iField, elemPos, QQ, nScalars,nSize, neigh ) )

      press = rho * cs2

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

      ! moments for compressible model
      !moments = (/ rho, rho*ux(iElem), rho*uy(iElem), press+rho*ux(iElem)**2, &
      !  & press+rho*uy(iElem)**2, rho*ux(iElem)*uy(iElem), cs2*rho*uy(iElem), &
      !  & cs2*rho*ux(iElem), cs2*(press+rho*ux(iElem)**2+rho*uy(iElem)**2) /)
      ! moments for in-compressible model
      moments = (/ rho, rho0*vel(1), rho0*vel(2), press+rho0*vel(1)**2, &
        & press+rho0*vel(2)**2, rho0*vel(1)*vel(2), cs2*rho0*vel(2), &
        & cs2*rho0*vel(1), cs2*(press+rho0*vel(1)**2+rho0*vel(2)**2) /)
      !write(*,*) 'Moments ', moments
      !write(*,*) 'knownMomPos ', me%elemLvl(iLevel)%moments(iElem)%knownMom_pos

      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

!!    do iElem = 1, globBC%nElems( iLevel )
!!      !ux = 0.0_rk
!!      do iDir = 1,QQ
!!        fTmp_2(iDir) = state(                                                  &
!!&?FETCH?( iDir, iField, globBC%elemLvl(iLevel)%elem%val( iElem ), QQ, nScalars, neigh,nElems ) &
!!        & )
!!      end do
!!      rhoTmp = sum(fTmp_2)
!!      call derVarPos%velFromState( state  = fTmp_2, &
!!        &                          iField = iField, &
!!        &                          nElems = 1,      &
!!        &                          varSys = varSys, &
!!        &                          layout = layout, &
!!        &                          res    = uxTmp   )
!!      write(*,*) 'iElem ', iElem
!!      write(*,*) 'density ', rhoTmp, cs2*rhoTmp
!!!      write(*,*) 'velocity ', uxTmp
!!!      write(*,*) 'velocity ', uxTmp*rhoTmp
!!    end do
!
!  stop

  end subroutine moments_outflow
! ****************************************************************************** !
?? ENDIF

! ****************************************************************************** !
  !> 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 pressure_expol_slow( 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, ii, ff
    logical :: axisNormal
    ! ---------------------------------------------------------------------------

    ! 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   = ff,                    &
        &                            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 ii = 1, QQ * 4
    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)     )
      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

      end if
      state( me%links(iLevel)%val(iLink) ) = &
        & state( me%links(iLevel)%val(iLink) ) + dble(QQ*ii - nScalars*ii)
    end do ! iLink, linkwise treatment ends here
end do

    ! 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

  end subroutine pressure_expol_slow
! ****************************************************************************** !

end module mus_bc_fluid_experimental_module
! ****************************************************************************** !
