! Copyright (c) 2012-2020 Kannan Masilamani <kannan.masilamani@uni-siegen.de>
! Copyright (c) 2012-2014 Simon Zimny <s.zimny@grs-sim.de>
! Copyright (c) 2012-2013 Manuel Hasert <m.hasert@grs-sim.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 Raphael Haupt <raphael.haupt@uni-siegen.de>
! Copyright (c) 2019-2020 Peter Vitt <peter.vitt2@uni-siegen.de>
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions are met:
!
! 1. Redistributions of source code must retain the above copyright notice,
! this list of conditions and the following disclaimer.
!
! 2. Redistributions in binary form must reproduce the above copyright notice,
! this list of conditions and the following disclaimer in the documentation
! and/or other materials provided with the distribution.
!
! THIS SOFTWARE IS PROVIDED BY THE UNIVERSITY OF SIEGEN “AS IS” AND ANY EXPRESS
! OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
! OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
! IN NO EVENT SHALL UNIVERSITY OF SIEGEN OR CONTRIBUTORS BE LIABLE FOR ANY
! DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
! (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
! ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
?? include 'header/lbm_macros.inc'
! ****************************************************************************** !
!> Boundary condition treatment routines for multispecies simulation
!!
module mus_bc_species_module
?? IF( .not. SOA) then
  use iso_c_binding, only: c_f_pointer

  ! include treelm modules
  use env_module,               only: rk, labelLen
  use tem_param_module,         only: div1_3, div1_36, div1_8, div3_4h,        &
    &                                 div1_4, div3_8, div3_2, div9_16, div3_16,&
    &                                 cs2inv, div1_6, cs4inv,                  &
    &                                 t2cs2inv, t2cs4inv, div1_18, cs2
  use tem_time_module,          only: tem_time_type
  use treelmesh_module,         only: treelmesh_type
  use tem_varSys_module,        only: tem_varSys_type
  use tem_math_module,          only: invert_matrix
  use tem_spatial_module,       only: tem_spatial_for
  use tem_spacetime_fun_module, only: tem_spacetime_for
  use tem_isNaN_module,         only: tem_isNaN
  use tem_debug_module,         only: dbgUnit
  use tem_property_module,      only: prp_solid
  use tem_construction_module,  only: tem_levelDesc_type
  use tem_stencil_module,       only: tem_stencilHeader_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_species_module,         only: mus_species_type
  use mus_eNRTL_module,           only: mus_calc_thermFactor,                  &
    &                                   mus_calc_MS_DiffMatrix
  use mus_varSys_module,          only: mus_varSys_data_type
  use mus_derQuanMSLiquid_module, only: momentumFromMacroLSE


  implicit none

  private

  public :: spc_inlet_eq, spc_inlet, spc_vel_bb

  public :: spc_outlet_zero_prsgrd
  public :: spc_outlet_eq, spc_outlet_vel, spc_outlet_expol
  public :: spc_moleFrac, spc_moleDiff_Flux
  public :: spc_moleFlux, spc_moleFlux_eq
  public :: spc_moleDens_eq
  public :: spc_blackbox_mem_ion, spc_blackbox_mem_solvent
  public :: spc_moleFrac_wtdf, spc_moleFrac_eq
  public :: spc_outflow, spc_inflow
  public :: spc_solvent_outflow, spc_solvent_inflow
  public :: spc_velocity_noneq_expol, spc_mole_fraction_noneq_expol
  ! moments BC
  public :: spc_moments_moleFrac
  public :: spc_moments_moleFlux
  public :: spc_moments_wall
  public :: spc_moments_vel

contains
! ****************************************************************************** !
  !> author: Kannan Masilamani
  !! Outlet boundary conditions with zero pressure gradient.
  !!
  !! These boundary conditions use the neighbor cells density and velocity in
  !! the equilibrium function. It is not necessary to specify density at
  !! boundary in the lua configuration file
  !! \( f = f^{eq}(\rho_{b-1},u_{b-1}) \)
  !!
  !! 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_outlet_zero_prsgrd( me, state, bcBuffer, globBC, levelDesc, &
    &                                tree, nSize, iLevel, sim_time, neigh,   &
    &                                layout, fieldProp, varPos, nScalars,    &
    &                                varSys, derVarPos, physics, iField,     &
    &                                mixture                                 )
    ! -------------------------------------------------------------------- !
    !> global boundary type
    class(boundary_type) :: me
    !> Current state vector of iLevel
    real(kind=rk), intent(inout) :: state(:)
    !> size of state array ( in terms of elements )
    integer, intent(in) :: nSize
    !> state values of boundary elements of all fields of iLevel
    real(kind=rk), intent(in) :: bcBuffer(:)
    !> iLevel descriptor
    type(tem_levelDesc_type), intent(in) :: levelDesc
    !> Treelm Mesh
    type(treelmesh_type), intent(in) :: tree
    !> fluid parameters and properties
    type(mus_field_prop_type), intent(in) :: fieldProp
    !> stencil layout information
    type(mus_scheme_layout_type), intent(in) :: layout
    !> the level On which this boundary was invoked
    integer, intent(in) :: iLevel
    !> connectivity array corresponding to state vector
    integer, intent(in) :: neigh(:)
    !> global time information
    type(tem_time_type), intent(in)  :: sim_time
    !> pointer to field variable in the state vector
    integer, intent(in) :: varPos(:)
    !> number of Scalars in the scheme var system
    integer, intent(in) :: nScalars
    !> scheme variable system
    type(tem_varSys_type), intent(in) :: varSys
    !> position of derived quantities in varsys
    type(mus_derVarPos_type), intent(in) :: derVarPos
    !> scheme global boundary type
    type(glob_boundary_type), intent(in) :: globBC
    !> scheme global boundary type
    type(mus_physics_type), intent(in) :: physics
    !> current field
    integer, intent(in) :: iField
    !> mixture info
    type(mus_mixture_type), intent(in) :: mixture
    ! -------------------------------------------------------------------- !
    real(kind=rk) :: rho, ux(3), usq, m_ratio
    integer :: iElem, iDir, iField_2, pos, QQ
    real(kind=rk) :: fTmp( layout%fStencil%QQ * varSys%nStateVars )
    integer :: elemPos, posInBuffer
    ! ---------------------------------------------------------------------------
    QQ = layout%fStencil%QQ

    m_ratio = fieldProp%species%molWeigRatio

    do iElem = 1, globBC%nElems( iLevel )
      posInBuffer = globBC%elemLvl( iLevel )%posInBcElemBuf%val(iElem)
      do iField_2 = 1, varSys%nStateVars
        do iDir = 1,layout%fStencil%QQ
          pos = varSys%method%val(iField_2)%state_varPos(iDir)
          fTmp(pos) = bcBuffer( &
          & ?IDX?( pos, posInBuffer, nScalars,globBC%nElems(iLevel) ) )
        end do
      end do

      !Calculate local density and velocity moments
      rho = 0._rk
      do iDir = 1,layout%fStencil%QQ
        rho = rho + fTmp(varPos(iDir))
      end do

      call derVarPos%velFromState( state  = fTmp,   &
        &                          iField = iField, &
        &                          nElems = 1,      &
        &                          varSys = varSys, &
        &                          layout = layout, &
        &                          res    = ux      )

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

      elemPos = globBC%elemLvl(iLevel)%elem%val( iElem )
      do iDir = 1,layout%fStencil%QQN
        ! Write the values
        if( globBC%elemLvl(iLevel)%bitmask%val(iDir, iElem )) then
          ! Only changed from save to fetch and added currT
          state(?FETCH?(iDir, iField, elemPos, QQ, nScalars, nSize, neigh)) =  &
            &       layout%weight( iDir )*rho*(2.0_rk*m_ratio+9._rk*( (        &
            &       layout%fStencil%cxDir( 1, layout%fStencil                  &
            &                                       %cxDirInv( iDir ))*ux(1)   &
            &    +  layout%fStencil%cxDir( 2, layout%fStencil                  &
            &                                       %cxDirInv( iDir ))*ux(2)   &
            &    +  layout%fStencil%cxDir( 3, layout%fStencil                  &
            &                                       %cxDirInv( iDir ))*ux(3)   &
            &    )**2 -  div1_3*usq))                                          &
            &    - state(                                                      &
            & ?FETCH?( iDir, iField, elemPos, QQ, nScalars, nSize, neigh ))
        end if
      end do
    end do

  end subroutine spc_outlet_zero_prsgrd
! ****************************************************************************** !

! ****************************************************************************** !
  !> Inlet boundary condition for defined species velocity and mole fraction
  !! Usage
  !! -----
  !!```lua
  !!boundary_condition = {
  !! { label = 'outlet',
  !!    kind = 'spc_inlet',
  !!    velocity = 'inlet_vel',
  !!    mole_fraction = 'inlet_mole'
  !! }
  !!}
  !!variable = {
  !! {
  !!   label = 'inlet_vel',
  !!   ncomponents = 3,
  !!   vartype = 'st_fun',
  !!   st_fun = {0.1,0.0,0.0}
  !! },
  !! {
  !!   label = 'inlet_mole',
  !!   ncomponents = 1,
  !!   vartype = 'st_fun',
  !!   st_fun = 0.1
  !! }
  !!}
  !!```
  !!
  !! This subroutine's interface must match the abstract interface definition
  !! [[boundaryRoutine]] in bc/[[mus_bc_header_module]].f90 in order to be
  !! callable via [[boundary_type:fnct]] function pointer.
  subroutine spc_inlet( 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) :: mass_dens
    real(kind=rk) :: vel_b(globBC%nElems(iLevel)*3), inv_vel
    real(kind=rk) :: moleFrac(globBC%nElems(iLevel))
    real(kind=rk) :: massFlux(3)
    integer :: iElem, iDir, QQ, posInBuffer, bcVel_pos, offset, elemPos
    integer :: bcMoleFrac_pos
    ! ------------------------------------------------------------------------
    QQ = layout%fStencil%QQ
    inv_vel = 1.0_rk / physics%fac( iLevel )%vel

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

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

    bcMoleFrac_pos = me%bc_states%moleFrac%varPos
    ! mole fraction
    call varSys%method%val(bcMoleFrac_pos)%get_valOfIndex( &
      & varSys  = varSys,                                  &
      & time    = sim_time,                                &
      & iLevel  = iLevel,                                  &
      & idx     = me%bc_states%moleFrac                    &
      &           %pntIndex%indexLvl(iLevel)               &
      &           %val(1:globBC%nElems(iLevel)),           &
      & nVals   = globBC%nElems(iLevel),                   &
      & res     = moleFrac(:)                              )

    do iElem = 1, globBC%nElems(iLevel)
      offset = (iElem-1)*3

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

      ! local density
      !mass_dens = sum(fTmp)
      mass_dens = mixture%moleDens0LB *  moleFrac(iElem) &
        &       * fieldProp%species%molWeight

      massFlux = mass_dens * vel_b(offset+1:offset+3)

      elemPos = globBC%elemLvl( iLevel )%elem%val( iElem )
      do iDir = 1, layout%fStencil%QQN
        if( globBC%elemLvl(iLevel)%bitmask%val( iDir, iElem )) then
          ! Depending on PUSH or pull, use + or - for cxDir, because directions
          ! are inverted
          state(?FETCH?( iDir, iField, elemPos, QQ, nScalars, nSize, neigh )) &
              & = fTmp(layout%fStencil%cxDirInv( iDir ))                      &
              & + layout%weight( iDir )*2._rk*cs2inv                          &
              &    * ( layout%fStencil%cxDirRK( 1, iDir )*massFlux(1)         &
              &    +   layout%fStencil%cxDirRK( 2, iDir )*massFlux(2)         &
              &    +   layout%fStencil%cxDirRK( 3, iDir )*massFlux(3)         )
        end if
      end do
    end do
  end subroutine spc_inlet
! ****************************************************************************** !


! ****************************************************************************** !
  !> Inlet species velocity equilibrium boundary with specified
  !! mixture averaged mass velocity and its molefraction
  !! mixture kinematic pressure is extrapolated here.
  !! Density and velocity of all fields are used to compute equilibrium
  !! Usage
  !! -----
  !!```lua
  !!boundary_condition = {
  !! { label = 'outlet',
  !!    kind = 'spc_inlet_eq',
  !!    velocity = 'inlet_vel',
  !!    mole_fraction = 'inlet_mole'
  !! }
  !!}
  !!variable = {
  !! {
  !!   label = 'inlet_vel',
  !!   ncomponents = 3,
  !!   vartype = 'st_fun',
  !!   st_fun = {0.1,0.0,0.0}
  !! },
  !! {
  !!   label = 'inlet_mole',
  !!   ncomponents = 1,
  !!   vartype = 'st_fun',
  !!   st_fun = 0.1
  !! }
  !!}
  !!```
  !!
  !! This subroutine's interface must match the abstract interface definition
  !! [[boundaryRoutine]] in bc/[[mus_bc_header_module]].f90 in order to be
  !! callable via [[boundary_type:fnct]] function pointer.
  subroutine spc_inlet_eq( me, state, bcBuffer, globBC, levelDesc, tree,      &
    &                      nSize, iLevel, sim_time, neigh, layout, fieldProp, &
    &                      varPos, nScalars, varSys, derVarPos, physics,      &
    &                      iField, mixture                                    )
    ! -------------------------------------------------------------------- !
    !> global boundary type
    class(boundary_type) :: me
    !> Current state vector of iLevel
    real(kind=rk), intent(inout) :: state(:)
    !> size of state array ( in terms of elements )
    integer, intent(in) :: nSize
    !> state values of boundary elements of all fields of iLevel
    real(kind=rk), intent(in) :: bcBuffer(:)
    !> iLevel descriptor
    type(tem_levelDesc_type), intent(in) :: levelDesc
    !> Treelm Mesh
    type(treelmesh_type), intent(in) :: tree
    !> fluid parameters and properties
    type(mus_field_prop_type), intent(in) :: fieldProp
    !> stencil layout information
    type(mus_scheme_layout_type), intent(in) :: layout
    !> the level On which this boundary was invoked
    integer, intent(in) :: iLevel
    !> connectivity array corresponding to state vector
    integer, intent(in) :: neigh(:)
    !> global time information
    type(tem_time_type), intent(in)  :: sim_time
    !> pointer to field variable in the state vector
    integer, intent(in) :: varPos(:)
    !> number of Scalars in the scheme var system
    integer, intent(in) :: nScalars
    !> scheme variable system
    type(tem_varSys_type), intent(in) :: varSys
    !> position of derived quantities in varsys
    type(mus_derVarPos_type), intent(in) :: derVarPos
    !> scheme global boundary type
    type(glob_boundary_type), intent(in) :: globBC
    !> scheme global boundary type
    type(mus_physics_type), intent(in) :: physics
    !> current field
    integer, intent(in) :: iField
    !> mixture info
    type(mus_mixture_type), intent(in) :: mixture
    ! -------------------------------------------------------------------- !
    real(kind=rk) :: fTmp( layout%fStencil%QQ*globBC%nElems(iLevel) &
      &              * varSys%nStateVars )
    ! fEq uses AOS layout
    real(kind=rk) :: fEq( layout%fStencil%QQ*globBC%nElems(iLevel) )
    real(kind=rk) :: mass_dens( globBC%nElems(iLevel)*varSys%nStateVars )
    !real(kind=rk) :: uxB(3*globBC%nElems(iLevel))   !< Velocity on boundary element
    !> Velocity on boundary element
    real(kind=rk) :: uxB(globBC%nElems(iLevel)*varSys%nStateVars*3)
    real(kind=rk) :: velocity(3,globBC%nElems(iLevel)*varSys%nStateVars)
    real(kind=rk) :: spc_vel(globBC%nElems(iLevel)*3), inv_vel
    real(kind=rk) :: moleFrac(globBC%nElems(iLevel))
    integer :: iElem, iDir, iFieldLoc, nFields, pos, QQ
    integer :: offset, bcVel_pos, bcMoleFrac_pos, elemPos, posInBuffer
    ! ------------------------------------------------------------------------
    QQ = layout%fStencil%QQ
    nFields = varSys%nStateVars
    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     = spc_vel                             )

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

    ! position of molefraction spacetime variable in varSys
    bcMoleFrac_pos = me%bc_states%moleFrac%varPos
    ! mole fraction
    call varSys%method%val(bcMoleFrac_pos)%get_valOfIndex( &
      & varSys  = varSys,                                  &
      & time    = sim_time,                                &
      & iLevel  = iLevel,                                  &
      & idx     = me%bc_states%moleFrac                    &
      &           %pntIndex%indexLvl(iLevel)               &
      &           %val(1:globBC%nElems(iLevel)),           &
      & nVals   = globBC%nElems(iLevel),                   &
      & res     = moleFrac                                 )

    ! copy state
    do iElem = 1, globBC%nElems(iLevel)
      posInBuffer = globBC%elemLvl( iLevel )%posInBcElemBuf%val(iElem)
      do iFieldLoc = 1, nFields
        do iDir = 1,layout%fStencil%QQ
          pos = varSys%method%val(iFieldLoc)%state_varPos(iDir)
          fTmp( ?IDX?( pos, iElem, nScalars,globBC%nElems(iLevel) ) ) = &
            & bcBuffer(                                                 &
            & ?IDX?( pos, posInBuffer, nScalars,globBC%nElems(iLevel) ))
        end do
      end do !iField
    end do !iElem

    ! Derive all species velocities at once since
    ! it is inefficient to derive each species velocity
    call derVarPos%velocitiesFromState( state  = fTmp,                  &
      &                                 iField = iField,                &
      &                                 nElems = globBC%nElems(iLevel), &
      &                                 varSys = varSys,                &
      &                                 layout = layout,                &
      &                                 res    = uxB                    )

    ! get density and velocity.
    ! if current field, use velocity and density defined in lua file.
    ! for other fields, derive density and velocity from state
    mass_dens = 0.0_rk
    do iFieldLoc = 1, nFields
      if (iFieldLoc == iField) then
        ! store velocity in input_loc array
        do iElem = 1, globBC%nElems(iLevel)
          offset = (iElem-1)*nFields + iFieldLoc
          velocity(:,offset) = spc_vel((iElem-1)*3+1 : iElem*3)

          ! compute current species mass density from specified molefraction at
          ! boundary
          ! rho = n_t * chi_i * m_i + massFrac_i * rho0 * KinePress/(cs2*phi_i)
          ! massFrac_i = rho_i/rho
          ! rho = n_t * chi_i * m_i + rho_i * KinePress/(cs2*phi_i)
          ! 1st term is zero order density
          ! 2nd term is second order density with kinematic mixture pressure,
          ! p = cs2*(sum(phi_k*rho_k) - min_a(m_a)*n0)/rho0
          !KM: using inital mixture number density rho0/mixtureMOlWeight
          !Using local tot_NuMdens increases density over time
          mass_dens( offset ) = mixture%moleDens0LB * moleFrac(iElem) &
            &                 * fieldProp%species%molWeight
        end do
      else
        do iElem = 1, globBC%nElems(iLevel)
          offset = (iElem-1)*nFields + iFieldLoc
          velocity(:,offset) =                             &
            & uxB((iElem-1)*nFields*3 + (iFieldLoc-1)*3 + 1 : &
            &     (iElem-1)*nFields*3 + iFieldLoc*3 )
          ! species density
          do iDir = 1, layout%fStencil%QQ
            pos = varSys%method%val(iFieldLoc)%state_varPos(iDir)
            mass_dens(offset) = mass_dens(offset) &
              & + fTmp(?IDX?( pos, iElem, nScalars,globBC%nElems(iLevel)) )
          end do !iDir
        end do !iElem
      end if
    end do !iField

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

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

  end subroutine spc_inlet_eq
! ****************************************************************************** !

! ****************************************************************************** !
  !> Inlet species velocity bounce back boundary with specified
  !! mixture averaged mass velocity and its molefraction
  !! mixture kinematic pressure is extrapolated here.
  !! Usage
  !! -----
  !!```lua
  !!boundary_condition = {
  !! { label = 'outlet',
  !!    kind = 'spc_vel_bb',
  !!    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_vel_bb( 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) :: mass_dens
    real(kind=rk) :: vel_b(globBC%nElems(iLevel)*3), inv_vel
    integer :: iElem, iDir, QQ, posInBuffer, bcVel_pos, offset, elemPos
    ! ------------------------------------------------------------------------
    QQ = layout%fStencil%QQ
    inv_vel = 1.0_rk / physics%fac( iLevel )%vel

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

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

    do iElem = 1, globBC%nElems(iLevel)
      offset = (iElem-1)*3

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

      elemPos = globBC%elemLvl( iLevel )%elem%val( iElem )
      do iDir = 1,layout%fStencil%QQN
        if( globBC%elemLvl(iLevel)%bitmask%val( iDir, iElem )) then
          ! Depending on PUSH or pull, use + or - for cxDir, because directions
          ! are inverted
          state(?FETCH?( iDir, iField, elemPos, QQ, nScalars, nSize, neigh )) &
              & = fTmp(layout%fStencil%cxDirInv( iDir ))                      &
              & + layout%weight( iDir )*2._rk*cs2inv*mass_dens                &
              &    * ( layout%fStencil%cxDirRK( 1, iDir )*vel_b(offset+1)     &
              &    +   layout%fStencil%cxDirRK( 2, iDir )*vel_b(offset+2)     &
              &    +   layout%fStencil%cxDirRK( 3, iDir )*vel_b(offset+3)     )
        end if
      end do
    end do
  end subroutine spc_vel_bb
! ****************************************************************************** !

  ! ************************************************************************ !
  !> species Outlet Pressure extrapolation boundary. NOT VERIFIED
  !!
  !! This is taken from the paper:
  !! M. Junk and Z. Yang,
  !! Asymptotic Analysis of Lattice Boltzmann Outflow Treatments
  !! Commun. Comput. Phys., pp. 1–11, 2011.
  !!
  !! * Pressure as defined in the configuration file
  !! Usage
  !! -----
  !!```lua
  !!boundary_condition = {
  !! { label = 'outlet',
  !!    kind = 'spc_outlet_expol',
  !!    pressure = 'p1'
  !!  }
  !!}
  !!variable = {
  !!  name = 'p1',
  !!  ncomponents = 1,
  !!  vartype = 'st_fun',
  !!  st_fun = 1.0
  !!}
  !!```
  !!
  !! This subroutine's interface must match the abstract interface definition
  !! [[boundaryRoutine]] in bc/[[mus_bc_header_module]].f90 in order to be
  !! callable via [[boundary_type:fnct]] function pointer.
  subroutine spc_outlet_expol( me, state, bcBuffer, globBC, levelDesc, tree,   &
    &                          nSize, iLevel, sim_time, neigh, layout,         &
    &                          fieldProp, varPos, nScalars, varSys, derVarPos, &
    &                          physics, iField, mixture                        )
    ! -------------------------------------------------------------------- !
    !> global boundary type
    class(boundary_type) :: me
    !> Current state vector of iLevel
    real(kind=rk), intent(inout) :: state(:)
    !> size of state array ( in terms of elements )
    integer, intent(in) :: nSize
    !> state values of boundary elements of all fields of iLevel
    real(kind=rk), intent(in) :: bcBuffer(:)
    !> iLevel descriptor
    type(tem_levelDesc_type), intent(in) :: levelDesc
    !> Treelm Mesh
    type(treelmesh_type), intent(in) :: tree
    !> fluid parameters and properties
    type(mus_field_prop_type), intent(in) :: fieldProp
    !> stencil layout information
    type(mus_scheme_layout_type), intent(in) :: layout
    !> the level On which this boundary was invoked
    integer, intent(in) :: iLevel
    !> connectivity array corresponding to state vector
    integer, intent(in) :: neigh(:)
    !> global time information
    type(tem_time_type), intent(in)  :: sim_time
    !> pointer to field variable in the state vector
    integer, intent(in) :: varPos(:)
    !> number of Scalars in the scheme var system
    integer, intent(in) :: nScalars
    !> scheme variable system
    type(tem_varSys_type), intent(in) :: varSys
    !> position of derived quantities in varsys
    type(mus_derVarPos_type), intent(in) :: derVarPos
    !> scheme global boundary type
    type(glob_boundary_type), intent(in) :: globBC
    !> scheme global boundary type
    type(mus_physics_type), intent(in) :: physics
    !> current field
    integer, intent(in) :: iField
    !> mixture info
    type(mus_mixture_type), intent(in) :: mixture
    ! -------------------------------------------------------------------- !
    real(kind=rk) :: fTmpAll( layout%fStencil%QQ*globBC%nElems(iLevel) &
      &              * varSys%nStateVars )
    real(kind=rk) :: fTmp_1
    real(kind=rk) :: fTmp_2
    ! fEq uses AOS layout
    real(kind=rk) :: fEq( layout%fStencil%QQ*globBC%nElems(iLevel) )
    integer :: iDir              !< Direction counter
    integer :: iElem             !< Element counter
    integer :: pos, ifieldloc, nFields, neighPos, QQ, elemPos
    ! ------------------------------------------------------------------------

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

    ! Run over all elements in list for this boundary
    do iElem = 1, globBC%nElems(iLevel)
      neighPos = me%neigh(iLevel)%posInState(1,iElem)
      do iFieldLoc = 1, nFields
        do iDir = 1,layout%fStencil%QQ
          pos = varSys%method%val(iFieldLoc)%state_varPos(iDir)
          !neighbor element
          fTmpAll( ?IDX?( pos, iElem, nScalars,globBC%nElems(iLevel) ) )     &
            & = state(?SAVE?( iDir,iFieldLoc,neighPos,QQ,nScalars,nSize,neigh ))
          !local element
          !fTmpAll( ?IDX?( pos, iElem, nScalars,globBC%nElems(iLevel) ) ) = bcBuffer(                    &
          !  & ?IDX?( pos, globBC%elemLvl( iLevel )%posInBcElemBuf%val( iElem ), nScalars,globBC%nElems(iLevel) ) )
        end do
      end do !iField
    end do !iElem

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


    ! Run over all elements in list for this boundary
    do iElem = 1, globBC%nElems(iLevel)
      elemPos = globBC%elemLvl( iLevel )%elem%val( iElem )
      ! Treat all directions, but actually we only want to treat
      ! non-normal directions
      ! Equation (3.2b)
      do iDir = 1, layout%fStencil%QQN
        if( globBC%elemLvl(iLevel)%bitmask%val( iDir, iElem )) then
        ! KM: pre-collision of current time step must be used for extrapolation
          fTmp_1 = me%neigh( iLevel )%neighBufferPre( 1,  iDir+(iElem-1)*QQ)
          fTmp_2 = me%neigh( iLevel )%neighBufferPre( 2,  iDir+(iElem-1)*QQ)

          ! update the incoming velocity direction
           state(?FETCH?( iDir, iField, elemPos, QQ, nScalars,nSize,neigh )) &
            & = 1.5_rk*fTmp_1 - 0.5_rk*fTmp_2
            !& = feq( iDir+(iElem-1)*QQ )
        end if
      end do

      ! then overwrite the normal direction with special treatment
      ! Equation (3.2a)
      iDir = globBC%elemLvl( iLevel )%normalInd%val( iElem )
state( ?FETCH?( iDir, iField, elemPos, QQ, nScalars, nSize,neigh ))  = &
        &  feq( ?IDX?(iDir, iElem, QQ) ) !&
        !& +(fTmpAll(?IDX?(varPos(layout%fStencil%cxDirInv(iDir)), iElem, nScalars))&
        !& -feq( ?IDX?(layout%fStencil%cxDirInv(iDir), iElem, layout%fStencil%QQ) ) )
        !& fTmpAll(?IDX?(varPos(iDir), iElem, nScalars))
    end do !iElem

  end subroutine spc_outlet_expol
! ****************************************************************************** !


! ****************************************************************************** !
  !> Outlet species velocity equilibrium boundary with specified
  !! mixture averaged mass velocity.
  !! molefraction is extrapolated here.
  !! Usage
  !! -----
  !!```lua
  !!boundary_condition = {
  !! { label = 'outlet',
  !!    kind = 'spc_outlet_vel',
  !!    velocity = 'inlet_vel',
  !! }
  !!}
  !!variable = {
  !! {
  !!   label = 'inlet_vel',
  !!   ncomponents = 3,
  !!   vartype = 'st_fun',
  !!   st_fun = {0.1,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_outlet_vel( 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*globBC%nElems(iLevel) &
      &              * varSys%nStateVars )
    ! fEq uses AOS layout
    real(kind=rk) :: fEq( layout%fStencil%QQ*globBC%nElems(iLevel) )
    real(kind=rk) :: mass_dens( globBC%nElems(iLevel)*varSys%nStateVars )
    !> Velocity on boundary element
    real(kind=rk) :: uxB(globBC%nElems(iLevel)*varSys%nStateVars*3)
    real(kind=rk) :: velocity(3,globBC%nElems(iLevel)*varSys%nStateVars)
    real(kind=rk) :: spc_vel(globBC%nElems(iLevel)*3), inv_vel
    integer :: iElem, iDir, iFieldLoc, nFields, pos, QQ
    integer :: offset, bcVel_pos, elemPos, posInBuffer
    ! ------------------------------------------------------------------------
    QQ = layout%fStencil%QQ
    nFields = varSys%nStateVars
    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     = spc_vel                             )

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

    ! copy state
    do iElem = 1, globBC%nElems(iLevel)
      posInBuffer = globBC%elemLvl( iLevel )%posInBcElemBuf%val( iElem )
      do iFieldLoc = 1, nFields
        do iDir = 1,layout%fStencil%QQ
          pos = varSys%method%val(iFieldLoc)%state_varPos(iDir)
          fTmp( ?IDX?( pos, iElem, nScalars,globBC%nElems(iLevel) ) )   &
            & = bcBuffer(                                               &
            & ?IDX?( pos, posInBuffer, nScalars,globBC%nElems(iLevel) ) )
        end do
      end do !iField
    end do !iElem

    ! Derive all species velocities at once since
    ! it is inefficient to derive each species velocity
    call derVarPos%velocitiesFromState( state  = fTmp,                  &
      &                                 iField = iField,                &
      &                                 nElems = globBC%nElems(iLevel), &
      &                                 varSys = varSys,                &
      &                                 layout = layout,                &
      &                                 res    = uxB                    )

    ! use same velocity for all species and extrapolate density
    do iFieldLoc = 1, nFields
        do iElem = 1, globBC%nElems(iLevel)
          offset = (iElem-1)*nFields + iFieldLoc
          velocity(:,offset) = spc_vel((iElem-1)*3+1 : iElem*3)
        end do
    end do !iField

    ! species density
    mass_dens = 0.0_rk
    do iElem = 1, globBC%nElems(iLevel)
      do iFieldLoc = 1, nFields
        offset = (iElem-1)*nFields + iFieldLoc
        do iDir = 1, layout%fStencil%QQ
          pos = varSys%method%val(iFieldLoc)%state_varPos(iDir)
          mass_dens(offset) = mass_dens(offset) &
            & + fTmp(?IDX?( pos, iElem, nScalars,globBC%nElems(iLevel)) )
        end do !iDir
      end do !iField
    end do !iElem

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

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

  end subroutine spc_outlet_vel
! ****************************************************************************** !


! ****************************************************************************** !
  !> Outlet mixture pressure species equilibrium boundary
  !! kinematic pressure is computed from pressure
  !! species density and velocity are extrapolated
  !! Usage
  !! -----
  !!```lua
  !!boundary_condition = {
  !! { label = 'outlet',
  !!    kind = 'spc_outlet_eq',
  !!    pressure = 'press0'
  !! }
  !!}
  !!variable = {
  !! {
  !!   label = 'press0',
  !!   ncomponents = 1,
  !!   vartype = 'st_fun',
  !!   st_fun = 1.0
  !! }
  !!}
  !!```
  !!
  !! This subroutine's interface must match the abstract interface definition
  !! [[boundaryRoutine]] in bc/[[mus_bc_header_module]].f90 in order to be
  !! callable via [[boundary_type:fnct]] function pointer.
  subroutine spc_outlet_eq( me, state, bcBuffer, globBC, levelDesc, tree,      &
    &                       nSize, iLevel, sim_time, neigh, layout, fieldProp, &
    &                       varPos, nScalars, varSys, derVarPos, physics,      &
    &                       iField, mixture                                    )
    ! -------------------------------------------------------------------- !
    !> global boundary type
    class(boundary_type) :: me
    !> Current state vector of iLevel
    real(kind=rk), intent(inout) :: state(:)
    !> size of state array ( in terms of elements )
    integer, intent(in) :: nSize
    !> state values of boundary elements of all fields of iLevel
    real(kind=rk), intent(in) :: bcBuffer(:)
    !> iLevel descriptor
    type(tem_levelDesc_type), intent(in) :: levelDesc
    !> Treelm Mesh
    type(treelmesh_type), intent(in) :: tree
    !> fluid parameters and properties
    type(mus_field_prop_type), intent(in) :: fieldProp
    !> stencil layout information
    type(mus_scheme_layout_type), intent(in) :: layout
    !> the level On which this boundary was invoked
    integer, intent(in) :: iLevel
    !> connectivity array corresponding to state vector
    integer, intent(in) :: neigh(:)
    !> global time information
    type(tem_time_type), intent(in)  :: sim_time
    !> pointer to field variable in the state vector
    integer, intent(in) :: varPos(:)
    !> number of Scalars in the scheme var system
    integer, intent(in) :: nScalars
    !> scheme variable system
    type(tem_varSys_type), intent(in) :: varSys
    !> position of derived quantities in varsys
    type(mus_derVarPos_type), intent(in) :: derVarPos
    !> scheme global boundary type
    type(glob_boundary_type), intent(in) :: globBC
    !> scheme global boundary type
    type(mus_physics_type), intent(in) :: physics
    !> current field
    integer, intent(in) :: iField
    !> mixture info
    type(mus_mixture_type), intent(in) :: mixture
    ! -------------------------------------------------------------------- !
    ! fEq uses AOS layout
    real(kind=rk) :: fEq( layout%fStencil%QQ*globBC%nElems(iLevel) )
    !store pdf of all fields for derive function pointer input
    real(kind=rk) :: fTmp( layout%fStencil%QQ*globBC%nElems(iLevel) &
      &              * varSys%nStateVars )
    integer :: iDir              ! Direction counter
    integer :: iElem             ! Element counter
    integer :: nFields, iFieldLoc, pos, offset, QQ!, bcPress_pos
!    real(kind=rk) :: press(globBC%nElems(iLevel))
    integer :: elemPos, posInBuffer, neighPos
    ! ---------------------------------------------------------------------------
    QQ = layout%fStencil%QQ
    nFields = varSys%nStateVars

    ! @todo KM: use pressure term to set species density at boundary
    ! get density value from the definition in the lua file
    ! position of boundary pressure in varSys
!KM!    bcPress_pos = me%bc_states%pressure%varPos
!KM!    ! get pressure variable from spacetime function
!KM!    call varSys%method%val(bcPress_pos)%get_valOfIndex( &
!KM!      & varSys  = varSys,                               &
!KM!      & time    = sim_time,                             &
!KM!      & iLevel  = iLevel,                               &
!KM!      & idx     = me%bc_states%pressure                 &
!KM!      &           %pntIndex%indexLvl(iLevel)            &
!KM!      &           %val(1:globBC%nElems(iLevel)),        &
!KM!      & nVals   = globBC%nElems(iLevel),                &
!KM!      & res     = press                                 )
!KM!
!KM!    ! If physical quantities are given, transform to lattice units by division
!KM!    ! with the conversion factor
!KM!    ! kinematic pressure = pressure/rho0
!KM!    press = (press/mixture%rho0)/(physics%fac( iLevel )%press/physics%rho0)

    ! Calculate the density of current element
    do iElem = 1, globBC%nElems(iLevel)
      posInBuffer = globBC%elemLvl( iLevel )%posInBcElemBuf%val( iElem )
      neighPos = me%neigh(iLevel)%posInState(1,iElem)
      do iFieldLoc = 1, nFields
        offset = (iElem-1)*nFields + iFieldLoc
        do iDir = 1,layout%fStencil%QQ
          pos = varSys%method%val(iFieldLoc)%state_varPos(iDir)
          !fTmp( ?IDX?( pos, iElem, nScalars,globBC%nElems(iLevel) ) )   &
          !  & = bcBuffer(                                               &
          !  & ?IDX?( pos, posInBuffer, nScalars,globBC%nElems(iLevel) ) )
          fTmp( ?IDX?( pos, iElem, nScalars,globBC%nElems(iLevel) ) )     &
            & = state(?SAVE?( iDir,iFieldLoc,neighPos,QQ,nScalars,nSize,neigh ))
        end do
      end do
    end do !iElem

    ! Extrapolate equilibrium from local
    ! Calculate the equilibrium distribution
    call derVarPos%equilFromState( state  = fTmp,                    &
      &                            iField = iField,                  &
      &                            nElems = globBC%nElems( iLevel ), &
      &                            varSys = varSys,                  &
      &                            layout = layout,                  &
      &                            res    = fEq                      )

    ! Run over all elements in list for this boundary
    do iElem = 1, globBC%nElems( iLevel )
      elemPos = globBC%elemLvl(iLevel)%elem%val( iElem )
      do iDir = 1,layout%fStencil%QQN
        if( globBC%elemLvl(iLevel)%bitmask%val( iDir, iElem )) then
          ! Now assign the values
          state(?FETCH?( iDir, iField, elemPos, QQ, nScalars, nSize,neigh ) ) &
            & = fEq( iDir+(iElem-1)*QQ )
        end if
      end do
    end do

  end subroutine spc_outlet_eq
! ****************************************************************************** !


! ****************************************************************************** !
  !> Mole fraction boundary condition
  !! Usage
  !! -----
  !!```lua
  !!boundary_condition = {
  !! { label = 'outlet',
  !!    kind = 'spc_moleFrac',
  !!    moleFraction = 0.0
  !!     }
  !!```
  !! Post collision pdf of incoming link is updated with equilibrium functions
  !! which is similar to initial condition.
  !! \f$ \bar{f^c_k} = f^{eq}_k(\rho_k, u_k) + \frac{\lambda}{2}
  !!             (f^{eq}_k(\rho_k, u_k) + f^{eq}_k(\rho_k, u^{eq}_k})) \f$
  !! Here,
  !! $u_k$ - velocity of species k from original pdf computed from LSE
  !! $u^{eq}_k$ - equilibrium velocity of species k
  !!
  !!
  !! 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_moleFrac( 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) :: moleFrac(globBC%nElems(iLevel))
    real(kind=rk) :: rho, press
    integer :: iElem, iDir, iFieldLoc, nFields, pos, QQ
    integer :: bcMoleFrac_pos, elemPos, posInBuffer
    real(kind=rk) :: fTmp_all( layout%fStencil%QQ*varSys%nStateVars )
    real(kind=rk) :: fEq, fEqStar, ucx, ucxStar, ucxQuad, ucxQuadStar
    real(kind=rk) :: usq, usqStar
    real(kind=rk) :: velAvg(3), velQuad(3), velQuadStar(3), eqVel(3)
    real(kind=rk) :: vel(3,varSys%nStateVars)
    real(kind=rk) :: mom(3,varSys%nStateVars)
    real(kind=rk) :: mass_dens(varSys%nStateVars), num_dens(varSys%nStateVars)
    real(kind=rk) :: totMassDens
    real(kind=rk) :: moleFrac_loc(varSys%nStateVars)
    real(kind=rk) :: molWeightInv(varSys%nStateVars), phi(varSys%nStateVars)
    real(kind=rk) :: resi_coeff( varSys%nStateVars, varSys%nStateVars )
    real(kind=rk) :: weight0_inv, paramBInv
    type(mus_varSys_data_type), pointer :: fPtr
    !---------------------------------------------------------------------------
    call C_F_POINTER( varSys%method%val(iField)%method_Data, fPtr )
!    write(*,*) 'boundary label ', trim(me%label)
    QQ = layout%fStencil%QQ
    nFields = varSys%nStateVars

    do iFieldLoc = 1, nFields
      ! species properties
      ! molecular weight inverse
      molWeightInv(iFieldLoc) = fPtr%solverData%scheme%field(iFieldLoc) &
        &                           %fieldProp%species%molWeightInv
      ! 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
    weight0_inv = 1.0_rk/layout%weight(layout%fStencil%restPosition)

    paramBInv = 1.0_rk / mixture%paramB

    ! position of molefraction spacetime variable in varSys
    bcMoleFrac_pos = me%bc_states%moleFrac%varPos
    ! mole fraction
    call varSys%method%val(bcMoleFrac_pos)%get_valOfIndex( &
      & varSys  = varSys,                                  &
      & time    = sim_time,                                &
      & iLevel  = iLevel,                                  &
      & idx     = me%bc_states%moleFrac                    &
      &           %pntIndex%indexLvl(iLevel)               &
      &           %val(1:globBC%nElems(iLevel)),           &
      & nVals   = globBC%nElems(iLevel),                   &
      & res     = moleFrac                                 )

    ! Calculate the density of current element
    do iElem = 1, globBC%nElems(iLevel)
      posInBuffer = globBC%elemLvl( iLevel )%posInBcElemBuf%val(iElem)
      ! 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 ) = bcBuffer(                                  &
            & ?IDX?( pos, posInBuffer, nScalars,globBC%nElems(iLevel) ))
        end do
      end do

      !KM: using inital mixture number density rho0/mixtureMOlWeight
      rho = mixture%moleDens0LB *  moleFrac(iElem) &
        & * fieldProp%species%molWeight
      press = cs2 * rho * fieldProp%species%molWeigRatio

      !local density and momentum
      mass_dens = 0.0_rk
      vel = 0.0_rk
      do iFieldLoc = 1, nFields
        mass_dens(iFieldLoc) = &
          & sum( fTmp_all( varSys%method%val(iFieldLoc)%state_varPos(:) ) )
        num_dens(iFieldLoc) = mass_dens(iFieldLoc) * molWeightInv(iFieldLoc)
        !velocity
        call derVarPos%momFromState( state  = fTmp_all,              &
          &                          iField = iFieldLoc,             &
          &                          nElems = 1,                     &
          &                          varSys = varSys,                &
          &                          layout = layout,                &
          &                          res    = mom(:, iFieldLoc)      )
      end do

      mass_dens(iField) = rho
      num_dens(iField) = rho / fieldProp%species%molWeight

      !mass flux
      do iFieldLoc = 1, nFields
        vel(:, iFieldLoc) = mom(:, iFieldLoc) / mass_dens(iFieldLoc)
      end do

      ! total mass density
      totMassDens = sum(mass_dens)

      ! mole fraction
      moleFrac_loc = num_dens/sum(num_dens)

      ! equilibrium velocity
      eqVel = vel(:, iField)
      do iFieldLoc = 1, nFields
        eqVel(:) = eqVel(:) + resi_coeff(iField, iFieldLoc)               &
          &      * phi(iField) * moleFrac_loc(iFieldLoc)                  &
          &      * ( vel(:, iFieldLoc) - vel( :, iField ) )               &
          &      / mixture%paramB
      end do

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

      velQuadStar(:) = mixture%theta_eq*velAvg(:)                       &
        &            + (1.0_rk-mixture%theta_eq) * eqVel(:)
      velQuad(:) = mixture%theta_eq*velAvg(:)                           &
        &            + (1.0_rk-mixture%theta_eq)*vel(:, iField)

      usqStar = dot_product(velQuadStar, velQuadStar)*t2cs2inv
      usq = dot_product(velQuad, velQuad)*t2cs2inv

      elemPos = globBC%elemLvl( iLevel )%elem%val( iElem )
      do iDir =1,layout%fStencil%QQN
        if( globBC%elemLvl(iLevel)%bitmask%val( iDir, iElem )) then
          ucx = dot_product(layout%fStencil%cxDir(:, iDir),             &
            &               vel(:,iField))
          ucxQuad = dot_product(layout%fStencil%cxDir(:, iDir),         &
            &               velQuad)

          ucxStar = dot_product(layout%fStencil%cxDir(:, iDir),         &
            &                   eqVel)
          ucxQuadStar = dot_product(layout%fStencil%cxDir(:, iDir),     &
            &                   velQuadStar)

          ! eqVel is actually is rho_i*eqVel so ucxStar is not multiplied
          ! with rho in below equation
          fEqStar = layout%weight(iDir) * rho                  &
            & * ( phi(iField) + ucxStar * cs2inv               &
            & + ucxQuadStar * ucxQuadStar * t2cs4inv - usqStar )

          fEq = layout%weight(iDir) * rho          &
            & * ( phi(iField) + ucx * cs2inv       &
            & + ucxQuad * ucxQuad * t2cs4inv - usq )

          if ( iDir == layout%fStencil%restPosition ) then
            ! equilibrium at rest
            fEq = layout%weight( iDir ) * rho * (                          &
              & ( weight0_inv + (1.0_rk-weight0_inv) * phi(iField) ) - usq )
            fEqStar = layout%weight( iDir ) * rho * (                        &
              & ( weight0_inv + (1.0_rk-weight0_inv) * phi(iField) ) - usqStar )
          end if

          ! set equilibrium
          ! setting transformed pdf in similar to initial condition i.e
          ! \bar{f^c} = f^{eq}(\rho_spc, u_spc) + \frac{\lambda}{2}
          !             (f^{eq}(\rho_spc, u_spc) + f^{eq}(\rho_spc, u^{eq}_spc}))
          state(?FETCH?( iDir, iField, elemPos, QQ, nScalars, nSize,neigh )) &
            & = fEq + mixture%omega_diff*0.5_rk*( fEq - fEqStar )
        end if
      end do! iDir
    end do

  end subroutine spc_moleFrac
! ****************************************************************************** !

! ****************************************************************************** !
  !> Mole fraction boundary condition
  !! Usage
  !! -----
  !!```lua
  !!boundary_condition = {
  !! { label = 'outlet',
  !!    kind = 'spc_molefrac_eq',
  !!    moleFraction = 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_moleFrac_eq( me, state, bcBuffer, globBC, levelDesc, tree,   &
    &                         nSize, iLevel, sim_time, neigh, layout,         &
    &                         fieldProp, varPos, nScalars, varSys, derVarPos, &
    &                         physics, iField, mixture                        )
    ! -------------------------------------------------------------------- !
    !> global boundary type
    class(boundary_type) :: me
    !> Current state vector of iLevel
    real(kind=rk), intent(inout) :: state(:)
    !> size of state array ( in terms of elements )
    integer, intent(in) :: nSize
    !> state values of boundary elements of all fields of iLevel
    real(kind=rk), intent(in) :: bcBuffer(:)
    !> iLevel descriptor
    type(tem_levelDesc_type), intent(in) :: levelDesc
    !> Treelm Mesh
    type(treelmesh_type), intent(in) :: tree
    !> fluid parameters and properties
    type(mus_field_prop_type), intent(in) :: fieldProp
    !> stencil layout information
    type(mus_scheme_layout_type), intent(in) :: layout
    !> the level On which this boundary was invoked
    integer, intent(in) :: iLevel
    !> connectivity array corresponding to state vector
    integer, intent(in) :: neigh(:)
    !> global time information
    type(tem_time_type), intent(in)  :: sim_time
    !> pointer to field variable in the state vector
    integer, intent(in) :: varPos(:)
    !> number of Scalars in the scheme var system
    integer, intent(in) :: nScalars
    !> scheme variable system
    type(tem_varSys_type), intent(in) :: varSys
    !> position of derived quantities in varsys
    type(mus_derVarPos_type), intent(in) :: derVarPos
    !> scheme global boundary type
    type(glob_boundary_type), intent(in) :: globBC
    !> scheme global boundary type
    type(mus_physics_type), intent(in) :: physics
    !> current field
    integer, intent(in) :: iField
    !> mixture info
    type(mus_mixture_type), intent(in) :: mixture
    ! -------------------------------------------------------------------- !
    ! fEq uses AOS layout
    real(kind=rk) :: fEq( layout%fStencil%QQ*globBC%nElems(iLevel) )
    !store pdf of all fields for derive function pointer input
    real(kind=rk) :: fTmp( layout%fStencil%QQ*globBC%nElems(iLevel) &
      &              * varSys%nStateVars )
    real(kind=rk) :: mass_dens( globBC%nElems(iLevel)*varSys%nStateVars )
    real(kind=rk) :: moleFrac(globBC%nElems(iLevel))
    integer :: iElem, iDir, iFieldLoc, nFields, pos
    !> Velocity on boundary element
    real(kind=rk) :: uxB(globBC%nElems(iLevel)*varSys%nStateVars*3)
    real(kind=rk) :: velocity(3,globBC%nElems(iLevel)*varSys%nStateVars)
    integer :: offset, QQ, bcMoleFrac_pos, elemPos, posInBuffer
    ! ---------------------------------------------------------------------------
    QQ = layout%fStencil%QQ
    nFields = varSys%nStateVars

    ! position of molefraction spacetime variable in varSys
    bcMoleFrac_pos = me%bc_states%moleFrac%varPos
    ! mole fraction
    call varSys%method%val(bcMoleFrac_pos)%get_valOfIndex( &
      & varSys  = varSys,                                  &
      & time    = sim_time,                                &
      & iLevel  = iLevel,                                  &
      & idx     = me%bc_states%moleFrac                    &
      &           %pntIndex%indexLvl(iLevel)               &
      &           %val(1:globBC%nElems(iLevel)),           &
      & nVals   = globBC%nElems(iLevel),                   &
      & res     = moleFrac                                 )

    ! copy state
    do iElem = 1, globBC%nElems(iLevel)
      posInBuffer = globBC%elemLvl( iLevel )%posInBcElemBuf%val( iElem )
      do iFieldLoc = 1, nFields
        do iDir = 1,layout%fStencil%QQ
          pos = varSys%method%val(iFieldLoc)%state_varPos(iDir)
          fTmp( ?IDX?( pos, iElem, nScalars,globBC%nElems(iLevel) ) )   &
            & = bcBuffer(                                               &
            & ?IDX?( pos, posInBuffer, nScalars,globBC%nElems(iLevel) ) )
        end do
      end do !iField
    end do !iElem

    ! get density and velocity.
    ! if current field, use velocity and density defined in lua file.
    ! for other fields, derive density and velocity from state
    mass_dens = 0.0_rk
    do iFieldLoc = 1, nFields
      if (iFieldLoc == iField) then
        ! store velocity in input_loc array
        do iElem = 1, globBC%nElems(iLevel)
          offset = (iElem-1)*nFields + iFieldLoc
          ! compute current species mass density from specified molefraction at
          ! boundary
          ! rho = n_t * chi_i * m_i + massFrac_i * rho0 * KinePress/(cs2*phi_i)
          ! massFrac_i = rho_i/rho
          ! rho = n_t * chi_i * m_i + rho_i * KinePress/(cs2*phi_i)
          ! 1st term is zero order density
          ! 2nd term is second order density with kinematic mixture pressure,
          ! p = cs2*(sum(phi_k*rho_k) - min_a(m_a)*n0)/rho0
          !KM: using inital mixture number density rho0/mixtureMOlWeight
          !Using local tot_NuMdens increases density over time
          !kinePress = cs2 * ( dot_product(phi, mass_dens)                     &
          !  &       - minval(molWeight)*mixture%moleDens0LB )/mixture%rho0LB
          mass_dens( offset ) = mixture%moleDens0LB * moleFrac(iElem) &
            &                 * fieldProp%species%molWeight
        end do
      else
        do iElem = 1, globBC%nElems(iLevel)
          offset = (iElem-1)*nFields + iFieldLoc
          ! species density
          do iDir = 1, layout%fStencil%QQ
            pos = varSys%method%val(iFieldLoc)%state_varPos(iDir)
            mass_dens(offset) = mass_dens(offset) &
              & + fTmp(?IDX?( pos, iElem, nScalars,globBC%nElems(iLevel)) )
          end do !iDir
        end do !iElem
      end if
    end do !iField

    ! Derive all species velocities at once since
    ! it is inefficient to derive each species velocity
    call derVarPos%velocitiesFromState( state  = fTmp,                  &
      &                                 iField = iField,                &
      &                                 nElems = globBC%nElems(iLevel), &
      &                                 varSys = varSys,                &
      &                                 layout = layout,                &
      &                                 res    = uxB                    )

    do iElem = 1, globBC%nElems(iLevel)
      offset = (iElem-1)*nFields
      do iFieldLoc = 1, nFields
        velocity(:, offset+iFieldLoc) =      &
          & uxB(offset*3 + (iFieldLoc-1)*3 + 1 : &
          &     offset*3 + iFieldLoc*3 )
      end do !iField
    end do !iElem

    fEq = 0.0_rk
    ! Calculate the equilibrium distribution
    call derVarPos%equilFromMacro( density  = mass_dens,             &
      &                            velocity = velocity,              &
      &                            iField   = iField,                &
      &                            nElems   = globBC%nElems(iLevel), &
      &                            varSys   = varSys,                &
      &                            layout   = layout,                &
      &                            res      = fEq                    )

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

  end subroutine spc_moleFrac_eq
! ****************************************************************************** !

! ****************************************************************************** !
  !> Mole density boundary condition
  !! Usage
  !! -----
  !!```lua
  !!boundary_condition = {
  !! { label = 'outlet',
  !!    kind = 'spc_moledens_eq',
  !!    moleFraction = 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_moleDens_eq( me, state, bcBuffer, globBC, levelDesc, tree,   &
    &                         nSize, iLevel, sim_time, neigh, layout,         &
    &                         fieldProp, varPos, nScalars, varSys, derVarPos, &
    &                         physics, iField, mixture                        )
    ! -------------------------------------------------------------------- !
    !> global boundary type
    class(boundary_type) :: me
    !> Current state vector of iLevel
    real(kind=rk), intent(inout) :: state(:)
    !> size of state array ( in terms of elements )
    integer, intent(in) :: nSize
    !> state values of boundary elements of all fields of iLevel
    real(kind=rk), intent(in) :: bcBuffer(:)
    !> iLevel descriptor
    type(tem_levelDesc_type), intent(in) :: levelDesc
    !> Treelm Mesh
    type(treelmesh_type), intent(in) :: tree
    !> fluid parameters and properties
    type(mus_field_prop_type), intent(in) :: fieldProp
    !> stencil layout information
    type(mus_scheme_layout_type), intent(in) :: layout
    !> the level On which this boundary was invoked
    integer, intent(in) :: iLevel
    !> connectivity array corresponding to state vector
    integer, intent(in) :: neigh(:)
    !> global time information
    type(tem_time_type), intent(in)  :: sim_time
    !> pointer to field variable in the state vector
    integer, intent(in) :: varPos(:)
    !> number of Scalars in the scheme var system
    integer, intent(in) :: nScalars
    !> scheme variable system
    type(tem_varSys_type), intent(in) :: varSys
    !> position of derived quantities in varsys
    type(mus_derVarPos_type), intent(in) :: derVarPos
    !> scheme global boundary type
    type(glob_boundary_type), intent(in) :: globBC
    !> scheme global boundary type
    type(mus_physics_type), intent(in) :: physics
    !> current field
    integer, intent(in) :: iField
    !> mixture info
    type(mus_mixture_type), intent(in) :: mixture
    ! -------------------------------------------------------------------- !
    ! ---------------------------------------------------------------------------
    ! fEq uses AOS layout
    real(kind=rk) :: fEq( layout%fStencil%QQ*globBC%nElems(iLevel) )
    !store pdf of all fields for derive function pointer input
    real(kind=rk) :: fTmp( layout%fStencil%QQ*globBC%nElems(iLevel) &
      &              * varSys%nStateVars )
    real(kind=rk) :: mass_dens( globBC%nElems(iLevel)*varSys%nStateVars )
    real(kind=rk) :: moleDens(globBC%nElems(iLevel))
    integer :: iElem, iDir, iFieldLoc, nFields, pos
    !> Velocity on boundary element
    real(kind=rk) :: uxB(globBC%nElems(iLevel)*varSys%nStateVars*3)
    real(kind=rk) :: velocity(3,globBC%nElems(iLevel)*varSys%nStateVars)
    integer :: offset, QQ, bcMoleDens_pos, elemPos, posInBuffer
    ! ---------------------------------------------------------------------------
    QQ = layout%fStencil%QQ
    nFields = varSys%nStateVars

    ! position of molefraction spacetime variable in varSys
    bcMoleDens_pos = me%bc_states%moleDens%varPos
    ! mole fraction
    call varSys%method%val(bcMoleDens_pos)%get_valOfIndex( &
      & varSys  = varSys,                                  &
      & time    = sim_time,                                &
      & iLevel  = iLevel,                                  &
      & idx     = me%bc_states%moleFrac                    &
      &           %pntIndex%indexLvl(iLevel)               &
      &           %val(1:globBC%nElems(iLevel)),           &
      & nVals   = globBC%nElems(iLevel),                   &
      & res     = moleDens                                 )

    ! Convert to lattice Unit
    moleDens = moleDens / physics%moleDens0

    ! copy state
    do iElem = 1, globBC%nElems(iLevel)
      posInBuffer = globBC%elemLvl( iLevel )%posInBcElemBuf%val( iElem )
      do iFieldLoc = 1, nFields
        do iDir = 1,layout%fStencil%QQ
          pos = varSys%method%val(iFieldLoc)%state_varPos(iDir)
          fTmp( ?IDX?( pos, iElem, nScalars,globBC%nElems(iLevel) ) )   &
            & = bcBuffer(                                               &
            & ?IDX?( pos, posInBuffer, nScalars,globBC%nElems(iLevel) ) )
        end do
      end do !iField
    end do !iElem

    ! get density and velocity.
    ! if current field, use velocity and density defined in lua file.
    ! for other fields, derive density and velocity from state
    mass_dens = 0.0_rk
    do iFieldLoc = 1, nFields
      if (iFieldLoc == iField) then
        ! store velocity in input_loc array
        do iElem = 1, globBC%nElems(iLevel)
          offset = (iElem-1)*nFields + iFieldLoc
          ! compute current species mass density from specified molefraction at
          ! boundary
          ! rho = n_t * chi_i * m_i + massFrac_i * rho0 * KinePress/(cs2*phi_i)
          ! massFrac_i = rho_i/rho
          ! rho = n_t * chi_i * m_i + rho_i * KinePress/(cs2*phi_i)
          ! 1st term is zero order density
          ! 2nd term is second order density with kinematic mixture pressure,
          ! p = cs2*(sum(phi_k*rho_k) - min_a(m_a)*n0)/rho0
          !KM: using inital mixture number density rho0/mixtureMOlWeight
          !Using local tot_NuMdens increases density over time
          !kinePress = cs2 * ( dot_product(phi, mass_dens)                     &
          !  &       - minval(molWeight)*mixture%moleDens0LB )/mixture%rho0LB
          mass_dens( offset ) = moleDens(iElem) * fieldProp%species%molWeight
        end do
      else
        do iElem = 1, globBC%nElems(iLevel)
          offset = (iElem-1)*nFields + iFieldLoc
          ! species density
          do iDir = 1, layout%fStencil%QQ
            pos = varSys%method%val(iFieldLoc)%state_varPos(iDir)
            mass_dens(offset) = mass_dens(offset) &
              & + fTmp(?IDX?( pos, iElem, nScalars,globBC%nElems(iLevel)) )
          end do !iDir
        end do !iElem
      end if
    end do !iField

    ! Derive all species velocities at once since
    ! it is inefficient to derive each species velocity
    call derVarPos%velocitiesFromState( state  = fTmp,                  &
      &                                 iField = iField,                &
      &                                 nElems = globBC%nElems(iLevel), &
      &                                 varSys = varSys,                &
      &                                 layout = layout,                &
      &                                 res    = uxB                    )

    do iElem = 1, globBC%nElems(iLevel)
      offset = (iElem-1)*nFields
      do iFieldLoc = 1, nFields
        velocity(:, offset+iFieldLoc) =      &
          & uxB(offset*3 + (iFieldLoc-1)*3 + 1 : &
          &     offset*3 + iFieldLoc*3 )
      end do !iField
    end do !iElem

    fEq = 0.0_rk
    ! Calculate the equilibrium distribution
    call derVarPos%equilFromMacro( density  = mass_dens,             &
      &                            velocity = velocity,              &
      &                            iField   = iField,                &
      &                            nElems   = globBC%nElems(iLevel), &
      &                            varSys   = varSys,                &
      &                            layout   = layout,                &
      &                            res      = fEq                    )

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

  end subroutine spc_moleDens_eq
! ****************************************************************************** !



! ****************************************************************************** !
  !> Mole fraction boundary condition with thermodynamic factor
  !! Usage
  !! -----
  !!```lua
  !!boundary_condition = {
  !! { label = 'outlet',
  !!    kind = 'spc_moleFrac_wtdf',
  !!    moleFraction = 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_moleFrac_wtdf( 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
    !store pdf of all fields for derive function pointer input
    real(kind=rk) :: fTmp( varSys%nScalars )
    real(kind=rk) :: mass_dens( varSys%nStateVars )
    real(kind=rk) :: num_dens( varSys%nStateVars )
    real(kind=rk) :: tot_numDens, tot_massDens
    real(kind=rk) :: moleFrac(globBC%nElems(iLevel))
    integer :: iElem, iDir, iFieldLoc, iField_2, iField_3, nFields, pos, QQ
    real(kind=rk), dimension(varSys%nStateVars) :: molWeight, phi
    real(kind=rk), dimension(varSys%nStateVars) :: moleFrac_loc
    real(kind=rk), dimension(3, varSys%nStateVars) :: first_moments, velocity
    real(kind=rk), &
      & dimension(varSys%nStateVars, varSys%nStateVars) :: matA, invA, &
      & resi_coeff, thermodynamic_fac, inv_thermodyn_fac, diff_coeff
    real(kind=rk) :: temp, press, phy_moleDens_fac
    real(kind=rk) :: omega, paramB, theta_eq
    real(kind=rk) :: eqVel(3), velAvg(3), velQuad(3), usqr, ucx, ucxQuadTerm
    type(mus_varSys_data_type), pointer :: fPtr
    integer :: bcMoleFrac_pos, posInBuffer, posInState
    ! ------------------------------------------------------------------------
    call C_F_POINTER( varSys%method%val(iField)%method_Data, fPtr )
    QQ = layout%fStencil%QQ
    nFields = varSys%nStateVars

    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

    !KM @todo check moleDens for multilevel
    phy_moleDens_fac = physics%moleDens0
    !fixed parameter
    paramB = mixture%paramB
    !equilibrium theta
    theta_eq = mixture%theta_eq

    omega = mixture%relaxLvl(iLevel)%omega_diff
    ! temperature
    temp = mixture%temp0
    ! atmospheric pressure
    press = mixture%atm_press
!write(*,*) 'iField ', iField
!write(*,*) 'resi_coeff ', resi_coeff
!write(*,*) 'molWeight moleFrac BC', molWeight
!write(*,*) 'phi ', phi
!stop

    ! position of molefraction spacetime variable in varSys
    bcMoleFrac_pos = me%bc_states%moleFrac%varPos
    ! mole fraction
    call varSys%method%val(bcMoleFrac_pos)%get_valOfIndex( &
      & varSys  = varSys,                                  &
      & time    = sim_time,                                &
      & iLevel  = iLevel,                                  &
      & idx     = me%bc_states%moleFrac                    &
      &           %pntIndex%indexLvl(iLevel)               &
      &           %val(1:globBC%nElems(iLevel)),           &
      & nVals   = globBC%nElems(iLevel),                   &
      & res     = moleFrac                                 )

    ! Calculate the density of current element
    do iElem = 1, globBC%nElems(iLevel)
!write(*,*) 'IElem ', iElem
      mass_dens = 0.0_rk
      first_moments = 0.0_rk
      posInBuffer = globBC%elemLvl( iLevel )%posInBcElemBuf%val( iElem )
      do iFieldLoc = 1, varSys%nStateVars
        do iDir = 1,layout%fStencil%QQ
          pos = varSys%method%val(iFieldLoc)%state_varPos(iDir)
          fTmp( pos ) = bcBuffer(                                   &
            & ?IDX?( pos, posInBuffer, nScalars,globBC%nElems(iLevel) ) )
          mass_dens(iFieldLoc) = mass_dens(iFieldLoc) + fTmp( pos )

          !field momentum (rho*u)
          first_moments( 1, iFieldLoc ) = first_moments( 1, iFieldLoc ) &
            & + fTmp( pos ) * layout%fStencil%cxDir( 1, iDir )

          first_moments( 2, iFieldLoc ) = first_moments( 2, iFieldLoc ) &
            & + fTmp( pos ) * layout%fStencil%cxDir( 2, iDir )

          first_moments( 3, iFieldLoc ) = first_moments( 3, iFieldLoc ) &
            & + fTmp( pos ) * layout%fStencil%cxDir( 3, iDir )

        end do
        num_dens(iFieldLoc) = mass_dens(iFieldLoc)/molWeight(iFieldLoc)
      end do

      ! update num_dens, massDens and moleFrac of current species with specified
      ! molefraction
      num_dens(iField) = moleFrac(iElem) * mixture%moleDens0LB
      mass_dens(iField) = num_dens(iField) * molWeight(iField)

      !total number density
      tot_NumDens = sum(num_dens)

      !total mass density
      tot_massDens = sum(mass_dens)

      !mole fraction
      moleFrac_loc(:) =  num_dens(:)/tot_NumDens
!write(*,*) 'num_dens ', num_dens
!write(*,*) 'num_dens*phy_moleDens_fac ', num_dens*phy_moleDens_fac

      ! MS-Diff coeff matrix from C++ code
      call mus_calc_MS_DiffMatrix( nFields, temp, press,                       &
        &                          num_dens*phy_moleDens_fac, diff_coeff )
!write(*,*) 'diff_coeff ', diff_coeff
      ! Convert to lattice unit
      resi_coeff = physics%fac(iLevel)%diffusivity/diff_coeff
!write(*,*) 'resi_coeff ', resi_coeff

      call mus_calc_thermFactor( nFields, temp, press, moleFrac_loc,           &
        &                        thermodynamic_fac )
!write(*,*) 'thermodyn_fac ', thermodynamic_fac

      inv_thermodyn_fac = invert_matrix( thermodynamic_fac )
!write(*,*) 'inv_thermodyn_fac ', inv_thermodyn_fac

      matA = 0.0_rk
      !build up matrix to solver LSE for actual velocity
      do iFieldLoc = 1, nFields
        !set diagonal part
        matA(iFieldLoc, iFieldLoc) = 1.0_rk
        do iField_2 = 1, nFields
          do iField_3 = 1, nFields
            matA(iFieldLoc, iField_2) = matA(iFieldLoc, iField_2)              &
              &                       + omega * 0.5_rk                         &
              &                       * inv_thermodyn_fac(iFieldLoc, iField_2) &
              &                       * resi_coeff(iField_2, iField_3)         &
              &                       * phi(iField_2) * moleFrac_loc(iField_3) &
              &                       / paramB
          end do
        end do
        !set non-diagonal part
        do iField_2 = 1, nFields
          do iField_3 = 1, nFields
            matA(iFieldLoc, iField_3) = matA(iFieldLoc, iField_3)              &
              &                       - omega * 0.5_rk                         &
              &                       * inv_thermodyn_fac(iFieldLoc, iField_2) &
              &                       * resi_coeff(iField_2, iField_3)         &
              &                       * phi(iField_3) * moleFrac_loc(iField_2) &
              &                       / paramB
          end do
        end do
      end do
!write(*,*) 'matA ', matA

      ! invert matrix
      invA = invert_matrix( matA )
!write(*,*) 'invA ', invA

      !actual velocity of all species
      velocity(1, :) = matmul( invA, first_moments(1,:) ) / mass_dens(:)
      velocity(2, :) = matmul( invA, first_moments(2,:) ) / mass_dens(:)
      velocity(3, :) = matmul( invA, first_moments(3,:) ) / mass_dens(:)
!write(*,*) 'velocity ', velocity

      ! equilibrium velocity with thermodynamic factor
      eqVel( : ) = mass_dens(iField)*velocity( :, iField )
      do iField_2 = 1, nFields
        do iField_3 = 1, nFields
          eqVel( : ) = eqVel( : )                                      &
            &        + inv_thermodyn_fac(iField, iField_2)             &
            &        * mass_dens(iField_2)                             &
            &        * resi_coeff( iField_2, iField_3 ) * phi(iField_2)&
            &        * moleFrac_loc(iField_3)                          &
            &        * (velocity(:, iField_3) - velocity(:,iField_2))  &
            &        / paramB
        end do
      end do
!write(*,*) 'eqVel ', eqVel
      !compute mass averaged mixture velocity
      velAvg(1) = dot_product( mass_dens, velocity(1,:) )/tot_massDens
      velAvg(2) = dot_product( mass_dens, velocity(2,:) )/tot_massDens
      velAvg(3) = dot_product( mass_dens, velocity(3,:) )/tot_massDens

      velQuad = theta_eq*velAvg                                              &
        &     + (1.0_rk - theta_eq) * eqVel(:) / mass_dens(iField)

      usqr = dot_product( velQuad, velQuad ) * t2cs2inv

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

          ucx = dot_product(                                                   &
            & dble(layout%fStencil%cxDir(:, iDir)), eqVel(:) )

          ucxQuadTerm = dot_product(                                           &
            & dble(layout%fStencil%cxDir(:, iDir)), velQuad )

          ! Equilibrium
          feq = layout%weight(iDir) * ( mass_dens(iField) * ( phi(iField)      &
            & + ucxQuadTerm * ucxQuadTerm * t2cs4inv - usqr ) + ucx * cs2inv )

          if ( iDir == layout%fStencil%restPosition ) then
          ! equilibrium at rest
            select case( trim(layout%fStencil%label) )
            case('d2q9')
              feq = layout%weight( iDir ) * mass_dens(iField) * ( &
                    & ( 9._rk - 5._rk * phi(iField) )/4._rk - usqr )
            case('d3q19')
              feq = layout%weight( iDir ) * mass_dens(iField) * ( &
                    & ( 3._rk - 2._rk * phi(iField) ) - usqr )
            end select
          end if

          ! set equilibrium
          state(                                                               &
            & ?FETCH?( iDir, iField, posInState, QQ, nScalars, nSize, neigh ))&
            & = fEq
!write(*,*) 'fEq ', fEq
        end if
      end do
    end do

  end subroutine spc_moleFrac_wtdf
! ****************************************************************************** !


! ****************************************************************************** !
  !> molar flux equilibrium boundary condition
  !! Usage
  !! -----
  !!```lua
  !!boundary_condition = {
  !! { label = 'outlet',
  !!    kind = 'spc_moleflux',
  !!    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_moleFlux_eq( me, state, bcBuffer, globBC, levelDesc, tree,   &
    &                         nSize, iLevel, sim_time, neigh, layout,         &
    &                         fieldProp, varPos, nScalars, varSys, derVarPos, &
    &                         physics, iField, mixture                        )
    ! -------------------------------------------------------------------- !
    !> global boundary type
    class(boundary_type) :: me
    !> Current state vector of iLevel
    real(kind=rk), intent(inout) :: state(:)
    !> size of state array ( in terms of elements )
    integer, intent(in) :: nSize
    !> state values of boundary elements of all fields of iLevel
    real(kind=rk), intent(in) :: bcBuffer(:)
    !> iLevel descriptor
    type(tem_levelDesc_type), intent(in) :: levelDesc
    !> Treelm Mesh
    type(treelmesh_type), intent(in) :: tree
    !> fluid parameters and properties
    type(mus_field_prop_type), intent(in) :: fieldProp
    !> stencil layout information
    type(mus_scheme_layout_type), intent(in) :: layout
    !> the level On which this boundary was invoked
    integer, intent(in) :: iLevel
    !> connectivity array corresponding to state vector
    integer, intent(in) :: neigh(:)
    !> global time information
    type(tem_time_type), intent(in)  :: sim_time
    !> pointer to field variable in the state vector
    integer, intent(in) :: varPos(:)
    !> number of Scalars in the scheme var system
    integer, intent(in) :: nScalars
    !> scheme variable system
    type(tem_varSys_type), intent(in) :: varSys
    !> position of derived quantities in varsys
    type(mus_derVarPos_type), intent(in) :: derVarPos
    !> scheme global boundary type
    type(glob_boundary_type), intent(in) :: globBC
    !> scheme global boundary type
    type(mus_physics_type), intent(in) :: physics
    !> current field
    integer, intent(in) :: iField
    !> mixture info
    type(mus_mixture_type), intent(in) :: mixture
    ! -------------------------------------------------------------------- !
    real(kind=rk) :: fEq( layout%fStencil%QQ*globBC%nElems(iLevel) )
    !store pdf of all fields for derive function pointer input
    real(kind=rk) :: fTmp( layout%fStencil%QQ*globBC%nElems(iLevel) &
      &              * varSys%nStateVars )
    real(kind=rk) :: moleFlux(globBC%nElems(iLevel)*3), inv_flux
    !> Velocity on boundary element
    real(kind=rk) :: uxB(globBC%nElems(iLevel)*varSys%nStateVars*3)
    real(kind=rk) :: mass_dens( globBC%nElems(iLevel)*varSys%nStateVars )
    real(kind=rk) :: velocity(3,globBC%nElems(iLevel)*varSys%nStateVars)
    real(kind=rk) :: num_dens
    integer :: iElem, iDir, iFieldLoc, nFields, pos, QQ
    integer :: offset, posInBuffer, bcMoleFlux_pos, elemPos
    ! ------------------------------------------------------------------------
    QQ = layout%fStencil%QQ
    nFields = varSys%nStateVars
    inv_flux = 1.0_rk / physics%fac( iLevel )%moleFlux

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

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

    ! copy state
    do iElem = 1, globBC%nElems(iLevel)
      posInBuffer = globBC%elemLvl( iLevel )%posInBcElemBuf%val( iElem )
      do iFieldLoc = 1, nFields
        do iDir = 1,layout%fStencil%QQ
          pos = varSys%method%val(iFieldLoc)%state_varPos(iDir)
          fTmp( ?IDX?( pos, iElem, nScalars,globBC%nElems(iLevel) ) ) = &
            & bcBuffer(                                                 &
            & ?IDX?( pos, posInBuffer, nScalars,globBC%nElems(iLevel) ) )
        end do
      end do !iField
    end do !iElem

    ! get density and velocity.
    ! if current field, use velocity and density defined in lua file.
    ! for other fields, derive density and velocity from state
    ! species density
    mass_dens = 0.0_rk
    do iElem = 1, globBC%nElems(iLevel)
      do iFieldLoc = 1, nFields
        offset = (iElem-1)*nFields + iFieldLoc
        do iDir = 1, layout%fStencil%QQ
          pos = varSys%method%val(iFieldLoc)%state_varPos(iDir)
          mass_dens(offset) = mass_dens(offset) &
            & + fTmp(?IDX?( pos, iElem, nScalars,globBC%nElems(iLevel)) )
        end do !iDir
      end do !iField
    end do !iElem

    ! Derive all species velocities at once since
    ! it is inefficient to derive each species velocity
    call derVarPos%velocitiesFromState( state  = fTmp,                  &
      &                                 iField = iField,                &
      &                                 nElems = globBC%nElems(iLevel), &
      &                                 varSys = varSys,                &
      &                                 layout = layout,                &
      &                                 res    = uxB                    )

    do iFieldLoc = 1, nFields
      if (iFieldLoc == iField) then
        ! store velocity in input_loc array
        do iElem = 1, globBC%nElems(iLevel)
          offset = (iElem-1)*nFields + iFieldLoc
          num_dens = mass_dens(offset) * fieldProp%species%molWeightInv
          velocity(:,offset) = moleFlux((iElem-1)*3+1:iElem*3) / num_dens
        end do
      else
        do iElem = 1, globBC%nElems(iLevel)
          offset = (iElem-1)*nFields + iFieldLoc
          velocity(:,offset) =                             &
            & uxB((iElem-1)*nFields*3 + (iFieldLoc-1)*3 + 1 : &
            &     (iElem-1)*nFields*3 + iFieldLoc*3 )
        end do !iElem
      end if
    end do !iField

    ! Calculate the equilibrium distribution
    call derVarPos%equilFromMacro( density  = mass_dens,             &
      &                            velocity = velocity,              &
      &                            iField   = iField,                &
      &                            nElems   = globBC%nElems(iLevel), &
      &                            varSys   = varSys,                &
      &                            layout   = layout,                &
      &                            res      = fEq                    )

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

  end subroutine spc_moleFlux_eq
! ****************************************************************************** !


! ****************************************************************************** !
  !> molar flux boundary condition like velocity bounce back bc type
  !! Usage
  !! -----
  !!```lua
  !!boundary_condition = {
  !! { label = 'outlet',
  !!    kind = 'spc_moleflux',
  !!    mole_flux = 'mole_flux'
  !! }
  !!}
  !!variable = {
  !!  name = 'mole_flux',
  !!  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_moleFlux( 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) :: molWeight
    real(kind=rk) :: massFlux(3)
    real(kind=rk) :: moleFlux(globBC%nElems(iLevel)*3), inv_flux
    integer :: iElem, iDir, QQ, posInBuffer, offset, bcMoleFlux_pos
    integer :: elemPos
    ! ------------------------------------------------------------------------
    QQ = layout%fStencil%QQ
    molWeight = fieldProp%species%molWeight
    inv_flux = 1.0_rk / physics%fac( iLevel )%moleFlux

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

    ! If physical quantities are given, transform to lattice units by division
    ! with the conversion factor
    moleFlux = moleFlux * inv_flux
    !moleFlux = moleFlux / physics%fac(iLevel)%flux

    ! Calculate the density of current element
    do iElem = 1, globBC%nElems(iLevel)
      posInBuffer = globBC%elemLvl( iLevel )%posInBcElemBuf%val( iElem )
      do iDir = 1,layout%fStencil%QQ
        fTmp( iDir ) = bcBuffer(                    &
       & ?IDX?( varpos(iDir), posInBuffer, nScalars,globBC%nElems(iLevel) ) )
      end do
      !caulate mass flux from moleFlux
      ! massFlux = moleflux * molWeight
      offset = (iElem-1)*3
      massFlux(1) = moleFlux(offset+1) * molWeight
      massFlux(2) = moleFlux(offset+2) * molWeight
      massFlux(3) = moleFlux(offset+3) * molWeight
      !write(dbgUnit(1),*) 'massFlux ', massFlux
      !write(*,*) iElem, 'mass_flux ', massFlux

      elemPos = globBC%elemLvl( iLevel )%elem%val( iElem )
      do iDir = 1,layout%fStencil%QQN
        if( globBC%elemLvl(iLevel)%bitmask%val( iDir, iElem )) then
          state(?FETCH?( iDir, iField, elemPos, QQ, nScalars, nSize,neigh )) &
              & = fTmp(layout%fStencil%cxDirInv( iDir ))            &
              & + layout%weight( iDir )*2._rk*cs2inv                &
              &    * ( layout%fStencil%cxDir( 1, iDir )*massFlux(1) &
              &    +   layout%fStencil%cxDir( 2, iDir )*massFlux(2) &
              &    +   layout%fStencil%cxDir( 3, iDir )*massFlux(3) )
        end if
      end do
    end do

  end subroutine spc_moleFlux
! ****************************************************************************** !

! **************************************************************************** !
  !> Inflow boundary condition for solvent based on non-Equilbrium
  !! extrapolation method. Similar to spc_velocity_noneq_expol except
  !! the mass density for solvent is enforced such that total moleDens0 is
  !! maintained.
  !! Default qVal=1.0.
  !! Usage
  !! -----
  !!```lua
  !!boundary_condition = {
  !! { label = 'inlet',
  !!    kind = 'spc_solvent_inflow',
  !!    velocity = {0.1,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_solvent_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) :: fTmpAll_1f(nScalars), fTmpAll_2f(nScalars)
    real(kind=rk) :: fEq_b, fEq_1f
    real(kind=rk) :: velAvg_1f(3), vel_b(3), eqVel_1f(3)
    real(kind=rk) :: usq_1f, ucx_1f, usq_b, ucx_b, ucxQuad_1f
    real(kind=rk), dimension(varSys%nStateVars) :: mass_dens_1f, mass_dens_2f
    real(kind=rk), dimension(varSys%nStateVars) :: num_dens_1f, num_dens_2f
    real(kind=rk) :: resi_coeff(varSys%nStateVars, varSys%nStateVars)
    real(kind=rk) :: mass_dens_b, num_dens_b, tot_NumDens_b
    real(kind=rk) :: omega_fac, paramBInv
    integer :: iElem, iDir, QQ, iFld, vPos
    integer :: bcVel_pos
    integer :: elemPos, elemOff, posInBuffer, neighPos
    real(kind=rk) :: vel_w(3*globBC%nElems(iLevel)), inv_vel
    type(mus_varSys_data_type), pointer :: fPtr
    integer :: nFields
    ! -------------------------------------------------------------------------
    nFields = varSys%nStateVars
    ! 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:globBC%nElems(iLevel)),      &
      & nVals   = globBC%nElems(iLevel),              &
      & res     = vel_w                               )

    ! convert physical velocity into LB velocity
    inv_vel = 1.0_rk / physics%fac( iLevel )%vel
    vel_w = vel_w * inv_vel

    call C_F_POINTER( varSys%method%val(iField)%method_Data, fPtr )
    ! resitivity coefficient of all species
    do iFld = 1, nFields
      !resi_coeff(iFld, :) = species(iFld)%resi_coeff(:)
      resi_coeff(iFld, :) = fPtr%solverData%scheme                          &
        &                                  %field(iFld)%fieldProp           &
        &                                              %species%resi_coeff(:)
    end do
    ! omega factor of LSE
    omega_fac = mixture%omega_diff * 0.5_rk
    paramBInv = 1.0_rk / mixture%paramB

    associate( auxField => fPtr%solverData%scheme%auxField(iLevel)%val,      &
      &        species => fPtr%solverData%scheme%field(:)%fieldProp%species, &
      &        neighBufferPre => me%neigh(iLevel)%neighBufferPre,            &
      &        stencil => layout%fStencil                                    )

      QQ = layout%fStencil%QQ

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

        ! Position of current boundary element in bcBuffer
        posInBuffer = globBC%elemLvl(iLevel)%posInBcElemBuf%val( iElem )
        elemOff = (posInBuffer-1)*nScalars
        ! Position of 2nd fluid neighbor in state array
        neighPos = me%neigh(iLevel)%posInState(1, iElem)
        mass_dens_2f = 0.0_rk
        do iFld = 1, nFields
          do iDir = 1, QQ
            ! Position of current field variable current dir in state array
            vPos = varSys%method%val(iFld)%state_varPos(iDir)
            ! State array of 1st fluid element
            fTmpAll_1f(vPos) = bcBuffer( elemOff + vPos )
            ! Pre-Collision values of second fluid element
            fTmpAll_2f(vPos) = state(                                       &
              & ?FETCH?(iDir, iFld, neighPos, QQ, nScalars, nSize, neigh) )
            ! mass density of second neighbor
            mass_dens_2f(iFld) = mass_dens_2f(iFld) + fTmpAll_2f(vPos)
          end do
          num_dens_2f(iFld) = mass_dens_2f(iFld) * species(iFld)%molWeightInv
        end do

        ! calculate spc density and velAvg from PDF of first fluid
        call calcDensAndVelsFromPDF( nFields    = nFields,         &
          &                          iField     = iField,          &
          &                          mass_dens  = mass_dens_1f,    &
          &                          num_dens   = num_dens_1f,     &
          &                          velAvg     = velAvg_1f,       &
          &                          eqVel      = eqVel_1f,        &
          &                          varSys     = varSys,          &
          &                          pdf        = fTmpAll_1f,      &
          &                          stencil    = layout%fStencil, &
          &                          species    = species,         &
          &                          resi_coeff = resi_coeff,      &
          &                          omega_fac  = omega_fac,       &
          &                          paramBInv  = paramBInv        )

        ! fluid node velocity square term of equilibrium function
        usq_1f = dot_product( velAvg_1f(1:3), velAvg_1f(1:3) ) * t2cs2inv

        ! Second-order extrapolation of number density of current species
        num_dens_b = ( 4.0_rk * num_dens_1f(iField) - num_dens_2f(iField) ) &
          &         / 3.0_rk
        ! Second-order extrapolation of numerical total number density
        tot_NumDens_b = ( 4.0_rk * sum(num_dens_1f) - sum(num_dens_2f) ) &
          &           / 3.0_rk
        ! Compute number density of current species at boundary element
        num_dens_b = mixture%moleDens0LB - ( tot_NumDens_b - num_dens_b )
        ! Mass density of current species at boundary element
        mass_dens_b = num_dens_b * species(iField)%molWeight
        !mass_dens_b = mixture%rho0LB - ( sum(mass_dens_1f) - mass_dens_1f(iField) )

        ! velocity at boundary node for qVal = 1.0
        vel_b = vel_w( (iElem-1)*3 + 1: iElem*3 )

        ! boundary node velocity square term of equilibrium function
        usq_b = dot_product( vel_b(1:3), vel_b(1:3) ) * t2cs2inv

        ! Position of this boundary element in state array
        elemPos = globBC%elemLvl(iLevel)%elem%val(iElem)

        do iDir = 1, layout%fStencil%QQN
          if ( globBC%elemLvl(iLevel)%bitmask%val(iDir, iElem) ) then
            ! Compute equilibrium on fluid node
            ucx_1f = dot_product( stencil%cxDirRK(1:3, iDir), velAvg_1f )
            ucxQuad_1f = dot_product( stencil%cxDirRK(1:3, iDir), velAvg_1f )

            fEq_1f = layout%weight(iDir) * mass_dens_1f(iField)      &
              &   * ( species(iField)%molWeigRatio + ucx_1f * cs2inv &
              &   + ucxQuad_1f * ucxQuad_1f * t2cs4inv - usq_1f      )

            ! Compute equilibrium on boundary node
            ucx_b = dot_product( stencil%cxDirRK(1:3, iDir), vel_b )

            fEq_b = layout%weight(iDir) * mass_dens_b               &
              &   * ( species(iField)%molWeigRatio + ucx_b * cs2inv &
              &   + ucx_b * ucx_b * t2cs4inv - usq_b                )

            state(?FETCH?(iDir, iField, elemPos, QQ, nScalars, nSize, neigh))  &
              & = fEq_b + ( fTmpAll_1f(varPos(iDir)) - fEq_1f )

          end if
        end do
      end do
    end associate

  end subroutine spc_solvent_inflow
! **************************************************************************** !


! **************************************************************************** !
  !> Velocity boundary condition based on non-Equilbrium extrapolation method.
  !! Default qVal=1.0.
  !! Usage
  !! -----
  !!```lua
  !!boundary_condition = {
  !! { label = 'inlet',
  !!    kind = 'spc_velocity_noneq_expol',
  !!    velocity = {0.1,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_velocity_noneq_expol( me, state, bcBuffer, globBC, levelDesc, &
    &                                  tree, nSize, iLevel, sim_time, neigh,   &
    &                                  layout, fieldProp, varPos, nScalars,    &
    &                                  varSys, derVarPos, physics, iField,     &
    &                                  mixture                                 )
    ! -------------------------------------------------------------------- !
    !> global boundary type
    class(boundary_type) :: me
    !> Current state vector of iLevel
    real(kind=rk), intent(inout) :: state(:)
    !> size of state array ( in terms of elements )
    integer, intent(in) :: nSize
    !> state values of boundary elements of all fields of iLevel
    real(kind=rk), intent(in) :: bcBuffer(:)
    !> iLevel descriptor
    type(tem_levelDesc_type), intent(in) :: levelDesc
    !> Treelm Mesh
    type(treelmesh_type), intent(in) :: tree
    !> fluid parameters and properties
    type(mus_field_prop_type), intent(in) :: fieldProp
    !> stencil layout information
    type(mus_scheme_layout_type), intent(in) :: layout
    !> the level On which this boundary was invoked
    integer, intent(in) :: iLevel
    !> connectivity array corresponding to state vector
    integer, intent(in) :: neigh(:)
    !> global time information
    type(tem_time_type), intent(in)  :: sim_time
    !> pointer to field variable in the state vector
    integer, intent(in) :: varPos(:)
    !> number of Scalars in the scheme var system
    integer, intent(in) :: nScalars
    !> scheme variable system
    type(tem_varSys_type), intent(in) :: varSys
    !> position of derived quantities in varsys
    type(mus_derVarPos_type), intent(in) :: derVarPos
    !> scheme global boundary type
    type(glob_boundary_type), intent(in) :: globBC
    !> scheme global boundary type
    type(mus_physics_type), intent(in) :: physics
    !> current field
    integer, intent(in) :: iField
    !> mixture info
    type(mus_mixture_type), intent(in) :: mixture
    ! -------------------------------------------------------------------- !
    real(kind=rk) :: fTmpAll_1f(nScalars), fTmpAll_2f(nScalars)
    real(kind=rk) :: fEq_b, fEq_1f
    real(kind=rk) :: velAvg_1f(3), vel_b(3), eqVel_1f(3)
    real(kind=rk) :: usq_1f, ucx_1f, usq_b, ucx_b
    real(kind=rk), dimension(varSys%nStateVars) :: mass_dens_1f, mass_dens_2f, &
      &                                            num_dens_1f
    real(kind=rk) :: resi_coeff(varSys%nStateVars, varSys%nStateVars)
    real(kind=rk) :: mass_dens_b
    real(kind=rk) :: omega_fac, paramBInv
    integer :: iElem, iDir, QQ, iFld, vPos
    integer :: bcVel_pos
    integer :: elemPos, elemOff, posInBuffer, neighPos
    real(kind=rk) :: vel_w(3*globBC%nElems(iLevel)), inv_vel
    type(mus_varSys_data_type), pointer :: fPtr
    integer :: nFields
    ! -------------------------------------------------------------------------
    nFields = varSys%nStateVars
    ! 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:globBC%nElems(iLevel)),      &
      & nVals   = globBC%nElems(iLevel),              &
      & res     = vel_w                               )

    ! convert physical velocity into LB velocity
    inv_vel = 1.0_rk / physics%fac( iLevel )%vel
    vel_w = vel_w * inv_vel

    call C_F_POINTER( varSys%method%val(iField)%method_Data, fPtr )
    ! resitivity coefficient of all species
    do iFld = 1, nFields
      !resi_coeff(iFld, :) = species(iFld)%resi_coeff(:)
      resi_coeff(iFld, :) = fPtr%solverData%scheme                          &
        &                                  %field(iFld)%fieldProp           &
        &                                              %species%resi_coeff(:)
    end do
    ! omega factor of LSE
    omega_fac = mixture%omega_diff * 0.5_rk
    paramBInv = 1.0_rk / mixture%paramB

    associate( auxField => fPtr%solverData%scheme%auxField(iLevel)%val,      &
      &        species => fPtr%solverData%scheme%field(:)%fieldProp%species, &
      &        neighBufferPre => me%neigh(iLevel)%neighBufferPre,            &
      &        stencil => layout%fStencil                                    )

      QQ = layout%fStencil%QQ

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

        ! Position of current boundary element in bcBuffer
        posInBuffer = globBC%elemLvl(iLevel)%posInBcElemBuf%val( iElem )
        elemOff = (posInBuffer-1)*nScalars
        ! Position of 2nd fluid neighbor in state array
        neighPos = me%neigh(iLevel)%posInState(1, iElem)
        mass_dens_2f = 0.0_rk
        do iFld = 1, nFields
          do iDir = 1, QQ
            ! Position of current field variable current dir in state array
            vPos = varSys%method%val(iFld)%state_varPos(iDir)
            ! State array of 1st fluid element
            fTmpAll_1f(vPos) = bcBuffer( elemOff + vPos )
            ! Pre-Collision values of second fluid element
            fTmpAll_2f(vPos) = state(                                       &
              & ?FETCH?(iDir, iFld, neighPos, QQ, nScalars, nSize, neigh) )
            ! mass density of second neighbor
            mass_dens_2f(iFld) = mass_dens_2f(iFld) + fTmpAll_2f(vPos)
          end do
        end do

        ! calculate spc density and velAvg from PDF of first fluid
        call calcDensAndVelsFromPDF( nFields    = nFields,         &
          &                          iField     = iField,          &
          &                          mass_dens  = mass_dens_1f,    &
          &                          num_dens   = num_dens_1f,     &
          &                          velAvg     = velAvg_1f,       &
          &                          eqVel      = eqVel_1f,        &
          &                          varSys     = varSys,          &
          &                          pdf        = fTmpAll_1f,      &
          &                          stencil    = layout%fStencil, &
          &                          species    = species,         &
          &                          resi_coeff = resi_coeff,      &
          &                          omega_fac  = omega_fac,       &
          &                          paramBInv  = paramBInv        )

        ! fluid node velocity square term of equilibrium function
        usq_1f = dot_product( velAvg_1f(1:3), velAvg_1f(1:3) ) * t2cs2inv

        ! density of current species at boundary node extrapolated from
        ! fluid elements
        mass_dens_b = ( 4.0_rk * mass_dens_1f(iField) - mass_dens_2f(iField) ) &
          &         / 3.0_rk

        ! velocity at boundary node for qVal = 1.0
        vel_b = vel_w( (iElem-1)*3 + 1: iElem*3 )

        ! boundary node velocity square term of equilibrium function
        usq_b = dot_product( vel_b(1:3), vel_b(1:3) ) * t2cs2inv

        ! Position of this boundary element in state array
        elemPos = globBC%elemLvl(iLevel)%elem%val(iElem)

        do iDir = 1, layout%fStencil%QQN
          if ( globBC%elemLvl(iLevel)%bitmask%val(iDir, iElem) ) then
            ! Compute equilibrium on fluid node
            ucx_1f = dot_product( stencil%cxDirRK(1:3, iDir), velAvg_1f )

            fEq_1f = layout%weight(iDir) * mass_dens_1f(iField)      &
              &   * ( species(iField)%molWeigRatio + ucx_1f * cs2inv &
              &   + ucx_1f * ucx_1f * t2cs4inv - usq_1f              )

            ! Compute equilibrium on boundary node
            ucx_b = dot_product( stencil%cxDirRK(1:3, iDir), vel_b )

            fEq_b = layout%weight(iDir) * mass_dens_b               &
              &   * ( species(iField)%molWeigRatio + ucx_b * cs2inv &
              &   + ucx_b * ucx_b * t2cs4inv - usq_b                )

            state(?FETCH?(iDir, iField, elemPos, QQ, nScalars, nSize, neigh))  &
              & = fEq_b + ( fTmpAll_1f(varPos(iDir)) - fEq_1f )

          end if
        end do
      end do
    end associate

  end subroutine spc_velocity_noneq_expol
! **************************************************************************** !


! **************************************************************************** !
  !> Mole fraction boundary condition for nonequilibrium extrapolation based.
  !! Default qVal=0.0.
  !! Usage
  !! -----
  !!```lua
  !!boundary_condition = {
  !! { label = 'spc_mole_fraction_noneq_expol',
  !!    kind = 'mole_fraction',
  !!    mole_fraction = 0.1
  !! }
  !!}
  !!```
  !!
  !! This subroutine's interface must match the abstract interface definition
  !! [[boundaryRoutine]] in bc/[[mus_bc_header_module]].f90 in order to be
  !! callable via [[boundary_type:fnct]] function pointer.
  subroutine spc_mole_fraction_noneq_expol( me, state, bcBuffer, globBC, &
    &          levelDesc, tree, nSize, iLevel, sim_time, neigh, layout,  &
    &          fieldProp, varPos, nScalars, varSys, derVarPos, physics,  &
    &          iField, mixture                                           )
    ! -------------------------------------------------------------------- !
    !> global boundary type
    class(boundary_type) :: me
    !> Current state vector of iLevel
    real(kind=rk), intent(inout) :: state(:)
    !> size of state array ( in terms of elements )
    integer, intent(in) :: nSize
    !> state values of boundary elements of all fields of iLevel
    real(kind=rk), intent(in) :: bcBuffer(:)
    !> iLevel descriptor
    type(tem_levelDesc_type), intent(in) :: levelDesc
    !> Treelm Mesh
    type(treelmesh_type), intent(in) :: tree
    !> fluid parameters and properties
    type(mus_field_prop_type), intent(in) :: fieldProp
    !> stencil layout information
    type(mus_scheme_layout_type), intent(in) :: layout
    !> the level On which this boundary was invoked
    integer, intent(in) :: iLevel
    !> connectivity array corresponding to state vector
    integer, intent(in) :: neigh(:)
    !> global time information
    type(tem_time_type), intent(in)  :: sim_time
    !> pointer to field variable in the state vector
    integer, intent(in) :: varPos(:)
    !> number of Scalars in the scheme var system
    integer, intent(in) :: nScalars
    !> scheme variable system
    type(tem_varSys_type), intent(in) :: varSys
    !> position of derived quantities in varsys
    type(mus_derVarPos_type), intent(in) :: derVarPos
    !> scheme global boundary type
    type(glob_boundary_type), intent(in) :: globBC
    !> scheme global boundary type
    type(mus_physics_type), intent(in) :: physics
    !> current field
    integer, intent(in) :: iField
    !> mixture info
    type(mus_mixture_type), intent(in) :: mixture
    ! -------------------------------------------------------------------- !
    real(kind=rk) :: fTmpAll_2f(nScalars)
    real(kind=rk) :: fEq_b, fEq_2f
    real(kind=rk) :: velAvg_2f(3), vel_b(3), eqVel_2f(3)
    real(kind=rk) :: usq_b, ucx_b, usq_2f, ucx_2f
    real(kind=rk), dimension(varSys%nStateVars) :: mass_dens_2f
    real(kind=rk), dimension(varSys%nStateVars) :: num_dens_2f
    real(kind=rk), dimension(3, varSys%nStateVars) :: velocity_2f
    real(kind=rk) :: resi_coeff(varSys%nStateVars, varSys%nStateVars)
    real(kind=rk) :: omega_fac, paramBInv
    integer :: iElem, iDir, QQ, iFld, vPos
    integer :: bcMoleFrac_pos, bcVel_pos
    integer :: elemPos, elemOff, neighPos, posInBuffer
    real(kind=rk) :: moleFrac_w(globBC%nElems(iLevel)), mass_dens_b
    type(mus_varSys_data_type), pointer :: fPtr
    integer :: nFields
    !> Inverse of lattice weight ar restPosition
    real(kind=rk) :: weight0Inv
    ! -------------------------------------------------------------------------
    weight0Inv = 1.0_rk / layout%weight(layout%fStencil%restPosition)
    nFields = varSys%nStateVars

    ! position of molefraction spacetime variable in varSys
    bcMoleFrac_pos = me%bc_states%moleFrac%varPos
    ! mole fraction
    call varSys%method%val(bcMoleFrac_pos)%get_valOfIndex( &
      & varSys  = varSys,                                  &
      & time    = sim_time,                                &
      & iLevel  = iLevel,                                  &
      & idx     = me%bc_states%moleFrac                    &
      &           %pntIndex%indexLvl(iLevel)               &
      &           %val(1:globBC%nElems(iLevel)),           &
      & nVals   = globBC%nElems(iLevel),                   &
      & res     = moleFrac_w                               )

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

    ! resitivity coefficient of all species
    do iFld = 1, nFields
      !resi_coeff(iFld, :) = species(iFld)%resi_coeff(:)
      resi_coeff(iFld, :) = fPtr%solverData%scheme                          &
        &                                  %field(iFld)%fieldProp           &
        &                                              %species%resi_coeff(:)
    end do
    ! omega factor of LSE
    omega_fac = mixture%omega_diff * 0.5_rk
    paramBInv = 1.0_rk / mixture%paramB

    associate( auxField => fPtr%solverData%scheme%auxField(iLevel)%val,      &
      &        species => fPtr%solverData%scheme%field(:)%fieldProp%species, &
      &        stencil => layout%fStencil                                    )

      QQ = layout%fStencil%QQ

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

        ! Position of current boundary element in bcBuffer
        posInBuffer = globBC%elemLvl(iLevel)%posInBcElemBuf%val( iElem )
        elemOff = (posInBuffer-1)*nScalars
        ! Position of neighbor element in state array
        neighPos = me%neigh(iLevel)%posInState(1, iElem)

        do iFld = 1, nFields
          do iDir = 1, QQ
            ! Position of current field variable current dir in state array
            vPos = varSys%method%val(iFld)%state_varPos(iDir)
            ! Pre-Collision values of neighbor fluid element
            fTmpAll_2f(vPos) = state(                                    &
              & ?FETCH?(iDir, iFld, neighPos, QQ, nScalars, nSize, neigh) )
            ! State array of 1st fluid element
            !fTmpAll_1f(vPos) = bcBuffer( elemOff + vPos )
          end do
        end do

        ! calculate spc density and velAvg from PDF of second fluid element
        call calcDensAndVelsFromPDF( nFields    = nFields,         &
          &                          iField     = iField,          &
          &                          mass_dens  = mass_dens_2f,    &
          &                          num_dens   = num_dens_2f,     &
          &                          velAvg     = velAvg_2f,       &
          &                          eqVel      = eqVel_2f,        &
          &                          varSys     = varSys,          &
          &                          pdf        = fTmpAll_2f,      &
          &                          stencil    = layout%fStencil, &
          &                          species    = species,         &
          &                          resi_coeff = resi_coeff,      &
          &                          omega_fac  = omega_fac,       &
          &                          paramBInv  = paramBInv        )

        ! fluid node velocity square term of equilibrium function
        usq_2f = dot_product( velAvg_2f(1:3), velAvg_2f(1:3) ) * t2cs2inv

        ! density of current species at boundary node
        mass_dens_b = mixture%moleDens0LB * moleFrac_w(iElem) &
          &         * species(iField)%molWeight

        ! velocity at boundary node extrapolated from 2nd fluid element
        vel_b = velAvg_2f

        ! boundary node velocity square term of equilibrium function
        usq_b = dot_product( vel_b(1:3), vel_b(1:3) ) * t2cs2inv

        ! Position of this boundary element in state array
        elemPos = globBC%elemLvl(iLevel)%elem%val(iElem)

        do iDir = 1, layout%fStencil%QQN
          if ( globBC%elemLvl(iLevel)%bitmask%val(iDir, iElem) ) then
            ! Compute equilibrium on fluid node
            ucx_2f = dot_product( stencil%cxDirRK(1:3, iDir), velAvg_2f )

            fEq_2f = layout%weight(iDir) * mass_dens_2f(iField)      &
              &   * ( species(iField)%molWeigRatio + ucx_2f * cs2inv &
              &   + ucx_2f * ucx_2f * t2cs4inv - usq_2f              )

            ! Compute equilibrium on boundary node
            ucx_b = dot_product( stencil%cxDirRK(1:3, iDir), vel_b )

            fEq_b = layout%weight(iDir) * mass_dens_b               &
              &   * ( species(iField)%molWeigRatio + ucx_b * cs2inv &
              &   + ucx_b * ucx_b * t2cs4inv - usq_b                )

            ! update incomping link
            state(?FETCH?(iDir, iField, elemPos, QQ, nScalars, nSize, neigh))  &
              & = fEq_b + ( fTmpAll_2f(varPos(iDir)) - fEq_2f )
          end if
        end do
      end do
    end associate

  end subroutine spc_mole_fraction_noneq_expol
! **************************************************************************** !

! **************************************************************************** !
  !> Open outflow boundary condition based on nonequilibrium extrapolation
  !! method. Default qVal = 0.0
  !! Usage
  !! -----
  !!```lua
  !!boundary_condition = {
  !! { label = 'outlet',
  !!    kind = 'spc_outflow',
  !! }
  !!}
  !!```
  !!
  !! 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_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) :: fTmpAll_2f(nScalars), fTmpAll_3f(nScalars)
    real(kind=rk) :: fEq_b, fEq_2f
    real(kind=rk) :: velAvg_2f(3), vel_b(3), eqVel_2f(3)
    real(kind=rk) :: velAvg_3f(3), eqVel_3f(3)
    real(kind=rk) :: usq_2f, ucx_2f, usq_b, ucx_b
    real(kind=rk), dimension(varSys%nStateVars) :: mass_dens_2f, mass_dens_3f
    real(kind=rk), dimension(varSys%nStateVars) :: num_dens_2f
    real(kind=rk) :: resi_coeff(varSys%nStateVars, varSys%nStateVars)
    real(kind=rk) :: omega_fac, paramBInv, mass_dens_b
    integer :: iElem, iDir, QQ, iFld, vPos
    integer :: bcMoleFrac_pos, bcVel_pos
    integer :: elemPos, elemOff, neighPos_2f, neighPos_3f
    type(mus_varSys_data_type), pointer :: fPtr
    integer :: nFields
    ! -------------------------------------------------------------------------
    nFields = varSys%nStateVars

    call C_F_POINTER( varSys%method%val(iField)%method_Data, fPtr )
    ! resitivity coefficient of all species
    do iFld = 1, nFields
      !resi_coeff(iFld, :) = species(iFld)%resi_coeff(:)
      resi_coeff(iFld, :) = fPtr%solverData%scheme                          &
        &                                  %field(iFld)%fieldProp           &
        &                                              %species%resi_coeff(:)
    end do
    ! omega factor of LSE
    omega_fac = mixture%omega_diff * 0.5_rk
    paramBInv = 1.0_rk / mixture%paramB

    associate( auxField => fPtr%solverData%scheme%auxField(iLevel)%val,      &
      &        species => fPtr%solverData%scheme%field(:)%fieldProp%species, &
      &        neighBufferPre => me%neigh(iLevel)%neighBufferPre,            &
      &        stencil => layout%fStencil                                    )

      QQ = layout%fStencil%QQ

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

        ! Position of second neighbor element in state array
        neighPos_2f = me%neigh(iLevel)%posInState(1, iElem)
        ! Position of third neighbor element in state array
        neighPos_3f = me%neigh(iLevel)%posInState(2, iElem)
        ! Compute only mass density from third neighbor
        mass_dens_3f = 0.0_rk
        do iFld = 1, nFields
          do iDir = 1, QQ
            ! Position of current field variable current dir in state array
            vPos = varSys%method%val(iFld)%state_varPos(iDir)
            ! Pre-Collision values of second fluid element
            fTmpAll_2f(vPos) = state(                                        &
              & ?FETCH?(iDir, iFld, neighPos_2f, QQ, nScalars, nSize, neigh) )
            ! Pre-Collision values of third fluid element
            fTmpAll_3f(vPos) = state(                                        &
              & ?FETCH?(iDir, iFld, neighPos_3f, QQ, nScalars, nSize, neigh) )
            ! mass density of third neighbor
            mass_dens_3f(iFld) = mass_dens_3f(iFld) + fTmpAll_3f(vPos)
          end do
        end do

        ! calculate spc density and velAvg from PDF of second fluid
        call calcDensAndVelsFromPDF( nFields    = nFields,         &
          &                          iField     = iField,          &
          &                          mass_dens  = mass_dens_2f,    &
          &                          num_dens   = num_dens_2f,     &
          &                          velAvg     = velAvg_2f,       &
          &                          eqVel      = eqVel_2f,        &
          &                          varSys     = varSys,          &
          &                          pdf        = fTmpAll_2f,      &
          &                          stencil    = layout%fStencil, &
          &                          species    = species,         &
          &                          resi_coeff = resi_coeff,      &
          &                          omega_fac  = omega_fac,       &
          &                          paramBInv  = paramBInv        )

        ! fluid node velocity square term of equilibrium function
        usq_2f = dot_product( velAvg_2f(1:3), velAvg_2f(1:3) ) * t2cs2inv

        ! second-order extrapolation of density of current species at boundary
        ! node from second and third fluid element
        mass_dens_b = ( 4.0_rk * mass_dens_2f(iField) - mass_dens_3f(iField) ) &
          &         / 3.0_rk

        ! first-order extrapolation of velocity at boundary node from second
        ! fluid element
        vel_b = velAvg_2f
        ! second-order extrapolation of velocity at boundary node
        ! vel_b =  ( 4.0_rk * velAvg_2f - velAvg_3f ) / 3.0_rk

        ! boundary node velocity square term of equilibrium function
        usq_b = dot_product( vel_b(1:3), vel_b(1:3) ) * t2cs2inv

        ! Position of this boundary element in state array
        elemPos = globBC%elemLvl(iLevel)%elem%val(iElem)

        do iDir = 1, layout%fStencil%QQN
          if ( globBC%elemLvl(iLevel)%bitmask%val(iDir, iElem) ) then
            ! Compute equilibrium on second fluid node
            ucx_2f = dot_product( stencil%cxDirRK(1:3, iDir), velAvg_2f )

            fEq_2f = layout%weight(iDir) * mass_dens_2f(iField)      &
              &   * ( species(iField)%molWeigRatio + ucx_2f * cs2inv &
              &   + ucx_2f * ucx_2f * t2cs4inv - usq_2f              )

            ! Compute equilibrium on boundary node
            ucx_b = dot_product( stencil%cxDirRK(1:3, iDir), vel_b )

            fEq_b = layout%weight(iDir) * mass_dens_b               &
              &   * ( species(iField)%molWeigRatio + ucx_b * cs2inv &
              &   + ucx_b * ucx_b * t2cs4inv - usq_b                )

            ! update incomping direction
            state(?FETCH?(iDir, iField, elemPos, QQ, nScalars, nSize, neigh))  &
              & = fEq_b + ( fTmpAll_2f(varPos(iDir)) - fEq_2f )
          end if
        end do
      end do
    end associate

  end subroutine spc_outflow
! **************************************************************************** !

! **************************************************************************** !
  !> Open outflow boundary condition for solvent based on nonequilibrium
  !! extrapolation. total moledens at boundary is enforced.
  !! method. Default qVal = 0.0
  !! Usage
  !! -----
  !!```lua
  !!boundary_condition = {
  !! { label = 'outlet',
  !!    kind = 'spc_solvent_outflow',
  !! }
  !!}
  !!```
  !!
  !! 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_solvent_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) :: fTmpAll_2f(nScalars), fTmpAll_3f(nScalars)
    real(kind=rk) :: fEq_b, fEq_2f
    real(kind=rk) :: velAvg_2f(3), vel_b(3), eqVel_2f(3)
    real(kind=rk) :: velAvg_3f(3), eqVel_3f(3)
    real(kind=rk) :: usq_2f, ucx_2f, usq_b, ucx_b
    real(kind=rk), dimension(varSys%nStateVars) :: mass_dens_2f, mass_dens_3f
    real(kind=rk), dimension(varSys%nStateVars) :: num_dens_2f, num_dens_3f
    real(kind=rk) :: resi_coeff(varSys%nStateVars, varSys%nStateVars)
    real(kind=rk) :: omega_fac, paramBInv, mass_dens_b, num_dens_b, &
      &              tot_NumDens_b
    integer :: iElem, iDir, QQ, iFld, vPos
    integer :: bcMoleFrac_pos, bcVel_pos
    integer :: elemPos, elemOff, neighPos_2f, neighPos_3f
    type(mus_varSys_data_type), pointer :: fPtr
    integer :: nFields
    ! -------------------------------------------------------------------------
    nFields = varSys%nStateVars

    call C_F_POINTER( varSys%method%val(iField)%method_Data, fPtr )
    ! resitivity coefficient of all species
    do iFld = 1, nFields
      !resi_coeff(iFld, :) = species(iFld)%resi_coeff(:)
      resi_coeff(iFld, :) = fPtr%solverData%scheme                          &
        &                                  %field(iFld)%fieldProp           &
        &                                              %species%resi_coeff(:)
    end do
    ! omega factor of LSE
    omega_fac = mixture%omega_diff * 0.5_rk
    paramBInv = 1.0_rk / mixture%paramB

    associate( auxField => fPtr%solverData%scheme%auxField(iLevel)%val,      &
      &        species => fPtr%solverData%scheme%field(:)%fieldProp%species, &
      &        neighBufferPre => me%neigh(iLevel)%neighBufferPre,            &
      &        stencil => layout%fStencil                                    )

      QQ = layout%fStencil%QQ

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

        ! Position of second neighbor element in state array
        neighPos_2f = me%neigh(iLevel)%posInState(1, iElem)
        ! Position of third neighbor element in state array
        neighPos_3f = me%neigh(iLevel)%posInState(2, iElem)
        ! Compute only mass density from third neighbor
        mass_dens_3f = 0.0_rk
        do iFld = 1, nFields
          do iDir = 1, QQ
            ! Position of current field variable current dir in state array
            vPos = varSys%method%val(iFld)%state_varPos(iDir)
            ! Pre-Collision values of second fluid element
            fTmpAll_2f(vPos) = state(                                       &
              & ?FETCH?(iDir, iFld, neighPos_2f, QQ, nScalars, nSize, neigh) )
            ! Pre-Collision values of third fluid element
            fTmpAll_3f(vPos) = state(                                       &
              & ?FETCH?(iDir, iFld, neighPos_3f, QQ, nScalars, nSize, neigh) )
            ! mass density of third neighbor
            mass_dens_3f(iFld) = mass_dens_3f(iFld) + fTmpAll_3f(vPos)
          end do
        end do

        ! calculate spc density and velAvg from PDF of second fluid
        call calcDensAndVelsFromPDF( nFields    = nFields,         &
          &                          iField     = iField,          &
          &                          mass_dens  = mass_dens_2f,    &
          &                          num_dens   = num_dens_2f,     &
          &                          velAvg     = velAvg_2f,       &
          &                          eqVel      = eqVel_2f,        &
          &                          varSys     = varSys,          &
          &                          pdf        = fTmpAll_2f,      &
          &                          stencil    = layout%fStencil, &
          &                          species    = species,         &
          &                          resi_coeff = resi_coeff,      &
          &                          omega_fac  = omega_fac,       &
          &                          paramBInv  = paramBInv        )

        ! fluid node velocity square term of equilibrium function
        usq_2f = dot_product( velAvg_2f(1:3), velAvg_2f(1:3) ) * t2cs2inv

        ! Second-order extrapolation of number density of current species
        num_dens_b = ( 4.0_rk * num_dens_2f(iField) - num_dens_3f(iField) ) &
          &         / 3.0_rk
        ! Second-order extrapolation of numerical total number density
        tot_NumDens_b = ( 4.0_rk * sum(num_dens_2f) - sum(num_dens_3f) ) &
          &           / 3.0_rk
        ! Compute number density of current species at boundary element
        num_dens_b = mixture%moleDens0LB - ( tot_NumDens_b - num_dens_b )
        ! Mass density of current species at boundary element
        mass_dens_b = num_dens_b * species(iField)%molWeight

        ! first-order extrapolation of velocity at boundary node from second
        ! fluid element
        vel_b = velAvg_2f
        ! second-order extrapolation of velocity at boundary node
        ! vel_b =  ( 4.0_rk * velAvg_2f - velAvg_3f ) / 3.0_rk

        ! boundary node velocity square term of equilibrium function
        usq_b = dot_product( vel_b(1:3), vel_b(1:3) ) * t2cs2inv

        ! Position of this boundary element in state array
        elemPos = globBC%elemLvl(iLevel)%elem%val(iElem)

        do iDir = 1, layout%fStencil%QQN
          if ( globBC%elemLvl(iLevel)%bitmask%val(iDir, iElem) ) then
            ! Compute equilibrium on second fluid node
            ucx_2f = dot_product( stencil%cxDirRK(1:3, iDir), velAvg_2f )

            fEq_2f = layout%weight(iDir) * mass_dens_2f(iField)      &
              &   * ( species(iField)%molWeigRatio + ucx_2f * cs2inv &
              &   + ucx_2f * ucx_2f * t2cs4inv - usq_2f              )

            ! Compute equilibrium on boundary node
            ucx_b = dot_product( stencil%cxDirRK(1:3, iDir), vel_b )

            fEq_b = layout%weight(iDir) * mass_dens_b               &
              &   * ( species(iField)%molWeigRatio + ucx_b * cs2inv &
              &   + ucx_b * ucx_b * t2cs4inv - usq_b                )

            ! update incomping direction
            state(?FETCH?(iDir, iField, elemPos, QQ, nScalars, nSize, neigh))  &
              & = fEq_b + ( fTmpAll_2f(varPos(iDir)) - fEq_2f )
          end if
        end do
      end do
    end associate

  end subroutine spc_solvent_outflow
! **************************************************************************** !

! ****************************************************************************** !
  !> Inflow boundary condition based on non-Equilbrium extrapolation method.
  !! Usage
  !! -----
  !!```lua
  !!boundary_condition = {
  !! { label = 'outlet',
  !!    kind = 'spc_inflow',
  !!    mole_fraction = 0.01,
  !!    velocity = {0.1,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_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) :: fTmpAll_f(nScalars), fTmpAll_ff(nScalars)
    real(kind=rk) :: pdf_f, fEq_b, fEq_f, fEq_ff, fnEq_f, fnEq_ff, fnEq_b
    real(kind=rk) :: eqVel_f(3), velAvg_f(3), eqVel_ff(3), velAvg_ff(3), &
      &              vel_b(3)
    real(kind=rk) :: usq_f, ucx_f, ucxQuad_f, usq_b, ucx_b, usq_ff, ucx_ff, &
      &              ucxQuad_ff
    real(kind=rk), dimension(varSys%nStateVars) :: mass_dens_f, mass_dens_ff
    real(kind=rk), dimension(3, varSys%nStateVars) :: velocity_f, velocity_ff
    real(kind=rk) :: resi_coeff(varSys%nStateVars, varSys%nStateVars)
    real(kind=rk) :: omega_fac, paramBInv
    integer :: iElem, iDir, QQ, iFld, vPos, posInBuffer
    integer :: bcMoleFrac_pos, bcVel_pos
    integer :: elemPos, elemOff, neighPos_1, neighPos_2
    real(kind=rk) :: vel_w(3*globBC%nElems(iLevel)), inv_vel
    real(kind=rk) :: moleFrac_w(globBC%nElems(iLevel)), mass_dens_b
    type(mus_varSys_data_type), pointer :: fPtr
    real(kind=rk) :: coeff_w, coeff_f, coeff_ff
    integer :: nFields
    !> Inverse of lattice weight ar restPosition
    real(kind=rk) :: weight0Inv
    ! ------------------------------------------------------------------------
    weight0Inv = 1.0_rk / layout%weight(layout%fStencil%restPosition)
    nFields = varSys%nStateVars
    ! coefficient for extrapolation of macroscopic quantities for qVal=0.5
    coeff_w = 5.0_rk/3.0_rk
    coeff_f = -0.5_rk
    coeff_ff = -1.0_rk/6.0_rk
    ! position of molefraction spacetime variable in varSys
    bcMoleFrac_pos = me%bc_states%moleFrac%varPos
    ! mole fraction
    call varSys%method%val(bcMoleFrac_pos)%get_valOfIndex( &
      & varSys  = varSys,                                  &
      & time    = sim_time,                                &
      & iLevel  = iLevel,                                  &
      & idx     = me%bc_states%moleFrac                    &
      &           %pntIndex%indexLvl(iLevel)               &
      &           %val(1:globBC%nElems(iLevel)),           &
      & nVals   = globBC%nElems(iLevel),                   &
      & res     = moleFrac_w                               )

    ! 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:globBC%nElems(iLevel)),      &
      & nVals   = globBC%nElems(iLevel),              &
      & res     = vel_w                               )

    ! convert physical velocity into LB velocity
    inv_vel = 1.0_rk / physics%fac( iLevel )%vel
    vel_w = vel_w * inv_vel

    call C_F_POINTER( varSys%method%val(iField)%method_Data, fPtr )
    ! resitivity coefficient of all species
    do iFld = 1, nFields
      !resi_coeff(iFld, :) = species(iFld)%resi_coeff(:)
      resi_coeff(iFld, :) = fPtr%solverData%scheme                          &
        &                                  %field(iFld)%fieldProp           &
        &                                              %species%resi_coeff(:)
    end do
    ! omega factor of LSE
    omega_fac = mixture%omega_diff * 0.5_rk
    paramBInv = 1.0_rk / mixture%paramB

    associate( auxField => fPtr%solverData%scheme%auxField(iLevel)%val,      &
      &        species => fPtr%solverData%scheme%field(:)%fieldProp%species, &
      &        stencil => layout%fStencil                                    )

      QQ = layout%fStencil%QQ

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

        ! Position of current boundary element in bcBuffer
        posInBuffer = globBC%elemLvl(iLevel)%posInBcElemBuf%val( iElem )
        ! Position of neighbor element in state array
        neighPos_1 = me%neigh(iLevel)%posInState(1, iElem)
        neighPos_2 = me%neigh(iLevel)%posInState(2, iElem)

        do iFld = 1, nFields
          do iDir = 1, QQ
            ! Position of current field variable current dir in state array
            vPos = varSys%method%val(iFld)%state_varPos(iDir)
            ! State array of current element
            fTmpAll_f(vPos) = bcBuffer( (posInBuffer-1)*nScalars + vPos )
            !fTmpAll_f(vPos) = state(                                    &
            !  & ?FETCH?(iDir, iFld, neighPos_1, QQ, nScalars, nSize, neigh) )
            ! Post-Collision values of neighbor element
            fTmpAll_ff(vPos) = state(                                    &
              & ?FETCH?(iDir, iFld, neighPos_1, QQ, nScalars, nSize, neigh) )
          end do
        end do

        ! calculate spc density, velocity, and velAvg from PDF of first fluid
        call calcMacrosFromPDF( varSys     = varSys,      &
          &                     pdf        = fTmpAll_f,   &
          &                     mass_dens  = mass_dens_f, &
          &                     velocity   = velocity_f,  &
          &                     velAvg     = velAvg_f,    &
          &                     eqVel      = eqVel_f,     &
          &                     species    = species      )


        ! calculate spc density, velocity, and velAvg from PDF of second fluid
        call calcMacrosFromPDF( varSys     = varSys,       &
          &                     pdf        = fTmpAll_ff,   &
          &                     mass_dens  = mass_dens_ff, &
          &                     velocity   = velocity_ff,  &
          &                     velAvg     = velAvg_ff,    &
          &                     eqVel      = eqVel_ff,     &
          &                     species    = species       )

        ! fluid node velocity square term of equilibrium function
        usq_f = dot_product( velAvg_f(1:3), velAvg_f(1:3) ) * t2cs2inv
        usq_ff = dot_product( velAvg_ff(1:3), velAvg_ff(1:3) ) * t2cs2inv

        ! density of current species at boundary node
        !mass_dens_b = (4.0_rk * mass_dens_f(iField) - mass_dens_ff(iField)) / 3.0_rk
        !if (mass_dens_b < 0.0_rk) then
        !  mass_dens_b = mixture%moleDens0LB * 0.001 &
        !    &         * species(iField)%molWeight
        !end if
        mass_dens_b = mixture%moleDens0LB * moleFrac_w(iElem) &
          &         * species(iField)%molWeight
        !mass_dens_b = ( 4.0_rk * mass_dens_f(iField) - mass_dens_ff(iField) ) &
        !  &         / 3.0_rk
        !mass_dens_b = mixture%rho0LB - (sum(mass_dens_f) - mass_dens_b)
        !mass_dens_b = 3.0_rk * mixture%moleDens0LB * moleFrac_w(iElem) &
        !  &         * species(iField)%molWeight                         &
        !  &         - mass_dens_f(iField) - mass_dens_ff(iField)
        !mass_dens_b = coeff_w * mixture%moleDens0LB * moleFrac_w(iElem) &
        !  &         * species(iField)%molWeight                         &
        !  &         + coeff_f * mass_dens_f(iField)                     &
        !  &         + coeff_ff * mass_dens_ff(iField)

        ! velocity at boundary node
        vel_b = vel_w( (iElem-1)*3 + 1: iElem*3 )
        !vel_b = velocity_ff(1:3, iField)
        !vel_b = coeff_w * vel_w( (iElem-1)*3 + 1: iElem*3 ) &
        !  &   + coeff_f * velocity_f(1:3, iField)           &
        !  &   + coeff_ff * velocity_ff(1:3, iField)
        !vel_b = coeff_w * vel_w( (iElem-1)*3 + 1: iElem*3 ) &
        !  &   + coeff_f * velAvg_f(1:3) + coeff_ff * velAvg_ff(1:3)

!write(dbgUnit(1),*) 'vel_w ', vel_w( (iElem-1)*3 + 1: iElem*3 )
!write(dbgUnit(1),*) 'mass_dens_b ', mass_dens_b, ' vel_b ', vel_b

        ! boundary node velocity square term of equilibrium function
        usq_b = dot_product( vel_b(1:3), vel_b(1:3) ) * t2cs2inv

        ! Position of this boundary element in state array
        elemPos = globBC%elemLvl(iLevel)%elem%val(iElem)

        do iDir = 1, layout%fStencil%QQN
          if ( globBC%elemLvl(iLevel)%bitmask%val(iDir, iElem) ) then
            ! Compute equilibrium on fluid node
            ucx_f = dot_product( stencil%cxDirRK(1:3, iDir), velAvg_f )
            ucxQuad_f = dot_product( stencil%cxDirRK(1:3, iDir), velAvg_f(1:3) )

            fEq_f = layout%weight(iDir) * mass_dens_f(iField)       &
              &   * ( species(iField)%molWeigRatio + ucx_f * cs2inv &
              &   + ucxQuad_f * ucxQuad_f * t2cs4inv - usq_f        )

            ! Compute equilibrium on fluid node
            ucx_ff = dot_product( stencil%cxDirRK(1:3, iDir), velAvg_ff )
            ucxQuad_ff = dot_product( stencil%cxDirRK(1:3, iDir), &
              &                       velAvg_ff(1:3) )

            fEq_ff = layout%weight(iDir) * mass_dens_ff(iField)      &
              &   * ( species(iField)%molWeigRatio + ucx_ff * cs2inv &
              &   + ucxQuad_ff * ucxQuad_ff * t2cs4inv - usq_ff      )

            ! Compute equilibrium on boundary node
            ucx_b = dot_product( stencil%cxDirRK(1:3, iDir), vel_b )

            fEq_b = layout%weight(iDir) * mass_dens_b               &
              &   * ( species(iField)%molWeigRatio + ucx_b * cs2inv &
              &   + ucx_b * ucx_b * t2cs4inv - usq_b                )

            ! Position of current field variable current dir in state array
            !pdf_f = bcBuffer( (posInBuffer-1)*nScalars + varPos(iDir) )

            fnEq_f = fTmpAll_f(varPos(iDir)) - fEq_f
            fnEq_ff = fTmpAll_ff(varPos(iDir)) - fEq_ff
            fnEq_b = 0.5_rk * fnEq_f + 0.5_rk * fnEq_ff

            state(?FETCH?(iDir, iField, elemPos, QQ, nScalars, nSize, neigh))  &
              !& = fEq_b + ( pdf_f - fEq_f )
              !& = fEq_b + fnEq_b
              & = fEq_b + fnEq_f
          end if
        end do

        !iDir = layout%fStencil%restPosition
        !fEq_b = layout%weight(iDir) * mass_dens_b * ( weight0Inv             &
        ! &    + (1.0_rk - weight0Inv) * species(iField)%molWeigRatio - usq_b )

        !state(?FETCH?(iDir, iField, elemPos, QQ, nScalars, nSize, neigh)) &
        !  & = fEq_b
      end do
    end associate
!    stop

  contains
    !> This routine computes macroscopic quantities like density, velocity
    !! equilibrium velocity, velocity average for given iField from PDF
    subroutine calcMacrosFromPDF( varSys, pdf, mass_dens, velocity, velAvg, &
      &                           eqVel, species                            )
      ! -----------------------------------------------------------------------
      type(tem_varSys_type), intent(in) :: varSys
      real(kind=rk), intent(in) :: pdf(varSys%nScalars)
      real(kind=rk), intent(out) :: mass_dens(nFields)
      real(kind=rk), intent(out) :: velocity(nFields, nFields)
      real(kind=rk), intent(out) :: velAvg(3)
      real(kind=rk), intent(out) :: eqVel(3)
      type(mus_species_type), intent(in) :: species(nFields)
      ! -----------------------------------------------------------------------
      integer :: iFld, iDir, vPos
      real(kind=rk) :: fTmp(layout%fStencil%QQ)
      real(kind=rk) :: tot_numDens, tot_massDens
      real(kind=rk) :: momentum(3, nFields), first_moments(3, nFields)
      real(kind=rk) :: num_dens(nFields), moleFrac(nFields)
      ! -----------------------------------------------------------------------
      do iFld = 1, nFields
        do iDir = 1, layout%fStencil%QQ
          ! Position of current field variable current dir in state array
          vPos = varSys%method%val(iFld)%state_varPos(iDir)
          ! State array of current element
          fTmp(iDir) = pdf(vPos)
        end do
        ! species density
        mass_dens(iFld) = sum(fTmp)
        ! number density
        num_dens(iFld) = mass_dens(iFld) * species(iFld)%molWeightInv
        ! first moments: momentum of transformed PDF
        first_moments(1, iFld) = sum( fTmp * layout%fStencil%cxDirRK(1, :) )
        first_moments(2, iFld) = sum( fTmp * layout%fStencil%cxDirRK(2, :) )
        first_moments(3, iFld) = sum( fTmp * layout%fStencil%cxDirRK(3, :) )
      end do

      !total number density
      tot_numDens = sum(num_dens)
      !total mass density
      tot_massDens = sum(mass_dens)
      ! mole fraction
      moleFrac(:) = num_dens(:) / tot_numDens

      ! actual momentum of all species obtained from solving
      ! Linear Equation system
      momentum = momentumFromMacroLSE(                      &
        &          moleFraction  = moleFrac,                &
        &          first_moments = first_moments,           &
        &          nFields       = nFields,                 &
        &          phi           = species(:)%molWeigRatio, &
        &          resi_coeff    = resi_coeff,              &
        &          omega_fac     = omega_fac,               &
        &          paramBInv     = paramBInv                )

      ! velocity of all species on current fluid node
      do iFld = 1, nFields
        velocity(:, iFld) = momentum(:, iFld) / mass_dens(iFld)
      end do

      ! equilibrium velocity for current species
      eqVel(:) = velocity(:, iField)
      do iFld = 1, nFields
        eqVel(:) = eqVel(:) + resi_coeff(iField, iFld)           &
          &      * species(iField)%molWeigRatio * moleFrac(iFld) &
          &      * ( velocity(:, iFld) - velocity(:, iField) )   &
          &      * paramBInv
      end do

      ! Mass averaged mixture velocity
      velAvg(1) = dot_product( mass_dens, velocity(1, :) ) / tot_massDens
      velAvg(2) = dot_product( mass_dens, velocity(2, :) ) / tot_massDens
      velAvg(3) = dot_product( mass_dens, velocity(3, :) ) / tot_massDens

    end subroutine calcMacrosFromPDF

  end subroutine spc_inflow
! ****************************************************************************** !



! ****************************************************************************** !
  !> molar diffusion flux boundary condition
  !! Usage
  !! -----
  !!```lua
  !!boundary_condition = {
  !! { label = 'outlet',
  !!    kind = 'spc_moleDiff_flux',
  !!    mole_diff_flux = 'mole_diff_flux'
  !! }
  !!}
  !!variable = {
  !!  name = 'mole_diff_flux',
  !!  ncomponents = 3,
  !!  vartype = 'st_fun',
  !!  st_fun = {0.06, 0.0, 0.0}
  !!}
  !!```
  !! @todo KM Adapt moleDiff flux boundary input_loc for equilibrium function
  !! pointer
  !!
  !! 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_moleDiff_Flux( 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 )
    !first moments of nSpecies
    real(kind=rk) :: first_moments( 3, varSys%nStateVars )
    real(kind=rk) :: momentum( 3, varSys%nStateVars )
    real(kind=rk) :: vel( 3, varSys%nStateVars )
    real(kind=rk) :: mixMolvel( 3 )
    real(kind=rk) :: mass_dens( varSys%nStateVars ), &
      &              num_dens( varSys%nStateVars )
    real(kind=rk) :: tot_NumDens, tot_MassDens
    real(kind=rk) :: moleDiff_Flux(globBC%nElems(iLevel)*3), inv_flux
    real(kind=rk) :: moleFrac(varSys%nStateVars)
    integer :: iElem, iDir, iFieldLoc, nFields, iComp, pos, QQ
    integer :: bcMoleDiffFlux_pos, posInBuffer, posInState
    real(kind=rk) :: molWeightInv(varSys%nStateVars), phi(varSys%nStateVars)
    real(kind=rk) :: resi_coeff( varSys%nStateVars, varSys%nStateVars )
    real(kind=rk) :: paramBInv, omega_fac
    type(mus_varSys_data_type), pointer :: fPtr
    !------------------------------------------------------------------------
    QQ = layout%fStencil%QQ
    nFields = varSys%nStateVars

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

    do iFieldLoc = 1, nFields
      ! species properties
      ! molecular weight inverse
      molWeightInv(iFieldLoc) = fPtr%solverData%scheme%field(iFieldLoc) &
        &                           %fieldProp%species%molWeightInv
      ! 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

    omega_fac = mixture%omega_diff * 0.5_rk
    paramBInv = 1.0_rk/mixture%paramB

    inv_flux = 1.0_rk / physics%fac( iLevel )%moleFlux
    ! position of boundary moleflux in varSys
    bcMoleDiffFlux_pos = me%bc_states%moleDiff_flux%varPos
    ! Get moleflux
    call varSys%method%val(bcMoleDiffFlux_pos)%get_valOfIndex( &
      & varSys  = varSys,                                      &
      & time    = sim_time,                                    &
      & iLevel  = iLevel,                                      &
      & idx     = me%bc_states%moleDiff_flux                   &
      &           %pntIndex%indexLvl(iLevel)                   &
      &           %val(1:globBC%nElems(iLevel)),               &
      & nVals   = globBC%nElems(iLevel),                       &
      & res     = moleDiff_Flux                                )

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

    ! Calculate the mole averaged mixture velocity w = sum(n_i v_i)/n
    do iElem = 1, globBC%nElems(iLevel)
      ! initialize the local velocity for every element
      first_moments = 0.0_rk
      momentum = 0.0_rk
      mass_dens = 0.0_rk
      num_dens = 0.0_rk
      posInBuffer = globBC%elemLvl( iLevel )%posInBcElemBuf%val( iElem )
      do iFieldLoc = 1, nFields
        ! all fields have some nComponents
        do iDir = 1,layout%fStencil%QQ
          ! get the position of the current direction of the depending variable
          ! in the global system (= position in the input)
          pos = varSys%method%val(iFieldLoc)%state_varPos(iDir)
          ! density
          mass_dens(iFieldLoc) = mass_dens(iFieldLoc) + bcBuffer(       &
            & ?IDX?( pos, posInBuffer, nScalars,globBC%nElems(iLevel) ) )

          first_moments( :, iFieldLoc ) = first_moments( :, iFieldLoc )   &
& + bcBuffer( ?IDX?( pos, posInBuffer, nScalars,globBC%nElems(iLevel) ) ) &
              & * layout%fStencil%cxDirRK( :, iDir )
        end do !iDir
        ! number density
        num_dens(iFieldLoc) = mass_dens(iFieldLoc) * molWeightInv(iFieldLoc)
      end do !iField

      !total number density
      tot_NumDens = sum(num_dens)
      !total mass density
      tot_massDens = sum(mass_dens)
      ! mole fraction
      moleFrac(:) = num_dens(:) / tot_NumDens

      ! momentum of all species
      momentum = momentumFromMacroLSE( moleFraction  = moleFrac,       &
        &                              first_moments = first_moments,  &
        &                              nFields       = nFields,        &
        &                              phi           = phi,            &
        &                              resi_coeff    = resi_coeff,     &
        &                              omega_fac     = omega_fac,      &
        &                              paramBInv     = paramBInv       )

      !velocity of all species
      do iFieldLoc = 1,nFields
        vel( :, iFieldLoc) = momentum( :, iFieldLoc) / mass_dens(iFieldLoc)
      end do

      !mole averaged mixture velocity, w = sum(n_i*v_i)/n_T
      do iComp = 1, 3
        mixMolVel(iComp) = dot_product( num_dens, vel(iComp,:) )/tot_numDens
      end do

      ! update current species velocity with defined mole diffusion flux
      ! species velocity at boundary element
      ! moleDiff_flux, J_i = n_i(v_i - w )
      ! v_i = J_i/n_i + w
      vel(:, iField) = moleDiff_Flux((iElem-1)*3+1:iElem*3) &
        &            / num_dens(iField)                     &
        &            + mixMolVel

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

      posInState = globBC%elemLvl( iLevel )%elem%val( iElem )
      do iDir = 1,layout%fStencil%QQN
        if( globBC%elemLvl(iLevel)%bitmask%val( iDir, iElem )) then
          state(?FETCH?( iDir, iField, posInState, QQ, nScalars, nSize,neigh )) &
            & = fEq( iDir )
        end if
      end do
    end do

  end subroutine spc_moleDiff_Flux
! ****************************************************************************** !


! ****************************************************************************** !
  !> This routine computes mole diffusion flux of the ionic species at the
  !! membrance using black box model and then mass density at the membrane
  !! boundary from mole diffusion flux. Then equilibrium is set at the boundary
  !! which is computed from mass density and velocity
  !!
  !! Black box model: \n
  !!  - for ionic species mole diffusion flux:
  !! $J^m_k (x^m,t) = \frac{t^m_k(x^m,t)}{z_k F} i^m(x^m,t)$
  !! $J^m_k$ - mole diffusion flux of ion in the membrane
  !! $t^m_k$ - transference number of ion of the membrane
  !!           (selective membrane property)
  !! $z_k$ - charge number of ion
  !! $F$ - Faraday constant
  !! $i^m(x^m,t) = n^l_x_1 \cdot i^l(x,t) = - n^r_x_1 \cdot i^r(x,t)$
  !!    - orthogonal projection of the current density at the left $i^l$ and
  !!      right $i^r$ side of the membrance.
  !! The normal vectors $n^l_1$ and $n^r_1$ point into the membrane surface
  !! $i(x,t) = F \sum_{j=1}^{n} z_j J^e_j(x,t)$
  !!    - local current density [A/m^3]
  !! $J^e_k = c^e_k(v^e_k - v) $ - mole diffusion flux of ion in the electrolyte
  !! In black box model: mixture averaged velocity at the membrane is assumed to
  !! be zero
  !!
  !! 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_blackbox_mem_ion( 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_all( layout%fStencil%QQ*varSys%nStateVars )
    real(kind=rk) :: mass_dens( varSys%nStateVars )
    real(kind=rk) :: num_dens( varSys%nStateVars )
    integer :: iElem, iDir, iFieldLoc, nFields, pos, QQ
    integer :: posInBuffer, posInState
    real(kind=rk) :: molWeight(varSys%nStateVars), chargeNr(varSys%nStateVars)
    real(kind=rk), dimension(3, varSys%nStateVars):: molar_flux, vel, momentum
    real(kind=rk) :: loc_curr_dens(3)
    real(kind=rk) :: mem_curr_dens, mem_molar_flux
    real(kind=rk) :: paramBInv, omega_fac
    real(kind=rk) :: molefrac(3), mem_mass_flux, elec_mass_flux(3)
    real(kind=rk) :: first_moments( 3, varSys%nStateVars )
    real(kind=rk) :: tot_NumDens, tot_MassDens
    real(kind=rk) :: molWeightInv(varSys%nStateVars), phi(varSys%nStateVars)
    real(kind=rk) :: resi_coeff( varSys%nStateVars, varSys%nStateVars )
    type(mus_varSys_data_type), pointer :: fPtr
    ! --------------------------------------------------------------------------
!write(dbgUnit(1),*) 'Boundary label ', trim(me%label)
!write(dbgUnit(1),*) 'iField ', iField
    QQ = layout%fStencil%QQ
    nFields = varSys%nStateVars

    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
      molWeightInv(iFieldLoc) =  fPtr%solverData%scheme%field(iFieldLoc)   &
        &                                    %fieldProp%species%molWeightInv
      ! charge number
      chargeNr(iFieldLoc) = fPtr%solverData%scheme%field(iFieldLoc) &
        &                                  %fieldProp%species%chargeNr
      ! 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

    omega_fac = mixture%omega_diff * 0.5_rk
    paramBInv = 1.0_rk/mixture%paramB

    ! Calculate the mole averaged mixture velocity w = sum(n_i v_i)/n
    do iElem = 1, globBC%nElems(iLevel)
!write(dbgUnit(1),*) 'iElem ', iElem
      ! initialize the local velocity for every element
      first_moments = 0.0_rk
      mass_dens = 0.0_rk
      posInBuffer = globBC%elemLvl( iLevel )%posInBcElemBuf%val(iElem)
      do iFieldLoc = 1, nFields
        ! all fields have some nComponents
        do iDir = 1,layout%fStencil%QQ
          ! get the position of the current direction of the depending variable
          ! in the global system (= position in the input)
          pos = varSys%method%val(iFieldLoc)%state_varPos(iDir)
          fTmp_all( pos ) = bcBuffer(                                   &
            & ?IDX?( pos, posInBuffer, nScalars,globBC%nElems(iLevel) ) )
          ! density
          mass_dens(iFieldLoc) = mass_dens(iFieldLoc) + fTmp_all(pos)
          first_moments( :, iFieldLoc ) = first_moments( :, iFieldLoc ) &
            & + fTmp_all(pos) * layout%fStencil%cxDirRK( :, iDir )
        end do !iDir
        ! number density
        num_dens(iFieldLoc) = mass_dens(iFieldLoc) * molWeightInv(iFieldLoc)
      end do !iField

      !total number density
      tot_NumDens = sum(num_dens)
      !total mass density
      tot_massDens = sum(mass_dens)
      ! mole fraction
      moleFrac(:) = num_dens(:) / tot_NumDens

      ! momentum of all species
      momentum = momentumFromMacroLSE( moleFraction  = moleFrac,      &
        &                              first_moments = first_moments, &
        &                              nFields       = nFields,       &
        &                              phi           = phi,           &
        &                              resi_coeff    = resi_coeff,    &
        &                              omega_fac     = omega_fac,     &
        &                              paramBInv     = paramBInv      )

      !velocity of all species
      do iFieldLoc = 1,nFields
        vel( :, iFieldLoc) = momentum( :, iFieldLoc) / mass_dens(iFieldLoc)
      end do

      ! molar diffusive flux = molar flux: Assuming mixture velocity is zero
      ! at membrane surface
      do iFieldLoc = 1, nFields
        molar_flux(:, iFieldLoc) = num_dens(iFieldLoc) * vel(:, iFieldLoc)
      end do

      ! KM: \todo 2016_08_05 Use current density variable
      ! local current density
      loc_curr_dens = 0.0_rk
      do iFieldLoc = 1, nFields
        loc_curr_dens(:) = loc_curr_dens(:) + chargeNr(iFieldLoc) &
          &              * molar_flux(:, iFieldLoc)
      end do

      loc_curr_dens(:) = loc_curr_dens * mixture%faradayLB

      ! project local current density on the membrane
      mem_curr_dens = dot_product(abs(globBC%elemLvl(iLevel)%normal%val(:, iElem)), &
        &                             loc_curr_dens)

      ! membrance molar flux of ionic species
      mem_molar_flux = me%blackbox_mem%transNr * mem_curr_dens                 &
        &            / ( chargeNr(iField) * mixture%faradayLB )

      ! membrance mass flux
      mem_mass_flux = mem_molar_flux*molWeight(iField)
      ! set normal direction of membrance flux
      elec_mass_flux = 0.0_rk
      if( abs(globBC%elemLvl(iLevel)%normal%val(1, iElem)) == 1 )                &
        & elec_mass_flux(1) = mem_mass_flux
      if( abs(globBC%elemLvl(iLevel)%normal%val(2, iElem)) == 1 )                &
        & elec_mass_flux(2) = mem_mass_flux
      if( abs(globBC%elemLvl(iLevel)%normal%val(3, iElem)) == 1 )                &
        & elec_mass_flux(3) = mem_mass_flux

!write(dbgUnit(1),*) 'transNr ', me%blackbox_mem%transNr
!write(dbgUnit(1),*) 'moleFlux ', molar_flux(:, iField)
!write(dbgUnit(1),*) 'moleFluxAll ', molar_flux(:, :)
!write(dbgUnit(1),*) 'loc_curr_dens ', loc_curr_dens
!write(dbgUnit(1),*) 'mem_curr_dens ', mem_curr_dens
!write(dbgUnit(1),*) 'mem_mass_flux ', mem_mass_flux
!write(dbgUnit(1),*) 'elec_mass_flux ', elec_mass_flux
      posInState = globBC%elemLvl(iLevel)%elem%val(iElem)
      do iDir = 1,layout%fStencil%QQN
        ! Write the values
        if( globBC%elemLvl(iLevel)%bitmask%val(iDir, iElem )) then
          state(                                                               &
            & ?FETCH?(iDir, iField, posInState, QQ, nScalars, nSize,neigh)) = &
          ! We need to get post-collision pdf in direction
          ! alpha- outgoing direction, which is the inverse direction of bitmask
          ! For PULL this means, get the outgoing one, as this is the one which
          ! will be bounced back
          ! For PUSH this means, get the already bounced back pdf back, so take
          ! the incoming
            & fTmp_all(varPos(layout%fStencil%cxDirInv(iDir)))               &
            &       + layout%weight( iDir )*2._rk*cs2inv                     &
            &       * ( layout%fStencil%cxDir( 1, iDir )*elec_mass_flux(1)   &
            &       +   layout%fStencil%cxDir( 2, iDir )*elec_mass_flux(2)   &
            &       +   layout%fStencil%cxDir( 3, iDir )*elec_mass_flux(3) )
        end if
      end do
    end do !iElem

  end subroutine spc_blackbox_mem_ion
! ****************************************************************************** !


! ****************************************************************************** !
  !> This routine computes mole diffusion flux of the solvent at the membrance
  !! using black box model and then mass density at the membrane boundary
  !! from mole diffusion flux. Then equilibrium is set at the boundary which is
  !! computed from mass density and velocity
  !! @todo KM: Not completeed. Requires ionic strength at left and right side of
  !! membrane which is not available at the moment
  !!
  !! Black box model: \n
  !!  - for solvent mole diffusion flux:
  !! $J^m_s (x^m,t) = \frac{t^m_s(x^m,t)}{F} i^m(x^m,t)
  !!                + LW(IS^l(x^m,t)-IS^r(x^m,t))$
  !!
  !! $J^m_s$ - mole diffusion flux of solvent through the membrane
  !! $t^m_k$ - transference number of solvent on the membrane
  !!           (selective membrane property)
  !! $F$ - Faraday constant
  !! $LW$ - osmotic solvent transport coefficient
  !! $IS^l$ and $IS^r$ are ionic strength of electrolyte solution at the left
  !! and right side of membrance
  !! $IS^l(x^m,t)=\frac{1}{2}\sum^n_{k=1} z^2_k c^m_k(x^m,t)$
  !! $c^m_k = \rho^m_k/m_k$ - molar concentration
  !! $i^m(x^m,t) = n^l_x_1 \cdot i^l(x,t) = - n^r_x_1 \cdot i^r(x,t)$
  !!    - orthogonal projection of the current density at the left $i^l$ and
  !!      right $i^r$ side of the membrance.
  !! The normal vectors $n^l_1$ and $n^r_1$ point into the membrane surface
  !! $i(x,t) = F \sum_{j=1}^{n} z_j J^e_j(x,t)$
  !!    - local current density [A/m^3]
  !! $J^e_k = c^e_k(v^e_k - v) $ - mole diffusion flux of ion in the electrolyte
  !! In black box model: mixture averaged velocity at the membrane is assumed to
  !! be zero
  !!
  !! 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_blackbox_mem_solvent( 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
    ! -------------------------------------------------------------------- !
    !store pdf of all fields for derive function pointer input
    real(kind=rk) :: fTmp_all( layout%fStencil%QQ*varSys%nStateVars )
    real(kind=rk) :: mass_dens( varSys%nStateVars ), &
      &              num_dens( varSys%nStateVars )
    integer :: iElem, iDir, iFieldLoc, nFields, pos, bndNormalDir, QQ
    real(kind=rk) :: molWeight(varSys%nStateVars), chargeNr(varSys%nStateVars)
    real(kind=rk), dimension(3, varSys%nStateVars):: molar_flux, vel, momentum
    real(kind=rk) :: loc_curr_dens(3)
    real(kind=rk) :: mem_curr_dens, mem_molar_flux
    real(kind=rk) :: molefrac(3), mem_mass_flux, elec_mass_flux(3)
    real(kind=rk) :: tot_NumDens, tot_MassDens
    real(kind=rk) :: first_moments( 3, varSys%nStateVars )
    real(kind=rk) :: molWeightInv(varSys%nStateVars), phi(varSys%nStateVars)
    real(kind=rk) :: resi_coeff( varSys%nStateVars, varSys%nStateVars )
    real(kind=rk) :: paramBInv, omega_fac
    type(mus_varSys_data_type), pointer :: fPtr
    integer :: posInBuffer, posInState
    ! ------------------------------------------------------------------------
    QQ = layout%fStencil%QQ
    nFields = varSys%nStateVars

    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
      molWeightInv(iFieldLoc) =  fPtr%solverData%scheme%field(iFieldLoc)   &
        &                                    %fieldProp%species%molWeightInv
      ! charge number
      chargeNr(iFieldLoc) = fPtr%solverData%scheme%field(iFieldLoc) &
        &                       %fieldProp%species%chargeNr
      ! 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

    omega_fac = mixture%omega_diff * 0.5_rk
    paramBInv = 1.0_rk/mixture%paramB

    ! Calculate the mole averaged mixture velocity w = sum(n_i v_i)/n
    do iElem = 1, globBC%nElems(iLevel)
      ! initialize the local velocity for every element
      first_moments = 0.0_rk
      mass_dens = 0.0_rk
      posInBuffer = globBC%elemLvl( iLevel )%posInBcElemBuf%val(iElem)
      do iFieldLoc = 1, nFields
        ! all fields have some nComponents
        do iDir = 1,layout%fStencil%QQ
          ! get the position of the current direction of the depending variable
          ! in the global system (= position in the input)
          pos = varSys%method%val(iFieldLoc)%state_varPos(iDir)
          fTmp_all( pos ) = bcBuffer(                                  &
            & ?IDX?( pos, posInBuffer, nScalars,globBC%nElems(iLevel) ))
          ! density
          mass_dens(iFieldLoc) = mass_dens(iFieldLoc) + fTmp_all(pos)
          first_moments( :, iFieldLoc ) = first_moments( :, iFieldLoc ) &
            & + fTmp_all(pos) * layout%fStencil%cxDirRK( :, iDir )
        end do !iDir
        ! number density
        num_dens(iFieldLoc) = mass_dens(iFieldLoc) * molWeightInv(iFieldLoc)
      end do !iField

      !total number density
      tot_NumDens = sum(num_dens)
      !total mass density
      tot_massDens = sum(mass_dens)
      ! mole fraction
      moleFrac(:) = num_dens(:) / tot_NumDens

      ! momentum of all species
      momentum = momentumFromMacroLSE( moleFraction  = moleFrac,      &
        &                              first_moments = first_moments, &
        &                              nFields       = nFields,       &
        &                              phi           = phi,           &
        &                              resi_coeff    = resi_coeff,    &
        &                              omega_fac     = omega_fac,     &
        &                              paramBInv     = paramBInv      )

      !velocity of all species
      do iFieldLoc = 1,nFields
        vel( :, iFieldLoc) = momentum( :, iFieldLoc) / mass_dens(iFieldLoc)
      end do

      ! molar diffusive flux = molar flux
      do iFieldLoc = 1, nFields
        molar_flux(:, iFieldLoc) = num_dens(iFieldLoc) &
          & * (vel(:, iFieldLoc))
      end do

      ! local current density
      loc_curr_dens = 0.0_rk
      do iFieldLoc = 1, nFields
        loc_curr_dens(:) = loc_curr_dens(:) + chargeNr(iFieldLoc) &
          &              * molar_flux(:, iFieldLoc)
      end do

      loc_curr_dens = loc_curr_dens * mixture%faradayLB

      ! normal direction point to boundary
      bndNormalDir = layout%fStencil%cxDirInv( globBC%elemLvl( iLevel )%       &
        &                                               normalInd%val( iElem ))
      ! project local current density on the membrane
      mem_curr_dens = dot_product(loc_curr_dens,                          &
        &                         layout%fStencil%cxDirRK(:, bndNormalDir))

      ! membrance molar flux of ionic species
      mem_molar_flux = me%blackbox_mem%transNr * mem_curr_dens                 &
        &            / ( mixture%faradayLB )

      ! membrance mass flux
      mem_mass_flux = mem_molar_flux*molWeight(iField)
      elec_mass_flux = 0.0_rk
      ! set normal direction of membrance flux
      if( abs(globBC%elemLvl(iLevel)%normal%val(1, iElem)) == 1 )                &
        & elec_mass_flux(1) = mem_mass_flux
      if( abs(globBC%elemLvl(iLevel)%normal%val(2, iElem)) == 1 )                &
        & elec_mass_flux(2) = mem_mass_flux
      if( abs(globBC%elemLvl(iLevel)%normal%val(3, iElem)) == 1 )                &
        & elec_mass_flux(3) = mem_mass_flux

      posInState = globBC%elemLvl(iLevel)%elem%val(iElem)
      do iDir = 1,layout%fStencil%QQN
        ! Write the values
        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, posInState, QQ, nScalars, nSize,neigh)) = &
          ! We need to get post-collision pdf in direction
          ! alpha- outgoing direction, which is the inverse direction of bitmask
          ! For PULL this means, get the outgoing one, as this is the one which
          ! will be bounced back
          ! For PUSH this means, get the already bounced back pdf back, so take
          ! the incoming
            & fTmp_all(varPos(layout%fStencil%cxDirInv(iDir)))             &
            &       + layout%weight( iDir )*2._rk*cs2inv                   &
            &       * ( layout%fStencil%cxDir( 1, iDir )*elec_mass_flux(1) &
            &       +   layout%fStencil%cxDir( 2, iDir )*elec_mass_flux(2) &
            &       +   layout%fStencil%cxDir( 3, iDir )*elec_mass_flux(3) )
        end if
      end do
    end do !iElem

  end subroutine spc_blackbox_mem_solvent
! ****************************************************************************** !


! ****************************************************************************** !
  !> author: Kannan Masilamani
  !! Moment based wall 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 spc_moments_wall( me, state, bcBuffer, globBC, levelDesc, tree,   &
    &                          nSize, iLevel, sim_time, neigh, layout,         &
    &                          fieldProp, varPos, nScalars, varSys, derVarPos, &
    &                          physics, iField, mixture                        )
    ! -------------------------------------------------------------------- !
    !> global boundary type
    class(boundary_type) :: me
    !> Current state vector of iLevel
    real(kind=rk), intent(inout) :: state(:)
    !> size of state array ( in terms of elements )
    integer, intent(in) :: nSize
    !> state values of boundary elements of all fields of iLevel
    real(kind=rk), intent(in) :: bcBuffer(:)
    !> iLevel descriptor
    type(tem_levelDesc_type), intent(in) :: levelDesc
    !> Treelm Mesh
    type(treelmesh_type), intent(in) :: tree
    !> fluid parameters and properties
    type(mus_field_prop_type), intent(in) :: fieldProp
    !> stencil layout information
    type(mus_scheme_layout_type), intent(in) :: layout
    !> the level On which this boundary was invoked
    integer, intent(in) :: iLevel
    !> connectivity array corresponding to state vector
    integer, intent(in) :: neigh(:)
    !> global time information
    type(tem_time_type), intent(in)  :: sim_time
    !> pointer to field variable in the state vector
    integer, intent(in) :: varPos(:)
    !> number of Scalars in the scheme var system
    integer, intent(in) :: nScalars
    !> scheme variable system
    type(tem_varSys_type), intent(in) :: varSys
    !> position of derived quantities in varsys
    type(mus_derVarPos_type), intent(in) :: derVarPos
    !> scheme global boundary type
    type(glob_boundary_type), intent(in) :: globBC
    !> scheme global boundary type
    type(mus_physics_type), intent(in) :: physics
    !> current field
    integer, intent(in) :: iField
    !> mixture info
    type(mus_mixture_type), intent(in) :: mixture
    ! -------------------------------------------------------------------- !
    real(kind=rk) :: rhoB, press
    integer :: iElem, iDir, tmpDir, 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 )
    real(kind=rk) :: first_moments(3), second_moments(6), third_moments(6), &
      &              fourth_moments(3)
    integer :: nMom(4)
    ! ---------------------------------------------------------------------------
    QQ = layout%fStencil%QQ
    nMom(1) = 1 + size(layout%moment%first_moments)
    nMom(2) = nMom(1) + size(layout%moment%second_moments)
    nMom(3) = nMom(2) + size(layout%moment%third_moments)
    nMom(4) = nMom(3) + size(layout%moment%fourth_moments)

!    write(*,*) 'Boundary label ', trim(me%label)
    do iElem = 1, globBC%nElems( iLevel )
      if( .not. btest( levelDesc%property(                        &
        &              globBC%elemLvl(iLevel)%elem%val(iElem)), prp_solid))then

      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 )
!      write(*,*) 'elemPos ', elemPos
!      write(*,*) 'treeID ', levelDesc%total(elemPos)

      iLink = 0
      fTmp = 0.0_rk
      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 ))
!        write(*,*) 'fTmp0 ',  fTmp

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

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

      ! moments for multispecies LBM model of transformed variable
      first_moments = (/ 0.0_rk, 0.0_rk, 0.0_rk /) !momX, momY, momZ
      second_moments = (/ press, & !momXX
        &                 press, & !momYY
        &                 press, & !momZZ
        &                 0.0_rk, & !momXY
        &                 0.0_rk, & !momYZ
        &                 0.0_rk /) !momXZ
      third_moments = 0.0_rk
      fourth_moments = cs2*press !momXXYY, momYYZZ, momZZXX

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

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

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

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

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

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

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

  end subroutine spc_moments_wall
! ****************************************************************************** !


! ****************************************************************************** !
  !> Mole fraction boundary condition
  !! Usage
  !! -----
  !!```lua
  !!boundary_condition = {
  !! { label = 'outlet',
  !!    kind = 'spc_moments_moleFrac',
  !!    moleFraction = 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_moleFrac( 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
    real(kind=rk) :: moleFrac(globBC%nElems(iLevel))
    integer :: iElem, iDir, iFieldLoc, nFields, pos
    real(kind=rk) :: molWeight(varSys%nStateVars), phi(varSys%nStateVars)
    real(kind=rk) :: resi_coeff( varSys%nStateVars, varSys%nStateVars )
    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) :: fTmp_all( layout%fStencil%QQ*varSys%nStateVars )
    real(kind=rk) :: vel_inG(3), velAvg(3), velQuad(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
    real(kind=rk) :: moleFrac_loc(varSys%nStateVars)
    real(kind=rk) :: moments( layout%fStencil%QQ )
    real(kind=rk) :: first_moments(3), second_moments(6), third_moments(6), &
      &              fourth_moments(3)
    integer :: nMom(4)
    type(mus_varSys_data_type), pointer :: fPtr
    integer :: bcMoleFrac_pos, posInBuffer
    ! ---------------------------------------------------------------------------
    QQ = layout%fStencil%QQ
    nFields = varSys%nStateVars

    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

    ! position of molefraction spacetime variable in varSys
    bcMoleFrac_pos = me%bc_states%moleFrac%varPos
    ! mole fraction
    call varSys%method%val(bcMoleFrac_pos)%get_valOfIndex( &
      & varSys  = varSys,                                  &
      & time    = sim_time,                                &
      & iLevel  = iLevel,                                  &
      & idx     = me%bc_states%moleFrac                    &
      &           %pntIndex%indexLvl(iLevel)               &
      &           %val(1:globBC%nElems(iLevel)),           &
      & nVals   = globBC%nElems(iLevel),                   &
      & res     = moleFrac                                 )

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

    ! Calculate the density of current element
    do iElem = 1, globBC%nElems(iLevel)
      fTmp = 0.0_rk
      nLinks = me%elemLvl(iLevel)%moments(iElem)%nUnKnownPdfs
      allocate(missing_links(nLinks))
      allocate(unKnown_fTmp(nLinks))
      allocate(rhs(nLinks))
      elemPos = globBC%elemLvl(iLevel)%elem%val( iElem )

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

      posInBuffer = globBC%elemLvl( iLevel )%posInBcElemBuf%val(iElem)
      ! local state vector to compute density and velocity of other species
      do iFieldLoc = 1, varSys%nStateVars
        do iDir = 1,layout%fStencil%QQ
          pos = varSys%method%val(iFieldLoc)%state_varPos(iDir)
          fTmp_all( pos ) = bcBuffer(                                  &
            & ?IDX?( pos, posInBuffer, nScalars,globBC%nElems(iLevel) ))
        end do
      end do

      !KM: using inital mixture number density rho0/mixtureMOlWeight
      rho = mixture%moleDens0LB *  moleFrac(iElem) &
        & * fieldProp%species%molWeight
      press = cs2 * rho * fieldProp%species%molWeigRatio

      !local density and momentum
      mass_dens = 0.0_rk
      do iFieldLoc = 1, nFields
        mass_dens(iFieldLoc) = &
          & sum(fTmp_all(varSys%method%val(iFieldLoc)%state_varPos(:)))
        num_dens(iFieldLoc) = mass_dens(iFieldLoc) &
          &                 / molWeight(iFieldLoc)
        !velocity
        call derVarPos%momFromState( state  = fTmp_all,         &
          &                          iField = iFieldLoc,        &
          &                          nElems = 1,                &
          &                          varSys = varSys,           &
          &                          layout = layout,           &
          &                          res    = mom(:, iFieldLoc) )
      end do

      mass_dens(iField) = rho
      num_dens(iField) = rho / fieldProp%species%molWeight

      ! total mass density
      totMassDens = sum(mass_dens)

      ! mole fraction
      moleFrac_loc = num_dens/sum(num_dens)

!do iFieldLoc = 1, nFields
!  write(dbgUnit(1),*) iFieldLoc, 'vel ', mom(:, iFieldLoc)/mass_dens(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_loc(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_loc(iField) / mixture%paramB )        &
          &                 * mom(:,iFieldLoc)
      end do

      vel_inG = mom_inG/rho
!write(dbgUnit(1),*) 'transformed momentum ', mom_inG
!write(dbgUnit(1),*) 'transformed velocity ', vel_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
      velQuad(:) = mixture%theta_eq*velAvg(:) &
        &        + (1.0_rk - mixture%theta_eq)*vel_inG

      ! moments for multispecies LBM model of transformed variable
      first_moments = mom_inG !momX, momY, momZ
      second_moments = (/ press + rho*velQuad(1)*velQuad(1), & !momXX
        &                 press + rho*velQuad(2)*velQuad(2), & !momYY
        &                 press + rho*velQuad(3)*velQuad(3), & !momZZ
        &                 rho*velQuad(1)*velQuad(2),         & !momXY
        &                 rho*velQuad(2)*velQuad(3),         & !momYZ
        &                 rho*velQuad(3)*velQuad(1) /)         !momXZ
      third_moments = (/ cs2*mom_inG(2), & !momXXY
        &                cs2*mom_inG(3), & !momXXZ
        &                cs2*mom_inG(1), & !momYYX
        &                cs2*mom_inG(3), & !momYYZ
        &                cs2*mom_inG(1), & !momZZX
        &                cs2*mom_inG(2) /) !momZZY
      fourth_moments = (/ &
        & cs2*(press+rho*(        velQuad(1)**2 +        velQuad(2)**2 &
        &                - 0.5_rk*velQuad(3)**2 )), &!momXXYY
        & cs2*(press+rho*(-0.5_rk*velQuad(1)**2 +        velQuad(2)**2 &
        &                + velQuad(3)**2 )), &!momYYZZ
        & cs2*(press+rho*(        velQuad(1)**2 - 0.5_rk*velQuad(2)**2 &
        &                + velQuad(3)**2 ))  &!momZZXX
        & /)

      ! moments for liquid mixture
      moments = 0.0_rk
      moments(1) = rho
      moments(2:nMom(1)) = &
        & first_moments(layout%moment%first_moments(:))
      moments(nMom(1)+1:nMom(2)) = &
        & second_moments(layout%moment%second_moments(:))
      moments(nMom(2)+1:nMom(3)) = &
        & third_moments(layout%moment%third_moments(:))
      moments(nMom(3)+1:nMom(4)) = &
        & fourth_moments(layout%moment%fourth_moments(:))

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

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

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

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

  end subroutine spc_moments_moleFrac
! ****************************************************************************** !


! ****************************************************************************** !
  !> molar flux boundary condition like moments velocity bc type
  !! Usage
  !! -----
  !!```lua
  !!boundary_condition = {
  !! { label = 'outlet',
  !!    kind = 'spc_moments_moleflux',
  !!    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_moleFlux( 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) :: moleFlux(globBC%nElems(iLevel)*3), inv_flux
    integer :: iElem, iDir, iFieldLoc, nFields, pos, tmpDir
    real(kind=rk) :: rho, press, dens_correction
    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) :: fTmp_all( layout%fStencil%QQ*varSys%nStateVars )
    real(kind=rk) :: moments( layout%fStencil%QQ )
    real(kind=rk) :: massFlux(3), 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), velQuad(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)
    real(kind=rk) :: molWeightInv
    real(kind=rk) :: first_moments(3), second_moments(6), third_moments(6), &
      &              fourth_moments(3)
    integer :: nMom(4)
    type(mus_varSys_data_type), pointer :: fPtr
    integer :: offset, bcMoleFlux_pos, posInBuffer
    ! ------------------------------------------------------------------------
!write(dbgUnit(1),*) 'Boundary label ', trim(me%label)
!write(dbgUnit(1),*) 'iField ', iField

    QQ = layout%fStencil%QQ
    nFields = varSys%nStateVars
    molWeightInv = fieldProp%species%molWeightInv

    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

    inv_flux = 1.0_rk / physics%fac( iLevel )%moleFlux
    ! position of boundary moleflux in varSys
    bcMoleFlux_pos = me%bc_states%moleFlux%varPos
    ! Get moleflux
    call varSys%method%val(bcMoleFlux_pos)%get_valOfIndex( &
      & varSys  = varSys,                                  &
      & time    = sim_time,                                &
      & iLevel  = iLevel,                                  &
      & idx     = me%bc_states%moleFlux                    &
      &           %pntIndex%indexLvl(iLevel)               &
      &           %val(1:globBC%nElems(iLevel)),           &
      & nVals   = globBC%nElems(iLevel),                   &
      & res     = moleFlux                                 )

    ! If physical quantities are given, transform to lattice units by division
    ! with the conversion factor
    moleFlux = moleFlux * inv_flux
!write(*,*) 'physical moleFlux ', physics%fac( iLevel )%moleFlux
!write(*,*) 'moleFlux ', moleFlux
!stop

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

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

      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
      fTmp = 0.0_rk
      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 ))
      !write(*,*) 'fTmp ', fTmp

      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(dbgUnit(1),*) 'density ', rho

      !caulate mass flux from moleFlux
      ! massFlux = moleflux * molWeight
      offset = (iElem-1)*3
      massFlux(1) = moleFlux(offset+1) * molWeightInv
      massFlux(2) = moleFlux(offset+2) * molWeightInv
      massFlux(3) = moleFlux(offset+3) * molWeightInv
!write(dbgUnit(1),*) 'massFlux ', massFlux

      ! 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(massFlux,                            &
        &                           globBC%elemLvl(iLevel)%normal%val(:, iElem))
      ! correction term for Multispecies when flux is defined
      ! rho_new = rho_old + dens_correction
      rho = rho + dens_correction
      !write(*,*) 'density after correction', rho

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

      posInBuffer = globBC%elemLvl( iLevel )%posInBcElemBuf%val(iElem)
      ! local state vector to compute density and velocity of other species
      do iFieldLoc = 1, varSys%nStateVars
        do iDir = 1,layout%fStencil%QQ
          pos = varSys%method%val(iFieldLoc)%state_varPos(iDir)
          fTmp_all( pos ) = bcBuffer(                                  &
            & ?IDX?( pos, posInBuffer, nScalars,globBC%nElems(iLevel) ))
        end do
      end do

      !local density and momentum
      mass_dens = 0.0_rk
      do iFieldLoc = 1, nFields
        if(iField == iFieldLoc) then
          ! massflux = moleflux * molecular_weight
          mom(:, iField) = massFlux
          mass_dens(iField) = rho
          num_dens(iField) = rho * molWeightInv
        else
          ! mass density
          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 if
      end do

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

      ! mole fraction
      moleFrac = num_dens/sum(num_dens)

      ! 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(dbgUnit(1),*) 'defined momentum ', massFlux
!write(dbgUnit(1),*) 'transformed momentum ', mom_inG
!write(dbgUnit(1),*) 'transformed velocity ', vel_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
      velQuad(:) = mixture%theta_eq*velAvg(:) &
        &        + (1.0_rk - mixture%theta_eq)*vel_inG

      ! moments for multispecies LBM model of transformed variable
      first_moments = mom_inG !momX, momY, momZ
      second_moments = (/ press + rho*velQuad(1)*velQuad(1), & !momXX
        &                 press + rho*velQuad(2)*velQuad(2), & !momYY
        &                 press + rho*velQuad(3)*velQuad(3), & !momZZ
        &                 rho*velQuad(1)*velQuad(2),         & !momXY
        &                 rho*velQuad(2)*velQuad(3),         & !momYZ
        &                 rho*velQuad(3)*velQuad(1) /)         !momXZ
      third_moments = (/ cs2*mom_inG(2), & !momXXY
        &                cs2*mom_inG(3), & !momXXZ
        &                cs2*mom_inG(1), & !momYYX
        &                cs2*mom_inG(3), & !momYYZ
        &                cs2*mom_inG(1), & !momZZX
        &                cs2*mom_inG(2) /) !momZZY
      fourth_moments = (/ &
        & cs2*(press+rho*(        velQuad(1)**2 +        velQuad(2)**2 &
        &                - 0.5_rk*velQuad(3)**2 )), &!momXXYY
        & cs2*(press+rho*(-0.5_rk*velQuad(1)**2 +        velQuad(2)**2 &
        &                + velQuad(3)**2 )), &!momYYZZ
        & cs2*(press+rho*(        velQuad(1)**2 - 0.5_rk*velQuad(2)**2 &
        &                + velQuad(3)**2 ))  &!momZZXX
        & /)

      ! moments for liquid mixture
      moments = 0.0_rk
      moments(1) = rho
      moments(2:nMom(1)) = &
        & first_moments(layout%moment%first_moments(:))
      moments(nMom(1)+1:nMom(2)) = &
        & second_moments(layout%moment%second_moments(:))
      moments(nMom(2)+1:nMom(3)) = &
        & third_moments(layout%moment%third_moments(:))
      moments(nMom(3)+1:nMom(4)) = &
        & fourth_moments(layout%moment%fourth_moments(:))

      !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

  end subroutine spc_moments_moleFlux
! ****************************************************************************** !


! ****************************************************************************** !
  !> velocity boundary condition like moments velocity bc type
  !! Usage
  !! -----
  !!```lua
  !!boundary_condition = {
  !! { label = 'outlet',
  !!    kind = 'spc_moments_vel',
  !!    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_vel( 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) :: vel_b(globBC%nElems(iLevel)*3), inv_vel
    integer :: iElem, iDir, iFieldLoc, nFields, pos, tmpDir
    real(kind=rk) :: rho, press, dens_correction
    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) :: fTmp_all( layout%fStencil%QQ*varSys%nStateVars )
    real(kind=rk) :: moments( layout%fStencil%QQ )
    real(kind=rk) :: vel(3), 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) :: velAvg(3), velQuad(3)
    real(kind=rk) :: totMassDens
    real(kind=rk) :: moleFrac(varSys%nStateVars)
    real(kind=rk) :: resi_coeff( varSys%nStateVars, varSys%nStateVars )
    real(kind=rk) :: molWeight(varSys%nStateVars), phi(varSys%nStateVars)
    real(kind=rk) :: first_moments(3), second_moments(6), third_moments(6), &
      &              fourth_moments(3)
    integer :: nMom(4)
    type(mus_varSys_data_type), pointer :: fPtr
    integer :: bcVel_pos, posInBuffer
    ! ------------------------------------------------------------------------
    !write(*,*) 'Boundary label ', trim(me%label)

    !write(*,*) 'iField ', iField

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

    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

    ! 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

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

!KM!write(dbgUnit(1),*) 'iField ', iField
    do iElem = 1, globBC%nElems(iLevel)
!KM!write(dbgUnit(1),*) 'iElem ', iElem
!KM!write(dbgUnit(1),*) 'vel ', ux(iElem), uy(iElem), uz(iElem)
      vel = vel_b( (iElem-1)*3+1 : iElem*3 )

      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
      fTmp = 0.0_rk
      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 ))

      ! Calculate the density of current element
      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
!KM!write(dbgUnit(1),*) 'density ', 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))
      ! correction term for Multispecies when vel is defined
      ! rho_new = rho_old/(1.0 - dens_correction)
      rho = rho/(1.0_rk-dens_correction)

!KM!write(dbgUnit(1),*) 'density after correction', rho
      !pressure
      press = cs2*rho*fieldProp%species%molWeigRatio
!KM!write(dbgUnit(1),*) 'press ', press

      ! local state vector to compute density and velocity of other species
      posInBuffer = globBC%elemLvl( iLevel )%posInBcElemBuf%val(iElem)
      do iFieldLoc = 1, varSys%nStateVars
        do iDir = 1,layout%fStencil%QQ
          pos = varSys%method%val(iFieldLoc)%state_varPos(iDir)
          fTmp_all( pos ) = bcBuffer(                                  &
            & ?IDX?( pos, posInBuffer, nScalars,globBC%nElems(iLevel) ))
        end do
      end do

      !local density and momentum
      mass_dens = 0.0_rk
      do iFieldLoc = 1, nFields
        if(iField == iFieldLoc) then
          ! massflux = moleflux * molecular_weight
          mom(:, iField) = rho * vel
          mass_dens(iField) = rho
          num_dens(iField) = rho / fieldProp%species%molWeight
        else
          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 if
!KM!write(dbgUnit(1),*) iFieldLoc, ' mom ', mom(:, iFieldLoc)
      end do

      ! total mass density
      totMassDens = sum(mass_dens)

      ! mole fraction
      moleFrac = num_dens/sum(num_dens)

!KM!do iFieldLoc = 1, nFields
!KM!  write(dbgUnit(1),*) iFieldLoc, 'vel ', mom(:, iFieldLoc)/mass_dens(iFieldLoc), &
!KM!    & 'dens ', mass_dens(iFieldLoc), 'moleFrac ', moleFrac(iFieldLoc)
!KM!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(dbgUnit(1),*) 'transformed momentum ', mom_inG
!write(dbgUnit(1),*) 'transformed velocity ', vel_inG

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

      velQuad(:) = mixture%theta_eq*velAvg(:) &
        &        + (1.0_rk - mixture%theta_eq)*vel_inG
!      write(*,*) 'velAvg ', velAvg

      ! moments for multispecies LBM model of transformed variable
      first_moments = mom_inG !momX, momY, momZ
      second_moments = (/ press + rho*velQuad(1)*velQuad(1), & !momXX
        &                 press + rho*velQuad(2)*velQuad(2), & !momYY
        &                 press + rho*velQuad(3)*velQuad(3), & !momZZ
        &                 rho*velQuad(1)*velQuad(2),         & !momXY
        &                 rho*velQuad(2)*velQuad(3),         & !momYZ
        &                 rho*velQuad(3)*velQuad(1) /)         !momXZ
      third_moments = (/ cs2*mom_inG(2), & !momXXY
        &                cs2*mom_inG(3), & !momXXZ
        &                cs2*mom_inG(1), & !momYYX
        &                cs2*mom_inG(3), & !momYYZ
        &                cs2*mom_inG(1), & !momZZX
        &                cs2*mom_inG(2) /) !momZZY
      fourth_moments = (/ &
        & cs2*(press+rho*(        velQuad(1)**2 +        velQuad(2)**2 &
        &                - 0.5_rk*velQuad(3)**2 )), &!momXXYY
        & cs2*(press+rho*(-0.5_rk*velQuad(1)**2 +        velQuad(2)**2 &
        &                + velQuad(3)**2 )), &!momYYZZ
        & cs2*(press+rho*(        velQuad(1)**2 - 0.5_rk*velQuad(2)**2 &
        &                + velQuad(3)**2 ))  &!momZZXX
        & /)

      ! moments for liquid mixture
      moments = 0.0_rk
      moments(1) = rho
      moments(2:nMom(1)) = &
        & first_moments(layout%moment%first_moments(:))
      moments(nMom(1)+1:nMom(2)) = &
        & second_moments(layout%moment%second_moments(:))
      moments(nMom(2)+1:nMom(3)) = &
        & third_moments(layout%moment%third_moments(:))
      moments(nMom(3)+1:nMom(4)) = &
        & fourth_moments(layout%moment%fourth_moments(:))

      !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

  end subroutine spc_moments_vel
! ****************************************************************************** !

! **************************************************************************** !
  !> This routine returns mass density of all species and mass averaged mixture
  !! velocity from given pdf of all species for single element.
  !! It is used in Nonequilbrium extrapolation based boundary conditions.
  subroutine calcDensAndVelsFromPDF( nFields, iField, mass_dens, num_dens,     &
    &                                velAvg, eqVel, varSys, pdf, stencil,      &
    &                                species, resi_coeff, omega_fac, paramBInv )
    ! --------------------------------------------------------------------------
    !> Number of fields
    integer, intent(in) :: nFields
    !> current field to compute eqVel
    integer, intent(in) :: iField
    !> Mass density of all species
    real(kind=rk), intent(out) :: mass_dens(nFields)
    !> Number density of all species
    real(kind=rk), intent(out) :: num_dens(nFields)
    !> Mass averaged velocity
    real(kind=rk), intent(out) :: velAvg(3)
    !> Variable system
    type(tem_varSys_type), intent(in) :: varSys
    !> State array of single element for all species
    real(kind=rk), intent(in) :: pdf(varSys%nScalars)
    !> Stencil header
    type(tem_stencilHeader_type), intent(in) :: stencil
    ! Equilibrium velocity of current species
    real(kind=rk), intent(out) :: eqVel(3)
    !> Species information
    type(mus_species_type), intent(in) :: species(nFields)
    !> Resisitivity coefficient matrix
    real(kind=rk), intent(in) :: resi_coeff(nFields, nFields)
    !> relaxation parameter, omega_diff*0.5_rk
    real(kind=rk), intent(in) :: omega_fac
    !> inverse of free parameter B
    real(kind=rk), intent(in) :: paramBInv
    ! --------------------------------------------------------------------------
    integer :: iFld, iDir
    real(kind=rk) :: fTmp(stencil%QQ)
    real(kind=rk) :: tot_numDensInv, tot_massDensInv
    real(kind=rk) :: momentum(3, nFields), first_moments(3, nFields)
    real(kind=rk) :: moleFrac(nFields), massFrac(nFields)
    real(kind=rk) :: velocity(nFields, nFields)
    ! --------------------------------------------------------------------------
    do iFld = 1, nFields
      do iDir = 1, stencil%QQ
        ! State array of iFld
        fTmp(iDir) = pdf( varSys%method%val(iFld)%state_varPos(iDir) )
      end do
      ! species density
      mass_dens(iFld) = sum(fTmp)
      ! number density
      num_dens(iFld) = mass_dens(iFld) * species(iFld)%molWeightInv
      ! first moments: momentum of transformed PDF
      first_moments(1, iFld) = sum( fTmp * stencil%cxDirRK(1, :) )
      first_moments(2, iFld) = sum( fTmp * stencil%cxDirRK(2, :) )
      first_moments(3, iFld) = sum( fTmp * stencil%cxDirRK(3, :) )
    end do

    !total mass density
    tot_massDensInv = 1.0_rk / sum(mass_dens)
    ! mass fraction
    massFrac(:) = mass_dens(:) * tot_massDensInv
    !total number density
    tot_numDensInv = 1.0_rk / sum(num_dens)
    ! mole fraction
    moleFrac(:) = num_dens(:) * tot_numDensInv

    ! actual momentum of all species obtained from solving
    ! Linear Equation system
    momentum = momentumFromMacroLSE(                      &
      &          moleFraction  = moleFrac,                &
      &          first_moments = first_moments,           &
      &          nFields       = nFields,                 &
      &          phi           = species(:)%molWeigRatio, &
      &          resi_coeff    = resi_coeff,              &
      &          omega_fac     = omega_fac,               &
      &          paramBInv     = paramBInv                )

    ! velocity of all species on current fluid node
    do iFld = 1, nFields
      velocity(:, iFld) = momentum(:, iFld) / mass_dens(iFld)
    end do

    ! equilibrium velocity for current species
    eqVel(:) = velocity(:, iField)
    do iFld = 1, nFields
      eqVel(:) = eqVel(:) + resi_coeff(iField, iFld)           &
        &      * species(iField)%molWeigRatio * moleFrac(iFld) &
        &      * ( velocity(:, iFld) - velocity(:, iField) )   &
        &      * paramBInv
    end do

    ! Mass averaged mixture velocity
    velAvg(1) = dot_product( massFrac, velocity(1, :) )
    velAvg(2) = dot_product( massFrac, velocity(2, :) )
    velAvg(3) = dot_product( massFrac, velocity(3, :) )

  end subroutine calcDensAndVelsFromPDF
! **************************************************************************** !

?? ENDIF
end module mus_bc_species_module
! ****************************************************************************** !
