! Copyright (c) 2012-2022 Kannan Masilamani <kannan.masilamani@dlr.de>
! Copyright (c) 2012-2013 Manuel Hasert <m.hasert@grs-sim.de>
! Copyright (c) 2012-2014 Simon Zimny <s.zimny@grs-sim.de>
! Copyright (c) 2012 Sathish Krishnan P S <s.krishnan@grs-sim.de>
! Copyright (c) 2012, 2014-2017 Jiaxing Qi <jiaxing.qi@uni-siegen.de>
! Copyright (c) 2013, 2015, 2019 Harald Klimach <harald.klimach@uni-siegen.de>
! Copyright (c) 2014-2015 Kartik Jain <kartik.jain@uni-siegen.de>
! Copyright (c) 2015-2017 Tobias Schneider <tobias1.schneider@student.uni-siegen.de>
! Copyright (c) 2016 Verena Krupp <verena.krupp@uni-siegen.de>
! Copyright (c) 2016-2017 Raphael Haupt <raphael.haupt@uni-siegen.de>
! Copyright (c) 2017 Sindhuja Budaraju <nagasai.budaraju@student.uni-siegen.de>
! Copyright (c) 2018, 2020 Jana Gericke <jana.gericke@uni-siegen.de>
! Copyright (c) 2019 Seyfettin Bilgi <seyfettin.bilgi@student.uni-siegen.de>
! Copyright (c) 2019-2020 Peter Vitt <peter.vitt2@uni-siegen.de>
! Copyright (c) 2021-2022 Gregorio Gerardo Spinelli <gregoriogerardo.spinelli@dlr.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.
!  See copyright notice in the COPYRIGHT file.
! ****************************************************************************** !
!> author: Manuel Hasert
!! This module contains general boundary routines
!!
!! like initializing boundary, setting boundary at every time step for each
!! field and other general boundary routines
!!
! 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.
module mus_bc_general_module

  ! include treelm modules
  use env_module,              only: rk, long_k, labelLen, pathLen, newUnit
  use tem_param_module,        only: cs, qOffset_inChar, q000
  use treelmesh_module,        only: treelmesh_type
  use tem_aux_module,          only: tem_abort
  use tem_float_module,        only: operator(.feq.)
  use tem_varSys_module,       only: tem_varSys_type
  use tem_varMap_module,       only: tem_varMap_type
  use tem_bc_prop_module,      only: tem_bc_prop_type
  use tem_topology_module,     only: tem_levelOF
  use tem_debug_module,        only: dbgUnit
  use tem_logging_module,      only: logUnit
  use tem_math_module,         only: invert_matrix
  use tem_geometry_module,     only: tem_BaryOfId,tem_ElemSizeLevel
  use tem_stencil_module,      only: tem_stencilHeader_type
  use tem_bc_module,           only: tem_bc_state_type
  use tem_grow_array_module,   only: init, append, destroy
  use tem_construction_module, only: tem_levelDesc_type
  use tem_timer_module,        only: tem_startTimer, tem_stopTimer

  ! include musubi modules
  use mus_field_module,                 only: mus_field_type
  use mus_field_prop_module,            only: mus_field_prop_type
  use mus_scheme_layout_module,         only: mus_scheme_layout_type
  use mus_scheme_header_module,         only: mus_scheme_header_type
  use mus_scheme_type_module,           only: array2D_type
  use mus_pdf_module,                   only: pdf_data_type
  use mus_bc_header_module,             only: boundary_type,                   &
    &                                         glob_boundary_type,              &
    &                                         mus_set_bcLinks, mus_set_bouzidi,&
    &                                         mus_alloc_bouzidi,               &
    &                                         mus_set_inletUbb,                &
    &                                         mus_set_inletBfl,                &
    &                                         mus_set_nonEqExpol,              &
    &                                         mus_set_outletExpol
  use mus_bc_fluid_module,              only: velocity_eq, vel_neq,            &
    &                                         outlet_nrbc, outlet_nrbc_eq,     &
    &                                         outlet_nrbc_incomp,              &
    &                                         inlet_nrbc, inlet_nrbc_incomp,   &
    &                                         outlet_dnt, pressure_eq,         &
    &                                         pressure_expol, press_neq,       &
    &                                         pressure_antiBounceBack,         &
    &                                         outlet_zero_prsgrd,              &
    &                                         mfr_bounceback, mfr_eq,          &
    &                                         velocity_bounceback,             &
    &                                         velocity_bounceback_incomp,      &
    &                                         velocity_bfl, bc_pdf,            &
    &                                         velocity_bfl_incomp
  use mus_bc_fluid_nonEqExpol_module,   only: velocity_nonEqExpol,             &
    &                                         pressure_nonEqExpol,             &
    &                                         velocity_nonEqExpol_curved,      &
    &                                         velocity_nonEqExpol_incomp,      &
    &                                         pressure_nonEqExpol_incomp,      &
    &                                         velocity_nonEqExpol_curved_incomp
  use mus_bc_fluid_wall_module,         only: slip_wall,                &
    &                                         wall_libb, do_nothing,    &
    &                                         spc_slip_wall
  use mus_bc_fluid_turbulent_module,    only: turbulent_wall,                  &
    &                                  turbulent_wall_libb,                    &
    &                                  turbulent_wall_eq,                      &
    &                                  turbulent_wall_eq_curved,               &
    &                                  turbulent_wall_noneq_expol,             &
    &                                  turbulent_wall_noneq_expol_curved,      &
    &                                  turbulent_wall_eq_incomp,               &
    &                                  turbulent_wall_eq_curved_incomp,        &
    &                                  turbulent_wall_noneq_expol_incomp,      &
    &                                  turbulent_wall_noneq_expol_curved_incomp
  use mus_bc_fluid_experimental_module, only: pressure_expol_slow
  use mus_bc_fluid_module,              only: pressure_momentsbased,           &
    &                                         moments_wall,                    &
    &                                         velocity_momentsbased,           &
    &                                         velocity_momentsbased_incomp,    &
    &                                         pressure_momentsbased_incomp
  use mus_bc_fluid_experimental_module, only: moments_inflow, moments_outflow, &
    &                                         spc_moments_outflow,             &
    &                                         spc_bb_wall, spc_bb_vel_test
  use mus_bc_passiveScalar_module,      only: outlet_pasScal, inlet_pasScal
  use mus_bc_species_module,            only: spc_outlet_zero_prsgrd,          &
    &                                         spc_moleFrac, spc_moleFlux,      &
    &                                         spc_moleFrac_wtdf,               &
    &                                         spc_moleFrac_eq, spc_inlet,      &
    &                                         spc_moleDens_eq,                 &
    &                                         spc_mole_fraction_noneq_expol,   &
    &                                         spc_velocity_noneq_expol,        &
    &                                         spc_moleDiff_Flux, spc_inlet_eq, &
    &                                         spc_outlet_eq, spc_outlet_vel,   &
    &                                         spc_outlet_expol,                &
    &                                         spc_outflow,                     &
    &                                         spc_inflow,                      &
    &                                         spc_solvent_inflow,              &
    &                                         spc_solvent_outflow,             &
    &                                         spc_moleFlux_eq, spc_vel_bb,     &
    &                                         spc_blackbox_mem_ion,            &
    &                                         spc_blackbox_mem_solvent,        &
    &                                         spc_moments_moleFrac,            &
    &                                         spc_moments_moleFlux,            &
    &                                         spc_moments_wall,                &
    &                                         spc_moments_vel
  use mus_bc_poisson_module,            only: potential_nonEqExpol,        &
    &                                         potential_nonEqExpol_curved, &
    &                                         potential_neumann,           &
    &                                         potential_neumann_curved
  use mus_bc_nernstPlanck_module,       only: moleDens_nonEqExpol,        &
    &                                         moleDens_nonEqExpol_curved, &
    &                                         moleDens_neumann,           &
    &                                         moleDens_neumann_curved

  use mus_relaxationParam_module,       only: mus_viscosity_type
  use mus_auxField_module,              only: mus_auxFieldVar_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_timer_module,                 only: mus_timerHandles

  implicit none

  private

  public :: set_boundary
  public :: mus_init_boundary
  public :: mus_get_points_fromBC

contains


! ****************************************************************************** !
  !> Call the functions associated with each boundary condition
  !!
  !! Loop over each field and Run over all the boundary conditions for
  !! current iLevel and call the function pointer.
  !! The function pointer was before assigned in init_boundary
  !! This routine is being called from the
  !! [[mus_control_module:do_recursive_multilevel]] "main control routine"
  !! before the compute (advection_relaxation) kernel call
  subroutine set_boundary( field, pdf, levelDesc, tree, iLevel, nBCs, params, &
    &                      layout, varSys, derVarPos, globBC, mixture,        &
    &                      physics, state )
    ! ---------------------------------------------------------------------------
    !> fluid parameters and properties
    type( mus_field_type ), intent(inout) :: field(:)
    !> contains global state vector
    type( pdf_data_type ), intent(inout) :: pdf
    !> state arrays fo current iLevel both now and next
    real(kind=rk), intent(inout) :: state(:,:)
    !> global type contains iLevel descriptor
    type( tem_levelDesc_type ), intent(in) :: levelDesc
    !> global treelm mesh
    type( treelmesh_type ), intent(in) ::tree
    !> the iLevel on which this boundary was invoked
    integer, intent(in) :: iLevel
    !> number of BC
    integer, intent(in) :: nBCs
    !> global parameters
    type(mus_param_type),intent(in) :: params
    !> stencil layout information
    type( mus_scheme_layout_type ), intent(in) ::layout
    !> 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
    !> mixture info
    type(mus_mixture_type), intent(in) :: mixture
    ! --------------------------------------------------------------------------
    integer :: iField, nFields, iBnd
    ! --------------------------------------------------------------------------

    nFields = size( field )

    if ( nBCs > 0 ) then
      call tem_startTimer( timerHandle = mus_timerHandles%bcBuffer(iLevel) )
      call fill_bcBuffer(                                 &
        &    currState = state( :, pdf%nNext ),           &
        &    neigh     = pdf%neigh,                       &
        &    bcBuffer  = pdf%bcBuffer,                    &
        &    nSize     = pdf%nSize,                       &
        &    nElems_bc = levelDesc%bc_elemBuffer%nVals,   &
        &    posInTotal= levelDesc%bc_elemBuffer%val,     &
        &    nFields   = nFields,                         &
        &    QQ        = layout%fStencil%QQ,              &
        &    varSys    = varSys                           )
      call tem_stopTimer( timerHandle = mus_timerHandles%bcBuffer(iLevel) )

      ! @todo: As neigh buffer is inside
      !        field( iField )%bc( iBnd )%neigh( iLevel ),
      !        we may fill it in the following fileds and BCs loop.
      call fill_neighBuffer(                       &
        &    prevstate  = state( :, pdf%nNow ),    &
        &    currstate  = state( :, pdf%nNext ),   &
        &    neigh  = pdf%neigh,                   &
        &    globBC = globBC,                      &
        &    nBCs   = nBCs,                        &
        &    field  = field,                       &
        &    varSys = varSys,                      &
        &    QQ     = layout%fStencil%QQ,          &
        &    nSize  = pdf%nSize,                   &
        &    iLevel = iLevel                       )

      ! loop over fields
      do iField = 1, nFields

        ! Treat all boundary conditions
        do iBnd = 1, nBCs

          ! write(dbgUnit(10),*) 'Do boundary condition: '// trim( globBC( iBnd )%label )

          call tem_startTimer( timerHandle =  mus_timerHandles%setBnd(iBnd) )
          call field(iField)%bc( iBnd )%fnct(                   &
            &      state       = state( :, pdf%nNext ),         &
            &      nSize       = pdf%nSize,                     &
            ! we want to write into the current time step
            ! the field pointers have been inverted, so the newest
            ! information is stored in nNow
            &      bcBuffer    = pdf%bcBuffer,                  &
            &      globBC      = globBC( iBnd ),                &
            &      levelDesc   = levelDesc,                     &
            &      tree        = tree,                          &
            &      iLevel      = iLevel,                        &
            &      sim_time    = params%general%simControl%now, &
            ! &      params      = params,                        &
            &      physics     = physics,                       &
            &      mixture     = mixture,                       &
            &      iField      = iField,                        &
            &      neigh       = pdf%neigh,                     &
            &      layout      = layout,                        &
            &      fieldProp   = field(iField)%fieldProp,       &
            &      varSys      = varSys,                        &
            &      derVarPos   = derVarPos(iField),             &
            ! each field may has multiple state variables
            &      varPos      = varSys%method%val(iField)      &
            &                      %state_varPos,               &
            &      nScalars    = varSys%nScalars                )

          call tem_stopTimer( timerHandle =  mus_timerHandles%setBnd(iBnd) )
        end do  ! iBnd
      end do  ! iField
    end if
  end subroutine set_boundary
! ****************************************************************************** !


! ****************************************************************************** !
  ! > This routine call init_boundary to initialize boundary for each field
  !!
  subroutine mus_init_boundary( field, pdf, tree, levelDesc, layout,     &
    &                           schemeHeader, varSys, derVarPos, globBC, &
    &                           bc_prop, state, auxField )
    !---------------------------------------------------------------------------
    !> fluid parameters and properties
    type( mus_field_type ), intent(inout) :: field(:)
    !> global treelm mesh
    type( treelmesh_type ), intent(in) :: tree
    !> contains global state vector
    type( pdf_data_type ), intent(inout) :: pdf(tree%global%minLevel &
      &                                         :tree%global%maxLevel)
    !> state array
    type( array2D_type ), intent(inout) :: state(tree%global%minLevel &
      &                                          :tree%global%maxLevel)
    !> AuxField array
    type( mus_auxFieldVar_type), intent(in) :: auxField(tree%global%minLevel &
      &                                                 :tree%global%maxLevel)
    !> scheme layout type
    type( mus_scheme_layout_type ), intent(in) ::layout
    !> scheme header info
    type( mus_scheme_header_type ), intent(in) :: schemeHeader
    !> Level Descriptor
    type( tem_leveldesc_type ), intent(in) :: levelDesc(tree%global%minLevel:tree%global%maxLevel)
    !> 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(inout) :: globBC(:)
    !> boundary property type
    type( tem_bc_prop_type ) :: bc_prop
    ! ---------------------------------------------------------------------------
    integer :: iField, iBC, iLevel, minLevel, maxLevel
    ! ---------------------------------------------------------------------------

    minLevel = tree%global%minLevel
    maxLevel = tree%global%maxLevel

    if ( bc_prop%nBCtypes > 0 ) then
      !> Check prerequisite for multi-species boundary conditions
      if (size(field)>1) then
        call check_BCs_preRequisite_MS( field, bc_prop%nBCtypes )
      end if

      do iField = 1, size(field)
        ! Set up links which are needed to be updated
        write(logUnit(7), *) 'Counting actual number of links that need '//&
          &                  'to be update:'
        do iBC = 1, bc_prop%nBCtypes
          if ( .not. globBC( iBC )%isWall ) then
            do iLevel = minLevel, maxLevel
              if ( globBC(iBC)%nElems(iLevel) >= 0 ) then
                write(logUnit(7), "(2(A,I0))") &
                  &  trim(globBC(iBC)%label)//', level: ', iLevel, &
                  &  ', local nElems: ', globBC(iBC)%nElems(iLevel)
                call mus_set_bcLinks( iField   = iField,                     &
                  &                   QQ       = layout%fStencil%QQ,         &
                  &                   QQN      = layout%fStencil%QQN,        &
                  &                   nScalars = varSys%nScalars,            &
                  &                   nElems   = globBC(iBC)%nElems(iLevel), &
                  &                   nSize    = pdf(iLevel)%nSize,          &
                  &                   elemLvl  = globBC(iBC)%elemLvl(iLevel),&
                  &                   neigh    = pdf(iLevel)%neigh,          &
                  &                   links    = field(iField)%bc(iBC)       &
                  &                                           %links(iLevel) )
                write(logUnit(7), "(A,I0)") &
                  &     ' nLinks: ', field(iField)%bc(iBC)%links(iLevel)%nVals
              end if
            end do ! iLevel
          end if
        end do ! iBC

        call init_boundary_single(                                     &
          &    bc           = field(iField)%bc,                        &
          &    pdf          = pdf(minLevel:maxLevel),                  &
          &    state        = state(minLevel:maxLevel),                &
          &    auxField     = auxField(minLevel:maxLevel),             &
          &    tree         = tree,                                    &
          &    leveldesc    = levelDesc,                               &
          &    layout       = layout,                                  &
          &    schemeHeader = schemeHeader,                            &
          &    varPos       = varSys%method%val(iField)%state_varPos,  &
          &    varSys       = varSys,                                  &
          &    dervarPos    = derVarPos(iField),                       &
          &    globBC       = globBC,                                  &
          &    bc_prop      = bc_prop,                                 &
          &    fieldProp    = field(iField)%fieldProp                  )
      end do ! iField
    end if

  end subroutine mus_init_boundary
! ****************************************************************************** !


! ****************************************************************************** !
  !> This subroutine sets the right boundary conditions for the different
  !! boundaries.
  !!
  subroutine init_boundary_single(bc, pdf, tree, levelDesc, layout,          &
    &                             schemeHeader, varPos, varSys, derVarPos,   &
    &                             globBC,bc_prop, state, auxField, fieldProp )
    ! ---------------------------------------------------------------------------
    !> global array boundary type
    type( boundary_type ) :: bc(:)
    !> global treelm mesh
    type( treelmesh_type ), intent(in) ::tree
    !> contains global state vector
    type( pdf_data_type ), intent(in) :: pdf( tree%global%minLevel &
      &                                     : tree%global%maxLevel )
    !> contains global state vector
    type( array2D_type ), intent(in) :: state( tree%global%minLevel &
      &                                      : tree%global%maxLevel )
    !> AuxField array
    type( mus_auxFieldVar_type), intent(in) :: auxField(tree%global%minLevel &
      &                                                 :tree%global%maxLevel)
    !> scheme layout type
    type( mus_scheme_layout_type ), intent(in) ::layout
    !> scheme header info
    type( mus_scheme_header_type ), intent(in) :: schemeHeader
    !> global pdf type
    type( tem_leveldesc_type ), intent(in) :: levelDesc( tree%global%minLevel &
      &                                                : tree%global%maxLevel )
    !> varPos of current field variable
    integer, intent(in) :: varPos(:)
    !> 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(inout) :: globBC(:)
    !> boundary property type
    type( tem_bc_prop_type ), intent(in) :: bc_prop
    !> fluid parameters and properties
    type( mus_field_prop_type ), intent(in) :: fieldProp
    ! ---------------------------------------------------------------------------
    integer :: iBnd, iLevel, minLevel, maxLevel, nBCs
    logical :: isWall
    ! ---------------------------------------------------------------------------

    minLevel = tree%global%minLevel
    maxLevel = tree%global%maxLevel
    nBCs     = bc_prop%nBCtypes

    !> The boundary conditions were read in the order of config.lua
    !! the order in the boundary condition description file however
    !! might be different.
    !! We have to match the labels.
    write(logUnit(1),'(A,I0,A)') 'Initialize ', nBCs, ' boundary conditions'

    do iBnd = 1, nBCs
      if( trim(bc( iBnd )%label ) /= trim( globBC( iBnd )%label )) then
        write(logUnit(1),*) 'Error: The boundary with label from the mesh:'    &
          &                 //trim( globBC(iBnd)%label )
        write(logUnit(1),*) 'does not match the one specified in the lua file:'&
          &                 //trim(bc(iBnd)%label)
        write(logUnit(1),*) 'Stopping ...'
        call tem_abort()
      end if

      write(logUnit(1),*) 'Initializing BC: '//trim(bc(iBnd)%label)

      isWall = .false.
      select case (trim(bc(iBnd)%BC_kind))
      case('wall', 'symmetry')
        isWall = .true.
        bc( iBnd )%fnct => do_nothing
      case('wall_libb', 'velocity_bounceback', 'velocity_bfl')
        bc( iBnd )%evalBcVar_link = .true.

        ! Here we have to allocate and set the q-values if it is not
        ! provided by seeder and set the qVal from musubi.lua
        if( .not. globBC(iBnd)%hasQVal .and.        &
          & .not. globBC(iBnd)%qValInitialized ) then
          call init_qVals( refQval   = bc( iBnd )%qVal, &
            &              minLevel  = minLevel,        &
            &              maxLevel  = maxLevel,        &
            &              layout    = layout,          &
            &              globBC    = globBC(iBnd)     )
        end if

        select case (trim(bc( iBnd )%BC_kind))
        case('wall_libb')
          isWall = .true.
          bc( iBnd )%fnct => wall_libb
          ! set link-wise data
          call mus_alloc_bouzidi( me       = bc(iBnd)%bouzidi,       &
            &    nVals    = bc(iBnd)%links(minLevel:maxLevel)%nVals, &
            &    minLevel = minLevel,                                &
            &    maxLevel = maxLevel                                 )
          do iLevel = minLevel, maxLevel
            call mus_set_bouzidi( iLevel   = iLevel,                    &
              &                   QQ       = layout%fStencil%QQ,        &
              &                   QQN      = layout%fStencil%QQN,       &
              &                   nScalars = varSys%nScalars,           &
              &                   globBC   = globBC(iBnd),              &
              &                   cxDirInv = layout%fStencil%cxDirInv,  &
              &                   varPos   = varPos,                    &
              &                   bouzidi  = bc( iBnd )%bouzidi(iLevel) )
          end do
        case('velocity_bounceback')
          select case(trim(schemeHeader%kind))
          case('fluid')
            bc( iBnd )%fnct => velocity_bounceback
          case('fluid_incompressible')
            bc( iBnd )%fnct => velocity_bounceback_incomp
          case default
            call tem_abort('Unknown scheme kind for velocity_bounceback')
          end select
          ! set link-wise data
          ! bc(iBnd)%links is built in mus_set_bcLinks
          call mus_set_inletUbb( me           = bc(iBnd)%inletUbbQVal, &
            &                    tree         = tree,                  &
            &                    stencil      = layout%fStencil,       &
            &                    nScalars     = varSys%nScalars,       &
            &                    globBC       = globBC(iBnd),          &
            &                    levelDesc    = levelDesc(minLevel:maxLevel), &
            &                    varPos       = varPos,                &
            &        nLinks = bc(iBnd)%links(minLevel:maxLevel)%nVals, &
            &                    minLevel     = minLevel,              &
            &                    maxLevel     = maxLevel               )
        case('velocity_bfl')
          select case(trim(schemeHeader%kind))
          case('fluid')
            bc( iBnd )%fnct => velocity_bfl
          case('fluid_incompressible')
            bc( iBnd )%fnct => velocity_bfl_incomp
          case default
            call tem_abort('Unknown scheme kind for velocity_bfl')
          end select
          ! set link-wise data
          call mus_set_inletBfl(                                       &
            &     me        = bc(iBnd)%inletBfl,                       &
            &     tree      = tree,                                    &
            &     stencil   = layout%fStencil,                         &
            &     nScalars  = varSys%nScalars,                         &
            &     globBC    = globBC(iBnd),                            &
            &     levelDesc = levelDesc(minLevel:maxLevel),            &
            &     varPos    = varPos,                                  &
            &     nLinks    = bc(iBnd)%links(minLevel:maxLevel)%nVals, &
            &     minLevel  = minLevel,                                &
            &     maxLevel  = maxLevel                                 )
        end select

      case('turbulent_wall', 'turbulent_wall_noneq_expol', 'turbulent_wall_eq')
        isWall = .true.
        ! Here we have to allocate and set the q-values if it is not
        ! provided by seeder and set the qVal from musubi.lua
        if (.not. globBC(iBnd)%hasQVal .and.       &
          & .not. globBC(iBnd)%qValInitialized) then
          call init_qVals( refQval   = bc( iBnd )%qVal, &
            &              minLevel  = minLevel,        &
            &              maxLevel  = maxLevel,        &
            &              layout    = layout,          &
            &              globBC    = globBC(iBnd)     )
        end if

        ! allocate turbulent viscosity on boundary elements,
        ! and friction velocity and normal distance to boundary on first neigbor
        ! of boundary elements.
        ! Initialize friction velocity from stream-wise
        ! velocity component computed on first neighbor and normal distance to
        ! to boundary.
        ! \todo KM: Calculate friction velocity from wall shear
        ! stress computed on the first neighbor using non-equilibrium pdf
        write(logUnit(10), "(A)") 'Initializing turbulent_wall ...'
        allocate(bc(iBnd)%turbwallFunc%dataOnLvl(minLevel: maxLevel))
        do iLevel = minLevel, maxLevel
          call mus_init_turb_wallFunc( bc        = bc(iBnd),                 &
            &                          globBC    = globBC(iBnd),             &
            &                          auxField  = auxField(iLevel)%val(:),  &
            &                          viscKine  = fieldProp%fluid%viscKine, &
            &                          derVarPos = derVarPos,                &
            &                          varSys    = varSys,                   &
            &                          stencil   = layout%fStencil,          &
            &                          iLevel    = iLevel                    )
        end do

        select case (trim(bc(iBnd)%BC_kind))
        case ('turbulent_wall_noneq_expol')
          select case(trim(schemeHeader%kind))
          case('fluid')
            if (bc(iBnd)%curved) then
              bc( iBnd )%fnct => turbulent_wall_noneq_expol_curved
            else
              bc( iBnd )%fnct => turbulent_wall_noneq_expol
            end if
          case('fluid_incompressible')
            if (bc(iBnd)%curved) then
              bc( iBnd )%fnct => turbulent_wall_noneq_expol_curved_incomp
            else
              bc( iBnd )%fnct => turbulent_wall_noneq_expol_incomp
            end if
          case default
            call tem_abort('Unknown scheme kind for '//trim(bc(iBnd)%BC_kind))
          end select
        case ('turbulent_wall_eq')
          select case(trim(schemeHeader%kind))
          case('fluid')
            if (bc(iBnd)%curved) then
              bc( iBnd )%fnct => turbulent_wall_eq_curved
            else
              bc( iBnd )%fnct => turbulent_wall_eq
            end if
          case('fluid_incompressible')
            if (bc(iBnd)%curved) then
              bc( iBnd )%fnct => turbulent_wall_eq_curved_incomp
            else
              bc( iBnd )%fnct => turbulent_wall_eq_incomp
            end if
          case default
            call tem_abort('Unknown scheme kind for '//trim(bc(iBnd)%BC_kind))
          end select
        case ('turbulent_wall')
          select case(trim(schemeHeader%kind))
          case('fluid')
            if (bc(iBnd)%curved) then
              bc( iBnd )%fnct => turbulent_wall_libb
              ! set link-wise data for bouzidi wall linear interpolation
              call mus_alloc_bouzidi( me       = bc(iBnd)%bouzidi,       &
                &    nVals    = bc(iBnd)%links(minLevel:maxLevel)%nVals, &
                &    minLevel = minLevel,                                &
                &    maxLevel = maxLevel                                 )
              do iLevel = minLevel, maxLevel
                call mus_set_bouzidi( iLevel   = iLevel,                    &
                  &                   QQ       = layout%fStencil%QQ,        &
                  &                   QQN      = layout%fStencil%QQN,       &
                  &                   nScalars = varSys%nScalars,           &
                  &                   globBC   = globBC(iBnd),              &
                  &                   cxDirInv = layout%fStencil%cxDirInv,  &
                  &                   varPos   = varPos,                    &
                  &                   bouzidi  = bc( iBnd )%bouzidi(iLevel) )
              end do
            else
              bc( iBnd )%fnct => turbulent_wall
            end if
          case default
            call tem_abort('Unknown scheme kind for '//trim(bc(iBnd)%BC_kind))
          end select
        end select

      case('velocity_eq')
        bc( iBnd )%fnct => velocity_eq
      case('mfr_bounceback')
        bc( iBnd )%fnct => mfr_bounceback
      case('mfr_eq')
        bc( iBnd )%fnct => mfr_eq
      case('vel_neq')
        bc( iBnd )%fnct => vel_neq
      case('bc_pdf')
        bc( iBnd )%fnct => bc_pdf
      case('outlet_nrbc','outlet_nrbc_eq', 'inlet_nrbc')
        select case (trim(bc( iBnd )%BC_kind))
          case( 'outlet_nrbc' )
            select case(trim(schemeHeader%kind))
            case ('fluid')
              bc( iBnd )%fnct => outlet_nrbc
            case('fluid_incompressible')
              bc( iBnd )%fnct => outlet_nrbc_incomp
            case default
              call tem_abort('Unknown scheme kind for outlet_nrbc')
            end select
          case( 'outlet_nrbc_eq' )
            bc( iBnd )%fnct => outlet_nrbc_eq
          case( 'inlet_nrbc' )
            select case(trim(schemeHeader%kind))
            case ('fluid')
              bc( iBnd )%fnct => inlet_nrbc
            case('fluid_incompressible')
              bc( iBnd )%fnct => inlet_nrbc_incomp
            case default
              call tem_abort('Unknown scheme kind for inlet_nrbc')
            end select
        end select

        bc( iBnd )%nrbc%cs_mod = cs * sqrt( bc( iBnd )%nrbc%kappa )

        write(logUnit(5),"(A)")       'Initializing NRBC with parameters:'
        do iLevel = minLevel, maxLevel
          call init_nrbc( bc        = bc(iBnd),                                 &
            &             state     = state(iLevel)%val(:,pdf(iLevel)%nNext),   &
            &             nSize     = pdf(iLevel)%nSize,                        &
            &             neigh     = pdf(iLevel)%neigh(:),                     &
            &             layout    = layout,                                   &
            &             level     = iLevel,                                   &
            &             nScalars  = varSys%nScalars,                          &
            &             varSys    = varSys,                                   &
            &             derVarPos = derVarPos,                                &
            &             elemPos   = globBC(iBnd)%elemLvl(iLevel)%elem%val(:), &
            &             nElems    = globBC(iBnd)%nElems(iLevel)               )
        end do
      case('outlet_dnt')
        bc( iBnd )%fnct => outlet_dnt
      case('pressure_expol', 'pressure_expol_slow')
        ! set link-wise data
        call mus_set_outletExpol(                                             &
          &            me          = bc(iBnd)%outletExpol,                    &
          &            stencil     = layout%fStencil,                         &
          &            globBC      = globBC(iBnd),                            &
          &            nLinks      = bc(iBnd)%links(minLevel:maxLevel)%nVals, &
          &            minLevel    = minLevel,                                &
          &            maxLevel    = maxLevel                                 )
        if( trim(bc(iBnd)%BC_kind) == 'pressure_expol' ) then
          bc( iBnd )%fnct => pressure_expol
        else
          bc( iBnd )%fnct => pressure_expol_slow
        end if
      case('press_neq')
        bc( iBnd )%fnct => press_neq
      case('pressure_eq')
        bc( iBnd )%fnct => pressure_eq
      case('pressure_antibounceback')
        bc( iBnd )%fnct => pressure_antiBounceBack
      case('outlet_zero_prsgrd')
        bc( iBnd )%fnct => outlet_zero_prsgrd
      ! passive scalar boundary conditions
      case('flekkoy_inlet')
        bc( iBnd )%fnct => inlet_pasScal
      case('flekkoy_outlet')
        bc( iBnd )%fnct => outlet_pasScal
      ! multispecies boundary conditions
      case('spc_slip_wall')
        isWall = .true.
        bc( iBnd )%fnct => spc_slip_wall
      case('slip_wall')
        isWall = .true.
        bc( iBnd )%fnct => slip_wall
      case('spc_bb_wall')
        bc( iBnd )%fnct => spc_bb_wall
      case('spc_bb_vel_test')
        bc( iBnd )%fnct => spc_bb_vel_test
      case('spc_outlet_zero_prsgrd')
        bc( iBnd )%fnct => spc_outlet_zero_prsgrd
      case('spc_inlet_eq')
        bc( iBnd )%fnct => spc_inlet_eq
      case('spc_inlet')
        bc( iBnd )%fnct => spc_inlet
      case('spc_inflow')
        bc( iBnd )%fnct => spc_inflow
      case('spc_solvent_inflow')
        bc( iBnd )%fnct => spc_solvent_inflow
      case('spc_velocity_noneq_expol')
        bc( iBnd )%fnct => spc_velocity_noneq_expol
      case('spc_mole_fraction_noneq_expol')
        bc( iBnd )%fnct => spc_mole_fraction_noneq_expol
      case('spc_vel_bb')
        bc( iBnd )%fnct => spc_vel_bb
      case('spc_outlet_eq')
        bc( iBnd )%fnct => spc_outlet_eq
      case('spc_outlet_expol')
        bc( iBnd )%fnct => spc_outlet_expol
      case('spc_outflow')
        bc( iBnd )%fnct => spc_outflow
      case('spc_solvent_outflow')
        bc( iBnd )%fnct => spc_solvent_outflow
      case('spc_outlet_vel')
        bc( iBnd )%fnct => spc_outlet_vel
      case('spc_molefrac')
        bc( iBnd )%fnct => spc_moleFrac
      case('spc_molefrac_eq')
        bc( iBnd )%fnct => spc_moleFrac_eq
      case('spc_moledens_eq')
        bc( iBnd )%fnct => spc_moleDens_eq
      case('spc_molefrac_wtdf')
        bc( iBnd )%fnct => spc_moleFrac_wtdf
      case('spc_moleflux')
        bc( iBnd )%fnct => spc_moleFlux
      case('spc_moleflux_eq')
        bc( iBnd )%fnct => spc_moleFlux_eq
      case('spc_molediff_flux')
        bc( iBnd )%fnct => spc_moleDiff_Flux
      case( 'spc_blackbox_mem_ion')
        bc( iBnd )%fnct => spc_blackbox_mem_ion
      case( 'spc_blackbox_mem_solvent')
        bc( iBnd )%fnct => spc_blackbox_mem_solvent

      ! cases which use the non-equilibrium extrapolation (nonEqExpol)
      case('potential_noneq_expol', 'potential_neumann',    &
        &  'velocity_noneq_expol', 'pressure_noneq_expol' )
        bc( iBnd )%evalBcVar_link = .true.

        ! Here we have to allocate and set the q-values if it is not
        ! provided by seeder and set the qVal from musubi.lua
        if( .not. globBC(iBnd)%hasQVal .and.        &
          & .not. globBC(iBnd)%qValInitialized ) then
          call init_qVals( refQval   = bc(iBnd)%qVal, &
            &              minLevel  = minLevel,      &
            &              maxLevel  = maxLevel,      &
            &              layout    = layout,        &
            &              globBC    = globBC(iBnd)   )
        end if

        ! set link-wise data
        call mus_set_nonEqExpol(                                     &
          &     me        = bc(iBnd)%nonEqExpol,                     &
          &     curved    = bc(iBnd)%curved,                         &
          &     tree      = tree,                                    &
          &     stencil   = layout%fStencil,                         &
          &     nScalars  = varSys%nScalars,                         &
          &     globBC    = globBC(iBnd),                            &
          &     bc_neigh  = bc(iBnd)%neigh(minLevel:maxLevel),       &
          &     pdf       = pdf(minLevel:maxLevel),                  &
          &     levelDesc = levelDesc(minLevel:maxLevel),            &
          &     varPos    = varPos,                                  &
          &     nLinks    = bc(iBnd)%links(minLevel:maxLevel)%nVals, &
          &     minLevel  = minLevel,                                &
          &     maxLevel  = maxLevel                                 )

        select case ( trim(bc(iBnd)%BC_kind ) )
        case('potential_noneq_expol')
          if (bc(iBnd)%curved) then
            bc( iBnd )%fnct => potential_nonEqExpol_curved
          else
            bc( iBnd )%fnct => potential_nonEqExpol
          end if
        case('potential_neumann')
          if (bc(iBnd)%curved) then
            bc( iBnd )%fnct => potential_neumann_curved
          else
            bc( iBnd )%fnct => potential_neumann
          end if
        case('velocity_noneq_expol')
          if (trim(schemeHeader%relaxation(1:3)) == 'trt') then
            call tem_abort('velocity_noneq_expol is not supported for trt!')
          end if
          select case(trim(schemeHeader%kind))
          case('fluid')
            if (bc(iBnd)%curved) then
              bc( iBnd )%fnct => velocity_nonEqExpol_curved
            else
              bc( iBnd )%fnct => velocity_nonEqExpol
            end if
          case('fluid_incompressible')
            if (bc(iBnd)%curved) then
              bc( iBnd )%fnct => velocity_nonEqExpol_curved_incomp
            else
              bc( iBnd )%fnct => velocity_nonEqExpol_incomp
            end if
          case default
            call tem_abort('Unknown scheme kind for velocity_noneq_expol')
          end select
        case('pressure_noneq_expol')
          if (trim(schemeHeader%relaxation(1:3)) == 'trt') then
            call tem_abort('velocity_noneq_expol is not supported for trt!')
          end if
          select case(trim(schemeHeader%kind))
          case('fluid')
            bc( iBnd )%fnct => pressure_nonEqExpol
          case('fluid_incompressible')
            bc( iBnd )%fnct => pressure_nonEqExpol_incomp
          case default
            call tem_abort('Unknown scheme kind for pressure_noneq_expol')
          end select
        end select

      ! cases which use the nernst_planck
      case('moledens_noneq_expol', 'moledens_neumann')
        bc( iBnd )%evalBcVar_link = .true.

        ! Here we have to allocate and set the q-values if it is not
        ! provided by seeder and set the qVal from musubi.lua
        if( .not. globBC(iBnd)%hasQVal .and.        &
          & .not. globBC(iBnd)%qValInitialized ) then
          call init_qVals( refQval   = bc( iBnd )%qVal, &
            &              minLevel  = minLevel,        &
            &              maxLevel  = maxLevel,        &
            &              layout    = layout,          &
            &              globBC    = globBC(iBnd)     )
        end if

        ! set link-wise data
        call mus_set_nonEqExpol(                                     &
          &     me        = bc(iBnd)%nonEqExpol,                     &
          &     curved    = bc(iBnd)%curved,                         &
          &     tree      = tree,                                    &
          &     stencil   = layout%fStencil,                         &
          &     nScalars  = varSys%nScalars,                         &
          &     globBC    = globBC(iBnd),                            &
          &     bc_neigh  = bc(iBnd)%neigh(minLevel:maxLevel),       &
          &     pdf       = pdf(minLevel:maxLevel),                  &
          &     levelDesc = levelDesc(minLevel:maxLevel),            &
          &     varPos    = varPos,                                  &
          &     nLinks    = bc(iBnd)%links(minLevel:maxLevel)%nVals, &
          &     minLevel  = minLevel,                                &
          &     maxLevel  = maxLevel                                 )

        select case ( trim(bc(iBnd)%BC_kind ) )
        case('moleDens_noneq_expol')
          if (bc(iBnd)%curved) then
            bc( iBnd )%fnct => moleDens_nonEqExpol_curved
          else
            bc( iBnd )%fnct => moleDens_nonEqExpol
          end if
        case('moledens_neumann')
          if (bc(iBnd)%curved) then
            bc( iBnd )%fnct => moleDens_neumann_curved
          else
            bc( iBnd )%fnct => moleDens_neumann
          end if
        end select

      case( 'pressure_momentsbased', 'moments_wall', 'velocity_momentsbased', &
        &   'moments_inflow', 'moments_outflow',                              &
        &   'spc_moments_molefrac', 'spc_moments_moleflux',                   &
        &   'spc_moments_wall', 'spc_moments_vel', 'spc_moments_outflow'      )
        select case (trim(bc( iBnd  )%BC_kind))
        case( 'pressure_momentsbased')
          select case (trim(schemeHeader%kind))
          case ('fluid')
            bc( iBnd )%fnct => pressure_momentsbased
          case ('fluid_incompressible')
            bc( iBnd )%fnct => pressure_momentsbased_incomp
          case default
            write(logUnit(1),*) 'Chosen boundary kind '                        &
              &                 //trim(bc(iBnd)%BC_kind)//' is not supported'  &
              &                //'scheme kind'//trim(schemeHeader%kind)
            call tem_abort()
          end select
        case( 'moments_wall')
          bc( iBnd )%fnct => moments_wall
        case( 'velocity_momentsbased')
          select case (trim(schemeHeader%kind))
          case ('fluid')
            bc( iBnd )%fnct => velocity_momentsbased
          case('fluid_incompressible')
            bc( iBnd )%fnct => velocity_momentsbased_incomp
          case default
            write(logUnit(1),*) 'Chosen boundary kind '                        &
              &                 //trim(bc(iBnd)%BC_kind)//' is not supported'  &
              &                //'scheme kind'//trim(schemeHeader%kind)
            call tem_abort()
          end select
        case( 'moments_inflow')
          bc( iBnd )%fnct => moments_inflow
        case( 'moments_outflow')
          bc( iBnd )%fnct => moments_outflow
        case( 'spc_moments_molefrac')
          bc( iBnd )%fnct => spc_moments_moleFrac
        case( 'spc_moments_moleflux')
          bc( iBnd )%fnct => spc_moments_moleFlux
        case( 'spc_moments_outflow')
          bc( iBnd )%fnct => spc_moments_outflow
        case( 'spc_moments_wall')
          bc( iBnd )%fnct => spc_moments_wall
        case( 'spc_moments_vel')
          bc( iBnd )%fnct => spc_moments_vel
        end select
        call init_momentsBC( bc        = bc(iBnd),     &
          &                  leveldesc = levelDesc,    &
          &                  layout    = layout,       &
          &                  globBC    = globBC(iBnd), &
          &                  minLevel  = minLevel,     &
          &                  maxLevel  = maxLevel      )
      case default
        write(logUnit(1),*)'BC type '//trim(bc( iBnd  )%bc_kind)// &
          & ' for label '//trim(bc( iBnd  )%label)// 'not found'
        call tem_abort()
      end select

      if ( .not. isWall ) then
        ! setup indices for each boundaries
        call mus_setupIndices_forBC( bc       = bc(iBnd),        &
          &                          globBC   = globBC(iBnd),    &
          &                          tree     = tree,            &
          &                          stencil  = layout%fStencil, &
          &                          levelDesc= levelDesc,       &
          &                          varSys   = varSys,          &
          &                          minLevel = minLevel,        &
          &                          maxLevel = maxLevel         )
      end if

    end do ! iBnd

  end subroutine init_boundary_single
! ****************************************************************************** !

  ! ************************************************************************** !
  !> Check prerequisite for some boundary conditions
  subroutine check_BCs_preRequisite_MS(field, nBCs)
    ! --------------------------------------------------------------------------
    !> fluid parameters and properties
    type( mus_field_type ), intent(inout) :: field(:)
    !> Number of boundary types
    integer, intent(in) :: nBCs
    ! --------------------------------------------------------------------------
    integer :: iField, iBnd, counter(nBCs)
    logical :: checkBnd(nBCs)
    character(len=labelLen) :: checkKind(nBCs)
    ! --------------------------------------------------------------------------
    counter = 0
    checkBnd = .false.
    do iField = 1, size(field)
      do iBnd = 1, nBCs
        select case (trim(field(iField)%bc(iBnd)%BC_kind))
        case ('spc_inlet_eq', 'spc_inlet', 'spc_outlet_vel')
          ! this boundaries can be applied only if all species has this kind
          counter(iBnd) = counter(iBnd) + 1
          checkBnd(iBnd) = .true.
          checkKind(iBnd) = field(iField)%bc(iBnd)%BC_kind
        end select
      end do !iBnd
    end do !iField

    do iBnd = 1, nBCs
      if (checkBnd(iBnd)) then
        if (counter(iBnd) /= size(field)) then
           call tem_abort( 'Error: Not all species has kind: ' &
             &            //trim(checkKind(iBnd)) )
        end if
      end if
    end do
  end subroutine check_BCs_preRequisite_MS
  ! **************************************************************************** !

  ! **************************************************************************** !
  !> Initialize the values required for the moments BC
  subroutine init_momentsBC( bc, leveldesc, layout, globBC, minLevel, maxLevel )
    ! ---------------------------------------------------------------------------
    !> Level range
    integer, intent(in) :: minLevel, maxLevel
    !> global array boundary type
    type( boundary_type ), intent(inout) :: bc
    !> Level descriptor
    type( tem_leveldesc_type ), intent(in) :: levelDesc(minLevel:maxLevel)
    !> Layout
    type( mus_scheme_layout_type), intent(in) :: layout
    !> scheme global boundary type
    type( glob_boundary_type ), intent(inout) :: globBC
    ! ---------------------------------------------------------------------------
    integer :: iElem, iLevel, iDir
    integer :: d2q9_xNormal(3), d2q9_yNormal(3)
    integer :: d3q19_xNormal(5), d3q19_yNormal(5), d3q19_zNormal(5)
    integer, allocatable, dimension(:) :: xNormal_mom, yNormal_mom,            &
      &                                   zNormal_mom,                         &
      &                                   xyNormal_mom, yzNormal_mom,          &
      &                                   xzNormal_mom,                        &
      &                                   xyzNormal_mom
    integer, allocatable, dimension(:,:) :: xNorm_links, yNorm_links,          &
      &                                     zNorm_links, xyNorm_links,         &
      &                                     yzNorm_links, xzNorm_links,        &
      &                                     xyzNorm_links
    character(len=labelLen) :: normalIndex
    real(kind=rk) :: normal(3)
    integer :: nLinks, iLink, normal_nLinks, edge_nLinks, corner_nLinks
    integer, allocatable :: missing_links(:)
    real(kind=rk), allocatable :: unKnown_Mat(:,:)
    integer :: elemPos
    logical :: bitmask( layout%fStencil%QQN )
    integer(kind=long_k), allocatable :: corner_elems(:)
    integer(kind=long_k) :: treeID
    logical :: corner_node, update_allMoments
    integer :: iCorner
    ! ---------------------------------------------------------------------------
!    write(dbgUnit(1),*) 'Boundary label ', trim(bc%label)

    ! @todo KM: move this generic info to tem_stencil_module or tem_param_module
    ! known moments positions in moments array
    select case (trim(bc%BC_kind))
    case('pressure_momentsbased', 'moments_outflow', 'spc_moments_molefrac')
      d2q9_xNormal = (/ 1, 5, 3/) !m0, mY, mYY
      d2q9_yNormal = (/ 1, 4, 2/) !m0, mX, mXX
      d3q19_xNormal = (/ 1, 3, 4, 6, 7/) !m0, mY, mZ, mYY, mZZ
      d3q19_yNormal = (/ 1, 2, 4, 5, 7/) !m0, mX, mZ, mXX, mZZ
      d3q19_zNormal = (/ 1, 2, 3, 5, 6/) !m0, mX, mY, mXX, mYY
    case('moments_wall','velocity_momentsbased','moments_inflow', &
      &  'spc_moments_moleflux','spc_moments_wall','spc_moments_vel', &
      & 'spc_moments_outflow' )
      ! ordering is based on such that 1st knownMoments pos is
      ! moment in normal direction and rest follows
      d2q9_xNormal = (/ 2, 3, 5/) !mX, mY, mYY
      d2q9_yNormal = (/ 3, 2, 4/) !mX, mY, mXX
      d3q19_xNormal = (/ 2, 3, 4, 6, 7/) !mX, mY, mZ, mYY, mZZ
      d3q19_yNormal = (/ 3, 2, 4, 5, 7/) !mX, mY, mZ, mXX, mZZ
      d3q19_zNormal = (/ 4, 2, 3, 5, 6/) !mX, mY, mZ, mXX, mYY
    case default
      d2q9_xNormal = 0
      d2q9_yNormal = 0
      d3q19_xNormal = 0
      d3q19_yNormal = 0
      d3q19_zNormal = 0
      write(logUnit(1),*)'This boundary kind is not supported'
      call tem_abort()
    end select

    select case (trim(layout%fStencil%label))
    case('d2q9')
      normal_nLinks = 3
      ! in 2d edge and corner or same
      edge_nLinks = 5
      corner_nLinks = 5
      !axis-normal
      allocate(xNormal_mom(normal_nLinks))
      allocate(yNormal_mom(normal_nLinks))
      !edge normal
      allocate(xyNormal_mom(edge_nLinks))
      xNormal_mom = d2q9_xNormal
      yNormal_mom = d2q9_yNormal
      xyNormal_mom = (/ 1, 2, 3, 4, 5/) !m0, mX, mY, mXX, mYY

      ! normal links
      allocate(xNorm_links(normal_nLinks,2))
      xNorm_links = reshape((/ 1, 5, 6,    & !x-
        &                      3, 7, 8 /), & !x+
        &                  (/normal_nLinks,2/))

      allocate(yNorm_links(normal_nLinks,2))
      yNorm_links = reshape((/ 2, 5, 7,    & !y-
        &                      4, 6, 8 /), & !y+
        &                  (/normal_nLinks,2/))

      allocate(xyNorm_links(edge_nLinks,4))
      xyNorm_links = reshape((/ 1, 2, 5, 6, 7,   & !x-,y-
        &                       1, 4, 5, 6, 8,   & !x-.y+
        &                       2, 3, 5, 7, 8,   & !x+,y-
        &                       3, 4, 6, 7, 8 /),& !x+,y+
        &                   (/edge_nLinks,4/))
      allocate(zNorm_links(normal_nLinks,2))
      zNorm_links = 0
      allocate(yzNorm_links(edge_nLinks,4))
      yzNorm_links = 0
      allocate(xzNorm_links(edge_nLinks,4))
      xzNorm_links = 0
      allocate(xyzNorm_links(corner_nLinks,8))
      xyzNorm_links = 0
    case('d3q19')
      normal_nLinks = 5
      edge_nLinks = 9
      corner_nLinks = 12
      !axis-normal
      allocate(xNormal_mom(normal_nLinks))
      allocate(yNormal_mom(normal_nLinks))
      allocate(zNormal_mom(normal_nLinks))
      !edge normal
      allocate(xyNormal_mom(edge_nLinks))
      allocate(yzNormal_mom(edge_nLinks))
      allocate(xzNormal_mom(edge_nLinks))
      !corner normal
      allocate(xyzNormal_mom(corner_nLinks))
      xNormal_mom = d3q19_xNormal
      yNormal_mom = d3q19_yNormal
      zNormal_mom = d3q19_zNormal
      ! m0, mX, mY, mZ, mXX, mYY, mZZ, mYZ, mZZX
      xyNormal_mom = (/ 1, 2, 3, 4, 5, 6, 7, 9, 15/)
      ! m0, mX, mY, mZ, mXX, mYY, mZZ, mXZ, mXXY
      yzNormal_mom = (/ 1, 2, 3, 4, 5, 6, 7, 10, 11/)
      ! m0, mX, mY, mZ, mXX, mYY, mZZ, mXY, mYYZ
      xzNormal_mom = (/ 1, 2, 3, 4, 5, 6, 7, 8, 14/)
      ! m0, mX, mY, mZ, mXX, mYY, mZZ, mXY, mYZ, mXXY, mYYX, mZZY
      xyzNormal_mom = (/ 1, 2, 3, 4, 5, 6, 7, 8, 9, 11, 13, 16/)

      ! normal links
      allocate(xNorm_links(normal_nLinks,2))
      xNorm_links = reshape((/ 1, 11, 13, 15, 16,    & !x-
        &                      4, 12, 14, 17, 18 /), & !x+
        &                  (/normal_nLinks,2/))

      allocate(yNorm_links(normal_nLinks,2))
      yNorm_links = reshape((/ 2,  7,  8, 15, 17,    & !y-
        &                      5,  9, 10, 16, 18 /), & !y+
        &                  (/normal_nLinks,2/) )

      allocate(zNorm_links(normal_nLinks,2))
      zNorm_links = reshape((/ 3,  7,  9, 11, 12,    & !z-
        &                      6,  8, 10, 13, 14 /), & !z+
        &                  (/normal_nLinks,2/))

      allocate(xyNorm_links(edge_nLinks,4))
      xyNorm_links = reshape((/ 1,  2,  7,  8, 11, 13, 15, 16, 17,   & !x-,y-
        &                       1,  5,  9, 10, 11, 13, 15, 16, 18,   & !x-,y+
        &                       2,  4,  7,  8, 12, 14, 15, 17, 18,   & !x+,y-
        &                       4,  5,  9, 10, 12, 14, 16, 17, 18 /),& !x+,y+
        &                   (/edge_nLinks,4/))

      allocate(yzNorm_links(edge_nLinks,4))
      yzNorm_links = reshape((/ 2,  3,  7,  8,  9, 11, 12, 15, 17,   & !y-,z-
        &                       2,  6,  7,  8, 10, 13, 14, 15, 17,   & !y-,z+
        &                       3,  5,  7,  9, 10, 11, 12, 16, 18,   & !y+,z-
        &                       5,  6,  8,  9, 10, 13, 14, 16, 18 /),& !y+,z+
        &                   (/edge_nLinks,4/))

      allocate(xzNorm_links(edge_nLinks,4))
      xzNorm_links = reshape((/ 1,  3,  7,  9, 11, 12, 13, 15, 16,   & !x-,z-
        &                       1,  6,  8, 10, 11, 13, 14, 15, 16,   & !x-,z+
        &                       3,  4,  7,  9, 11, 12, 14, 17, 18,   & !x+,z-
        &                       4,  6,  8, 10, 12, 13, 14, 17, 18 /),& !x+,z+
        &                   (/edge_nLinks,4/))

      allocate(xyzNorm_links(corner_nLinks,8))
      xyzNorm_links = reshape((/ &
        &        1,  2,  3,  7,  8,  9, 11, 12, 13, 15, 16, 17,   & !x-,y-,z-
        &        1,  2,  6,  7,  8, 10, 11, 13, 14, 15, 16, 17,   & !x-,y-,z+
        &        1,  3,  5,  7,  9, 10, 11, 12, 13, 15, 16, 18,   & !x-,y+,z-
        &        1,  5,  6,  8,  9, 10, 11, 13, 14, 15, 16, 18,   & !x-,y+,z+
        &        2,  3,  4,  7,  8,  9, 11, 12, 14, 15, 17, 18,   & !x+,y-,z-
        &        2,  4,  6,  7,  8, 10, 12, 13, 14, 15, 17, 18,   & !x+,y-,z+
        &        3,  4,  5,  7,  9, 10, 11, 12, 14, 16, 17, 18,   & !x+,y+,z-
        &        4,  5,  6,  8,  9, 10, 12, 13, 14, 16, 17, 18 /),& !x+,y+,z+
        &                    (/corner_nLinks,8/))
    case default
      normal_nLinks = 0
      edge_nLinks   = 0
      corner_nLinks = 0
      allocate(xNormal_mom(0))
      allocate(yNormal_mom(0))
      allocate(zNormal_mom(0))
      allocate(xyNormal_mom(0))
      allocate(yzNormal_mom(0))
      allocate(xzNormal_mom(0))
      write(logUnit(1),*) trim(layout%fStencil%label)//&
        &' layout is not supported for this boundary'
      call tem_abort()
    end select


    do iLevel = minLevel, maxLevel
      allocate( bc%elemLvl( iLevel )%moments( globBC%nElems( iLevel )))

      ! list of corner elements treeIDs
      allocate(corner_elems(globBC%cornerBC%nElems(iLevel)))
      corner_elems =                                                           &
       & levelDesc(iLevel)%total(globBC%cornerBC%elemLvl(iLevel)%elem%val(:))

      iCorner = 0
      do iElem = 1, globBC%nElems( iLevel )
        ! when moments combination are not defined properly
        ! update all moments
        update_allMoments = .false.
        ! if corrent node is corner node i.e intersected by multiple boundaries
        corner_node = .false.

!write(dbgUnit(1),*) 'iElem ', iElem
        treeID = levelDesc(iLevel)%total(globBC%elemLvl(iLevel)%elem%val(iElem))
        if ( any(corner_elems == treeID) ) corner_node = .true.
!write(dbgUnit(1),*) 'corner_node ', corner_node

        ! update bitmask, normal and normalInd with cornerBC information
        ! if corrent node is corner node
        if (corner_node) then
          iCorner = iCorner + 1
          globBC%elemLvl(iLevel)%bitmask%val(:, iElem) = &
            & globBC%cornerBC%elemLvl(iLevel)%bitmask%val(:, iCorner)
          globBC%elemLvl(iLevel)%normal%val(:, iElem) =  &
            & globBC%cornerBC%elemLvl(iLevel)%normal%val(:, iCorner)
          globBC%elemLvl(iLevel)%normalInd%val(iElem) =  &
            & globBC%cornerBC%elemLvl(iLevel)%normalInd%val(iCorner)
        end if ! corner node

        ! element position in state array and treeID list
        elemPos = globBC%elemLvl(iLevel)%elem%val(iElem)
        ! normal of this corrent node
        normal = globBC%elemLvl(iLevel)%normal%val(:,iElem)
        bitmask = globBC%elemLvl(iLevel)%bitmask%val(:, iElem)

        nLinks = count(globBC%elemLvl(iLevel)%bitmask%val(:, iElem))
        bc%elemLvl(iLevel)%moments(iElem)%nUnKnownPdfs = nLinks

        allocate(bc%elemLvl(iLevel)%moments(iElem)%knownMom_pos(nLinks))
        allocate(missing_links(nLinks))
        iLink = 0
        do iDir = 1,layout%fStencil%QQN
          if (globBC%elemLvl(iLevel)%bitmask%val(iDir, iElem)) then
            iLink = iLink+1
            missing_links(iLink) = iDir
          endif
        end do

        normalIndex = 'default'

        ! axis normals
        if (nLinks == normal_nLinks) then
          do iDir = 1, 2
            if ( all(missing_links == xNorm_links(:, iDir)) ) &
              & normalIndex = 'xNormal'
            if ( all(missing_links == yNorm_links(:, iDir)) ) &
              & normalIndex = 'yNormal'
            if ( all(missing_links == zNorm_links(:, iDir)) ) &
              & normalIndex = 'zNormal'
          end do
        else if (nLinks == edge_nLinks) then
        ! edge direction
          do iDir = 1, 4
            if ( all(missing_links == xyNorm_links(:, iDir)) ) &
              & normalIndex = 'xyNormal'

            if ( all(missing_links == yzNorm_links(:, iDir)) ) &
              & normalIndex = 'yzNormal'

            if ( all(missing_links == xzNorm_links(:, iDir)) ) &
              & normalIndex = 'xzNormal'
          end do
        else if (nLinks == corner_nLinks) then
        ! corner direction
          do iDir = 1, 8
            if ( all(missing_links == xyzNorm_links(:,iDir)) ) &
              & normalIndex = 'xyzNormal'
          end do
        end if

!write(dbgUnit(1),*) 'normal ', globBC%elemLvl(iLevel)%normal%val(:, iElem)
!write(dbgUnit(1),*) 'normalIndex ', normalIndex
!write(dbgUnit(1),*) 'missing_links ', missing_links
        select case(trim(normalIndex))
        case ('xNormal')
          bc%elemLvl(iLevel)%moments(iElem)%knownMom_pos = xNormal_mom
        case ('yNormal')
          bc%elemLvl(iLevel)%moments(iElem)%knownMom_pos = yNormal_mom
        case ('zNormal')
          bc%elemLvl(iLevel)%moments(iElem)%knownMom_pos = zNormal_mom
        case ('xyNormal')
          bc%elemLvl(iLevel)%moments(iElem)%knownMom_pos = xyNormal_mom
        case ('yzNormal')
          bc%elemLvl(iLevel)%moments(iElem)%knownMom_pos = yzNormal_mom
        case ('xzNormal')
          bc%elemLvl(iLevel)%moments(iElem)%knownMom_pos = xzNormal_mom
        case ('xyzNormal')
          bc%elemLvl(iLevel)%moments(iElem)%knownMom_pos = xyzNormal_mom
        case default
          ! undefined normal so update all links for this boundary
          globBC%elemLvl(iLevel)%bitmask%val(:, iElem) = .true.
          nLinks = layout%fStencil%QQ
          deallocate(bc%elemLvl(iLevel)%moments(iElem)%knownMom_pos)
          allocate(bc%elemLvl(iLevel)%moments(iElem)%knownMom_pos(nLinks))
          do iLink = 1,nLinks
            bc%elemLvl(iLevel)%moments(iElem)%knownMom_pos(iLink) = iLink
          enddo

          deallocate(missing_links)
          allocate(missing_links(nLinks))
          iLink = 0
          do iDir = 1,layout%fStencil%QQN
            if (globBC%elemLvl(iLevel)%bitmask%val(iDir, iElem)) then
              iLink = iLink+1
              missing_links(iLink) = iDir
            endif
          end do
        end select

        bc%elemLvl(iLevel)%moments(iElem)%nUnKnownPdfs = nLinks
        allocate(bc%elemLvl(iLevel)%moments(iElem)                   &
          &                      %unKnownPdfs_MatInv(nLinks,nLinks))
        allocate(unKnown_Mat(nLinks,nLinks))

        do iLink = 1, nLinks
          unKnown_Mat(:,iLink) = layout%moment%toMoments%A( &
            & bc%elemLvl(iLevel)%moments(iElem)%knownMom_pos, &
            & missing_links(iLink) )
        end do

        bc%elemLvl(iLevel)%moments(iElem)%unKnownPdfs_MatInv &
          & = invert_matrix(unKnown_Mat)

        deallocate(unKnown_Mat)
        deallocate(missing_links)

      end do ! iElem
    end do ! level

  end subroutine init_momentsBC
  ! **************************************************************************** !


! ****************************************************************************** !
  !> assign qVal to corresponding BC and level-wise if qVal not provided by
  !! seeder. qVal from seeder are assigned in assignBCList in
  !! mus_construction_module, So set qVal from config only when it is not
  !! provided by seeder.
  subroutine init_qVals(refQVal, layout, globBC, minLevel, maxLevel)
    ! ---------------------------------------------------------------------------
    real(kind=rk),                 intent(in) :: refQVal
    type( mus_scheme_layout_type), intent(in) :: layout
    type( glob_boundary_type ), intent(inout) :: globBC
    integer,                       intent(in) :: minLevel, maxLevel
    ! ---------------------------------------------------------------------------
    integer :: iElem, iLevel, iDir, invDir, QQN
    ! ---------------------------------------------------------------------------

    globBC%qValInitialized = .true.
    QQN = layout%fStencil%QQN

    ! 2. allocate qVal array level-wise
    lvlLoop: do iLevel = minLevel, maxLevel

      ! if qVal is not available from seeder then initialize qVal
      call init( me     = globBC%elemLvl(iLevel)%qVal, &
        &        width  = QQN,                         &
        &        length = globBC%nElems(iLevel)        )
        globBC%elemLvl( iLevel )%qVal%val = 0.0_rk

      ! loop over elements
      do iElem = 1, globBC%nElems(iLevel)
        ! loop over directions
        do iDir = 1, QQN
          ! Bitmask is true for incoming direction.
          ! so use invDir to access qVal and cxDirRK
          if ( globBC%elemLvl(iLevel)%bitmask%val(iDir, iElem) ) then
            invDir  = layout%fStencil%cxDirInv(iDir)
            ! if no qVal from seeder then use refQVal
            globBC%elemLvl(iLevel)%qVal%val(invDir, iElem) = refQVal
          end if
        end do !iDir
      end do !iElem
    end do lvlLoop

  end subroutine init_qVals
! ****************************************************************************** !


! ****************************************************************************** !
  !> Initialize the values required for the characteristic boundary conditions
  !!
  !! Simply calculate the macroscopic values from the current pdf distributions
  !! in the boundary element and its neighbor elements and store it to the
  !! %nrbc%lodi field
  !! NOTE:
  !! To be consistent, the same calculation procedure of the LODI values need to
  !! be performed as in the
  !! [[mus_bc_fluid_module:outlet_nrbc]] "NRBC routine" itself
  !!
  subroutine init_nrbc( bc, state, nSize, neigh, layout, level, &
    &                   nScalars, varSys, derVarPos, elemPos, nElems )
    ! ---------------------------------------------------------------------------
    !> global array boundary type
    type( boundary_type ) :: bc
    !> level
    integer, intent(in) :: level
    !> State array
    real(kind=rk), intent(in) :: state(:)
    !> nSize
    integer, intent(in) :: nSize
    !> neighbor array
    integer, intent(in) :: neigh(:)
    !> fluid parameters
    type( mus_scheme_layout_type), intent(in) :: layout
    !> number of scalars in global 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
    !> BC elements positions in total list
    integer, intent(in) :: elemPos(:)
    !> number of BC elements
    integer, intent(in) :: nElems
    ! ---------------------------------------------------------------------------
    integer :: iElem, iDir, iField, QQ
    integer :: iNeigh, neighPos, posInTotal
    real(kind=rk) :: ff( layout%fStencil%QQ )
    ! ---------------------------------------------------------------------------

    QQ = layout%fStencil%QQ
    iField = 1

    ! Assign the initial values to the lodi variables
    allocate( bc%elemLvl( level )%lodi( 4, 3, nElems))
    bc%elemLvl( level )%lodi = 0._rk

    ! reset
    do iElem = 1, nElems
      posInTotal = elemPos( iElem )

      ! adopt initial density from the fluid cells
      do iDir = 1, QQ
        ff(iDir) = state( (posintotal-1)*nscalars+idir+(ifield-1)*qq )
      end do !iDir
      bc%elemLvl( level )%lodi( 1, 1, iElem ) = sum( ff )

      ! adapt initial velocities from the fluid cells
      call derVarPos%velFromState(                                &
        &           state  = ff,                                  &
        &           iField = iField,                              &
        &           nElems = 1,                                   &
        &           varSys = varSys,                              &
        &           layout = layout,                              &
        &           res    = bc%elemLvl(level)%lodi(2:4, 1, iElem) )

      do iNeigh = 1,2
        ! Get the values for the neighbors as well
        neighPos = bc%neigh( level )%posInState( iNeigh, iElem )
        if( neighPos > 0 ) then
          do iDir = 1,QQ
            ff(iDir) = state( (neighpos-1)*nscalars+idir+(ifield-1)*qq)
          end do
          bc%elemLvl( level )%lodi( 1, 1+iNeigh, iElem ) = sum( ff )

          ! velocities
          call derVarPos%velFromState(                                  &
            &       state  = ff,                                        &
            &       iField = iField,                                    &
            &       nElems = 1,                                         &
            &       varSys = varSys,                                    &
            &       layout = layout,                                    &
            &       res    = bc%elemLvl(level)%lodi(2:4, 1+iNeigh, iElem) )
        end if
      end do ! iNeigh
    end do ! iElem

  end subroutine init_nrbc
! ****************************************************************************** !


  ! -------------------------------------------------------------------------- !
  !> This routine allocates turbulent viscosity and friction velocity on
  !! boundary elements. It also initialize friction velocity from stream-wise
  !! velocity component on first neighbor in normal direction.
  subroutine mus_init_turb_wallFunc(bc, globBC, auxField, viscKine, derVarPos,&
    &                                varSys, stencil, iLevel)
    ! -------------------------------------------------------------------- !
    !> Field bc which contains turbwallFunc type and neighbor info
    type(boundary_type), intent(inout) :: bc
    !> global bc of current boundary with elemPos and normal info
    type(glob_boundary_type), intent(inout) :: globBC
    !> auxField array
    real(kind=rk), intent(in) :: auxField(:)
    !> Kinematic viscosity
    type(mus_viscosity_type) :: viscKine
    !> position of derived quantities in varsys
    type(mus_derVarPos_type), intent(in) :: derVarPos
    !> scheme variable system
    type( tem_varSys_type ), intent(in) :: varSys
    !> stencil info
    type(tem_stencilHeader_type), intent(in) :: stencil
    !> current level
    integer, intent(in) :: iLevel
    ! -------------------------------------------------------------------- !
    integer :: iElem, neighPos, elemOff, elemPos, normalInd_inv, vel_pos(3)
    integer :: nElems
    real(kind=rk) :: velNeigh(3), normal(3), streamwise(3), unitSW(3)
    real(kind=rk) :: velSW, normalMag, qVal, distToBnd, streamwise_mag
    ! -------------------------------------------------------------------- !
    nElems = globBC%nElems(iLevel)
    allocate( bc%turbwallFunc%dataOnLvl(iLevel)%tVisc(nElems) )
    bc%turbwallFunc%dataOnLvl(iLevel)%tVisc = 0.0_rk

    allocate( bc%turbwallFunc%dataOnLvl(iLevel)%velTau(nElems) )
    bc%turbwallFunc%dataOnLvl(iLevel)%velTau = 0.0_rk

    allocate( bc%turbwallFunc%dataOnLvl(iLevel)%distToBnd(nElems) )
    bc%turbwallFunc%dataOnLvl(iLevel)%distToBnd = 0.0_rk

    allocate( bc%turbwallFunc%dataOnLvl(iLevel)%neighDistToBnd(nElems) )
    bc%turbwallFunc%dataOnLvl(iLevel)%neighDistToBnd = 0.0_rk

    ! Position of velocity in auxField array
    vel_pos = varSys%method%val(derVarPos%velocity)%auxField_varPos(1:3)

    do iElem = 1, nElems
      ! Get velocity on the first neighbor element from auxField array
      neighPos = bc%neigh(iLevel)%posInState(1, iElem)
      elemOff = (neighPos-1)*varSys%nAuxScalars
      velNeigh(1) = auxField(elemOff + vel_pos(1))
      velNeigh(2) = auxField(elemOff + vel_pos(2))
      velNeigh(3) = auxField(elemOff + vel_pos(3))

      ! Calculate local stream-wise unit vector
      ! e' = (u - (u . n) . n) / |(u - (u . n) . n|
      !
      ! NormalInd in bc elem is pointing into the domain so we invert it
      normalInd_inv = stencil%cxDirInv( globBC%elemLvl(iLevel)  &
        &                                 %normalInd%val(iElem) )
      normal = stencil%cxDirRK(:, normalInd_inv)
      streamwise = velNeigh - dot_product( velNeigh, normal ) * normal
      streamwise_mag = sqrt(dot_product(streamwise, streamwise))
      ! if velocity is zero then compute streamwise unit vector from bc normal
      if (streamwise_mag .feq. 0.0_rk)then
        call tem_abort("Error: In Init_turb_wallFunc. Stream-wise velocity &
          &             mag is zero. Initialite flow with non-zero velocity")
      end if
      ! Unit vector
      unitSW = streamwise / streamwise_mag

      ! stream-wise velocity component
      velSW = dot_product(velNeigh, unitSW)

      ! use qValue in normal direction as distance of firode from boundary
      normalMag = sqrt(dot_product(normal, normal))
      qVal = globBC%elemLvl(iLevel)%qVal%val(normalInd_inv, iElem)
      distToBnd = qVal + normalMag
      bc%turbwallFunc%dataOnLvl(iLevel)%distToBnd(iElem) = qVal
      bc%turbwallFunc%dataOnLvl(iLevel)%neighDistToBnd(iElem) = distToBnd

      ! Initialize friction velocity on first neighbor with werner and wengle
      ! wall model in linear profile i.e. u+=y+ -> u_tau = sqrt(u * nu / y)
      elemPos = globBC%elemLvl(iLevel)%elem%val(iElem)
      bc%turbwallFunc%dataOnLvl(iLevel)%velTau(iElem)              &
        & = sqrt( velSW * viscKine%dataOnLvl(iLevel)%val( elemPos ) &
        &        / distToBnd )
    end do
  end subroutine mus_init_turb_wallFunc
  ! -------------------------------------------------------------------------- !


! ****************************************************************************** !
  !> Transfer pre- and post-collision PDF of neighbors of boundary elements
  !! into neighBufferPre and neighBufferPost.
  !! Access to state array
  !! ---------------------
  !!
  !! * NeighBufferPre: FETCH
  !! * NeighBufferPost: SAVE
  !! @todo: for PUSH, post-collision is not correct yet.
  !!
  subroutine fill_neighBuffer( prevstate, currstate, neigh, globBC, nBCs, &
    &                          field, varSys, QQ, nSize, iLevel )
    !---------------------------------------------------------------------------
    !> Previous state vector of iLevel
    real(kind=rk),intent(in) :: prevstate(:)
    !> Current state vector of iLevel
    real(kind=rk),intent(in) :: currstate(:)
    !> connectivity array corresponding to state vector
    integer,      intent(in) :: neigh(:)
    !> scheme global boundary type
    type( glob_boundary_type ), intent(in) :: globBC(:)
    !> number of total BC
    integer, intent(in) :: nBCs
    !> fluid parameters and properties
    type( mus_field_type ), intent(inout) :: field(:)
    !> scheme variable system
    type( tem_varSys_type ), intent(in) :: varSys
    !> number of total links
    integer, intent(in) :: QQ
    !> number of elements in state vector
    integer, intent(in) :: nSize
    !> the iLevel on which this boundary was invoked
    integer, intent(in) :: iLevel
    ! ---------------------------------------------------------------------------
    integer :: iField, iElem, iDir, iBnd, iNeigh
    integer :: pdfVarPos(QQ), neighPos, myPos
    ! ---------------------------------------------------------------------------
    ! write(dbgUnit(5),*) 'Fill neighBufferPre and neighBufferPost'

    do iField = 1, size( field )

      pdfVarPos(:) = varSys%method%val(iField)%state_varPos(1:QQ)

      ! fill boundary neighbor pre- and post-collision state values
      do iBnd = 1, nBCs
        call tem_startTimer( timerHandle = mus_timerHandles%setBnd(iBnd) )
        if (field(iField)%bc( iBnd )%requireNeighBufPre) then
! write(dbgUnit(5),"(A,I0)") 'iBC: ', iBnd
          do iNeigh = 1, field(iField)%bc( iBnd )%nNeighs
! write(dbgUnit(5),"(A,I0)") 'iNeigh: ', iNeigh
            do iElem = 1, globBC( iBnd )%nElems( iLevel )
              ! neighbor position in total list
              neighPos =  field(iField)%bc( iBnd )%neigh( iLevel )          &
                &                                 %posInState( iNeigh, iElem )
! write(dbgUnit(5),"(2(A,I0))") 'iElem: ', iElem, ', neighPos: ', neighPos
              do iDir = 1, QQ
                ! pre-collision neigh
                field(iField)%bc( iBnd )%neigh( iLevel )%neighBufferPre( &
  & iNeigh, (iElem-1)*QQ+iDir ) = prevstate(                             &
  &  neigh (( idir-1)* nsize+ neighpos)+( ifield-1)* qq+ varsys%nscalars*0 )
! write(dbgUnit(5),"(A,I0)") 'iDir: ', iDir
! write(dbgUnit(5),"(3(A,I0),A)") &
!   & 'neighBufferPre(', iNeigh, ',',(iElem-1)*QQ+iDir, ') = state(', &
!   &  neigh (( idir-1)* nsize+ neighpos)+( ifield-1)* qq+ varsys%nscalars*0, &
!   & ')'
              end do ! iDir
            end do ! iElem
          end do ! iNeigh
        end if ! neighBufferPre

        if (field(iField)%bc( iBnd )%requireNeighBufPre_nNext) then
! write(dbgUnit(5),"(A,I0)") 'iBC: ', iBnd
          do iNeigh = 1, field(iField)%bc( iBnd )%nNeighs
! write(dbgUnit(5),"(A,I0)") 'iNeigh: ', iNeigh
            do iElem = 1, globBC( iBnd )%nElems( iLevel )
              ! neighbor position in total list
              neighPos =  field(iField)%bc( iBnd )%neigh( iLevel )          &
                &                                 %posInState( iNeigh, iElem )
! write(dbgUnit(5),"(2(A,I0))") 'iElem: ', iElem, ', neighPos: ', neighPos
              do iDir = 1, QQ
                ! pre-collision neigh
                field(iField)%bc( iBnd )%neigh( iLevel )%neighBufferPre_nNext( &
  & iNeigh, (iElem-1)*QQ+iDir ) = currstate(                                   &
  &  neigh (( idir-1)* nsize+ neighpos)+( ifield-1)* qq+ varsys%nscalars*0 )
! write(dbgUnit(5),"(A,I0)") 'iDir: ', iDir
! write(dbgUnit(5),"(3(A,I0),A)") &
!   & 'neighBufferPre_nNext(', iNeigh, ',',(iElem-1)*QQ+iDir, ') = state(', &
!   &  neigh (( idir-1)* nsize+ neighpos)+( ifield-1)* qq+ varsys%nscalars*0, &
!   & ')'
              end do ! iDir
            end do ! iElem
          end do ! iNeigh
        end if ! neighBufferPre_nNext

        if (field(iField)%bc( iBnd )%requireNeighBufPost) then
! write(dbgUnit(5),"(A,I0)") 'iBC: ', iBnd
          do iNeigh = 1, field(iField)%bc( iBnd )%nNeighs
! write(dbgUnit(5),"(A,I0)") 'iNeigh: ', iNeigh
            do iElem = 1, globBC( iBnd )%nElems( iLevel )
              ! neighbor position in total list
              neighPos =  field(iField)%bc( iBnd )%neigh( iLevel )          &
                &                                 %posInState( iNeigh, iElem )
! write(dbgUnit(5),"(2(A,I0))") 'iElem: ', iElem, ', neighPos: ', neighPos
              do iDir = 1, QQ
                ! for PULL, this is post-collision value
                ! for PUSH, this is  pre-collision value
                field(iField)%bc( iBnd )%neigh( iLevel )%neighBufferPost( &
  & iNeigh, (iElem-1)*QQ+iDir ) = currstate(          &
  & ( neighpos-1)* varsys%nscalars+ idir+( ifield-1)* qq )
!  & ( neighpos-1)* varsys%nscalars+ pdfvarpos(idir))

! write(dbgUnit(5),"(A,I0)") 'iDir: ', iDir
! write(dbgUnit(5),"(3(A,I0),A)") &
!   & 'neighBufferPost(', iNeigh, ',',(iElem-1)*QQ+iDir, ') = state(', &
!   & ( neighpos-1)* varsys%nscalars+ pdfvarpos(idir), &
!   & ')'
              end do ! iDir
            end do ! iElem
          end do ! iNeigh
        end if !neighBufferPost

        if ( field(iField)%bc( iBnd )%useComputeNeigh ) then
          do iElem = 1, globBC( iBnd )%nElems( iLevel )
            ! my (bnd element) position in total list
            myPos = globBC( iBnd )%elemLvl( iLevel )%elem%val( iElem )
            do iDir = 1, QQ
              ! get post-collision state values of my neighbors
              ! computeNeighBuf always uses AOS layout
field(iField)%bc(iBnd)%neigh(iLevel)%computeNeighBuf( (iElem-1)*QQ+iDir ) =    &
& currstate(  neigh (( idir-1)* nsize+ mypos)+( ifield-1)* qq+ varsys%nscalars*0 )
            end do ! iDir
          end do ! iElem
        end if ! useComputeNeigh

        call tem_stopTimer( timerHandle = mus_timerHandles%setBnd(iBnd) )
      end do ! iBnd
    end do ! iField

  end subroutine fill_neighBuffer
! ****************************************************************************** !


! ****************************************************************************** !
  !> Transfer pdf of boundary elements into bcBuffer which is used by all
  !! boundary routines.
  ! for PULL, this is post-collision value
  ! for PUSH, this is  pre-collision value
  subroutine fill_bcBuffer( bcBuffer, currState, neigh, nSize, nElems_bc, &
    &                       posInTotal, nFields, QQ, varSys )
    ! ---------------------------------------------------------------------------
    !> state values of all boundary elements
    real(kind=rk),intent(out) :: bcBuffer(:)
    !> Current state vector
    real(kind=rk),intent(in) :: currState(:)
    !> connectivity array corresponding to state vector
    integer,      intent(in) :: neigh(:)
    !> nSize
    integer, intent(in) :: nSize
    !> number of boundary elements
    integer, intent(in) :: nElems_bc
    !> positions in total list of boundary elements
    integer, intent(in) :: posInTotal(nElems_bc)
    !> Number of fields
    integer, intent(in) :: nFields
    !> number of total links
    integer, intent(in) :: QQ
    !> scheme variable system
    type( tem_varSys_type ), intent(in) :: varSys
    ! ---------------------------------------------------------------------------
    integer :: iElem, iField, varPos, iDir
    ! ---------------------------------------------------------------------------

    ! write(dbgUnit(10), "(A)") ' Fill bcBuffer. '

    ! Loop over fields
    do iField = 1, nFields
      !NEC$ ivdep
      !dir$ ivdep
      !ibm* independent
      do iElem = 1, nElems_bc
        do iDir = 1, QQ
          ! bcBuffer always uses AOS data structure
          varPos = varSys%method%val(iField)%state_varPos(iDir)
          bcBuffer( varPos+(iElem-1)*varSys%nScalars ) = currState( &
& ( posintotal(ielem)-1)* varsys%nscalars+ idir+( ifield-1)* qq )
        end do
      end do ! iElem
    end do

  end subroutine fill_bcBuffer
! ****************************************************************************** !

! ***************************************************************************** !
  !> Get Surface points on boundary elements.
  !! For boundary state variable which are evaluated linkwise, extract surface
  !! points for each link and for non-link based variables project barycenter
  !! on the boundary surface.
  !! Return real coordinates on boundary surface and offset bit which encodes
  !! direction.
  subroutine mus_get_points_fromBC( bc, globBC, tree, stencil, total, iLevel, &
    &                               nPoints, points, offset_bit )
    !---------------------------------------------------------------------------
    !> Field boundary type
    type(boundary_type), intent(in) :: bc
    !> for number of elements in boundary and position in buffer
    type(glob_boundary_type), intent(in)     :: globBC
    !> global treelm mesh
    type(treelmesh_type), intent(in) ::tree
    !> global pdf type
    ! type(tem_levelDesc_type), intent(in) :: levelDesc
    integer(kind=long_k), intent(in) :: total(:)
    !> for directions
    type(tem_stencilHeader_type), intent(in) :: stencil
    !> Current level
    integer, intent(in) :: iLevel
    !> Number of points
    integer, intent(out) :: nPoints
    !> 3-d real coordinates on which boundary variables are evaluated
    real(kind=rk), allocatable, intent(out) :: points(:,:)
    !> Offset bit encodes direction of boundary.
    !! used by apesmate to translate space coordinate in the offset direction
    !! to determine the treeID in remote domain
    character, allocatable, intent(out) :: offset_bit(:)
    !---------------------------------------------------------------------------
    real(kind=rk) :: dx, bary(3)
    real(kind=rk), allocatable :: qVal(:)
    integer :: iElem, iDir, iPnt, invDir
    logical :: qValIsZero
    !---------------------------------------------------------------------------
    write(logUnit(10),*) ' Get points from BC '
    ! Element size on current level
    dx = tem_ElemSizeLevel(tree, iLevel)
    allocate(qVal(stencil%QQN))

    ! if qVal is not provided by seeder and default is zero then no need to
    ! extract point on boundary surface. Just extract on barycenter
    if ((bc%qVal .feq. 0.0_rk) .and. .not. globBC%hasQVal &
      & .and. .not. bc%curved) then
      qValIsZero = .true.
    else
      qValIsZero = .false.
    end if

    ! Get qValue if boundary variables are evaulated link wise barycenter
    ! projection on boundary surface only if qVal is not zero
    if (bc%evalBcVar_link .and. .not. qValIsZero) then
      nPoints = bc%links(iLevel)%nVals
    else
      ! boundary variables are evaluated at bary center projection on
      ! boundary surface along boundary normal
      nPoints = globBC%nElems(iLevel)
    end if

    ! if no points on this level for this boundary, leave this routine
    ! if (nPoints==0) return

    allocate(points(nPoints,3))
    allocate(offset_bit(nPoints))

    iPnt = 0
    ! loop over elements
    do iElem = 1, globBC%nElems(ilevel)
      bary(:) = tem_BaryOfId( tree, total(globbc%elemLvl(iLevel)  &
        &                                       %elem%val(iElem)) )

      ! Set qValue to 0.5 if qValue is not defined.
      ! Mostly qVal is already set for link based boundary conditions
      if( .not. globBC%hasQVal .and.        &
        & .not. globBC%qValInitialized ) then
        qVal = bc%qVal
      else
        qVal = globBC%elemLvl(iLevel)%qVal%val(:, iElem)
      end if

      ! boundary variables are evaluated at link wise projection of barycenter
      ! on boundary surface
      if (bc%evalBcVar_link .and. .not. qValIsZero) then

        ! loop over directions
        do iDir = 1, stencil%QQN
          ! Bitmask is true for incoming direction.
          ! so use invDir to access qVal and cxDirRK
          if ( globBC%elemLvl(iLevel)%bitmask%val(iDir, iElem) ) then
            iPnt   = iPnt + 1
            invDir  = stencil%cxDirInv(iDir)
            ! position of boundary surface
            points(iPnt, :) = bary(:) + dx * qVal(invDir)  &
              &             * stencil%cxDirRK(:, invDir)

            offset_bit(iPnt) = qOffset_inChar( stencil%map( invDir ) )
          end if
        end do !iDir
      elseif (bc%BC_kind == 'bc_pdf') then
        ! For this boundary, boundary elements must overlap between two domains
        ! so return barycenter and offset 0
        iPnt = iPnt + 1
        points(iPnt, :) = bary(:)
        offset_bit(iPnt) = qOffset_inChar(q000)
      else

        ! boundary variables are evaluated at bary center projection on
        ! boundary surface along boundary normal
        iPnt   = iPnt + 1
        invDir  = stencil%cxDirInv(globBC%elemLvl(iLevel)%normalInd%val(iElem))
        if (invDir<=stencil%QQN) then
          points(iPnt, :) = bary(:) + dx * qVal(invDir)   &
            &             * stencil%cxDirRK(:, invDir)
          offset_bit(iPnt) = qOffset_inChar( stencil%map( invDir ) )
        else
          points(iPnt, :) = bary(:)
          offset_bit(iPnt) = qOffset_inChar(q000)
        end if

      end if !bc%evalBcVar_link
    end do !iElem

  end subroutine mus_get_points_fromBC
! ***************************************************************************** !


! ***************************************************************************** !
  !> This routine setup indices for boundary variables in bc_State_type
  !! pntIndex for the points on which boundaries are treated.
  subroutine mus_setupIndices_forBC( bc, globBC, tree, stencil, levelDesc, &
    &                                varSys, minLevel, maxLevel)
    ! --------------------------------------------------------------------------
    !> Field boundary type
    type(boundary_type), target, intent(inout) :: bc
    !> for number of elements in boundary and position in buffer
    type(glob_boundary_type), intent(in)     :: globBC
    !> global treelm mesh
    type(treelmesh_type), intent(in) ::tree
    !> Min and Max level in mesh
    integer, intent(in) :: minLevel, maxLevel
    !> global pdf type
    type(tem_levelDesc_type), intent(in) :: levelDesc(minLevel:maxLevel)
    !> for directions
    type(tem_stencilHeader_type), intent(in) :: stencil
    !> Global variable system
    type(tem_varSys_type), intent(in) :: varSys
    ! --------------------------------------------------------------------------
    !> Number of points
    integer :: nPoints
    !> 3-d real coordinates on which boundary variables are evaluated
    real(kind=rk), allocatable :: points(:,:)
    !> Offset bit encodes direction of boundary.
    !! used by apesmate to translate space coordinate in the offset direction
    !! to determine the treeID in remote domain
    character, allocatable :: offset_bit(:)
    integer, allocatable :: idx(:)
    integer :: iLevel, iVar
    character(len=labelLen) :: bc_varName
    type(tem_bc_state_type), pointer :: bc_state => NULL()
    character(len=pathLen) :: isSurface
    ! --------------------------------------------------------------------------


    write(logUnit(10),*) ' Setup indices for BC: ', trim(bc%label)
    do iLevel = minLevel, maxLevel
      write(logUnit(10),*) 'iLevel: ', iLevel
      call mus_get_points_fromBC( bc         = bc,                      &
        &                         globBC     = globBC,                  &
        &                         tree       = tree,                    &
        &                         stencil    = stencil,                 &
        &                         total      = levelDesc(iLevel)%total, &
        &                         iLevel     = iLevel,                  &
        &                         nPoints    = nPoints,                 &
        &                         points     = points,                  &
        &                         offset_bit = offset_bit               )

      ! if no points on this level for this boundary, goto next level
      ! KM: Skipping append idx results in not allocated indexLvl%val array
      ! which causes problems in calling get_valOfIndex for intel compiler
      !if (nPoints==0) cycle

      allocate(idx(nPoints))
      ! store points in spacetime function and store indices of points
      ! in bc_state
      do iVar = 1, bc%varDict%nVals
        bc_varName = trim(bc%varDict%val(iVar)%key)
        write(logUnit(10),*) ' Storing index for variable: ', trim(bc_varName)
        isSurface = 'isSurface = true'

        bc_state => NULL()
        select case (trim(bc_varName))
        case ('velocity')
          bc_state => bc%bc_states%velocity
        case ('pdf')
          bc_state => bc%bc_states%pdf
          ! For pdf boundary, both domain boundary layer must overlap
          isSurface = 'isSurface = false'
        case ('pressure')
          bc_state => bc%bc_states%pressure
        case ('mass_flowrate')
          bc_state => bc%bc_states%massFlowRate
        case ('mole_fraction')
          bc_state => bc%bc_states%moleFrac
        case ('mole_density')
          bc_state => bc%bc_states%moleDens
        case ('mole_flux')
          bc_state => bc%bc_states%moleFlux
        case ('mole_diff_flux')
          bc_state => bc%bc_states%moleDiff_flux
        case ('potential')
          bc_state => bc%bc_states%potential
        case ('surface_charge_density')
          bc_state => bc%bc_states%surChargeDens
        case default
          write(logUnit(1),*) 'Error: Unknown boundary variable: "'// &
            & trim(bc_varName)//'"'
          call tem_abort()
        end select

        if (.not. associated(bc_state)) then
          call tem_abort('Error: bc_state is not assosiated')
        end if

        ! set params
        write(logUnit(10),*) ' Set params in bc variable: '          &
          &                //trim(varSys%varName%val(bc_state%varPos))
        call varSys%method%val(bc_state%varPos)%set_params( &
          & varSys   = varSys,                              &
          & instring = isSurface                            )

        call varSys%method%val(bc_state%varPos)%setup_indices( &
          & varSys     = varSys,                               &
          & point      = points,                               &
          & offset_bit = offset_bit,                           &
          & iLevel     = iLevel,                               &
          & tree       = tree,                                 &
          & nPnts      = nPoints,                              &
          & idx        = idx                                   )

        if (nPoints == 0) then
          ! KM: Intel compiler fails if indexLvl is not allocated and accessed
          ! in get_valOfIndex
          ! initialize array with size zero
          call init(bc_state%pntIndex%indexLvl(iLevel))
        else
          call append(bc_state%pntIndex%indexLvl(iLevel), idx)
        end if
      end do !iVar
      deallocate(idx)
    end do !iLevel

  end subroutine mus_setupIndices_forBC
! ***************************************************************************** !


end module mus_bc_general_module
! ****************************************************************************** !
