! Copyright (c) 2012-2013 Manuel Hasert <m.hasert@grs-sim.de>
! Copyright (c) 2012-2014 Simon Zimny <s.zimny@grs-sim.de>
! Copyright (c) 2012-2017, 2019-2020 Kannan Masilamani <kannan.masilamani@uni-siegen.de>
! Copyright (c) 2012-2016 Jiaxing Qi <jiaxing.qi@uni-siegen.de>
! Copyright (c) 2013, 2019 Harald Klimach <harald.klimach@uni-siegen.de>
! Copyright (c) 2015-2016 Tobias Schneider <tobias1.schneider@student.uni-siegen.de>
! Copyright (c) 2017-2018, 2020 Raphael Haupt <raphael.haupt@uni-siegen.de>
! Copyright (c) 2020 Jana Gericke <jana.gericke@uni-siegen.de>
! Copyright (c) 2020 Peter Vitt <peter.vitt2@uni-siegen.de>
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions are met:
!
! 1. Redistributions of source code must retain the above copyright notice,
! this list of conditions and the following disclaimer.
!
! 2. Redistributions in binary form must reproduce the above copyright notice,
! this list of conditions and the following disclaimer in the documentation
! and/or other materials provided with the distribution.
!
! THIS SOFTWARE IS PROVIDED BY THE UNIVERSITY OF SIEGEN “AS IS” AND ANY EXPRESS
! OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
! OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
! IN NO EVENT SHALL UNIVERSITY OF SIEGEN OR CONTRIBUTORS BE LIABLE FOR ANY
! DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
! (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
! ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
?? include 'header/lbm_macros.inc'
?? include 'header/lbm_deriveMacros.inc'
?? include 'treelm/source/logMacros.inc'
! ****************************************************************************** !
!> Boundary condition wall treatment routines
!!
!! This module contains higher order wall treatments
!! A detailed description on the implementation details are given
!! in [[tem_bc_module]].
!!
module mus_bc_fluid_wall_module

  ! include treelm modules
  use env_module,               only: rk
  use tem_param_module,         only: cs2inv, rho0
  use tem_time_module,          only: tem_time_type
  use treelmesh_module,         only: treelmesh_type
  use tem_varSys_module,        only: tem_varSys_type
  use tem_debug_module,         only: dbgUnit
  use tem_geometry_module,      only: tem_ElemSizeLevel
  use tem_property_module,      only: prp_solid
  use tem_construction_module,  only: tem_levelDesc_type

  ! include musubi modules
  use mus_bc_header_module,      only: boundary_type, glob_boundary_type
  use mus_scheme_layout_module,  only: mus_scheme_layout_type
  use mus_field_prop_module,     only: mus_field_prop_type
  use mus_derVarPos_module,      only: mus_derVarPos_type
  use mus_param_module,          only: mus_param_type
  use mus_physics_module,        only: mus_physics_type
  use mus_mixture_module,        only: mus_mixture_type

  implicit none

  private

  ! public :: wall_multiReflection
  public :: slip_wall, spc_slip_wall
  public :: wall_libb
  public :: do_nothing
  public :: turbulent_wall_powerLaw_incomp

contains

! ****************************************************************************** !
  !> slip-wall boundary condition. Slip defined by a slip factor
  !!
  !! \li Normal velocity,\f$ u_n = 0 \f$
  !! \li Tangential velocity, \f$ \frac{\partial u_t}{\partial n} = 0 \f$
  !! \li Pressure, \f$ \frac{\partial P}{\partial n} = 0 \f$
  !! For slip-wall boundary, the slip factor will be multiplied by the velocity
  !! if slip factor = 1, then it is full/free-slip and if slip factor = 0, then
  !! it is no-slip
  !!
  !! @todo KM: Currently, free-slip boundary works only for axis-parallel planes.
  !!           Need to extend it for arbitrary geometries
  !!
  !! This subroutine's interface must match the abstract interface definition
  !! [[boundaryRoutine]] in bc/[[mus_bc_header_module]].f90 in order to be
  !! callable via [[boundary_type:fnct]] function pointer.
  subroutine slip_wall( me, state, bcBuffer, globBC, levelDesc, tree, nSize,  &
    &                   iLevel, sim_time, neigh, layout, fieldProp, varPos,   &
    &                   nScalars, varSys, derVarPos, physics, iField, mixture )
    ! -------------------------------------------------------------------- !
    !> global boundary type
    class(boundary_type) :: me
    !> Current state vector of iLevel
    real(kind=rk), intent(inout) :: state(:)
    !> size of state array ( in terms of elements )
    integer, intent(in) :: nSize
    !> state values of boundary elements of all fields of iLevel
    real(kind=rk), intent(in) :: bcBuffer(:)
    !> iLevel descriptor
    type(tem_levelDesc_type), intent(in) :: levelDesc
    !> Treelm Mesh
    type(treelmesh_type), intent(in) :: tree
    !> fluid parameters and properties
    type(mus_field_prop_type), intent(in) :: fieldProp
    !> stencil layout information
    type(mus_scheme_layout_type), intent(in) :: layout
    !> the level On which this boundary was invoked
    integer, intent(in) :: iLevel
    !> connectivity array corresponding to state vector
    integer, intent(in) :: neigh(:)
    !> global time information
    type(tem_time_type), intent(in)  :: sim_time
    !> pointer to field variable in the state vector
    integer, intent(in) :: varPos(:)
    !> number of Scalars in the scheme var system
    integer, intent(in) :: nScalars
    !> scheme variable system
    type(tem_varSys_type), intent(in) :: varSys
    !> position of derived quantities in varsys
    type(mus_derVarPos_type), intent(in) :: derVarPos
    !> scheme global boundary type
    type(glob_boundary_type), intent(in) :: globBC
    !> scheme global boundary type
    type(mus_physics_type), intent(in) :: physics
    !> current field
    integer, intent(in) :: iField
    !> mixture info
    type(mus_mixture_type), intent(in) :: mixture
    ! -------------------------------------------------------------------- !
    ! defining local variables
    real(kind=rk) :: fTmp( layout%fStencil%QQ*globBC%nElems(iLevel) )
    real(kind=rk) :: vel(3*globBC%nElems(iLevel)) ! Velocity on boundary element
    real(kind=rk) :: velTmp(3), rho
    integer :: iELem, iDir, bndNormalDir, QQ, posInBuffer
    ! ---------------------------------------------------------------------------

    QQ = layout%fStencil%QQ

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

    ! Get local velocity
    call derVarPos%velFromState( state  = fTmp ,                 &
      &                          iField = iField,                &
      &                          nElems = globBC%nElems(iLevel), &
      &                          varSys = varSys,                &
      &                          layout = layout,                &
      &                          res    = vel                    )

    do iElem = 1, globBC%nElems(iLevel)
      velTmp = vel((iElem-1)*3+1 : iELem*3) * me%slip_fac
      rho = sum(fTmp( (iElem-1)*QQ+1: (iElem-1)*QQ+QQ ))
      bndNormalDir = layout%fStencil%cxDirInv( globBC%elemLvl( iLevel )%       &
        &                                                normalInd%val( iElem ))
      !write(dbgUnit(1),*) 'bndNormalDir ',  bndNormalDir
      if( abs(layout%fStencil%cxDir( 1, bndNormalDir )) == 1) velTmp(1) = 0.0_rk
      if( abs(layout%fStencil%cxDir( 2, bndNormalDir )) == 1) velTmp(2) = 0.0_rk
      if( abs(layout%fStencil%cxDir( 3, bndNormalDir )) == 1) velTmp(3) = 0.0_rk
      !write(dbgUnit(1),*) 'velTmp ', velTmp

      do iDir = 1, layout%fStencil%QQN
        ! Write the values
        if( globBC%elemLvl(iLevel)%bitmask%val(iDir, iElem )) then
          ! Depending on PUSH or pull, use + or - for cxDir, because directions
          ! are inverted
          state(                                                               &
& ?FETCH?(iDir, iField, globBC%elemLvl(iLevel)%elem%val(iElem), QQ, nScalars, nSize,neigh))=&
          ! We need to get post-collision pdf in direction
          ! alpha- outgoing direction, which is the inverse direction of bitmask
          ! For PULL this means, get the outgoing one, as this is the one which
          ! will be bounced back
          ! For PUSH this means, get the already bounced back pdf back, so take
          ! the incoming
            & fTmp(?IDX?(layout%fStencil%cxDirInv(iDir),iElem,QQ,globBC%nElems(iLevel)))&
            &       - layout%weight( iDir )*6._rk*rho                          &
            &       * ( layout%fStencil%cxDir( 1, layout%fStencil%             &
            &                                       cxDirInv( iDir ))*velTmp(1)&
            &       +   layout%fStencil%cxDir( 2, layout%fStencil%             &
            &                                       cxDirInv( iDir ))*velTmp(2)&
            &       +   layout%fStencil%cxDir( 3, layout%fStencil%             &
            &                                       cxDirInv( iDir ))*velTmp(3))
        end if
      end do
    end do !iElem

  end subroutine slip_wall
! ****************************************************************************** !


! ****************************************************************************** !
  !> slip-wall boundary condition. Slip defined by a slip factor
  !!
  !! * Normal velocity,\( u_n = 0 \)
  !! * Tangential velocity, \( \frac{\partial u_t}{\partial n} = 0 \)
  !! * Pressure, \( \frac{\partial P}{\partial n} = 0 \)
  !!
  !! For slip-wall boundary, the slip factor will be multiplied by the velocity
  !! if slip factor = 1, then it is full/free-slip and if slip factor = 0, then
  !! it is no-slip
  !!
  !! @todo KM: Currently, free-slip boundary works only for axis-parallel planes.
  !!           Need to extend it for arbitrary geometries
  !!
  !! This subroutine's interface must match the abstract interface definition
  !! [[boundaryRoutine]] in bc/[[mus_bc_header_module]].f90 in order to be
  !! callable via [[boundary_type:fnct]] function pointer.
  subroutine spc_slip_wall( me, state, bcBuffer, globBC, levelDesc, tree,      &
    &                       nSize, iLevel, sim_time, neigh, layout, fieldProp, &
    &                       varPos, nScalars, varSys, derVarPos, physics,      &
    &                       iField, mixture                                    )
    ! -------------------------------------------------------------------- !
    !> global boundary type
    class(boundary_type) :: me
    !> Current state vector of iLevel
    real(kind=rk), intent(inout) :: state(:)
    !> size of state array ( in terms of elements )
    integer, intent(in) :: nSize
    !> state values of boundary elements of all fields of iLevel
    real(kind=rk), intent(in) :: bcBuffer(:)
    !> iLevel descriptor
    type(tem_levelDesc_type), intent(in) :: levelDesc
    !> Treelm Mesh
    type(treelmesh_type), intent(in) :: tree
    !> fluid parameters and properties
    type(mus_field_prop_type), intent(in) :: fieldProp
    !> stencil layout information
    type(mus_scheme_layout_type), intent(in) :: layout
    !> the level On which this boundary was invoked
    integer, intent(in) :: iLevel
    !> connectivity array corresponding to state vector
    integer, intent(in) :: neigh(:)
    !> global time information
    type(tem_time_type), intent(in)  :: sim_time
    !> pointer to field variable in the state vector
    integer, intent(in) :: varPos(:)
    !> number of Scalars in the scheme var system
    integer, intent(in) :: nScalars
    !> scheme variable system
    type(tem_varSys_type), intent(in) :: varSys
    !> position of derived quantities in varsys
    type(mus_derVarPos_type), intent(in) :: derVarPos
    !> scheme global boundary type
    type(glob_boundary_type), intent(in) :: globBC
    !> scheme global boundary type
    type(mus_physics_type), intent(in) :: physics
    !> current field
    integer, intent(in) :: iField
    !> mixture info
    type(mus_mixture_type), intent(in) :: mixture
    ! -------------------------------------------------------------------- !
    real(kind=rk) :: fTmp_all( layout%fStencil%QQ*globBC%nElems(iLevel) &
      &              * varSys%nStateVars )
    real(kind=rk) :: fTmp(layout%fStencil%QQ)
    real(kind=rk) :: mom(3*globBC%nElems(iLevel)), momTmp(3)
    integer :: iELem, iDir, bndNormalDir, pos, iFieldLoc, QQ, posInBuffer
    ! ------------------------------------------------------------------------
    QQ = layout%fStencil%QQ

    do iElem = 1, globBC%nElems(iLevel)
      posInBuffer = globBC%elemLvl( iLevel )%posInBcElemBuf%val( iElem )
      do iFieldLoc = 1, varSys%nStateVars
        do iDir = 1, QQ
          pos = varSys%method%val(iFieldLoc)%state_varPos(iDir)
          fTmp_all( pos+(iElem-1)*nScalars ) &
            &  = bcBuffer( pos+(posInBuffer-1)*nScalars )
        end do
      end do
    end do

    call derVarPos%momFromState( state  = fTmp_all,                &
      &                          iField = iField,                  &
      &                          nElems = globBC%nElems( iLevel ), &
      &                          varSys = varSys,                  &
      &                          layout = layout,                  &
      &                          res    = mom                      )


    do iElem = 1, globBC%nElems(iLevel)
      if( .not. btest( levelDesc%property(                        &
        &              globBC%elemLvl(iLevel)%elem%val(iElem)), prp_solid))then

      momTmp(1) = mom((iElem-1)*3+1) * me%slip_fac
      momTmp(2) = mom((iElem-1)*3+2) * me%slip_fac
      momTmp(3) = mom((iElem-1)*3+3) * me%slip_fac

      bndNormalDir = layout%fStencil%cxDirInv( globBC%elemLvl( iLevel )%       &
        &                                                normalInd%val( iElem ))
      if( abs(layout%fStencil%cxDir( 1, bndNormalDir )) == 1) momTmp(1) = 0.0_rk
      if( abs(layout%fStencil%cxDir( 2, bndNormalDir )) == 1) momTmp(2) = 0.0_rk
      if( abs(layout%fStencil%cxDir( 3, bndNormalDir )) == 1) momTmp(3) = 0.0_rk

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

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

  end subroutine spc_slip_wall
! ****************************************************************************** !

! ****************************************************************************** !
  !> No comment yet!
  !!
  !! @TODO add comment
  !!
  !! This subroutine's interface must match the abstract interface definition
  !! [[boundaryRoutine]] in bc/[[mus_bc_header_module]].f90 in order to be
  !! callable via [[boundary_type:fnct]] function pointer.
  subroutine wall_libb( me, state, bcBuffer, globBC, levelDesc, tree, nSize,  &
    &                   iLevel, sim_time, neigh, layout, fieldProp, varPos,   &
    &                   nScalars, varSys, derVarPos, physics, iField, mixture )
    ! -------------------------------------------------------------------- !
    !> global boundary type
    class(boundary_type) :: me
    !> Current state vector of iLevel
    real(kind=rk), intent(inout) :: state(:)
    !> size of state array ( in terms of elements )
    integer, intent(in) :: nSize
    !> state values of boundary elements of all fields of iLevel
    real(kind=rk), intent(in) :: bcBuffer(:)
    !> iLevel descriptor
    type(tem_levelDesc_type), intent(in) :: levelDesc
    !> Treelm Mesh
    type(treelmesh_type), intent(in) :: tree
    !> fluid parameters and properties
    type(mus_field_prop_type), intent(in) :: fieldProp
    !> stencil layout information
    type(mus_scheme_layout_type), intent(in) :: layout
    !> the level On which this boundary was invoked
    integer, intent(in) :: iLevel
    !> connectivity array corresponding to state vector
    integer, intent(in) :: neigh(:)
    !> global time information
    type(tem_time_type), intent(in)  :: sim_time
    !> pointer to field variable in the state vector
    integer, intent(in) :: varPos(:)
    !> number of Scalars in the scheme var system
    integer, intent(in) :: nScalars
    !> scheme variable system
    type(tem_varSys_type), intent(in) :: varSys
    !> position of derived quantities in varsys
    type(mus_derVarPos_type), intent(in) :: derVarPos
    !> scheme global boundary type
    type(glob_boundary_type), intent(in) :: globBC
    !> scheme global boundary type
    type(mus_physics_type), intent(in) :: physics
    !> current field
    integer, intent(in) :: iField
    !> mixture info
    type(mus_mixture_type), intent(in) :: mixture
    ! -------------------------------------------------------------------- !
    real(kind=rk) :: fIn, fOut, fNgh
    real(kind=rk) :: cIn, cOut, cNgh
    integer :: iLink
    ! ---------------------------------------------------------------------------

    !NEC$ ivdep
    !DIR$ ivdep
    !IBM* independent
    do iLink = 1, me%links(iLevel)%nVals

      cIn  = me%bouzidi(iLevel)% cIn( iLink )
      cOut = me%bouzidi(iLevel)%cOut( iLink )
      cNgh = me%bouzidi(iLevel)%cNgh( iLink )

      fIn  = bcBuffer( me%bouzidi(iLevel)% inPos(iLink) )
      fOut = bcBuffer( me%bouzidi(iLevel)%outPos(iLink) )
      fNgh = me%neigh(iLevel)%computeNeighBuf(me%bouzidi(iLevel)%nghPos(iLink))

      state( me%links(iLevel)%val(iLink) ) = cIn*fIn + cOut*fOut + cNgh*fNgh

    end do ! iLink

  end subroutine wall_libb
! ****************************************************************************** !

! ****************************************************************************** !
  !> No comment yet!
  !!
  !! @TODO add comment
  !!
  !! This subroutine's interface must match the abstract interface definition
  !! [[boundaryRoutine]] in bc/[[mus_bc_header_module]].f90 in order to be
  !! callable via [[boundary_type:fnct]] function pointer.
  subroutine do_nothing( me, state, bcBuffer, globBC, levelDesc, tree, nSize,  &
    &                    iLevel, sim_time, neigh, layout, fieldProp, varPos,   &
    &                    nScalars, varSys, derVarPos, physics, iField, mixture )
    ! -------------------------------------------------------------------- !
    !> global boundary type
    class(boundary_type) :: me
    !> Current state vector of iLevel
    real(kind=rk), intent(inout) :: state(:)
    !> size of state array ( in terms of elements )
    integer, intent(in) :: nSize
    !> state values of boundary elements of all fields of iLevel
    real(kind=rk), intent(in) :: bcBuffer(:)
    !> iLevel descriptor
    type(tem_levelDesc_type), intent(in) :: levelDesc
    !> Treelm Mesh
    type(treelmesh_type), intent(in) :: tree
    !> fluid parameters and properties
    type(mus_field_prop_type), intent(in) :: fieldProp
    !> stencil layout information
    type(mus_scheme_layout_type), intent(in) :: layout
    !> the level On which this boundary was invoked
    integer, intent(in) :: iLevel
    !> connectivity array corresponding to state vector
    integer, intent(in) :: neigh(:)
    !> global time information
    type(tem_time_type), intent(in)  :: sim_time
    !> pointer to field variable in the state vector
    integer, intent(in) :: varPos(:)
    !> number of Scalars in the scheme var system
    integer, intent(in) :: nScalars
    !> scheme variable system
    type(tem_varSys_type), intent(in) :: varSys
    !> position of derived quantities in varsys
    type(mus_derVarPos_type), intent(in) :: derVarPos
    !> scheme global boundary type
    type(glob_boundary_type), intent(in) :: globBC
    !> scheme global boundary type
    type(mus_physics_type), intent(in) :: physics
    !> current field
    integer, intent(in) :: iField
    !> mixture info
    type(mus_mixture_type), intent(in) :: mixture
    ! -------------------------------------------------------------------- !
  end subroutine do_nothing
! ****************************************************************************** !

  ! ************************************************************************** !
  !> Power-law wall function for turbulent flow.
  !! Calculate velocity and along with that the corresponding PFD at the wall
  !! using a wall function based on the power-law.
  !! The implementation is based on the following paper:
  !! Haussmann, M. et al. (2019) ‘Large-eddy simulation coupled with wall models
  !! for turbulent channel flows at high Reynolds numbers with a lattice
  !! Boltzmann method — Application to Coriolis mass flowmeter’, Computers &
  !! Mathematics with Applications. Elsevier Ltd, 78(10), pp. 3285–3302.
  !! The main formulas used from it are mentioned in the code below.
  !!
  !! The model constants are chosen according to Werner and Wengle:
  !! Wengle, H. and Werner, H. (1993) ‘Large-eddy Simulation of Turbulent Flow
  !! Over Sharp-edged Obstacles in a Plate Channel’, (1985), pp. 192–199.
  !!
  !! Note, that we use information of the second fluid element FF, to specify
  !! quantities at the first fluid element F which is adjacent to the wall W.
  !!
  !! Usage
  !! -----
  !!```lua
  !!boundary_condition = {
  !!  {
  !!    label = 'wall',
  !!    kind = 'turbulent_wall_powerlaw'
  !!  }
  !!}
  !!
  !! This subroutine's interface must match the abstract interface definition
  !! [[boundaryRoutine]] in bc/[[mus_bc_header_module]].f90 in order to be
  !! callable via [[boundary_type:fnct]] function pointer.
  subroutine turbulent_wall_powerLaw_incomp( me, state, bcBuffer, globBC,    &
    &                                        levelDesc, tree, nSize, iLevel, &
    &                                        sim_time, neigh, layout,        &
    &                                        fieldProp, varPos, nScalars,    &
    &                                        varSys, derVarPos, physics,     &
    &                                        iField, mixture                 )
    ! -------------------------------------------------------------------- !
    !> global boundary type
    class(boundary_type) :: me
    !> Current state vector of iLevel
    real(kind=rk), intent(inout) :: state(:)
    !> size of state array ( in terms of elements )
    integer, intent(in) :: nSize
    !> state values of boundary elements of all fields of iLevel
    real(kind=rk), intent(in) :: bcBuffer(:)
    !> iLevel descriptor
    type(tem_levelDesc_type), intent(in) :: levelDesc
    !> Treelm Mesh
    type(treelmesh_type), intent(in) :: tree
    !> fluid parameters and properties
    type(mus_field_prop_type), intent(in) :: fieldProp
    !> stencil layout information
    type(mus_scheme_layout_type), intent(in) :: layout
    !> the level On which this boundary was invoked
    integer, intent(in) :: iLevel
    !> connectivity array corresponding to state vector
    integer, intent(in) :: neigh(:)
    !> global time information
    type(tem_time_type), intent(in)  :: sim_time
    !> pointer to field variable in the state vector
    integer, intent(in) :: varPos(:)
    !> number of Scalars in the scheme var system
    integer, intent(in) :: nScalars
    !> scheme variable system
    type(tem_varSys_type), intent(in) :: varSys
    !> position of derived quantities in varsys
    type(mus_derVarPos_type), intent(in) :: derVarPos
    !> scheme global boundary type
    type(glob_boundary_type), intent(in) :: globBC
    !> scheme global boundary type
    type(mus_physics_type), intent(in) :: physics
    !> current field
    integer, intent(in) :: iField
    !> mixture info
    type(mus_mixture_type), intent(in) :: mixture
    ! -------------------------------------------------------------------- !
    ! --------------------------------------------------------------------------
    ! Constants for the power-law model according to Werner & Wengle (1993)
    real(kind=rk) :: m, C_m, yPlus_v
    ! Auxiliary constants in order to calculate them only once
    real(kind=rk) :: onePlusM, oneMinusM, expC, velLimiter
    ! Number of links corresponding to DdQq with w q=QQ
    integer :: QQ
    ! Self-describing variables, loop indices, etc.
    integer :: iDir, iElem, normalIn, normalOut
    real(kind=rk) :: normal(3)
    ! Quantities related to first fluid element F (the one next to the wall)
    ! Its position in the state array
    integer :: elemPos
    ! Density on first fluid element F (the one next to the wall)
    real(kind=rk) :: rhoF
    ! Velocity on first fluid element F and its stream-wise component (_SW)
    real(kind=rk) :: velF(3), velF_SW
    ! Auxiliary term used for calculation of velF_SW and for the if condition
    real(kind=rk) :: auxVelF
    ! Its equilibrium part (Eq)
    real(kind=rk) :: feqF(layout%fStencil%QQ)
    ! Quantities related to second fluid element FF
    ! Its position in the state array
    integer :: neighPos
    ! Density on FF and its inverse
    real(kind=rk) :: rhoFF, rhoFF_inv
    ! Velocity on second fluid element FF and its stream-wise component (_SW)
    real(kind=rk) :: velFF(3), velFF_SW
    ! Temporary pdf of second fluid element FF
    real(kind=rk) :: pdfFF( layout%fStencil%QQ )
    ! Its equilibrium (eq) and non-equilibrium (nonEq) part
    real(kind=rk) :: feqFF(layout%fStencil%QQ), fneqFF(layout%fStencil%QQ)
    ! Unit vector in stream-wise direction (SW) and aux-variable for its calc.
    real(kind=rk) :: unitSW(3), numerator(3)
    ! Kinematic viscosity and corresponding omega of FF
    real(kind=rk) :: viscKineFF, omegaKineFF
    ! Element size (used for the following distance calculation: yWF, yWFF)
    real(kind=rk) :: dx
    ! Distance from physical wall W to bary center of first fluid element F
    real(kind=rk) :: yWF
    ! Distance from physical wall W to bary center of second fluid element FF
    real(kind=rk) :: yWFF
    ! Wall-shear stress
    real(kind=rk) :: tau_w
    ! qVal used for workaround at the beginning. (Later one it should be
    ! distinguished between curved and straight boundaries.
    real(kind=rk) :: qVal
    ! --------------------------------------------------------------------------
    ! Set model constants to suggested values proposed by Werner & Wengle
    m = 1.0_rk / 7.0_rk
    C_m = 8.3_rk
    ! Height of normalized viscous sublayer
    yPlus_v = 11.81_rk

    ! Auxiliary constants, which are used multiple times in the element-loop to
    ! calculate tau_w
    onePlusM = 1 + m
    oneMinusM = 1 - m
    expC = onePlusM / oneMinusM
    ! Auxiliary variable used for the condition in the calculation of tau_w
    velLimiter = 0.5_rk * yPlus_v ** 2

    ! Load stencil for velocity space (DdQq) with q=QQ
    QQ = layout%fStencil%QQ

    ! Get element size on current level (dx). It's needed for distance calculat.
    dx = tem_ElemSizeLevel(tree, iLevel)

    ! Start of the main calculations by looping over elements with boundaries.
    do iElem = 1, globBC%nElems(iLevel)

      ! 0) Preliminary steps.
      ! Current element position in state array. (Used for state in step 10.)
      elemPos = globBC%elemLvl(iLevel)%elem%val(iElem)
      ! Get kinematic viscosity on second fluid element (viscKineFF) and
      ! corresponding omega (omegaKineFF). Access it via the position of the
      ! neighbor element (FF) in state array.
      neighPos =  me%neigh(iLevel)%posInState(1,iElem)
      viscKineFF = fieldProp%fluid%viscKine%dataOnLvl(iLevel)%val(neighPos)
      omegaKineFF = fieldProp%fluid%viscKine%omLvl(iLevel)%val(neighPos)

      ! 1) Get density on second fluid element (rhoFF) as sum of pdf.
      ! First, get PDF for FF.
      pdfFF = me%neigh(iLevel)%neighBufferPost(1,(iElem-1)*QQ+1:(iElem-1)*QQ+QQ)
      ! Second, calculate density as sum of pdf.
      rhoFF = sum(pdfFF)
      ! Third, calculate its inverse, which is needed as input for velMacro.
      rhoFF_inv = 1.0_rk / rhoFF

      ! 2) Compute velocity on second fluid element (velFF) using a macro with
      !    pdf, inverse of density and stencil as input.
?? copy :: velocity(velFF, pdfFF, rhoFF_inv, layout%fStencil)

      ! 3) Calculate unit vector for stream-wise direction (unitSW) (eq.28).
      ! First, get normal for that. In Musubi it's pointing into the domain.
      normalIn = globBC%elemLvl(iLevel)%normalInd%val(iElem)
      ! But we need the one poiting out, so invert it.
      normalOut = layout%fStencil%cxDirInv(normalIn)
      normal = layout%fStencil%cxDirRK(:, normalOut)
      ! Second, compute unitSW using velFF and normal via intermediate step.
      numerator = velFF - dot_product( velFF, normal ) * normal
      unitSW = numerator / sqrt( dot_product( numerator, numerator ) )

      ! 4) Calculate stream-wise velocity on fluid element (velFF_SW) (eq.29).
      velFF_SW = dot_product( velFF(:), unitSW(:) )

      ! 5) Calculate distance from physical wall W to second fluid element FF
      !    using qVal
      ! TODO: Enable qVals for straight and curved. At the moment, they are set:
      qVal = 0.5_rk
      ! Use qVal to calculate distance from wall W to first fluid node F: yWF
      yWF = qVal * dx * sqrt( dot_product( normal, normal) )
      ! Use yWF to calculate distance from wall W to second fluid node FF: yWFF
      yWFF = yWF + dx * sqrt( dot_product( normal, normal) )

      ! 6) Use stream-wise vel. of FF (velFF_SW) and distance of FF to the wall
      !    (yWFF) to calculate the wall-shear stress (eq.22).
      ! If yPlus is in range of viscous sublayer, use linear profile.
      if ( velFF_SW <= ( velLimiter * viscKineFF  / yWFF ) ) then
        tau_w = ( 2.0_rk * rho0 * viscKineFF * velFF_SW ) / yWFF
      ! If it's greater, use logarithmic profile.
      else
        tau_w = rho0 * ( ( oneMinusM / 2.0_rk ) * C_m ** expC                  &
          &                 * ( viscKineFF / yWFF ) ** onePlusM                &
          &                 + ( onePlusM / C_m ) * ( viscKineFF / yWFF ) ** m  &
          &                 * velFF_SW  )                                      &
          &               ** ( 2.0_rk / onePlusM )
      end if

      ! 7) Calculate stream-wise velocity on first fluid F (velF_SW) using yWF
      !    (step 5) and tau_w (step 6), (eq.21).
      ! Before, calculate the term they have in common (auxVelF):
      auxVelF = yWF / viscKineFF * sqrt( tau_w / rho0 )
      ! If yPlus is in range of viscous sublayer, use linear profile.
      if ( auxVelF <= yPlus_v ) then
        velF_SW = auxVelF
      ! If it's greater, use logarithmic profile.
      else
        velF_SW = C_m * auxVelF ** m
      end if

      ! 8) Compute velocity vector of first fluid (velF) using velF_SW & unitSW.
      velF(:) = velF_SW * unitSW(:)

      ! 9) Set density on first fluid F (rhoF) to the one of second (rhoFF).
      rhoF = rhoFF

      ! 10) Reconstruct the PDF on first fluid F as sum of Eq and nonEq part.
      !     Eq is calculated on F, while nonEq is calculated on FF (based on
      !     extrapolation scheme proposed by Guo et al. (see nonEqExpol)).
      ! First, get eq parts on F (fEqF) and FF (fEqFF) using a macro with
      ! density and velocity as input
      do iDir = 1, layout%fStencil%QQ
?? copy :: EqDir_incomp(fEqf(iDir), rhoF, velF, layout, iDir, rho0, cs2inv, rk)
      end do
      do iDir = 1, layout%fStencil%QQ
?? copy :: EqDir_incomp(fEqFF(iDir), rhoFF, velFF, layout, iDir, rho0, cs2inv, rk)
      end do
      ! Second, calculate nonEq part for FF (fneqFF). As it is calculated from
      ! post-collision PDF, we have to divide it by (1-omega).
      fneqFF = ( pdfFF - feqFF ) / ( 1.0_rk - omegaKineFF )
      ! Finally, reconstruct the PDF on first fluid F as sum of Eq and nonEq and
      ! write it into the state for first fluid F using the FETCH-macro. This
      ! has to be done for all directions of QQ.
      do iDir = 1, layout%fStencil%QQ
        state( ?FETCH?( iDir, iField, elemPos, QQ, nScalars,nSize, neigh)) &
          & = feqF(iDir) + fneqFF(iDir)
      end do

    end do ! iElem

  end subroutine turbulent_wall_powerLaw_incomp
! ****************************************************************************** !

end module mus_bc_fluid_wall_module
! ****************************************************************************** !
