! Copyright (c) 2018-2020 Kannan Masilamani <kannan.masilamani@uni-siegen.de>
! Copyright (c) 2019 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.
! ************************************************************************** !
!> author: Kannan Masilamani
!! Average Interpolation of flow quantities between different grid levels
!!
?? include 'header/lbm_macros.inc'
?? include 'header/lbm_deriveMacros.inc'
?? include 'header/lbm_d2q9Macros.inc'
?? include 'treelm/source/logMacros.inc'
module mus_interpolate_average_module
  use iso_c_binding, only: c_loc, c_ptr, c_f_pointer

  ! include treelm modules
  use env_module,              only: rk
  use tem_aux_module,          only: tem_abort
  use tem_element_module,      only: eT_GhostFromCoarser, &
    &                                eT_ghostFromFiner
  use tem_param_module,        only: cs2inv, cs2, PI, div1_2, div1_9, div4_9,&
    &                                div1_36, div2_3, rho0, rho0Inv
  use tem_comm_env_module,     only: tem_comm_env_type
  use tem_debug_module,        only: dbgUnit
  use tem_construction_module, only: tem_levelDesc_type
  use tem_matrix_module,       only: tem_matrix_type
  use tem_stencil_module,      only: tem_stencilHeader_type
  use tem_logging_module,      only: logUnit
  use tem_varSys_module,       only: tem_varSys_type
  use tem_time_module,         only: tem_time_type

?? IF (loglvl >= 1) THEN
  use env_module,            only: long_k
  use tem_topology_module,   only: tem_coordOfID
  use tem_param_module,      only: PI
?? ENDIF

  ! include musubi modules
?? IF( .not. SOA) then
  use mus_interpolate_debug_module,  only: TGV_2D
?? ENDIF
  use mus_scheme_layout_module,      only: mus_scheme_layout_type
  use mus_physics_module,            only: mus_physics_type
  use mus_interpolate_header_module, only: mus_interpolation_method_type
  use mus_interpolate_tools_module,  only: mus_intp_convertMomToPDF3D,        &
    &                                      mus_intp_convertMomToPDF3D_incomp, &
    &                                      mus_intp_convertMomToPDF2D,        &
    &                                      mus_intp_convertMomToPDF2D_incomp, &
    &                                      mus_intp_getMoments
  use mus_field_prop_module,         only: mus_field_prop_type
  use mus_fluid_module,              only: mus_fluid_type
  use mus_relaxationParam_module,    only: mus_calcOmegaFromVisc
  use mus_derVarPos_module,          only: mus_derVarPos_type

  implicit none

  private

  public :: fillArbiMyGhostsFromFiner_avg
  public :: fillArbiFinerGhostsFromMe_weighAvg
  public :: fillMyGhostsFromFiner_avg
  public :: fillMyGhostsFromFiner_avgLES
  public :: fillMyGhostsFromFiner_avgIncomp
  public :: fillMyGhostsFromFiner_avgIncompLES
  public :: fillMyGhostsFromFiner_avg2D
  public :: fillMyGhostsFromFiner_avg2DIncomp
  public :: fillFinerGhostsFromMe_weighAvg
  public :: fillFinerGhostsFromMe_weighAvgLES
  public :: fillFinerGhostsFromMe_weighAvgIncomp
  public :: fillFinerGhostsFromMe_weighAvgIncompLES
  public :: fillFinerGhostsFromMe_weighAvg2D
  public :: fillFinerGhostsFromMe_weighAvg2DIncomp

  contains


! **************************************************************************** !
  !> Interpolate auxiliary field from fine source to coarse target
  !!
  !! This subroutine's interface must match the abstract interface definition
  !! [[intpRoutine_arbitraryVal]] in intp/[[mus_interpolate_header_module]].f90
  !! in order to be callable via
  !! [[mus_interpolation_method_type:do_intpArbiVal]] function pointer.
  subroutine fillArbiMyGhostsFromFiner_avg( method, tLevelDesc, level,    &
    &                                       stencil, sVal, snSize, tVal,  &
    &                                       tnSize, nTargets, targetList, &
    &                                       nScalars                      )
    ! ------------------------------------------------------------------ !
    class(mus_interpolation_method_type), intent(in) :: method

    !> my refinement level
    integer, intent(in) :: level

    !> stencil header
    type(tem_stencilHeader_type), intent(in) :: stencil

    !> State vector of SOURCE FLUID elements
    real(kind=rk), intent(in) :: sVal(:)
    integer, intent(in) :: snSize

    !> State vector of TARGET GHOST elements
    real(kind=rk), intent(inout) :: tVal(:)
    integer, intent(in) :: tnSize

    !> level descriptor on target level
    type( tem_levelDesc_type ), intent(in) :: tLevelDesc

    !> List of target elements ( their position in depSource list )
    integer, intent(in) :: nTargets
    integer, intent(in) :: targetList(nTargets)

    !> number of scalars to interpolate
    integer, intent(in) :: nScalars
    ! ------------------------------------------------------------------ !
    ! --------------------------------------------------------------------------
    integer :: sourceLevel    ! level of source elements
    integer :: sourceElem     ! treeId of current source element
    integer :: targetLevel    ! level of target elements
    integer :: targetElem     ! treeId of current source element
    integer :: iElem          ! current target element (for outer loop)
    integer :: indElem        ! element counter for indirection list
    integer :: iSourceElem    ! current source element (for inner loop)
    integer :: nSourceElems   ! number of source elements for the current target
    real(kind=rk) :: inv_nSourceElems
    real(kind=rk) :: tArbi(nScalars)  ! target auxField
    real(kind=rk) :: sArbi(nScalars)  ! temp source ArbiField
    ! --------------------------------------------------------------------------
    sourceLevel = level + 1
    targetLevel = level

    ! Treat all coarse target elements
    do indElem = 1, nTargets

      iElem = targetList( indElem )

      ! Read the target element
      targetElem = iElem + tLevelDesc%offset( 1, eT_ghostFromFiner)

      ! Get how many fine source elements we have for interpolation.
      nSourceElems = tLevelDesc%depFromFiner( iElem )%elem%nVals
      inv_nSourceElems = 1.0_rk / dble(nSourceElems)

      tArbi = 0.0_rk
      ! Now loop over all fine source elements for this target:
      do iSourceElem = 1, nSourceElems

        ! Get the source element
        sourceElem = tLevelDesc%depFromFiner( iElem ) &
          &                    %elem%val( iSourceElem )

        ! Get souce auxilary variables
        sArbi(1:nScalars) = sVal((sourceElem-1)*nScalars+1:sourceElem*nScalars)

        ! target aux
        tArbi = sArbi + tArbi

      end do  ! iSourceElem

      ! interpolate all auxiliary variables by average
      tVal((targetElem-1)*nScalars+1: targetElem*nScalars) &
        & = tArbi(:) * inv_nSourceElems

    enddo

  end subroutine fillArbiMyGhostsFromFiner_avg
! **************************************************************************** !

! **************************************************************************** !
  !> [Average interpolation](../page/features/intp_methods.html) of ghostFromFiner
  !! The interpolation procedure used in this routine is:\n
  !! 1. Calculate moments from source
  !! 2. Interpolate moments on target by simple average
  !! 3. Store target auxilary field
  !! 4. Compute viscosity on target element and compute source and target omega
  !! 5. Get nonEq scaling factor depeding on scheme layout and relaxation
  !! 6. Calculate Equilibrium and nonEquilibrium
  !! 7. calculate target: Eq + Scale * nonEquilibrium
  !!
  !! This subroutine's interface must match the abstract interface definition
  !! [[intpRoutine]] in intp/[[mus_interpolate_header_module]].f90 in order to
  !! be callable via [[mus_interpolation_method_type:do_intp]] function pointer.
  subroutine fillMyGhostsFromFiner_avg( method, fieldProp, tLevelDesc, level, &
    &                                   sState, snSize, tState, tnSize,       &
    &                                   tAuxField, layout, nTargets,          &
    &                                   targetList, physics, time, varSys,    &
    &                                   derVarPos                             )
    ! -------------------------------------------------------------------- !
    class(mus_interpolation_method_type), intent(in) :: method

    !> Array of field properties (fluid or species)
    type(mus_field_prop_type), target, intent(in) :: fieldProp(:)

    !> level descriptor on target level
    type( tem_levelDesc_type ), intent(in) :: tLevelDesc

    !> my refinement level
    integer, intent(in) :: level

    !> State vector of SOURCE FLUID elements
    real(kind=rk), intent(in) :: sState(:)
    ! integer, intent(in) :: sNeigh(:)
    integer, intent(in) :: snSize

    !> State vector of TARGET GHOST elements
    real(kind=rk), intent(inout) :: tState(:)
    ! integer, intent(in) :: tNeigh(:)
    integer, intent(in) :: tnSize

    !> AuxField variable to fill on target GHOST elements
    real(kind=rk), intent(inout) :: tAuxField(:)

    !> the layout used
    type( mus_scheme_layout_type ), intent(in) :: layout

    !> List of target elements ( their position in depSource list )
    integer, intent(in) :: nTargets
    integer, intent(in) :: targetList(nTargets)

    !> physics type to convert lattice to physics SI unit and vice versa
    !! @todo: This can be replaced by scale factor by level
    type( mus_physics_type ), intent(in) :: physics

    !> time required to compute viscosity on target element barycenter
    type(tem_time_type), intent(in) :: time

    !> scheme variable system
    type( tem_varSys_type ), intent(in) :: varSys

    !> position of all derive variable in varSys for all fields
    type(mus_derVarPos_type), intent(in) :: derVarPos(:)
    ! -------------------------------------------------------------------- !
    integer :: sourceLevel    ! level of source elements
    integer :: sourceElem     ! treeId of current source element
    integer :: targetLevel    ! level of target elements
    integer :: targetElem     ! treeId of current source element
    integer :: iDir           ! current direction (discrete velocity) for loop
    integer :: iElem          ! current target element (for outer loop)
    integer :: indElem        ! element counter for indirection list
    integer :: iSourceElem    ! current source element (for inner loop)
    integer :: nSourceElems   ! number of source elements for the current target
    real(kind=rk) :: tPDF( layout%fStencil%QQ )
    integer :: QQ
    real(kind=rk) :: inv_nSourceElems
    real(kind=rk) :: tMom(layout%fStencil%QQ)  ! target moment calculation
    real(kind=rk) :: srcMom(layout%fStencil%QQ)  ! temp moment calculation
    type(mus_fluid_type), pointer :: fluid
    real(kind=rk) :: nonEqScalingFacs(layout%fStencil%QQ)
    real(kind=rk) :: sOmegaKine, tOmegaKine, tVisc
    real(kind=rk) :: invRho
    integer :: nScalars, tOffset
    integer :: dens_pos, vel_pos(3)
    ! --------------------------------------------------------------------------
    fluid => fieldProp(1)%fluid
    nScalars = varSys%nScalars
    QQ = layout%fStencil%QQ
    dens_pos = varSys%method%val(derVarPos(1)%density)%auxField_varPos(1)
    vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:3)

    sourceLevel = level + 1
    targetLevel = level

    ! Treat all coarse target elements
    do indElem = 1, nTargets

      iElem = targetList( indElem )

      ! Read the target element
      targetElem = iElem + tLevelDesc%offset( 1, eT_ghostFromFiner)


      ! Get how many fine source elements we have for interpolation.
      nSourceElems = tLevelDesc%depFromFiner( iElem )%elem%nVals
      inv_nSourceElems = 1.0_rk / dble(nSourceElems)

      tMom = 0.0_rk
      ! Now loop over all fine source elements for this target:
      do iSourceElem = 1, nSourceElems

        ! Get the source element
        sourceElem = tLevelDesc%depFromFiner( iElem ) &
          &                    %elem%val( iSourceElem )

        ! Get macroscopic moments
        srcMom = mus_intp_getMoments(state     = sstate,                 &
          &                          elem      = sourceElem,             &
          &                          QQ        = QQ,                     &
          &                          nScalars  = nScalars,               &
          &                          nSize     = snSize,                 &
          &                          toMoments = layout%moment%toMoments )

        tMom = srcMom + tMom

      end do  ! iSourceElem

      ! interpolate all moments by average
      tMom(:) = tMom(:) * inv_nSourceElems

      ! store interpolated target auxField
      invRho = 1.0_rk/tMom(1)
      tOffset = (targetElem-1)*varSys%nAuxScalars
      tAuxField(tOffset+dens_pos)   = tMom(1)
      tAuxField(tOffset+vel_pos(1)) = tMom(layout%moment%first_moments(1))*invRho
      tAuxField(tOffset+vel_pos(2)) = tMom(layout%moment%first_moments(2))*invRho
      tAuxField(tOffset+vel_pos(3)) = tMom(layout%moment%first_moments(3))*invRho

      ! get normalized kinematic viscosity on target element
      tVisc = fluid%viscKine%dataOnLvl(targetLevel)%val(targetElem)

      ! relation between coarse and fine grid kinematic viscosity:
      ! v^s_f = 2 v^s_c
      ! calculate omega on source and target level
      sOmegaKine = mus_calcOmegaFromVisc(2.0_rk*tVisc)
      tOmegaKine = mus_calcOmegaFromVisc(tVisc)

      ! scaling factors for nonEquilibrium moments
      nonEqScalingFacs = fluid%nonEqScalingFacs(           &
        & omegaKine_SRC = sOmegaKine,                      &
        & omegaKine_TGT = tOmegaKine,                      &
        & omegaBulk_SRC = fluid%omegaBulkLvl(sourceLevel), &
        & omegaBulk_TGT = fluid%omegaBulkLvl(targetLevel), &
        & scaleFac      = 2.0_rk,                          &
        & QQ            = QQ                               )
      !write(dbgUnit(1),*) 'avg nonEqFacs', nonEqScalingFacs

      ! Convert moment to PDF
      tPDF = mus_intp_convertMomToPDF3D( moments          = tMom,             &
        &                                nonEqScalingFacs = nonEqScalingFacs, &
        &                                layout           = layout            )

      ! Now write the resulting pdf in the current direction to the target
      ! Element position
      do iDir = 1, QQ
        tstate(?IDX?( iDir, targetElem, nscalars, tnSize ) ) &
          & = tPDF(iDir)
      enddo ! iDir

    enddo

  end subroutine fillMyGhostsFromFiner_avg
! **************************************************************************** !

! **************************************************************************** !
  !> [Average interpolation](../page/features/intp_methods.html) of ghostFromFiner
  !! The interpolation procedure used in this routine is:\n
  !! 1. calculate moments from source
  !! 2. Interpolate moments on target by simple average
  !! 3. Calculate Equilibrium and nonEquilibrium
  !! 4. calculate target: Eq + Scale * nonEquilibrium
  !!
  !! This subroutine's interface must match the abstract interface definition
  !! [[intpRoutine]] in intp/[[mus_interpolate_header_module]].f90 in order to
  !! be callable via [[mus_interpolation_method_type:do_intp]] function pointer.
  subroutine fillMyGhostsFromFiner_avgLES( method, fieldProp, tLevelDesc, &
    &                                      level, sState, snSize, tState, &
    &                                      tnSize, tAuxField, layout,     &
    &                                      nTargets, targetList, physics, &
    &                                      time, varSys, derVarPos        )
    ! -------------------------------------------------------------------- !
    class(mus_interpolation_method_type), intent(in) :: method

    !> Array of field properties (fluid or species)
    type(mus_field_prop_type), target, intent(in) :: fieldProp(:)

    !> level descriptor on target level
    type( tem_levelDesc_type ), intent(in) :: tLevelDesc

    !> my refinement level
    integer, intent(in) :: level

    !> State vector of SOURCE FLUID elements
    real(kind=rk), intent(in) :: sState(:)
    ! integer, intent(in) :: sNeigh(:)
    integer, intent(in) :: snSize

    !> State vector of TARGET GHOST elements
    real(kind=rk), intent(inout) :: tState(:)
    ! integer, intent(in) :: tNeigh(:)
    integer, intent(in) :: tnSize

    !> AuxField variable to fill on target GHOST elements
    real(kind=rk), intent(inout) :: tAuxField(:)

    !> the layout used
    type( mus_scheme_layout_type ), intent(in) :: layout

    !> List of target elements ( their position in depSource list )
    integer, intent(in) :: nTargets
    integer, intent(in) :: targetList(nTargets)

    !> physics type to convert lattice to physics SI unit and vice versa
    !! @todo: This can be replaced by scale factor by level
    type( mus_physics_type ), intent(in) :: physics

    !> time required to compute viscosity on target element barycenter
    type(tem_time_type), intent(in) :: time

    !> scheme variable system
    type( tem_varSys_type ), intent(in) :: varSys

    !> position of all derive variable in varSys for all fields
    type(mus_derVarPos_type), intent(in) :: derVarPos(:)
    ! -------------------------------------------------------------------- !
    integer :: sourceLevel    ! level of source elements
    integer :: sourceElem     ! treeId of current source element
    integer :: targetLevel    ! level of target elements
    integer :: targetElem     ! treeId of current source element
    integer :: iDir           ! current direction (discrete velocity) for loop
    integer :: iElem          ! current target element (for outer loop)
    integer :: indElem        ! element counter for indirection list
    integer :: iSourceElem    ! current source element (for inner loop)
    integer :: nSourceElems   ! number of source elements for the current target
    real(kind=rk) :: tPDF( layout%fStencil%QQ )
    integer :: QQ
    real(kind=rk) :: inv_nSourceElems
    real(kind=rk) :: tMom(layout%fStencil%QQ)  ! target moment calculation
    real(kind=rk) :: srcMom(layout%fStencil%QQ)  ! temp moment calculation
    type(mus_fluid_type), pointer :: fluid
    real(kind=rk) :: nonEqScalingFacs(layout%fStencil%QQ)
    real(kind=rk) :: sOmegaKine, tOmegaKine, sVisc, tVisc, invRho
    real(kind=rk) :: tTurbVisc
    integer :: nScalars, tOffset
    integer :: dens_pos, vel_pos(3)
    ! --------------------------------------------------------------------------
    fluid => fieldProp(1)%fluid
    nScalars = varSys%nScalars
    QQ = layout%fStencil%QQ
    dens_pos = varSys%method%val(derVarPos(1)%density)%auxField_varPos(1)
    vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:3)

    sourceLevel = level + 1
    targetLevel = level

    ! Treat all coarse target elements
    do indElem = 1, nTargets

      iElem = targetList( indElem )

      ! Read the target element
      targetElem = iElem + tLevelDesc%offset( 1, eT_ghostFromFiner)

      ! Get how many fine source elements we have for interpolation.
      nSourceElems = tLevelDesc%depFromFiner( iElem )%elem%nVals
      inv_nSourceElems = 1.0_rk / dble(nSourceElems)

      tMom = 0.0_rk
      tTurbVisc = 0.0_rk
      ! Now loop over all fine source elements for this target:
      do iSourceElem = 1, nSourceElems

        ! Get the source element
        sourceElem = tLevelDesc%depFromFiner( iElem ) &
          &                    %elem%val( iSourceElem )

        ! Get macroscopic moments
        srcMom = mus_intp_getMoments(state     = sstate,                 &
          &                          elem      = sourceElem,             &
          &                          QQ        = QQ,                     &
          &                          nScalars  = nScalars,               &
          &                          nSize     = snSize,                 &
          &                          toMoments = layout%moment%toMoments )

        tMom = srcMom + tMom

        tTurbVisc = tTurbVisc + fluid%turbulence%dataOnLvl(sourceLevel) &
          &                                      %visc(sourceElem)

      end do  ! iSourceElem

      ! interpolate all moments by average
      tMom(:) = tMom(:) * inv_nSourceElems

      ! store interpolated target auxField
      invRho = 1.0_rk/tMom(1)
      tOffset = (targetElem-1)*varSys%nAuxScalars
      tAuxField(tOffset+dens_pos)   = tMom(1)
      tAuxField(tOffset+vel_pos(1)) = tMom(layout%moment%first_moments(1))*invRho
      tAuxField(tOffset+vel_pos(2)) = tMom(layout%moment%first_moments(2))*invRho
      tAuxField(tOffset+vel_pos(3)) = tMom(layout%moment%first_moments(3))*invRho

      ! interpolate turbulent viscosity on target element
      tTurbVisc = tTurbVisc * inv_nSourceElems
      ! scale interpolated turbulent viscosity to target element
      fluid%turbulence%dataOnLvl(targetLevel)%visc(targetElem) &
        & = fluid%turbulence%fac_f2c*tTurbVisc

      ! get normalized kinematic viscosity on target element
      tVisc = fluid%viscKine%dataOnLvl(targetLevel)%val(targetElem)

      ! relation between coarse and fine grid kinematic viscosity:
      ! v^s_f = 2 v^s_c
      ! total viscosity on source element
      sVisc = 2.0_rk * tVisc + tTurbVisc
      ! total viscosity on target element
      tVisc = tVisc + fluid%turbulence%fac_f2c*tTurbVisc

      ! calculate omega on source and target level
      sOmegaKine = mus_calcOmegaFromVisc(sVisc)
      tOmegaKine = mus_calcOmegaFromVisc(tVisc)

      ! scaling factors for nonEquilibrium moments
      nonEqScalingFacs = fluid%nonEqScalingFacs(           &
        & omegaKine_SRC = sOmegaKine,                      &
        & omegaKine_TGT = tOmegaKine,                      &
        & omegaBulk_SRC = fluid%omegaBulkLvl(sourceLevel), &
        & omegaBulk_TGT = fluid%omegaBulkLvl(targetLevel), &
        & scaleFac      = 2.0_rk,                          &
        & QQ            = QQ                               )
      !write(dbgUnit(1),*) 'avg nonEqFacs', nonEqScalingFacs

      ! Convert moment to PDF
      tPDF = mus_intp_convertMomToPDF3D( moments          = tMom,             &
        &                                nonEqScalingFacs = nonEqScalingFacs, &
        &                                layout           = layout            )

      ! Now write the resulting pdf in the current direction to the target
      ! Element position
      do iDir = 1, QQ
        tstate(?IDX?( iDir, targetElem, nscalars, tnSize ) ) &
          & = tPDF(iDir)
      enddo ! iDir

    enddo

  end subroutine fillMyGhostsFromFiner_avgLES
! **************************************************************************** !

! **************************************************************************** !
  !> [Average interpolation](../page/features/intp_methods.html) of ghostFromFiner
  !! The interpolation procedure used in this routine is:\n
  !! 1. Calculate moments from source
  !! 2. Interpolate moments on target by simple average
  !! 3. Store target auxilary field
  !! 4. Get viscosity on target element and compute source and target omega
  !! 5. Get nonEq scaling factor depeding on scheme layout and relaxation
  !! 6. Calculate Equilibrium and nonEquilibrium
  !! 7. calculate target: Eq + Scale * nonEquilibrium
  !! This routine is used by 3D acoustic average interpolation for
  !! incompressible model.
  !!
  !! This subroutine's interface must match the abstract interface definition
  !! [[intpRoutine]] in intp/[[mus_interpolate_header_module]].f90 in order to
  !! be callable via [[mus_interpolation_method_type:do_intp]] function pointer.
  subroutine fillMyGhostsFromFiner_avgIncomp( method, fieldProp, tLevelDesc, &
    &                                         level, sState, snSize, tState, &
    &                                         tnSize, tAuxField, layout,     &
    &                                         nTargets, targetList, physics, &
    &                                         time, varSys, derVarPos        )
    ! -------------------------------------------------------------------- !
    class(mus_interpolation_method_type), intent(in) :: method

    !> Array of field properties (fluid or species)
    type(mus_field_prop_type), target, intent(in) :: fieldProp(:)

    !> level descriptor on target level
    type( tem_levelDesc_type ), intent(in) :: tLevelDesc

    !> my refinement level
    integer, intent(in) :: level

    !> State vector of SOURCE FLUID elements
    real(kind=rk), intent(in) :: sState(:)
    ! integer, intent(in) :: sNeigh(:)
    integer, intent(in) :: snSize

    !> State vector of TARGET GHOST elements
    real(kind=rk), intent(inout) :: tState(:)
    ! integer, intent(in) :: tNeigh(:)
    integer, intent(in) :: tnSize

    !> AuxField variable to fill on target GHOST elements
    real(kind=rk), intent(inout) :: tAuxField(:)

    !> the layout used
    type( mus_scheme_layout_type ), intent(in) :: layout

    !> List of target elements ( their position in depSource list )
    integer, intent(in) :: nTargets
    integer, intent(in) :: targetList(nTargets)

    !> physics type to convert lattice to physics SI unit and vice versa
    !! @todo: This can be replaced by scale factor by level
    type( mus_physics_type ), intent(in) :: physics

    !> time required to compute viscosity on target element barycenter
    type(tem_time_type), intent(in) :: time

    !> scheme variable system
    type( tem_varSys_type ), intent(in) :: varSys

    !> position of all derive variable in varSys for all fields
    type(mus_derVarPos_type), intent(in) :: derVarPos(:)
    ! -------------------------------------------------------------------- !
    integer :: sourceLevel    ! level of source elements
    integer :: sourceElem     ! treeId of current source element
    integer :: targetLevel    ! level of target elements
    integer :: targetElem     ! treeId of current source element
    integer :: iDir           ! current direction (discrete velocity) for loop
    integer :: iElem          ! current target element (for outer loop)
    integer :: indElem        ! element counter for indirection list
    integer :: iSourceElem    ! current source element (for inner loop)
    integer :: nSourceElems   ! number of source elements for the current target
    real(kind=rk) :: tPDF( layout%fStencil%QQ ) ! pdf to reconstruct from
    ! moments of the source elements' pdf
    real(kind=rk) :: srcMom(layout%fStencil%QQ)  ! temp moment calculation
    real(kind=rk) :: tMom(layout%fStencil%QQ)  ! target moment calculation
    integer :: QQ
    real(kind=rk) :: inv_nSourceElems
    type(mus_fluid_type), pointer :: fluid
    real(kind=rk) :: nonEqScalingFacs(layout%fStencil%QQ)
    real(kind=rk) :: sOmegaKine, tOmegaKine, tVisc
    integer :: nScalars, tOffset
    integer :: dens_pos, vel_pos(3)
    ! --------------------------------------------------------------------------
    fluid => fieldProp(1)%fluid
    QQ = layout%fStencil%QQ
    nScalars = varSys%nScalars
    dens_pos = varSys%method%val(derVarPos(1)%density)%auxField_varPos(1)
    vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:3)

    sourceLevel = level + 1
    targetLevel = level

    ! Treat all coarse target elements
    do indElem = 1, nTargets

      iElem = targetList( indElem )

      ! Read the target element treeId
      targetElem = iElem + tLevelDesc%offset( 1, eT_ghostFromFiner)

      ! Find out how many fine source elements we have for interpolation.
      nSourceELems = tLevelDesc%depFromFiner( iElem )%elem%nVals
      inv_nSourceElems = 1.0_rk / dble(nSourceElems)

      tMom = 0.0_rk
      ! Now loop over all fine source elements for this target:
      do iSourceElem = 1, nSourceElems

        ! Get the source element position
        sourceElem = tLevelDesc%depFromFiner( iElem ) &
          &                    %elem%val( iSourceElem )

        ! Get macroscopic moments
        srcMom = mus_intp_getMoments(state     = sstate,                 &
          &                          elem      = sourceElem,             &
          &                          QQ        = QQ,                     &
          &                          nScalars  = nScalars,               &
          &                          nSize     = snSize,                 &
          &                          toMoments = layout%moment%toMoments )

        tMom = srcMom + tMom

      end do

      !interpolate moments to target element
      tMom(:) = tMom(:) * inv_nSourceElems

      ! store interpolated target auxField
      tOffset = (targetElem-1)*varSys%nAuxScalars
      tAuxField(tOffset+dens_pos)   = tMom(1)
      tAuxField(tOffset+vel_pos(1)) = tMom(layout%moment%first_moments(1))*rho0Inv
      tAuxField(tOffset+vel_pos(2)) = tMom(layout%moment%first_moments(2))*rho0Inv
      tAuxField(tOffset+vel_pos(3)) = tMom(layout%moment%first_moments(3))*rho0Inv

      ! get normalized kinematic viscosity on target element
      tVisc = fluid%viscKine%dataOnLvl(targetLevel)%val(targetElem)

      ! relation between coarse and fine grid kinematic viscosity:
      ! v^s_f = 2 v^s_c
      ! calculate omega on source (fine) and target level (coarse)
      sOmegaKine = mus_calcOmegaFromVisc(2.0_rk * tVisc)
      tOmegaKine = mus_calcOmegaFromVisc(tVisc)

      ! Get scaling factors for nonequilibrium moments
      nonEqScalingFacs = fluid%nonEqScalingFacs(           &
        & omegaKine_SRC = sOmegaKine,                      &
        & omegaKine_TGT = tOmegaKine,                      &
        & omegaBulk_SRC = fluid%omegaBulkLvl(sourceLevel), &
        & omegaBulk_TGT = fluid%omegaBulkLvl(targetLevel), &
        & scaleFac      = 2.0_rk,                          &
        & QQ            = QQ                               )

      ! Convert moment to PDF
      tPDF = mus_intp_convertMomToPDF3D_incomp(                              &
        &                               moments          = tMom,             &
        &                               nonEqScalingFacs = nonEqScalingFacs, &
        &                               layout           = layout            )

      ! Now write the resulting pdf in the current direction to the target
      ! Element position
      do iDir = 1,QQ
        tstate(?IDX?( iDir, targetElem, nScalars, tnSize )) &
          & = tPDF(iDir)
      end do
    end do

  end subroutine fillMyGhostsFromFiner_avgIncomp
! **************************************************************************** !

! **************************************************************************** !
  !> [Average interpolation](../page/features/intp_methods.html) of ghostFromFiner
  !! The interpolation procedure used in this routine is:\n
  !! 1. Calculate moments from source
  !! 2. Interpolate moments on target by simple average
  !! 3. Store target auxilary field
  !! 4. Get viscosity on target element and compute source and target omega
  !! 5. Get nonEq scaling factor depeding on scheme layout and relaxation
  !! 6. Calculate Equilibrium and nonEquilibrium
  !! 7. calculate target: Eq + Scale * nonEquilibrium
  !! This routine is used by 3D acoustic average interpolation for
  !! incompressible LES model.
  !!
  !! This subroutine's interface must match the abstract interface definition
  !! [[intpRoutine]] in intp/[[mus_interpolate_header_module]].f90 in order to
  !! be callable via [[mus_interpolation_method_type:do_intp]] function pointer.
  subroutine fillMyGhostsFromFiner_avgIncompLES( method, fieldProp,         &
    & tLevelDesc, level, sState, snSize, tState, tnSize, tAuxField, layout, &
    & nTargets, targetList, physics, time, varSys, derVarPos                )
    ! -------------------------------------------------------------------- !
    class(mus_interpolation_method_type), intent(in) :: method

    !> Array of field properties (fluid or species)
    type(mus_field_prop_type), target, intent(in) :: fieldProp(:)

    !> level descriptor on target level
    type( tem_levelDesc_type ), intent(in) :: tLevelDesc

    !> my refinement level
    integer, intent(in) :: level

    !> State vector of SOURCE FLUID elements
    real(kind=rk), intent(in) :: sState(:)
    ! integer, intent(in) :: sNeigh(:)
    integer, intent(in) :: snSize

    !> State vector of TARGET GHOST elements
    real(kind=rk), intent(inout) :: tState(:)
    ! integer, intent(in) :: tNeigh(:)
    integer, intent(in) :: tnSize

    !> AuxField variable to fill on target GHOST elements
    real(kind=rk), intent(inout) :: tAuxField(:)

    !> the layout used
    type( mus_scheme_layout_type ), intent(in) :: layout

    !> List of target elements ( their position in depSource list )
    integer, intent(in) :: nTargets
    integer, intent(in) :: targetList(nTargets)

    !> physics type to convert lattice to physics SI unit and vice versa
    !! @todo: This can be replaced by scale factor by level
    type( mus_physics_type ), intent(in) :: physics

    !> time required to compute viscosity on target element barycenter
    type(tem_time_type), intent(in) :: time

    !> scheme variable system
    type( tem_varSys_type ), intent(in) :: varSys

    !> position of all derive variable in varSys for all fields
    type(mus_derVarPos_type), intent(in) :: derVarPos(:)
    ! -------------------------------------------------------------------- !
    integer :: sourceLevel    ! level of source elements
    integer :: sourceElem     ! treeId of current source element
    integer :: targetLevel    ! level of target elements
    integer :: targetElem     ! treeId of current source element
    integer :: iDir           ! current direction (discrete velocity) for loop
    integer :: iElem          ! current target element (for outer loop)
    integer :: indElem        ! element counter for indirection list
    integer :: iSourceElem    ! current source element (for inner loop)
    integer :: nSourceElems   ! number of source elements for the current target
    real(kind=rk) :: tPDF( layout%fStencil%QQ ) ! pdf to reconstruct from
    ! moments of the source elements' pdf
    real(kind=rk) :: srcMom(layout%fStencil%QQ)  ! temp moment calculation
    real(kind=rk) :: tMom(layout%fStencil%QQ)  ! target moment calculation
    integer :: QQ
    real(kind=rk) :: inv_nSourceElems
    type(mus_fluid_type), pointer :: fluid
    real(kind=rk) :: nonEqScalingFacs(layout%fStencil%QQ)
    real(kind=rk) :: sOmegaKine, tOmegaKine, sVisc, tVisc
    real(kind=rk) ::  tTurbVisc
    integer :: nScalars, tOffset
    integer :: dens_pos, vel_pos(3)
    ! --------------------------------------------------------------------------
    fluid => fieldProp(1)%fluid
    QQ = layout%fStencil%QQ
    nScalars = varSys%nScalars
    dens_pos = varSys%method%val(derVarPos(1)%density)%auxField_varPos(1)
    vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:3)

    sourceLevel = level + 1
    targetLevel = level

    ! Treat all coarse target elements
    do indElem = 1, nTargets

      iElem = targetList( indElem )

      ! Read the target element treeId
      targetElem = iElem + tLevelDesc%offset( 1, eT_ghostFromFiner)

      ! Find out how many fine source elements we have for interpolation.
      nSourceELems = tLevelDesc%depFromFiner( iElem )%elem%nVals
      inv_nSourceElems = 1.0_rk / dble(nSourceElems)

      tMom = 0.0_rk
      tTurbVisc = 0.0_rk
      ! Now loop over all fine source elements for this target:
      do iSourceElem = 1, nSourceElems

        ! Get the source element position
        sourceElem = tLevelDesc%depFromFiner( iElem ) &
          &                    %elem%val( iSourceElem )

        ! Get macroscopic moments
        srcMom = mus_intp_getMoments(state     = sstate,                 &
          &                          elem      = sourceElem,             &
          &                          QQ        = QQ,                     &
          &                          nScalars  = nScalars,               &
          &                          nSize     = snSize,                 &
          &                          toMoments = layout%moment%toMoments )

        tMom = srcMom + tMom

        !interpolate turbulent viscosity to target element
        tTurbVisc = tTurbVisc + fluid%turbulence%dataOnLvl(sourceLevel) &
          &                                      %visc(sourceElem)

      end do

      !interpolate moments to target element
      tMom(:) = tMom(:) * inv_nSourceElems

      ! store interpolated target auxField
      tOffset = (targetElem-1)*varSys%nAuxScalars
      tAuxField(tOffset+dens_pos)   = tMom(1)
      tAuxField(tOffset+vel_pos(1)) = tMom(layout%moment%first_moments(1))*rho0Inv
      tAuxField(tOffset+vel_pos(2)) = tMom(layout%moment%first_moments(2))*rho0Inv
      tAuxField(tOffset+vel_pos(3)) = tMom(layout%moment%first_moments(3))*rho0Inv

      ! interpolate turbulent viscosity on target element
      tTurbVisc = tTurbVisc * inv_nSourceElems
      ! scale interpolated turbulent viscosity to target element
      fluid%turbulence%dataOnLvl(targetLevel)%visc(targetElem) &
        & = fluid%turbulence%fac_f2c*tTurbVisc

      ! get normalized kinematic viscosity on target element
      tVisc = fluid%viscKine%dataOnLvl(targetLevel)%val(targetElem)

      ! relation between coarse and fine grid kinematic viscosity:
      ! v^s_f = 2 v^s_c
      ! total viscosity on source element
      sVisc = 2.0_rk * tVisc + tTurbVisc
      ! total viscosity on target element
      tVisc = tVisc + fluid%turbulence%fac_f2c*tTurbVisc

      ! calculate omega on source (fine) and target level (coarse)
      sOmegaKine = mus_calcOmegaFromVisc(sVisc)
      tOmegaKine = mus_calcOmegaFromVisc(tVisc)

      ! Get scaling factors for nonequilibrium moments
      nonEqScalingFacs = fluid%nonEqScalingFacs(           &
        & omegaKine_SRC = sOmegaKine,                      &
        & omegaKine_TGT = tOmegaKine,                      &
        & omegaBulk_SRC = fluid%omegaBulkLvl(sourceLevel), &
        & omegaBulk_TGT = fluid%omegaBulkLvl(targetLevel), &
        & scaleFac      = 2.0_rk,                          &
        & QQ            = QQ                               )

      ! Convert moment to PDF
      tPDF = mus_intp_convertMomToPDF3D_incomp(                              &
        &                               moments          = tMom,             &
        &                               nonEqScalingFacs = nonEqScalingFacs, &
        &                               layout           = layout            )

      ! Now write the resulting pdf in the current direction to the target
      ! Element position
      do iDir = 1,QQ
        tstate(?IDX?( iDir, targetElem, nScalars, tnSize )) &
          & = tPDF(iDir)
      end do
    end do

  end subroutine fillMyGhostsFromFiner_avgIncompLES
! **************************************************************************** !


! **************************************************************************** !
  !> Fill coarse target ghost from fine source fluid by average interpolation.
  !! 1. For each target, calculate rho, vel and fNeq of its sources.
  !! 2. Average rho and vel, then calculate fEq.
  !! 3. Average fNeq and scale.
  !! 4. set target by f = fEq + fNeq
  !! This routine is used by 2D acoustic average interpolation.
  !!
  !! This subroutine's interface must match the abstract interface definition
  !! [[intpRoutine]] in intp/[[mus_interpolate_header_module]].f90 in order to
  !! be callable via [[mus_interpolation_method_type:do_intp]] function pointer.
  subroutine fillMyGhostsFromFiner_avg2D( method, fieldProp, tLevelDesc,       &
    &                                     level, sState, snSize, tState,       &
    &                                     tnSize, tAuxField, layout, nTargets, &
    &                                     targetList, physics, time, varSys,   &
    &                                     derVarPos                            )
    ! -------------------------------------------------------------------- !
    class(mus_interpolation_method_type), intent(in) :: method

    !> Array of field properties (fluid or species)
    type(mus_field_prop_type), target, intent(in) :: fieldProp(:)

    !> level descriptor on target level
    type( tem_levelDesc_type ), intent(in) :: tLevelDesc

    !> my refinement level
    integer, intent(in) :: level

    !> State vector of SOURCE FLUID elements
    real(kind=rk), intent(in) :: sState(:)
    ! integer, intent(in) :: sNeigh(:)
    integer, intent(in) :: snSize

    !> State vector of TARGET GHOST elements
    real(kind=rk), intent(inout) :: tState(:)
    ! integer, intent(in) :: tNeigh(:)
    integer, intent(in) :: tnSize

    !> AuxField variable to fill on target GHOST elements
    real(kind=rk), intent(inout) :: tAuxField(:)

    !> the layout used
    type( mus_scheme_layout_type ), intent(in) :: layout

    !> List of target elements ( their position in depSource list )
    integer, intent(in) :: nTargets
    integer, intent(in) :: targetList(nTargets)

    !> physics type to convert lattice to physics SI unit and vice versa
    !! @todo: This can be replaced by scale factor by level
    type( mus_physics_type ), intent(in) :: physics

    !> time required to compute viscosity on target element barycenter
    type(tem_time_type), intent(in) :: time

    !> scheme variable system
    type( tem_varSys_type ), intent(in) :: varSys

    !> position of all derive variable in varSys for all fields
    type(mus_derVarPos_type), intent(in) :: derVarPos(:)
    ! -------------------------------------------------------------------- !
    integer :: sourceLevel    ! level of source elements
    integer :: sourceElem     ! treeId of current source element
    integer :: targetLevel    ! level of target elements
    integer :: targetElem     ! treeId of current source element
    integer :: iDir           ! current direction (discrete velocity) for loop
    integer :: iElem          ! current target element (for outer loop)
    integer :: indElem        ! element counter for indirection list
    integer :: iSourceElem    ! current source element (for inner loop)
    integer :: nSourceElems   ! number of source elements for the current target
    real(kind=rk) :: tPDF( layout%fStencil%QQ )
    integer :: QQ
    real(kind=rk) :: inv_nSourceElems
    real(kind=rk) :: tMom(layout%fStencil%QQ)  ! target moment calculation
    real(kind=rk) :: srcMom(layout%fStencil%QQ)  ! temp moment calculation
    type(mus_fluid_type), pointer :: fluid
    real(kind=rk) :: nonEqScalingFacs(layout%fStencil%QQ)
    real(kind=rk) :: sOmegaKine, tOmegaKine, tVisc, invRho
    integer :: nScalars, tOffset
    integer :: dens_pos, vel_pos(3)
    ! --------------------------------------------------------------------------
!write(dbgUnit(1),*) 'Entering fillCoarser average'
    fluid => fieldProp(1)%fluid
    nScalars = varSys%nScalars
    QQ = layout%fStencil%QQ
    dens_pos = varSys%method%val(derVarPos(1)%density)%auxField_varPos(1)
    vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:3)

    sourceLevel = level + 1
    targetLevel = level


    ! Treat all coarse target elements
    do indElem = 1, nTargets

      iElem = targetList( indElem )

      ! Read the target element treeId
      targetElem = iElem + tLevelDesc%offset( 1, eT_ghostFromFiner)
!write(dbgUnit(1),*) 'Target iElem ', indElem, ' ID ', tLevelDesc%total(targetElem)

      ! Find out how many fine source elements we have for interpolation.
      nSourceELems = tLevelDesc%depFromFiner( iElem )%elem%nVals
      inv_nSourceElems = 1.0_rk / dble(nSourceElems)

      ! Now loop over all fine source elements for this target:
      tMom = 0.0_rk
      do iSourceElem = 1, nSourceElems

        ! Get the source element position
        sourceElem = tLevelDesc%depFromFiner( iElem ) &
          &                    %elem%val( iSourceElem )

        ! Get macroscopic moments
        srcMom = mus_intp_getMoments(state     = sstate,                 &
          &                          elem      = sourceElem,             &
          &                          QQ        = QQ,                     &
          &                          nScalars  = nScalars,               &
          &                          nSize     = snSize,                 &
          &                          toMoments = layout%moment%toMoments )

        tMom = srcMom + tMom
      end do ! iSourceElem

      ! interpolate all moments by average
      tMom(:) = tMom(:) * inv_nSourceElems

      ! store interpolated target auxField
      invRho = 1.0_rk/tMom(1)
      tOffset = (targetElem-1)*varSys%nAuxScalars
      tAuxField(tOffset+dens_pos)   = tMom(1)
      tAuxField(tOffset+vel_pos(1)) = tMom(layout%moment%first_moments(1))*invRho
      tAuxField(tOffset+vel_pos(2)) = tMom(layout%moment%first_moments(2))*invRho
      tAuxField(tOffset+vel_pos(3)) = 0.0_rk

      ! get normalized kinematic viscosity on target element
      tVisc = fluid%viscKine%dataOnLvl(targetLevel)%val(targetElem)

      ! relation between coarse and fine grid kinematic viscosity:
      ! v^s_f = 2 v^s_c
      ! calculate omega on source and target level
      sOmegaKine = mus_calcOmegaFromVisc(2.0_rk*tVisc)
      tOmegaKine = mus_calcOmegaFromVisc(tVisc)

      ! scaling factors for nonEquilibrium moments
      nonEqScalingFacs = fluid%nonEqScalingFacs(           &
        & omegaKine_SRC = sOmegaKine,                      &
        & omegaKine_TGT = tOmegaKine,                      &
        & omegaBulk_SRC = fluid%omegaBulkLvl(sourceLevel), &
        & omegaBulk_TGT = fluid%omegaBulkLvl(targetLevel), &
        & scaleFac      = 2.0_rk,                          &
        & QQ            = QQ                               )

      ! Convert moment to PDF
      tPDF = mus_intp_convertMomToPDF2D( moments          = tMom,             &
        &                                nonEqScalingFacs = nonEqScalingFacs, &
        &                                layout           = layout            )

      do iDir = 1,QQ
        tstate(?IDX?( iDir, targetElem, nScalars, tnSize )) &
          & = tPDF(iDir)
      end do
!write(dbgUnit(1),*) 'targetRho-rho0 ', rho - rho0
    end do

  end subroutine fillMyGhostsFromFiner_avg2D
! **************************************************************************** !

! **************************************************************************** !
  !> Fill coarse target ghost from fine source fluid by average interpolation.
  !! 1. For each target, calculate rho, vel and fNeq of its sources.
  !! 2. Average rho and vel, then calculate fEq.
  !! 3. Average fNeq and scale.
  !! 4. set target by f = fEq + fNeq
  !! This routine is used by 2D incomp acoustic average interpolation.
  !!
  !! This subroutine's interface must match the abstract interface definition
  !! [[intpRoutine]] in intp/[[mus_interpolate_header_module]].f90 in order to
  !! be callable via [[mus_interpolation_method_type:do_intp]] function pointer.
  subroutine fillMyGhostsFromFiner_avg2DIncomp( method, fieldProp, tLevelDesc, &
    &                                           level, sState, snSize, tState, &
    &                                           tnSize, tAuxField, layout,     &
    &                                           nTargets, targetList, physics, &
    &                                           time, varSys, derVarPos        )
    ! -------------------------------------------------------------------- !
    class(mus_interpolation_method_type), intent(in) :: method

    !> Array of field properties (fluid or species)
    type(mus_field_prop_type), target, intent(in) :: fieldProp(:)

    !> level descriptor on target level
    type( tem_levelDesc_type ), intent(in) :: tLevelDesc

    !> my refinement level
    integer, intent(in) :: level

    !> State vector of SOURCE FLUID elements
    real(kind=rk), intent(in) :: sState(:)
    ! integer, intent(in) :: sNeigh(:)
    integer, intent(in) :: snSize

    !> State vector of TARGET GHOST elements
    real(kind=rk), intent(inout) :: tState(:)
    ! integer, intent(in) :: tNeigh(:)
    integer, intent(in) :: tnSize

    !> AuxField variable to fill on target GHOST elements
    real(kind=rk), intent(inout) :: tAuxField(:)

    !> the layout used
    type( mus_scheme_layout_type ), intent(in) :: layout

    !> List of target elements ( their position in depSource list )
    integer, intent(in) :: nTargets
    integer, intent(in) :: targetList(nTargets)

    !> physics type to convert lattice to physics SI unit and vice versa
    !! @todo: This can be replaced by scale factor by level
    type( mus_physics_type ), intent(in) :: physics

    !> time required to compute viscosity on target element barycenter
    type(tem_time_type), intent(in) :: time

    !> scheme variable system
    type( tem_varSys_type ), intent(in) :: varSys

    !> position of all derive variable in varSys for all fields
    type(mus_derVarPos_type), intent(in) :: derVarPos(:)
    ! -------------------------------------------------------------------- !
    integer :: sourceLevel    ! level of source elements
    integer :: sourceElem     ! treeId of current source element
    integer :: targetLevel    ! level of target elements
    integer :: targetElem     ! treeId of current source element
    integer :: iDir           ! current direction (discrete velocity) for loop
    integer :: iElem          ! current target element (for outer loop)
    integer :: indElem        ! element counter for indirection list
    integer :: iSourceElem    ! current source element (for inner loop)
    integer :: nSourceElems   ! number of source elements for the current target
    ! current non-equilibrium part
    real(kind=rk) :: tPDF( layout%fStencil%QQ )
    integer :: QQ
    real(kind=rk) :: inv_nSourceElems
    real(kind=rk) :: tMom(layout%fStencil%QQ)  ! target moment calculation
    real(kind=rk) :: srcMom(layout%fStencil%QQ)  ! temp moment calculation
    type(mus_fluid_type), pointer :: fluid
    real(kind=rk) :: nonEqScalingFacs(layout%fStencil%QQ)
    real(kind=rk) :: sOmegaKine, tOmegaKine, tVisc
    integer :: nScalars, tOffset
    integer :: dens_pos, vel_pos(3)
    ! --------------------------------------------------------------------------
!write(dbgUnit(1),*) 'Entering fillCoarser average'
    fluid => fieldProp(1)%fluid
    nScalars = varSys%nScalars
    QQ = layout%fStencil%QQ
    dens_pos = varSys%method%val(derVarPos(1)%density)%auxField_varPos(1)
    vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:3)

    sourceLevel = level + 1
    targetLevel = level

    ! Treat all coarse target elements
    do indElem = 1, nTargets

      iElem = targetList( indElem )

      ! Read the target element treeId
      targetElem = iElem + tLevelDesc%offset( 1, eT_ghostFromFiner)
!write(dbgUnit(1),*) 'Target iElem ', indElem, ' ID ', tLevelDesc%total(targetElem)

      ! Find out how many fine source elements we have for interpolation.
      nSourceELems = tLevelDesc%depFromFiner( iElem )%elem%nVals
      inv_nSourceElems = 1.0_rk / dble(nSourceElems)

      ! Now loop over all fine source elements for this target:
      tMom = 0.0_rk
      do iSourceElem = 1, nSourceElems

        ! Get the source element position
        sourceElem = tLevelDesc%depFromFiner( iElem ) &
          &                    %elem%val( iSourceElem )

        ! Get macroscopic moments
        srcMom = mus_intp_getMoments(state     = sstate,                 &
          &                          elem      = sourceElem,             &
          &                          QQ        = QQ,                     &
          &                          nScalars  = nScalars,               &
          &                          nSize     = snSize,                 &
          &                          toMoments = layout%moment%toMoments )

        tMom = srcMom + tMom
      end do ! iSourceElem

      ! interpolate all moments by average
      tMom(:) = tMom(:) * inv_nSourceElems
!write(dbgUnit(1),*) 'before scaling tMom ', tMom

      ! store interpolated target auxField
      tOffset = (targetElem-1)*varSys%nAuxScalars
      tAuxField(tOffset+dens_pos)   = tMom(1)
      tAuxField(tOffset+vel_pos(1)) = tMom(layout%moment%first_moments(1))*rho0Inv
      tAuxField(tOffset+vel_pos(2)) = tMom(layout%moment%first_moments(2))*rho0Inv
      tAuxField(tOffset+vel_pos(3)) = 0.0_rk

      ! get normalized kinematic viscosity on target element
      tVisc = fluid%viscKine%dataOnLvl(targetLevel)%val(targetElem)

      ! relation between coarse and fine grid kinematic viscosity:
      ! v^s_f = 2 v^s_c
      ! calculate omega on source (fine) and target level (coarse)
      sOmegaKine = mus_calcOmegaFromVisc(2.0_rk * tVisc)
      tOmegaKine = mus_calcOmegaFromVisc(tVisc)

      ! Get scaling factors for nonequilibrium moments
      nonEqScalingFacs = fluid%nonEqScalingFacs(           &
        & omegaKine_SRC = sOmegaKine,                      &
        & omegaKine_TGT = tOmegaKine,                      &
        & omegaBulk_SRC = fluid%omegaBulkLvl(sourceLevel), &
        & omegaBulk_TGT = fluid%omegaBulkLvl(targetLevel), &
        & scaleFac      = 2.0_rk,                          &
        & QQ            = QQ                               )

      ! Convert moment to PDF
      tPDF = mus_intp_convertMomToPDF2D_incomp(                              &
        &                               moments          = tMom,             &
        &                               nonEqScalingFacs = nonEqScalingFacs, &
        &                               layout           = layout            )

      do iDir = 1,QQ
        tstate(?IDX?( iDir, targetElem, nScalars, tnSize )) &
          & = tPDF(iDir)
      end do
!write(dbgUnit(1),*) 'tPDF ', tPDF
    end do

  end subroutine fillMyGhostsFromFiner_avg2DIncomp
! **************************************************************************** !

! **************************************************************************** !
  !> Interpolate auxiliary field from coarse source to fine target
  !!
  !! This subroutine's interface must match the abstract interface definition
  !! [[intpRoutine_arbitraryVal]] in intp/[[mus_interpolate_header_module]].f90
  !! in order to be callable via
  !! [[mus_interpolation_method_type:do_intpArbiVal]] function pointer.
  subroutine fillArbiFinerGhostsFromMe_weighAvg( method, tLevelDesc, level,    &
    &                                            stencil, sVal, snSize, tVal,  &
    &                                            tnSize, nTargets, targetList, &
    &                                            nScalars                      )
    ! ------------------------------------------------------------------ !
    class(mus_interpolation_method_type), intent(in) :: method

    !> my refinement level
    integer, intent(in) :: level

    !> stencil header
    type(tem_stencilHeader_type), intent(in) :: stencil

    !> State vector of SOURCE FLUID elements
    real(kind=rk), intent(in) :: sVal(:)
    integer, intent(in) :: snSize

    !> State vector of TARGET GHOST elements
    real(kind=rk), intent(inout) :: tVal(:)
    integer, intent(in) :: tnSize

    !> level descriptor on target level
    type( tem_levelDesc_type ), intent(in) :: tLevelDesc

    !> List of target elements ( their position in depSource list )
    integer, intent(in) :: nTargets
    integer, intent(in) :: targetList(nTargets)

    !> number of scalars to interpolate
    integer, intent(in) :: nScalars
    ! ------------------------------------------------------------------ !
    ! --------------------------------------------------------------------------
    integer :: sourceLevel    ! level of source elements
    integer :: sourceElem     ! treeId of current source element
    integer :: targetLevel    ! level of target elements
    integer :: targetElem     ! treeId of current source element
    integer :: iElem          ! current target element (for outer loop)
    integer :: indElem        ! element counter for indirection list
    integer :: iSourceElem    ! current source element (for inner loop)
    integer :: nSourceElems   ! number of source elements for the current target
    real(kind=rk) :: weight( stencil%QQ )
    real(kind=rk) :: sArbi(nScalars, stencil%QQ)  ! temp source ArbiField
    integer :: iVar
    ! --------------------------------------------------------------------------
    sourceLevel = level
    targetLevel = level + 1

    ! Treat all coarse target elements
    do indElem = 1, nTargets

      iElem = targetList( indElem )

      ! Read the target element
      targetElem = iElem + tLevelDesc%offset( 1, eT_ghostFromCoarser)

      ! Get how many fine source elements we have for interpolation.
      nSourceElems = tLevelDesc%depFromCoarser( iElem )%elem%nVals

      ! Read the pre-calculated weights (like in Average intp.)
      weight(1:nSourceElems) = tLevelDesc%depFromCoarser( iElem ) &
        &                                %weight(1:nSourceElems)

      ! Now loop over all fine source elements for this target:
      do iSourceElem = 1, nSourceElems

        ! Get the source element
        sourceElem = tLevelDesc%depFromCoarser( iElem ) &
          &                    %elem%val( iSourceElem )

        ! Get souce auxilary variables
        sArbi(:, iSourceElem) = sVal( (sourceElem-1)*nScalars+1 &
          &                           : sourceElem*nScalars       )

      end do  ! iSourceElem

      ! interpolate all auxiliary variables by average
      do iVar = 1, nScalars
        tVal((targetElem-1)*nScalars+iVar)                    &
          & = sum( weight(1:nSourceElems)*sArbi(iVar,1:nSourceElems) )
      end do

    enddo

  end subroutine fillArbiFinerGhostsFromMe_weighAvg
! **************************************************************************** !

! **************************************************************************** !
  !> [Linear interpolation](../page/features/intp_methods.html) of ghostFromFiner
  !!
  !! This subroutine's interface must match the abstract interface definition
  !! [[intpRoutine]] in intp/[[mus_interpolate_header_module]].f90 in order to
  !! be callable via [[mus_interpolation_method_type:do_intp]] function pointer.
  subroutine fillFinerGhostsFromMe_weighAvg( method, fieldProp, tLevelDesc, &
    &                                        level, sState, snSize, tState, &
    &                                        tnSize, tAuxField, layout,     &
    &                                        nTargets, targetList, physics, &
    &                                        time, varSys, derVarPos        )
    ! -------------------------------------------------------------------- !
    class(mus_interpolation_method_type), intent(in) :: method

    !> Array of field properties (fluid or species)
    type(mus_field_prop_type), target, intent(in) :: fieldProp(:)

    !> level descriptor on target level
    type( tem_levelDesc_type ), intent(in) :: tLevelDesc

    !> my refinement level
    integer, intent(in) :: level

    !> State vector of SOURCE FLUID elements
    real(kind=rk), intent(in) :: sState(:)
    ! integer, intent(in) :: sNeigh(:)
    integer, intent(in) :: snSize

    !> State vector of TARGET GHOST elements
    real(kind=rk), intent(inout) :: tState(:)
    ! integer, intent(in) :: tNeigh(:)
    integer, intent(in) :: tnSize

    !> AuxField variable to fill on target GHOST elements
    real(kind=rk), intent(inout) :: tAuxField(:)

    !> the layout used
    type( mus_scheme_layout_type ), intent(in) :: layout

    !> List of target elements ( their position in depSource list )
    integer, intent(in) :: nTargets
    integer, intent(in) :: targetList(nTargets)

    !> physics type to convert lattice to physics SI unit and vice versa
    !! @todo: This can be replaced by scale factor by level
    type( mus_physics_type ), intent(in) :: physics

    !> time required to compute viscosity on target element barycenter
    type(tem_time_type), intent(in) :: time

    !> scheme variable system
    type( tem_varSys_type ), intent(in) :: varSys

    !> position of all derive variable in varSys for all fields
    type(mus_derVarPos_type), intent(in) :: derVarPos(:)
    ! -------------------------------------------------------------------- !
    integer :: sourceLevel    ! level of source elements
    integer :: sourceElem     ! treeId of current source element
    integer :: targetLevel    ! level of target elements
    integer :: targetElem     ! treeId of current source element
    integer :: iDir           ! current direction (discrete velocity) for loop
    integer :: iElem          ! current target element (for outer loop)
    integer :: indElem        ! element counter for indirection list
    integer :: iSourceElem    ! current source element (for inner loop)
    integer :: nSourceElems   ! number of source elements for the current target
    real(kind=rk) :: tPDF( layout%fStencil%QQ )
    real(kind=rk) :: weight( layout%fStencil%QQ )
    real(kind=rk) :: srcMom( layout%fStencil%QQ, layout%fStencil%QQ )
    real(kind=rk) :: tMom(layout%fStencil%QQ)  ! target moment calculation
    real(kind=rk) :: nonEqScalingFacs(layout%fStencil%QQ)
    integer :: iMom, QQ
    type(mus_fluid_type), pointer :: fluid
    real(kind=rk) :: sOmegaKine, tOmegaKine, tVisc, invRho
    integer :: nScalars, tOffset
    integer :: dens_pos, vel_pos(3)
    ! --------------------------------------------------------------------------
    fluid => fieldProp(1)%fluid
    nScalars = varSys%nScalars
    QQ = layout%fStencil%QQ
    dens_pos = varSys%method%val(derVarPos(1)%density)%auxField_varPos(1)
    vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:3)

    sourceLevel = level
    targetLevel = level + 1

    ! Treat all fine target elements:
    do indElem = 1, nTargets

      iElem = targetList( indElem )

      ! Read the target element treeId
      targetElem = iElem + tLevelDesc%offset( 1, eT_ghostFromCoarser)

      ! Find out how many fine source elements we have for interpolation.
      nSourceElems = tLevelDesc%depFromCoarser( iElem )%elem%nVals

      ! Read the pre-calculated weights (like in Average intp.)
      weight(1:nSourceElems) = tLevelDesc%depFromCoarser( iElem ) &
        &                                %weight(1:nSourceElems)

      ! Now loop over all fine source elements for this target:
      do iSourceElem = 1, nSourceElems

        ! Get the source element position
        sourceElem = tLevelDesc%depFromCoarser( iElem ) &
          &                    %elem%val( iSourceElem )

        ! Get macroscopic moments
        srcMom(:, iSourceElem) = mus_intp_getMoments(                         &
          &                               state     = sstate,                 &
          &                               elem      = sourceElem,             &
          &                               QQ        = QQ,                     &
          &                               nScalars  = nScalars,               &
          &                               nSize     = snSize,                 &
          &                               toMoments = layout%moment%toMoments )

      end do ! iSourceElem

      ! interpolate all moments by weighted average
      do iMom = 1,QQ
        tMom(iMom) = sum(weight(1:nSourceElems)*srcMom(iMom,1:nSourceElems))
      end do

      ! store interpolated target auxField
      invRho = 1.0_rk/tMom(1)
      tOffset = (targetElem-1)*varSys%nAuxScalars
      tAuxField(tOffset+dens_pos)   = tMom(1)
      tAuxField(tOffset+vel_pos(1)) = tMom(layout%moment%first_moments(1))*invRho
      tAuxField(tOffset+vel_pos(2)) = tMom(layout%moment%first_moments(2))*invRho
      tAuxField(tOffset+vel_pos(3)) = tMom(layout%moment%first_moments(3))*invRho

      ! get normalized kinematic viscosity on target element
      tVisc = fluid%viscKine%dataOnLvl(targetLevel)%val(targetElem)
      ! relation between coarse and fine grid kinematic viscosity:
      ! v^s_f = 2 v^s_c
      ! calculate omega on source and target level
      sOmegaKine = mus_calcOmegaFromVisc(0.5_rk*tVisc)
      tOmegaKine = mus_calcOmegaFromVisc(tVisc)

      ! Get scaling factors for nonequilibrium moments
      nonEqScalingFacs = fluid%nonEqScalingFacs(           &
        & omegaKine_SRC = sOmegaKine,                      &
        & omegaKine_TGT = tOmegaKine,                      &
        & omegaBulk_SRC = fluid%omegaBulkLvl(sourceLevel), &
        & omegaBulk_TGT = fluid%omegaBulkLvl(targetLevel), &
        & scaleFac      = 0.5_rk, QQ = QQ                  )

      ! Convert moment to PDF
      tPDF = mus_intp_convertMomToPDF3D( moments          = tMom,             &
        &                                nonEqScalingFacs = nonEqScalingFacs, &
        &                                layout           = layout            )

      ! Now write the resulting pdf to the target element
      do iDir = 1, QQ
        tstate(?IDX?(iDir, targetElem, nScalars, tnSize ) ) &
          & = tPDF( iDir )
      end do
    enddo ! element loop

  end subroutine fillFinerGhostsFromMe_weighAvg
! **************************************************************************** !


! **************************************************************************** !
  !> [Linear interpolation](../page/features/intp_methods.html) of ghostFromFiner
  !!
  !! This subroutine's interface must match the abstract interface definition
  !! [[intpRoutine]] in intp/[[mus_interpolate_header_module]].f90 in order to
  !! be callable via [[mus_interpolation_method_type:do_intp]] function pointer.
  subroutine fillFinerGhostsFromMe_weighAvgLES( method, fieldProp, tLevelDesc, &
    &                                           level, sState, snSize, tState, &
    &                                           tnSize, tAuxField, layout,     &
    &                                           nTargets, targetList, physics, &
    &                                           time, varSys, derVarPos        )
    ! -------------------------------------------------------------------- !
    class(mus_interpolation_method_type), intent(in) :: method

    !> Array of field properties (fluid or species)
    type(mus_field_prop_type), target, intent(in) :: fieldProp(:)

    !> level descriptor on target level
    type( tem_levelDesc_type ), intent(in) :: tLevelDesc

    !> my refinement level
    integer, intent(in) :: level

    !> State vector of SOURCE FLUID elements
    real(kind=rk), intent(in) :: sState(:)
    ! integer, intent(in) :: sNeigh(:)
    integer, intent(in) :: snSize

    !> State vector of TARGET GHOST elements
    real(kind=rk), intent(inout) :: tState(:)
    ! integer, intent(in) :: tNeigh(:)
    integer, intent(in) :: tnSize

    !> AuxField variable to fill on target GHOST elements
    real(kind=rk), intent(inout) :: tAuxField(:)

    !> the layout used
    type( mus_scheme_layout_type ), intent(in) :: layout

    !> List of target elements ( their position in depSource list )
    integer, intent(in) :: nTargets
    integer, intent(in) :: targetList(nTargets)

    !> physics type to convert lattice to physics SI unit and vice versa
    !! @todo: This can be replaced by scale factor by level
    type( mus_physics_type ), intent(in) :: physics

    !> time required to compute viscosity on target element barycenter
    type(tem_time_type), intent(in) :: time

    !> scheme variable system
    type( tem_varSys_type ), intent(in) :: varSys

    !> position of all derive variable in varSys for all fields
    type(mus_derVarPos_type), intent(in) :: derVarPos(:)
    ! -------------------------------------------------------------------- !
    integer :: sourceLevel    ! level of source elements
    integer :: sourceElem     ! treeId of current source element
    integer :: targetLevel    ! level of target elements
    integer :: targetElem     ! treeId of current source element
    integer :: iDir           ! current direction (discrete velocity) for loop
    integer :: iElem          ! current target element (for outer loop)
    integer :: indElem        ! element counter for indirection list
    integer :: iSourceElem    ! current source element (for inner loop)
    integer :: nSourceElems   ! number of source elements for the current target
    real(kind=rk) :: tPDF( layout%fStencil%QQ )
    real(kind=rk) :: weight( layout%fStencil%QQ )
    real(kind=rk) :: srcMom( layout%fStencil%QQ, layout%fStencil%QQ )
    real(kind=rk) :: tMom(layout%fStencil%QQ)  ! target moment calculation
    integer :: iMom, QQ
    type(mus_fluid_type), pointer :: fluid
    real(kind=rk) :: nonEqScalingFacs(layout%fStencil%QQ)
    real(kind=rk) :: sOmegaKine, tOmegaKine, sVisc, tVisc, invRho
    real(kind=rk) :: sTurbVisc(layout%fStencil%QQ), tTurbVisc
    integer :: nScalars, tOffset
    integer :: dens_pos, vel_pos(3)
    ! --------------------------------------------------------------------------
    fluid => fieldProp(1)%fluid
    nScalars = varSys%nScalars
    QQ = layout%fStencil%QQ
    dens_pos = varSys%method%val(derVarPos(1)%density)%auxField_varPos(1)
    vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:3)

    sourceLevel = level
    targetLevel = level + 1

    ! Treat all fine target elements:
    do indElem = 1, nTargets

      iElem = targetList( indElem )

      ! Read the target element treeId
      targetElem = iElem + tLevelDesc%offset( 1, eT_ghostFromCoarser)

      ! Find out how many fine source elements we have for interpolation.
      nSourceElems = tLevelDesc%depFromCoarser( iElem )%elem%nVals

      ! Read the pre-calculated weights (like in Average intp.)
      weight(1:nSourceElems) = tLevelDesc%depFromCoarser( iElem ) &
        &                                %weight(1:nSourceElems)

      ! Now loop over all fine source elements for this target:
      do iSourceElem = 1, nSourceElems

        ! Get the source element position
        sourceElem = tLevelDesc%depFromCoarser( iElem ) &
          &                    %elem%val( iSourceElem )

        ! Get macroscopic moments
        srcMom(:, iSourceElem) = mus_intp_getMoments(                         &
          &                               state     = sstate,                 &
          &                               elem      = sourceElem,             &
          &                               QQ        = QQ,                     &
          &                               nScalars  = nScalars,               &
          &                               nSize     = snSize,                 &
          &                               toMoments = layout%moment%toMoments )

        ! get turbulent viscosity
        sTurbVisc(iSourceElem) = fluid%turbulence%dataOnLvl(sourceLevel) &
          &                                       %visc(sourceElem)
      end do ! iSourceElem

      ! interpolate all moments by weighted average
      do iMom = 1,QQ
        tMom(iMom) = sum(weight(1:nSourceElems)*srcMom(iMom,1:nSourceElems))
      end do

      ! store interpolated target auxField
      invRho = 1.0_rk/tMom(1)
      tOffset = (targetElem-1)*varSys%nAuxScalars
      tAuxField(tOffset+dens_pos)   = tMom(1)
      tAuxField(tOffset+vel_pos(1)) = tMom(layout%moment%first_moments(1))*invRho
      tAuxField(tOffset+vel_pos(2)) = tMom(layout%moment%first_moments(2))*invRho
      tAuxField(tOffset+vel_pos(3)) = tMom(layout%moment%first_moments(3))*invRho

      ! interpolate turbulent viscosity
      tTurbVisc = sum(weight(1:nSourceElems)*sTurbVisc(1:nSourceElems))

      ! scale interpolated turbulent viscosity to target element
      fluid%turbulence%dataOnLvl(targetLevel)%visc(targetElem) &
        & = fluid%turbulence%fac_c2f*tTurbVisc

      ! get normalized kinematic viscosity on target element
      tVisc = fluid%viscKine%dataOnLvl(targetLevel)%val(targetElem)

      ! relation between coarse and fine grid kinematic viscosity:
      ! v^s_f = 2 v^s_c
      ! total viscosity on source element
      sVisc = 0.5_rk * tVisc + tTurbVisc
      ! total viscosity on target element
      tVisc = tVisc + fluid%turbulence%fac_c2f*tTurbVisc

      ! calculate omega on source and target level
      sOmegaKine = mus_calcOmegaFromVisc(sVisc)
      tOmegaKine = mus_calcOmegaFromVisc(tVisc)

      ! Get scaling factors for nonequilibrium moments
      nonEqScalingFacs = fluid%nonEqScalingFacs(           &
        & omegaKine_SRC = sOmegaKine,                      &
        & omegaKine_TGT = tOmegaKine,                      &
        & omegaBulk_SRC = fluid%omegaBulkLvl(sourceLevel), &
        & omegaBulk_TGT = fluid%omegaBulkLvl(targetLevel), &
        & scaleFac      = 0.5_rk, QQ = QQ                  )

      ! Convert moment to PDF
      tPDF = mus_intp_convertMomToPDF3D( moments          = tMom,             &
        &                                nonEqScalingFacs = nonEqScalingFacs, &
        &                                layout           = layout            )

      ! Now write the resulting pdf to the target element
      do iDir = 1, QQ
        tstate(?IDX?(iDir, targetElem, nScalars, tnSize ) ) &
          & = tPDF( iDir )
      end do
    enddo ! element loop

  end subroutine fillFinerGhostsFromMe_weighAvgLES
! **************************************************************************** !


! **************************************************************************** !
  !> Fill finer ghost from coarser fluid:
  !! 1. Compute moments for all source elements, save in momBuf
  !! 2. For each target, interpolate all moments
  !! 3. Store target auxilary field
  !! 4. Get viscosity on target element and compute source and target omega
  !! 5. Get nonEq scaling factor depeding on scheme layout and relaxation
  !! 6. Convert moments to PDF
  !! 7. Calculate target: Eq + Scale * nonEquilibrium
  !! This routine is used by 2D and 3D acoustic weighted average interpolation
  !! for incompressible model.
  !!
  !! This subroutine's interface must match the abstract interface definition
  !! [[intpRoutine]] in intp/[[mus_interpolate_header_module]].f90 in order to
  !! be callable via [[mus_interpolation_method_type:do_intp]] function pointer.
  subroutine fillFinerGhostsFromMe_weighAvgIncomp( method, fieldProp,       &
    & tLevelDesc, level, sState, snSize, tState, tnSize, tAuxField, layout, &
    & nTargets, targetList, physics, time, varSys, derVarPos                )
    ! -------------------------------------------------------------------- !
    class(mus_interpolation_method_type), intent(in) :: method

    !> Array of field properties (fluid or species)
    type(mus_field_prop_type), target, intent(in) :: fieldProp(:)

    !> level descriptor on target level
    type( tem_levelDesc_type ), intent(in) :: tLevelDesc

    !> my refinement level
    integer, intent(in) :: level

    !> State vector of SOURCE FLUID elements
    real(kind=rk), intent(in) :: sState(:)
    ! integer, intent(in) :: sNeigh(:)
    integer, intent(in) :: snSize

    !> State vector of TARGET GHOST elements
    real(kind=rk), intent(inout) :: tState(:)
    ! integer, intent(in) :: tNeigh(:)
    integer, intent(in) :: tnSize

    !> AuxField variable to fill on target GHOST elements
    real(kind=rk), intent(inout) :: tAuxField(:)

    !> the layout used
    type( mus_scheme_layout_type ), intent(in) :: layout

    !> List of target elements ( their position in depSource list )
    integer, intent(in) :: nTargets
    integer, intent(in) :: targetList(nTargets)

    !> physics type to convert lattice to physics SI unit and vice versa
    !! @todo: This can be replaced by scale factor by level
    type( mus_physics_type ), intent(in) :: physics

    !> time required to compute viscosity on target element barycenter
    type(tem_time_type), intent(in) :: time

    !> scheme variable system
    type( tem_varSys_type ), intent(in) :: varSys

    !> position of all derive variable in varSys for all fields
    type(mus_derVarPos_type), intent(in) :: derVarPos(:)
    ! -------------------------------------------------------------------- !
    integer :: sourceLevel    ! level of source elements
    integer :: sourceElem     ! treeId of current source element
    integer :: targetLevel    ! level of target elements
    integer :: targetElem     ! treeId of current source element
    integer :: iDir           ! current direction (discrete velocity) for loop
    integer :: iElem          ! current target element (for outer loop)
    integer :: indElem        ! element counter for indirection list
    integer :: iSourceElem    ! current source element (for inner loop)
    integer :: nSourceElems   ! number of source elements for the current target
    real(kind=rk) :: tPDF( layout%fStencil%QQ )
    real(kind=rk) :: weight( layout%fStencil%QQ )
    ! moments of the source elements' pdf
    real(kind=rk) :: srcMom( layout%fStencil%QQ, layout%fStencil%QQ )
    real(kind=rk) :: tMom(layout%fStencil%QQ)  ! target moment calculation
    integer :: iMom, QQ
    type(mus_fluid_type), pointer :: fluid
    real(kind=rk) :: nonEqScalingFacs(layout%fStencil%QQ)
    real(kind=rk) :: sOmegaKine, tOmegaKine, tVisc
    integer :: nScalars, tOffset
    integer :: dens_pos, vel_pos(3)
    ! --------------------------------------------------------------------------
    fluid => fieldProp(1)%fluid
    nScalars = varSys%nScalars
    QQ = layout%fStencil%QQ
    dens_pos = varSys%method%val(derVarPos(1)%density)%auxField_varPos(1)
    vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:3)

    sourceLevel = level
    targetLevel = level + 1

    ! Treat all fine target elements:
    do indElem = 1, nTargets

      iElem = targetList( indElem )

      ! Read the target element treeId
      targetElem = iElem + tLevelDesc%offset( 1, eT_ghostFromCoarser)

      ! Find out how many fine source elements we have for interpolation.
      ! Usually 8, but can differ at corners, obstacles, boundaries...
      nSourceElems = tLevelDesc%depFromCoarser( iElem )%elem%nVals

      ! Read the pre-calculated weights (like in Average intp.)
      weight(1:nSourceElems) = &
        &            tLevelDesc%depFromCoarser( iElem )%weight(1:nSourceElems)

      ! Now loop over all fine source elements for this target:
      do iSourceElem = 1, nSourceElems

        ! Get the source element position
        sourceElem = tLevelDesc%depFromCoarser( iElem ) &
          &                    %elem%val( iSourceElem )

        ! Get macroscopic moments
        srcMom(:, iSourceElem) = mus_intp_getMoments(                         &
          &                               state     = sstate,                 &
          &                               elem      = sourceElem,             &
          &                               QQ        = QQ,                     &
          &                               nScalars  = nScalars,               &
          &                               nSize     = snSize,                 &
          &                               toMoments = layout%moment%toMoments )

      end do ! iSourceElem

      ! interpolate moments (rho, momentum and shear stress) to tatrget element
      do iMom = 1,QQ
        tMom(iMom) = sum(weight(1:nSourceElems)*srcMom(iMom,1:nSourceElems))
      end do

      ! store interpolated target auxField
      tOffset = (targetElem-1)*varSys%nAuxScalars
      tAuxField(tOffset+dens_pos)   = tMom(1)
      tAuxField(tOffset+vel_pos(1)) = tMom(layout%moment%first_moments(1))*rho0Inv
      tAuxField(tOffset+vel_pos(2)) = tMom(layout%moment%first_moments(2))*rho0Inv
      tAuxField(tOffset+vel_pos(3)) = tMom(layout%moment%first_moments(3))*rho0Inv

      ! get normalized kinematic viscosity on target element
      tVisc = fluid%viscKine%dataOnLvl(targetLevel)%val(targetElem)

      ! relation between coarse and fine grid kinematic viscosity:
      ! v^s_f = 2 v^s_c
      ! calculate omega on source and target level
      sOmegaKine = mus_calcOmegaFromVisc(0.5_rk * tVisc)
      tOmegaKine = mus_calcOmegaFromVisc(tVisc)

      ! Get scaling factors for nonequilibrium moments
      nonEqScalingFacs = fluid%nonEqScalingFacs(           &
        & omegaKine_SRC = sOmegaKine,                      &
        & omegaKine_TGT = tOmegaKine,                      &
        & omegaBulk_SRC = fluid%omegaBulkLvl(sourceLevel), &
        & omegaBulk_TGT = fluid%omegaBulkLvl(targetLevel), &
        & scaleFac      = 0.5_rk, QQ = QQ                  )

      ! Convert moment to PDF
      tPDF = mus_intp_convertMomToPDF3D_incomp(                               &
        &                                moments          = tMom,             &
        &                                nonEqScalingFacs = nonEqScalingFacs, &
        &                                layout           = layout            )

      ! Now write the resulting pdf in the current direction to the target
      ! element position
      do iDir = 1,QQ
        tstate(?IDX?(iDir, targetElem, nScalars, tnSize )) &
          &  = tPDF( iDir )
      end do

    end do

  end subroutine fillFinerGhostsFromMe_weighAvgIncomp
! **************************************************************************** !


! **************************************************************************** !
  !> Fill finer ghost from coarser fluid:
  !! 1. Compute moments for all source elements, save in momBuf
  !! 2. For each target, interpolate all moments
  !! 3. Store target auxilary field
  !! 4. Compute viscosity on target element and compute source and target omega
  !! 5. Get nonEq scaling factor depeding on scheme layout and relaxation
  !! 6. Convert moments to PDF
  !! 7. Calculate target: Eq + Scale * nonEquilibrium
  !! This routine is used by 2D and 3D acoustic weighted average interpolation
  !! for incompressible LES model.
  !!
  !! This subroutine's interface must match the abstract interface definition
  !! [[intpRoutine]] in intp/[[mus_interpolate_header_module]].f90 in order to
  !! be callable via [[mus_interpolation_method_type:do_intp]] function pointer.
  subroutine fillFinerGhostsFromMe_weighAvgIncompLES( method, fieldProp,    &
    & tLevelDesc, level, sState, snSize, tState, tnSize, tAuxField, layout, &
    & nTargets, targetList, physics, time, varSys, derVarPos                )
    ! -------------------------------------------------------------------- !
    class(mus_interpolation_method_type), intent(in) :: method

    !> Array of field properties (fluid or species)
    type(mus_field_prop_type), target, intent(in) :: fieldProp(:)

    !> level descriptor on target level
    type( tem_levelDesc_type ), intent(in) :: tLevelDesc

    !> my refinement level
    integer, intent(in) :: level

    !> State vector of SOURCE FLUID elements
    real(kind=rk), intent(in) :: sState(:)
    ! integer, intent(in) :: sNeigh(:)
    integer, intent(in) :: snSize

    !> State vector of TARGET GHOST elements
    real(kind=rk), intent(inout) :: tState(:)
    ! integer, intent(in) :: tNeigh(:)
    integer, intent(in) :: tnSize

    !> AuxField variable to fill on target GHOST elements
    real(kind=rk), intent(inout) :: tAuxField(:)

    !> the layout used
    type( mus_scheme_layout_type ), intent(in) :: layout

    !> List of target elements ( their position in depSource list )
    integer, intent(in) :: nTargets
    integer, intent(in) :: targetList(nTargets)

    !> physics type to convert lattice to physics SI unit and vice versa
    !! @todo: This can be replaced by scale factor by level
    type( mus_physics_type ), intent(in) :: physics

    !> time required to compute viscosity on target element barycenter
    type(tem_time_type), intent(in) :: time

    !> scheme variable system
    type( tem_varSys_type ), intent(in) :: varSys

    !> position of all derive variable in varSys for all fields
    type(mus_derVarPos_type), intent(in) :: derVarPos(:)
    ! -------------------------------------------------------------------- !
    integer :: sourceLevel    ! level of source elements
    integer :: sourceElem     ! treeId of current source element
    integer :: targetLevel    ! level of target elements
    integer :: targetElem     ! treeId of current source element
    integer :: iDir           ! current direction (discrete velocity) for loop
    integer :: iElem          ! current target element (for outer loop)
    integer :: indElem        ! element counter for indirection list
    integer :: iSourceElem    ! current source element (for inner loop)
    integer :: nSourceElems   ! number of source elements for the current target
    real(kind=rk) :: tPDF( layout%fStencil%QQ )
    real(kind=rk) :: weight( layout%fStencil%QQ )
    ! moments of the source elements' pdf
    real(kind=rk) :: srcMom( layout%fStencil%QQ, layout%fStencil%QQ )
    real(kind=rk) :: tMom(layout%fStencil%QQ)  ! target moment calculation
    integer :: iMom, QQ
    type(mus_fluid_type), pointer :: fluid
    real(kind=rk) :: nonEqScalingFacs(layout%fStencil%QQ)
    real(kind=rk) :: sOmegaKine, tOmegaKine, sVisc, tVisc
    real(kind=rk) :: sTurbVisc(layout%fStencil%QQ), tTurbVisc
    integer :: nScalars, tOffset
    integer :: dens_pos, vel_pos(3)
    ! --------------------------------------------------------------------------
    fluid => fieldProp(1)%fluid
    nScalars = varSys%nScalars
    QQ = layout%fStencil%QQ
    dens_pos = varSys%method%val(derVarPos(1)%density)%auxField_varPos(1)
    vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:3)

    sourceLevel = level
    targetLevel = level + 1

    ! Treat all fine target elements:
    do indElem = 1, nTargets

      iElem = targetList( indElem )

      ! Read the target element treeId
      targetElem = iElem + tLevelDesc%offset( 1, eT_ghostFromCoarser)

      ! Find out how many fine source elements we have for interpolation.
      ! Usually 8, but can differ at corners, obstacles, boundaries...
      nSourceElems = tLevelDesc%depFromCoarser( iElem )%elem%nVals

      ! Read the pre-calculated weights (like in Average intp.)
      weight(1:nSourceElems) = &
        &            tLevelDesc%depFromCoarser( iElem )%weight(1:nSourceElems)

      ! Now loop over all fine source elements for this target:
      do iSourceElem = 1, nSourceElems

        ! Get the source element position
        sourceElem = tLevelDesc%depFromCoarser( iElem ) &
          &                    %elem%val( iSourceElem )

        ! Get macroscopic moments
        srcMom(:, iSourceElem) = mus_intp_getMoments(                         &
          &                               state     = sstate,                 &
          &                               elem      = sourceElem,             &
          &                               QQ        = QQ,                     &
          &                               nScalars  = nScalars,               &
          &                               nSize     = snSize,                 &
          &                               toMoments = layout%moment%toMoments )

        ! get turbulent viscosity
        sTurbVisc(iSourceElem) = fluid%turbulence%dataOnLvl(sourceLevel) &
          &                                       %visc(sourceElem)

      end do ! iSourceElem

      ! interpolate moments (rho, momentum and shear stress) to tatrget element
      do iMom = 1,QQ
        tMom(iMom) = sum(weight(1:nSourceElems)*srcMom(iMom,1:nSourceElems))
      end do

      ! store interpolated target auxField
      tOffset = (targetElem-1)*varSys%nAuxScalars
      tAuxField(tOffset+dens_pos)   = tMom(1)
      tAuxField(tOffset+vel_pos(1)) = tMom(layout%moment%first_moments(1))*rho0Inv
      tAuxField(tOffset+vel_pos(2)) = tMom(layout%moment%first_moments(2))*rho0Inv
      tAuxField(tOffset+vel_pos(3)) = tMom(layout%moment%first_moments(3))*rho0Inv

      ! interpolate turbulent viscosity to target element
      tTurbVisc = sum(weight(1:nSourceElems)*sTurbVisc(1:nSourceElems))

      ! scale interpolated turbulent viscosity to target element
      fluid%turbulence%dataOnLvl(targetLevel)%visc(targetElem) &
        & = fluid%turbulence%fac_c2f*tTurbVisc

      ! get normalized kinematic viscosity on target element
      tVisc = fluid%viscKine%dataOnLvl(targetLevel)%val(targetElem)

      ! relation between coarse and fine grid kinematic viscosity:
      ! v^s_f = 2 v^s_c
      ! total viscosity on source element
      sVisc = 0.5_rk * tVisc + tTurbVisc
      ! total viscosity on target element
      tVisc = tVisc + fluid%turbulence%fac_c2f*tTurbVisc

      ! calculate omega on source and target level
      sOmegaKine = mus_calcOmegaFromVisc(sVisc)
      tOmegaKine = mus_calcOmegaFromVisc(tVisc)

      ! Get scaling factors for nonequilibrium moments
      nonEqScalingFacs = fluid%nonEqScalingFacs(           &
        & omegaKine_SRC = sOmegaKine,                      &
        & omegaKine_TGT = tOmegaKine,                      &
        & omegaBulk_SRC = fluid%omegaBulkLvl(sourceLevel), &
        & omegaBulk_TGT = fluid%omegaBulkLvl(targetLevel), &
        & scaleFac      = 0.5_rk, QQ = QQ                  )

      ! Convert moment to PDF
      tPDF = mus_intp_convertMomToPDF3D_incomp(                               &
        &                                moments          = tMom,             &
        &                                nonEqScalingFacs = nonEqScalingFacs, &
        &                                layout           = layout            )

      ! Now write the resulting pdf in the current direction to the target
      ! element position
      do iDir = 1,QQ
        tstate(?IDX?(iDir, targetElem, nScalars, tnSize )) &
          &  = tPDF( iDir )
      end do

    end do

  end subroutine fillFinerGhostsFromMe_weighAvgIncompLES
! **************************************************************************** !



! **************************************************************************** !
  !> [Linear interpolation](../page/features/intp_methods.html) of ghostFromFiner
  !!
  !! This subroutine's interface must match the abstract interface definition
  !! [[intpRoutine]] in intp/[[mus_interpolate_header_module]].f90 in order to
  !! be callable via [[mus_interpolation_method_type:do_intp]] function pointer.
  subroutine fillFinerGhostsFromMe_weighAvg2D( method, fieldProp, tLevelDesc, &
    &                                          level, sState, snSize, tState, &
    &                                          tnSize, tAuxField, layout,     &
    &                                          nTargets, targetList, physics, &
    &                                          time, varSys, derVarPos        )
    ! -------------------------------------------------------------------- !
    class(mus_interpolation_method_type), intent(in) :: method

    !> Array of field properties (fluid or species)
    type(mus_field_prop_type), target, intent(in) :: fieldProp(:)

    !> level descriptor on target level
    type( tem_levelDesc_type ), intent(in) :: tLevelDesc

    !> my refinement level
    integer, intent(in) :: level

    !> State vector of SOURCE FLUID elements
    real(kind=rk), intent(in) :: sState(:)
    ! integer, intent(in) :: sNeigh(:)
    integer, intent(in) :: snSize

    !> State vector of TARGET GHOST elements
    real(kind=rk), intent(inout) :: tState(:)
    ! integer, intent(in) :: tNeigh(:)
    integer, intent(in) :: tnSize

    !> AuxField variable to fill on target GHOST elements
    real(kind=rk), intent(inout) :: tAuxField(:)

    !> the layout used
    type( mus_scheme_layout_type ), intent(in) :: layout

    !> List of target elements ( their position in depSource list )
    integer, intent(in) :: nTargets
    integer, intent(in) :: targetList(nTargets)

    !> physics type to convert lattice to physics SI unit and vice versa
    !! @todo: This can be replaced by scale factor by level
    type( mus_physics_type ), intent(in) :: physics

    !> time required to compute viscosity on target element barycenter
    type(tem_time_type), intent(in) :: time

    !> scheme variable system
    type( tem_varSys_type ), intent(in) :: varSys

    !> position of all derive variable in varSys for all fields
    type(mus_derVarPos_type), intent(in) :: derVarPos(:)
    ! -------------------------------------------------------------------- !
    integer :: sourceLevel    ! level of source elements
    integer :: sourceElem     ! treeId of current source element
    integer :: targetLevel    ! level of target elements
    integer :: targetElem     ! treeId of current source element
    integer :: iDir           ! current direction (discrete velocity) for loop
    integer :: iElem          ! current target element (for outer loop)
    integer :: indElem        ! element counter for indirection list
    integer :: iSourceElem    ! current source element (for inner loop)
    integer :: nSourceElems   ! number of source elements for the current target
    real(kind=rk) :: tPDF( layout%fStencil%QQ )
    real(kind=rk) :: weight( layout%fStencil%QQ )
    ! moments of the source elements' pdf
    real(kind=rk) :: srcMom( layout%fStencil%QQ, layout%fStencil%QQ )
    real(kind=rk) :: tMom(layout%fStencil%QQ)  ! target moment calculation
    integer :: iMom, QQ
    type(mus_fluid_type), pointer :: fluid
    real(kind=rk) :: nonEqScalingFacs(layout%fStencil%QQ)
    real(kind=rk) :: sOmegaKine, tOmegaKine, tVisc, invRho
    integer :: nScalars, tOffset
    integer :: dens_pos, vel_pos(3)
    ! --------------------------------------------------------------------------
    fluid => fieldProp(1)%fluid
    nScalars = varSys%nScalars
    QQ = layout%fStencil%QQ
    dens_pos = varSys%method%val(derVarPos(1)%density)%auxField_varPos(1)
    vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:3)

    sourceLevel = level
    targetLevel = level + 1

    ! Treat all fine target elements:
    do indElem = 1, nTargets

      iElem = targetList( indElem )

      ! Read the target element treeId
      targetElem = iElem + tLevelDesc%offset( 1, eT_ghostFromCoarser)

      ! Find out how many fine source elements we have for interpolation.
      nSourceElems = tLevelDesc%depFromCoarser( iElem )%elem%nVals

      ! Read the pre-calculated weights (like in Average intp.)
      weight(1:nSourceElems) = tLevelDesc%depFromCoarser( iElem ) &
        &                                %weight(1:nSourceElems)

      ! Now loop over all fine source elements for this target:
      do iSourceElem = 1, nSourceElems

        ! Get the source element position
        sourceElem = tLevelDesc%depFromCoarser( iElem ) &
          &                    %elem%val( iSourceElem )

        ! Get macroscopic moments
        srcMom(:, iSourceElem) = mus_intp_getMoments(                         &
          &                               state     = sstate,                 &
          &                               elem      = sourceElem,             &
          &                               QQ        = QQ,                     &
          &                               nScalars  = nScalars,               &
          &                               nSize     = snSize,                 &
          &                               toMoments = layout%moment%toMoments )

      end do ! iSourceElem

      ! interpolate all moments by weighted average
      do iMom = 1,QQ
        tMom(iMom) = sum(weight(1:nSourceElems)*srcMom(iMom,1:nSourceElems))
      end do

      ! store interpolated target auxField
      invRho = 1.0_rk/tMom(1)
      tOffset = (targetElem-1)*varSys%nAuxScalars
      tAuxField(tOffset+dens_pos)   = tMom(1)
      tAuxField(tOffset+vel_pos(1)) = tMom(layout%moment%first_moments(1))*invRho
      tAuxField(tOffset+vel_pos(2)) = tMom(layout%moment%first_moments(2))*invRho
      tAuxField(tOffset+vel_pos(3)) = 0.0_rk

      ! get normalized kinematic viscosity on target element
      tVisc = fluid%viscKine%dataOnLvl(targetLevel)%val(targetElem)
      ! relation between coarse and fine grid kinematic viscosity:
      ! v^s_f = 2 v^s_c
      ! calculate omega on source and target level
      sOmegaKine = mus_calcOmegaFromVisc(0.5_rk*tVisc)
      tOmegaKine = mus_calcOmegaFromVisc(tVisc)

      ! Get scaling factors for nonequilibrium moments
      nonEqScalingFacs = fluid%nonEqScalingFacs(           &
        & omegaKine_SRC = sOmegaKine,                      &
        & omegaKine_TGT = tOmegaKine,                      &
        & omegaBulk_SRC = fluid%omegaBulkLvl(sourceLevel), &
        & omegaBulk_TGT = fluid%omegaBulkLvl(targetLevel), &
        & scaleFac      = 0.5_rk, QQ = QQ                  )

      ! Convert moment to PDF
      tPDF = mus_intp_convertMomToPDF2D( moments          = tMom,             &
        &                                nonEqScalingFacs = nonEqScalingFacs, &
        &                                layout           = layout            )

      ! Now write the resulting pdf to the target element
      do iDir = 1, QQ
        tstate(?IDX?(iDir, targetElem, nScalars, tnSize ) ) &
          & = tPDF( iDir )
      end do
    enddo ! element loop

  end subroutine fillFinerGhostsFromMe_weighAvg2D
! **************************************************************************** !

! **************************************************************************** !
  !> Fill finer ghost from coarser fluid.
  !! This routine is used by 2D and 3D incomp weighted average interpolation.
  !!
  !! This subroutine's interface must match the abstract interface definition
  !! [[intpRoutine]] in intp/[[mus_interpolate_header_module]].f90 in order to
  !! be callable via [[mus_interpolation_method_type:do_intp]] function pointer.
  subroutine fillFinerGhostsFromMe_weighAvg2DIncomp( method, fieldProp,     &
    & tLevelDesc, level, sState, snSize, tState, tnSize, tAuxField, layout, &
    & nTargets, targetList, physics, time, varSys, derVarPos                )
    ! -------------------------------------------------------------------- !
    class(mus_interpolation_method_type), intent(in) :: method

    !> Array of field properties (fluid or species)
    type(mus_field_prop_type), target, intent(in) :: fieldProp(:)

    !> level descriptor on target level
    type( tem_levelDesc_type ), intent(in) :: tLevelDesc

    !> my refinement level
    integer, intent(in) :: level

    !> State vector of SOURCE FLUID elements
    real(kind=rk), intent(in) :: sState(:)
    ! integer, intent(in) :: sNeigh(:)
    integer, intent(in) :: snSize

    !> State vector of TARGET GHOST elements
    real(kind=rk), intent(inout) :: tState(:)
    ! integer, intent(in) :: tNeigh(:)
    integer, intent(in) :: tnSize

    !> AuxField variable to fill on target GHOST elements
    real(kind=rk), intent(inout) :: tAuxField(:)

    !> the layout used
    type( mus_scheme_layout_type ), intent(in) :: layout

    !> List of target elements ( their position in depSource list )
    integer, intent(in) :: nTargets
    integer, intent(in) :: targetList(nTargets)

    !> physics type to convert lattice to physics SI unit and vice versa
    !! @todo: This can be replaced by scale factor by level
    type( mus_physics_type ), intent(in) :: physics

    !> time required to compute viscosity on target element barycenter
    type(tem_time_type), intent(in) :: time

    !> scheme variable system
    type( tem_varSys_type ), intent(in) :: varSys

    !> position of all derive variable in varSys for all fields
    type(mus_derVarPos_type), intent(in) :: derVarPos(:)
    ! -------------------------------------------------------------------- !
    integer :: sourceLevel    ! level of source elements
    integer :: sourceElem     ! treeId of current source element
    integer :: targetLevel    ! level of target elements
    integer :: targetElem     ! treeId of current source element
    integer :: iDir           ! current direction (discrete velocity) for loop
    integer :: iElem          ! current target element (for outer loop)
    integer :: indElem        ! element counter for indirection list
    integer :: iSourceElem    ! current source element (for inner loop)
    integer :: nSourceElems   ! number of source elements for the current target
    ! macroscopic velocities for all source elements
    real(kind=rk) :: tPDF( layout%fStencil%QQ )
    real(kind=rk) :: weight( layout%fStencil%QQ )
    ! moments of the source elements' pdf
    real(kind=rk) :: srcMom( layout%fStencil%QQ, layout%fStencil%QQ )
    real(kind=rk) :: tMom(layout%fStencil%QQ)  ! target moment calculation
    integer :: iMom, QQ
    type(mus_fluid_type), pointer :: fluid
    real(kind=rk) :: nonEqScalingFacs(layout%fStencil%QQ)
    real(kind=rk) :: sOmegaKine, tOmegaKine, tVisc
    integer :: nScalars, tOffset
    integer :: dens_pos, vel_pos(3)
    ! --------------------------------------------------------------------------
!write(dbgUnit(1),*) 'Entering fillFiner average'
    fluid => fieldProp(1)%fluid
    nScalars = varSys%nScalars
    QQ = layout%fStencil%QQ
    dens_pos = varSys%method%val(derVarPos(1)%density)%auxField_varPos(1)
    vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:3)

    sourceLevel = level
    targetLevel = level + 1

    ! Treat all fine target elements:
    do indElem = 1, nTargets

      iElem = targetList( indElem )

      ! Read the target element treeId
      targetElem = iElem + tLevelDesc%offset( 1, eT_ghostFromCoarser)
!write(dbgUnit(1),*) 'Target iElem ', indElem, ' ID ', tLevelDesc%total(targetElem)

      ! Find out how many fine source elements we have for interpolation.
      ! Usually 8, but can differ at corners, obstacles, boundaries...
      nSourceElems = tLevelDesc%depFromCoarser( iElem )%elem%nVals

      ! Read the pre-calculated weights (like in Average intp.)
      weight(1:nSourceElems) = tLevelDesc%depFromCoarser( iElem ) &
        &                                %weight(1:nSourceElems)

      ! Now loop over all fine source elements for this target:
      do iSourceElem = 1, nSourceElems
        ! Get the source element position
        sourceElem = tLevelDesc%depFromCoarser( iElem ) &
          &                    %elem%val( iSourceElem )

        ! Get macroscopic moments
        srcMom(:, iSourceElem) = mus_intp_getMoments(                         &
          &                               state     = sstate,                 &
          &                               elem      = sourceElem,             &
          &                               QQ        = QQ,                     &
          &                               nScalars  = nScalars,               &
          &                               nSize     = snSize,                 &
          &                               toMoments = layout%moment%toMoments )

      end do ! iSourceElem

      ! interpolate rho, momentum and shear stress by weighted average
      do iMom = 1, QQ
        tMom(iMom) = sum(weight(1:nSourceElems)*srcMom(iMom,1:nSourceElems))
      end do
!write(dbgUnit(1),*) 'before scaling tMom ', tMom

      ! store interpolated target auxField
      tOffset = (targetElem-1)*varSys%nAuxScalars
      tAuxField(tOffset+dens_pos)   = tMom(1)
      tAuxField(tOffset+vel_pos(1)) = tMom(layout%moment%first_moments(1))*rho0Inv
      tAuxField(tOffset+vel_pos(2)) = tMom(layout%moment%first_moments(2))*rho0Inv
      tAuxField(tOffset+vel_pos(3)) = 0.0_rk

      ! get normalized kinematic viscosity on target element
      tVisc = fluid%viscKine%dataOnLvl(targetLevel)%val(targetElem)

      ! relation between coarse and fine grid kinematic viscosity:
      ! v^s_f = 2 v^s_c
      ! calculate omega on source and target level
      sOmegaKine = mus_calcOmegaFromVisc(0.5_rk * tVisc)
      tOmegaKine = mus_calcOmegaFromVisc(tVisc)

      ! Get scaling factors for nonequilibrium moments
      nonEqScalingFacs = fluid%nonEqScalingFacs(           &
        & omegaKine_SRC = sOmegaKine,                      &
        & omegaKine_TGT = tOmegaKine,                      &
        & omegaBulk_SRC = fluid%omegaBulkLvl(sourceLevel), &
        & omegaBulk_TGT = fluid%omegaBulkLvl(targetLevel), &
        & scaleFac      = 0.5_rk, QQ = QQ                  )

      ! Convert moment to PDF
      tPDF = mus_intp_convertMomToPDF2D_incomp(                               &
        &                                moments          = tMom,             &
        &                                nonEqScalingFacs = nonEqScalingFacs, &
        &                                layout           = layout            )
!write(dbgUnit(1),*) 'tPDF ', tPDF
      do iDir = 1,QQ
        tstate(?IDX?(iDir, targetElem, nScalars, tnSize )) &
          &  = tPDF( iDir )
      end do

    end do

  end subroutine fillFinerGhostsFromMe_weighAvg2DIncomp
! **************************************************************************** !


end module mus_interpolate_average_module
! **************************************************************************** !
