mus_bc_species_module.f90 Source File


This file depends on

sourcefile~~mus_bc_species_module.f90~~EfferentGraph sourcefile~mus_bc_species_module.f90 mus_bc_species_module.f90 sourcefile~mus_field_prop_module.f90 mus_field_prop_module.f90 sourcefile~mus_bc_species_module.f90->sourcefile~mus_field_prop_module.f90 sourcefile~mus_enrtl_module.f90 mus_eNRTL_module.f90 sourcefile~mus_bc_species_module.f90->sourcefile~mus_enrtl_module.f90 sourcefile~mus_scheme_layout_module.f90 mus_scheme_layout_module.f90 sourcefile~mus_bc_species_module.f90->sourcefile~mus_scheme_layout_module.f90 sourcefile~mus_bc_header_module.f90 mus_bc_header_module.f90 sourcefile~mus_bc_species_module.f90->sourcefile~mus_bc_header_module.f90 sourcefile~mus_param_module.f90 mus_param_module.f90 sourcefile~mus_bc_species_module.f90->sourcefile~mus_param_module.f90 sourcefile~mus_species_module.f90 mus_species_module.f90 sourcefile~mus_bc_species_module.f90->sourcefile~mus_species_module.f90 sourcefile~mus_varsys_module.f90 mus_varSys_module.f90 sourcefile~mus_bc_species_module.f90->sourcefile~mus_varsys_module.f90 sourcefile~mus_derquanmsliquid_module.f90 mus_derQuanMSLiquid_module.f90 sourcefile~mus_bc_species_module.f90->sourcefile~mus_derquanmsliquid_module.f90 sourcefile~mus_dervarpos_module.f90 mus_derVarPos_module.f90 sourcefile~mus_bc_species_module.f90->sourcefile~mus_dervarpos_module.f90 sourcefile~mus_mixture_module.f90 mus_mixture_module.f90 sourcefile~mus_bc_species_module.f90->sourcefile~mus_mixture_module.f90 sourcefile~mus_physics_module.f90 mus_physics_module.f90 sourcefile~mus_bc_species_module.f90->sourcefile~mus_physics_module.f90 sourcefile~mus_field_prop_module.f90->sourcefile~mus_species_module.f90 sourcefile~mus_field_prop_module.f90->sourcefile~mus_physics_module.f90 sourcefile~mus_fluid_module.f90 mus_fluid_module.f90 sourcefile~mus_field_prop_module.f90->sourcefile~mus_fluid_module.f90 sourcefile~mus_scheme_header_module.f90 mus_scheme_header_module.f90 sourcefile~mus_field_prop_module.f90->sourcefile~mus_scheme_header_module.f90 sourcefile~mus_poisson_module.f90 mus_poisson_module.f90 sourcefile~mus_field_prop_module.f90->sourcefile~mus_poisson_module.f90 sourcefile~mus_moments_type_module.f90 mus_moments_type_module.f90 sourcefile~mus_scheme_layout_module.f90->sourcefile~mus_moments_type_module.f90 sourcefile~mus_bc_header_module.f90->sourcefile~mus_field_prop_module.f90 sourcefile~mus_bc_header_module.f90->sourcefile~mus_scheme_layout_module.f90 sourcefile~mus_bc_header_module.f90->sourcefile~mus_dervarpos_module.f90 sourcefile~mus_bc_header_module.f90->sourcefile~mus_mixture_module.f90 sourcefile~mus_bc_header_module.f90->sourcefile~mus_physics_module.f90 sourcefile~mus_pdf_module.f90 mus_pdf_module.f90 sourcefile~mus_bc_header_module.f90->sourcefile~mus_pdf_module.f90 sourcefile~mus_param_module.f90->sourcefile~mus_physics_module.f90 sourcefile~mus_abortcriteria_module.f90 mus_abortCriteria_module.f90 sourcefile~mus_param_module.f90->sourcefile~mus_abortcriteria_module.f90 sourcefile~mus_species_module.f90->sourcefile~mus_physics_module.f90 sourcefile~mus_varsys_module.f90->sourcefile~mus_physics_module.f90 sourcefile~mus_geom_module.f90 mus_geom_module.f90 sourcefile~mus_varsys_module.f90->sourcefile~mus_geom_module.f90 sourcefile~mus_scheme_type_module.f90 mus_scheme_type_module.f90 sourcefile~mus_varsys_module.f90->sourcefile~mus_scheme_type_module.f90 sourcefile~mus_derquanmsliquid_module.f90->sourcefile~mus_enrtl_module.f90 sourcefile~mus_derquanmsliquid_module.f90->sourcefile~mus_scheme_layout_module.f90 sourcefile~mus_derquanmsliquid_module.f90->sourcefile~mus_varsys_module.f90 sourcefile~mus_derquanmsliquid_module.f90->sourcefile~mus_dervarpos_module.f90 sourcefile~mus_derquanmsliquid_module.f90->sourcefile~mus_physics_module.f90 sourcefile~mus_source_var_module.f90 mus_source_var_module.f90 sourcefile~mus_derquanmsliquid_module.f90->sourcefile~mus_source_var_module.f90 sourcefile~mus_operation_var_module.f90 mus_operation_var_module.f90 sourcefile~mus_derquanmsliquid_module.f90->sourcefile~mus_operation_var_module.f90 sourcefile~mus_statevar_module.f90 mus_stateVar_module.f90 sourcefile~mus_derquanmsliquid_module.f90->sourcefile~mus_statevar_module.f90 sourcefile~mus_derquanmsliquid_module.f90->sourcefile~mus_scheme_type_module.f90 sourcefile~mus_derquanmsliquid_module.f90->sourcefile~mus_scheme_header_module.f90 sourcefile~mus_derquanmsliquid_module.f90->sourcefile~mus_pdf_module.f90 sourcefile~mus_dervarpos_module.f90->sourcefile~mus_scheme_layout_module.f90 sourcefile~mus_mixture_module.f90->sourcefile~mus_enrtl_module.f90 sourcefile~mus_mixture_module.f90->sourcefile~mus_physics_module.f90 sourcefile~mus_mixture_module.f90->sourcefile~mus_scheme_header_module.f90

Files dependent on this one

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

Contents


Source Code

! Copyright (c) 2012-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.
! Copyright (c) 2011-2013 Manuel Hasert <m.hasert@grs-sim.de>
! Copyright (c) 2011 Harald Klimach <harald.klimach@uni-siegen.de>
! Copyright (c) 2011 Konstantin Kleinheinz <k.kleinheinz@grs-sim.de>
! Copyright (c) 2011-2012 Simon Zimny <s.zimny@grs-sim.de>
! Copyright (c) 2012, 2014-2016 Jiaxing Qi <jiaxing.qi@uni-siegen.de>
! Copyright (c) 2012 Kartik Jain <kartik.jain@uni-siegen.de>
! Copyright (c) 2013-2015, 2019 Kannan Masilamani <kannan.masilamani@uni-siegen.de>
! Copyright (c) 2016 Tobias Schneider <tobias1.schneider@student.uni-siegen.de>
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions are met:
!
! 1. Redistributions of source code must retain the above copyright notice,
! this list of conditions and the following disclaimer.
!
! 2. Redistributions in binary form must reproduce the above copyright notice,
! this list of conditions and the following disclaimer in the documentation
! and/or other materials provided with the distribution.
!
! THIS SOFTWARE IS PROVIDED BY THE UNIVERSITY OF SIEGEN “AS IS” AND ANY EXPRESS
! OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
! OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
! IN NO EVENT SHALL UNIVERSITY OF SIEGEN OR CONTRIBUTORS BE LIABLE FOR ANY
! DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
! (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
! ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
! ****************************************************************************** !
!> Boundary condition treatment routines for multispecies simulation
!!
module mus_bc_species_module
  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( &
          & ( posinbuffer-1)* nscalars+ pos )
        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( neigh((idir-1)* nsize+ elempos)+( ifield-1)* qq+ nscalars*0) =  &
            &       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(                                                      &
            &  neigh (( idir-1)* nsize+ elempos)+( ifield-1)* qq+ nscalars*0)
        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( neigh (( idir-1)* nsize+ elempos)+( ifield-1)* qq+ nscalars*0) &
              & = 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( ( ielem-1)* nscalars+ pos ) = &
            & bcBuffer(                                                 &
            & ( posinbuffer-1)* nscalars+ pos)
        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(( ielem-1)* nscalars+ pos )
          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(neigh (( idir-1)* nsize+ elempos)+( ifield-1)* qq+ nscalars*0) &
            & = 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( neigh (( idir-1)* nsize+ elempos)+( ifield-1)* qq+ nscalars*0) &
              & = 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( ( ielem-1)* nscalars+ pos )     &
            & = state((neighpos-1)*nscalars+ idir+(ifieldloc-1)*qq)
          !local element
          !fTmpAll( ( ielem-1)* nscalars+ pos ) = bcBuffer(                    &
          !  & ( globbc%elemlvl( ilevel )%posinbcelembuf%val( ielem )-1)* nscalars+ pos )
        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(neigh (( idir-1)*nsize+ elempos)+( ifield-1)* qq+ nscalars*0) &
            & = 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( neigh (( idir-1)* nsize+ elempos)+( ifield-1)* qq+ nscalars*0)  = &
        &  feq( ( ielem-1)* qq+idir ) !&
        !& +(fTmpAll(( ielem-1)* nscalars+varpos(layout%fstencil%cxdirinv(idir)))&
        !& -feq( ( ielem-1)* layout%fstencil%qq+layout%fstencil%cxdirinv(idir) ) )
        !& fTmpAll(( ielem-1)* nscalars+varpos(idir))
    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( ( ielem-1)* nscalars+ pos )   &
            & = bcBuffer(                                               &
            & ( posinbuffer-1)* nscalars+ pos )
        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(( ielem-1)* nscalars+ pos )
        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(neigh (( idir-1)* nsize+ elempos)+( ifield-1)* qq+ nscalars*0) &
            & = 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( ( ielem-1)* nscalars+ pos )   &
          !  & = bcBuffer(                                               &
          !  & ( posinbuffer-1)* nscalars+ pos )
          fTmp( ( ielem-1)* nscalars+ pos )     &
            & = state((neighpos-1)*nscalars+ idir+(ifieldloc-1)*qq)
        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(neigh (( idir-1)* nsize+ elempos)+( ifield-1)* qq+ nscalars*0 ) &
            & = 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(                                  &
            & ( posinbuffer-1)* nscalars+ pos)
        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(neigh (( idir-1)* nsize+ elempos)+( ifield-1)* qq+ nscalars*0) &
            & = 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( ( ielem-1)* nscalars+ pos )   &
            & = bcBuffer(                                               &
            & ( posinbuffer-1)* nscalars+ pos )
        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(( ielem-1)* nscalars+ pos )
          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(neigh (( idir-1)* nsize+ elempos)+( ifield-1)* qq+ nscalars*0) &
            & = 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( ( ielem-1)* nscalars+ pos )   &
            & = bcBuffer(                                               &
            & ( posinbuffer-1)* nscalars+ pos )
        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(( ielem-1)* nscalars+ pos )
          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(neigh (( idir-1)* nsize+ elempos)+( ifield-1)* qq+ nscalars*0) &
            & = 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(                                   &
            & ( posinbuffer-1)* nscalars+ pos )
          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(                                                               &
            &  neigh (( idir-1)* nsize+ posinstate)+( ifield-1)* qq+ nscalars*0)&
            & = 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( ( ielem-1)* nscalars+ pos ) = &
            & bcBuffer(                                                 &
            & ( posinbuffer-1)* nscalars+ pos )
        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(( ielem-1)* nscalars+ pos )
        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(neigh (( idir-1)* nsize+ elempos)+( ifield-1)* qq+ nscalars*0) &
            & = 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(                    &
       & ( posinbuffer-1)* nscalars+ varpos(idir) )
      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(neigh (( idir-1)* nsize+ elempos)+( ifield-1)* qq+ nscalars*0) &
              & = 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(                                       &
              &  neigh((idir-1)* nsize+ neighpos)+( ifld-1)* qq+ nscalars*0 )
            ! 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( neigh((idir-1)* nsize+ elempos)+( ifield-1)* qq+ nscalars*0)  &
              & = 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(                                       &
              &  neigh((idir-1)* nsize+ neighpos)+( ifld-1)* qq+ nscalars*0 )
            ! 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( neigh((idir-1)* nsize+ elempos)+( ifield-1)* qq+ nscalars*0)  &
              & = 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) :: resi_coeff(varSys%nStateVars, varSys%nStateVars)
    real(kind=rk) :: omega_fac, paramBInv
    integer :: iElem, iDir, QQ, iFld, vPos
    integer :: bcMoleFrac_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(                                    &
              &  neigh((idir-1)* nsize+ neighpos)+( ifld-1)* qq+ nscalars*0 )
            ! 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( neigh((idir-1)* nsize+ elempos)+( ifield-1)* qq+ nscalars*0)  &
              & = 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) :: 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 :: elemPos, 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(                                        &
              &  neigh((idir-1)* nsize+ neighpos_2f)+( ifld-1)* qq+ nscalars*0 )
            ! Pre-Collision values of third fluid element
            fTmpAll_3f(vPos) = state(                                        &
              &  neigh((idir-1)* nsize+ neighpos_3f)+( ifld-1)* qq+ nscalars*0 )
            ! 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( neigh((idir-1)* nsize+ elempos)+( ifield-1)* qq+ nscalars*0)  &
              & = 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) :: 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 :: elemPos, 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(                                       &
              &  neigh((idir-1)* nsize+ neighpos_2f)+( ifld-1)* qq+ nscalars*0 )
            ! Pre-Collision values of third fluid element
            fTmpAll_3f(vPos) = state(                                       &
              &  neigh((idir-1)* nsize+ neighpos_3f)+( ifld-1)* qq+ nscalars*0 )
            ! 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