! Copyright (c) 2012-2013 Manuel Hasert <m.hasert@grs-sim.de>
! Copyright (c) 2012-2014 Simon Zimny <s.zimny@grs-sim.de>
! Copyright (c) 2013, 2020 Harald Klimach <harald.klimach@uni-siegen.de>
! Copyright (c) 2013-2015, 2018-2020 Kannan Masilamani <kannan.masilamani@uni-siegen.de>
! Copyright (c) 2014-2016 Jiaxing Qi <jiaxing.qi@uni-siegen.de>
! Copyright (c) 2014 Kartik Jain <kartik.jain@uni-siegen.de>
! Copyright (c) 2016 Tobias Schneider <tobias1.schneider@student.uni-siegen.de>
! Copyright (c) 2020 Peter Vitt <peter.vitt2@uni-siegen.de>
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions are met:
!
! 1. Redistributions of source code must retain the above copyright notice,
! this list of conditions and the following disclaimer.
!
! 2. Redistributions in binary form must reproduce the above copyright notice,
! this list of conditions and the following disclaimer in the documentation
! and/or other materials provided with the distribution.
!
! THIS SOFTWARE IS PROVIDED BY THE UNIVERSITY OF SIEGEN “AS IS” AND ANY EXPRESS
! OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
! OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
! IN NO EVENT SHALL UNIVERSITY OF SIEGEN OR CONTRIBUTORS BE LIABLE FOR ANY
! DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
! (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
! ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
! ****************************************************************************** !
?? include 'header/lbm_macros.inc'
?? include 'header/lbm_deriveMacros.inc'
!> summary: Interpolation of flow quantities between different grid levels
!!
!! Ghost elements are employed at grid level interfaces to provide valid
!! pdf values to the neighboring fluid elements. This way, the solvers can
!! act on elements of the same size only, treating the levels successively.
!! Target elements are the ghost elements, which have to be filled with
!! valid values.
!! Source elements are the fluid elements from other levels, from where to
!! take the input values for the interpolation.
!! The target ghost elements on the target level have corresponding source
!! fluid elements on the source level.
!!
!! [[tem_topology_module]] For a detailed description of the grid
!!
module mus_interpolate_compact_module
  use iso_c_binding, only: c_f_pointer

  ! include treelm modules
  use env_module,            only: rk
  use tem_element_module,    only: eT_GhostFromCoarser, &
    &                              eT_ghostFromFiner
  use tem_geometry_module,   only: tem_intp_bilinear
  use tem_logging_module,    only: logUnit
  use tem_param_module,      only: cs2inv, cs2, div1_4
  use tem_math_module,       only: tem_intp_trilinearReduced
  use tem_construction_module, only: tem_levelDesc_type
  use tem_varSys_module,       only: tem_varSys_type
  use tem_time_module,       only: tem_time_type

  ! include musubi modules
  use mus_moments_module,            only: get_moment
  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_field_prop_module,         only: mus_field_prop_type
  use mus_fluid_module,              only: mus_fluid_type
  use mus_derVarPos_module,          only: mus_derVarPos_type

  implicit none

  private

  public :: unitTest_interpolate_quadraticCompact
  public :: fillMyGhostsFromFiner_quadraticCompact2d
  public :: fillFinerGhostsFromMe_quadraticCompact2d
  public :: fillMyGhostsFromFiner_quadraticCompact3d
  public :: fillFinerGhostsFromMe_quadraticCompact3d

  integer, parameter :: QmomentumX = 1
  integer, parameter :: QmomentumY = 2
  integer, parameter :: QmomentumZ = 3
  integer, parameter :: QPxx       = 4
  integer, parameter :: QPyy       = 5
  integer, parameter :: QPzz       = 6
  integer, parameter :: QPxy       = 7
  integer, parameter :: QPyz       = 8
  integer, parameter :: QPzx       = 9
  integer, parameter :: c_x  = 1
  integer, parameter :: c_y  = 2
  integer, parameter :: c_z  = 3

 contains

! ****************************************************************************** !
  !> Interpolation with momentum gradients from the stress tensor
  !! source points to the target position located at targetCoord.
  !! The source points are arranged in a square from (0,0,0)x(1,1,1)
  !! The order of the source points are according to the morton curve
  !! See tem_param for the explicit ordering
  !!
  !!> Deactivated because actually not in use
  !pure function mus_interpolate_quadraticCompact3d_leastSq( srcMom,            &
  !  &                            targetCoord, Pab, nSourceElems, matrix )      &
  !  &                            result( phi )
  !  ! ---------------------------------------------------------------------------
  !  !>
  !  integer, intent(in) :: nSourceElems
  !  !> source values of the momentum on the square corners
  !  real(kind=rk), intent(in) :: srcMom(3,nSourceElems)
  !  !>
  !  real(kind=rk), intent(in) :: targetCoord(3)
  !  !> source values of the square corners
  !  real(kind=rk), intent(in) :: Pab(6,nSourceElems)
  !  !> relaxation time for the non-eq moments
  !  real(kind=rk), intent(in) :: matrix(:,:)
  !  !> result: interpolated value
  !  real(kind=rk) :: phi(3)
  !  ! ---------------------------------------------------------------------------
  !  ! Coefficients
  !  real(kind=rk) :: a(10), b(10), c(10)
  !  ! right and left hand side of the equation system
  !  real(kind=rk) :: rhs(9*nSourceElems), lhs(30)
  !  ! ---------------------------------------------------------------------------

  !  ! We extract momentum information completely on the view of the source
  !  ! coordinate system. Set the right hand side of the equation system
  !  rhs( 0*nSourceElems+1:1*nSourceElems)  = srcMom( 1, 1:nSourceElems)
  !  rhs( 1*nSourceElems+1:2*nSourceElems)  = srcMom( 2, 1:nSourceElems)
  !  rhs( 2*nSourceElems+1:3*nSourceElems)  = srcMom( 3, 1:nSourceElems)
  !  rhs( 3*nSourceElems+1:4*nSourceElems)  =    Pab( 1, 1:nSourceElems)
  !  rhs( 4*nSourceElems+1:5*nSourceElems)  =    Pab( 2, 1:nSourceElems)
  !  rhs( 5*nSourceElems+1:6*nSourceElems)  =    Pab( 3, 1:nSourceElems)
  !  rhs( 6*nSourceElems+1:7*nSourceElems)  =    Pab( 4, 1:nSourceElems)
  !  rhs( 7*nSourceElems+1:8*nSourceElems)  =    Pab( 5, 1:nSourceElems)
  !  rhs( 8*nSourceElems+1:9*nSourceElems)  =    Pab( 6, 1:nSourceElems)
  !  ! Solve the problem, where b = rhs, x = coefficients
  !  ! A*x = b
  !  ! overdetermined, solve the least Square fit problem
  !  ! (A^T)A*x = (A^T)b
  !  ! x = ((A^T)A)^-1*(A^T)b
  !  ! Solve linear system of equation with inverted matrix
  !  lhs = matmul( matrix, rhs )

  !  a(1:10) = lhs( 1:10)
  !  b(1:10) = lhs(11:20)
  !  c(1:10) = lhs(21:30)

  !  ! Evaluate the bubble function with the above calculated coefficients
  !  ! m_pi(x) = a1+a2*xcoord+a3*ycoord+a4*zcoord
  !  !             +a5*x^2 + a6*y^2 + a7*z^2
  !  !             +a8*x*y + a9*y*z + a10*x*z
  !  phi(1) = a(1) + a( 2)*targetCoord(c_x)                                     &
  !                + a( 3)*targetCoord(c_y)                                     &
  !                + a( 4)*targetCoord(c_z)                                     &
  !                + a( 5)*targetCoord(c_x)*targetCoord(c_x)                    &
  !                + a( 6)*targetCoord(c_y)*targetCoord(c_y)                    &
  !                + a( 7)*targetCoord(c_z)*targetCoord(c_z)                    &
  !                + a( 8)*targetCoord(c_x)*targetCoord(c_y)                    &
  !                + a( 9)*targetCoord(c_y)*targetCoord(c_z)                    &
  !                + a(10)*targetCoord(c_z)*targetCoord(c_x)
  !  phi(2) = b(1) + b( 2)*targetCoord(c_x)                                     &
  !                + b( 3)*targetCoord(c_y)                                     &
  !                + b( 4)*targetCoord(c_z)                                     &
  !                + b( 5)*targetCoord(c_x)*targetCoord(c_x)                    &
  !                + b( 6)*targetCoord(c_y)*targetCoord(c_y)                    &
  !                + b( 7)*targetCoord(c_z)*targetCoord(c_z)                    &
  !                + b( 8)*targetCoord(c_x)*targetCoord(c_y)                    &
  !                + b( 9)*targetCoord(c_y)*targetCoord(c_z)                    &
  !                + b(10)*targetCoord(c_z)*targetCoord(c_x)
  !  phi(3) = c(1) + c( 2)*targetCoord(c_x)                                     &
  !                + c( 3)*targetCoord(c_y)                                     &
  !                + c( 4)*targetCoord(c_z)                                     &
  !                + c( 5)*targetCoord(c_x)*targetCoord(c_x)                    &
  !                + c( 6)*targetCoord(c_y)*targetCoord(c_y)                    &
  !                + c( 7)*targetCoord(c_z)*targetCoord(c_z)                    &
  !                + c( 8)*targetCoord(c_x)*targetCoord(c_y)                    &
  !                + c( 9)*targetCoord(c_y)*targetCoord(c_z)                    &
  !                + c(10)*targetCoord(c_z)*targetCoord(c_x)

  !end function mus_interpolate_quadraticCompact3d_leastSq
! ****************************************************************************** !


! ****************************************************************************** !
  !> Determine the coefficients of quadratic polynomial
  !! from momentum and strain rate of source elements.
  !! Evaluate the result of polynomial at targetCoord.
  !! The source points are arranged in a square from (0,0)x(1,1)
  !! The order of the source points are according to the morton curve
  !!```
  !!  1      2      3      4
  !! (0,0); (1,0); (0,1); (1,1)
  !!```
  !! This routine solves the equation system:
  !!```
  !! A * x = b by x = A^(-1) * b
  !!```
  !!> Deactivated, because actually not in use
  !pure function mus_interpolate_quadraticCompact2d_leastSq( srcMom,          &
  !  &                                targetCoord, Pab, nSourceElems, invA )  &
  !  &                                result( phi )
  !  ! ---------------------------------------------------------------------------
  !  !> number of source elements
  !  integer, intent(in) :: nSourceElems
  !  !> Velocity or Momentum Ux, Uy
  !  real(kind=rk), intent(in) :: srcMom(2,nSourceElems)
  !  !> interpolation location within the square
  !  real(kind=rk), intent(in) :: targetCoord(2)
  !  !> Strain rate Sxx, Syy, Sxy
  !  real(kind=rk), intent(in) :: Pab(3,nSourceElems)
  !  !> inverted matrix
  !  real(kind=rk), intent(in) :: invA(:,:)
  !  !> interpolated value
  !  real(kind=rk) :: phi(2)
  !  ! ---------------------------------------------------------------------------
  !  real(kind=rk) :: a(6), b(6) ! Coefficients
  !  real(kind=rk) :: rhs(5*nSourceElems), lhs(12) ! Coefficients
  !  ! ---------------------------------------------------------------------------

  !  ! We extract momentum information completely on the view of the source
  !  ! coordinate system. Set the right hand side of the equation system
  !  rhs( 0*nSourceElems+1:1*nSourceElems)  = srcMom( 1, 1:nSourceElems)
  !  rhs( 1*nSourceElems+1:2*nSourceElems)  = srcMom( 2, 1:nSourceElems)
  !  rhs( 2*nSourceElems+1:3*nSourceElems)  =    Pab( 1, 1:nSourceElems)
  !  rhs( 3*nSourceElems+1:4*nSourceElems)  =    Pab( 2, 1:nSourceElems)
  !  rhs( 4*nSourceElems+1:5*nSourceElems)  =    Pab( 3, 1:nSourceElems)

  !  ! Solve the problem: A*x = b
  !  ! x = A^(-1) * b, where b = rhs, x = coefficients
  !  ! Solve linear system of equation with inverted matrix
  !  lhs = matmul( invA, rhs )

  !  a(1:6) = lhs(1:6)
  !  b(1:6) = lhs(7:12)
  !  ! Evaluate the bubble function with the above calculated coefficients
  !  ! m_pi(x) = a1+a2*xcoord+a3*ycoord+a4*xcoord*xcoord+a5*ycoord*ycoord+a6*xcoord*ycoord;
  !  phi(1) = a(1) + a(2)*targetCoord(c_x) &
  !    &           + a(3)*targetCoord(c_y) &
  !    &           + a(4)*targetCoord(c_x)*targetCoord(c_x) &
  !    &           + a(5)*targetCoord(c_y)*targetCoord(c_y) &
  !    &           + a(6)*targetCoord(c_x)*targetCoord(c_y)

  !  ! m_pi(y) = b1+b2*xcoord+b3*ycoord+b4*xcoord*xcoord+b5*ycoord*ycoord+b6*xcoord*ycoord;
  !  phi(2) = b(1) + b(2)*targetCoord(c_x) &
  !    &           + b(3)*targetCoord(c_y) &
  !    &           + b(4)*targetCoord(c_x)*targetCoord(c_x) &
  !    &           + b(5)*targetCoord(c_y)*targetCoord(c_y) &
  !    &           + b(6)*targetCoord(c_x)*targetCoord(c_y)

  !end function mus_interpolate_quadraticCompact2d_leastSq
! ****************************************************************************** !


! ****************************************************************************** !
  !> No comment yet!
  !!
  !! @TODO add comment
  !!
  !! 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_quadraticCompact3d( 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 :: iVal           ! value counter
    integer :: iElem          ! element counter 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) :: f( layout%fStencil%QQ ) ! pdf to reconstruct from
    real(kind=rk) :: fEq( layout%fStencil%QQ ) ! source equilirbium part
    ! moments of the source elements' pdf
    real(kind=rk) :: sourceMom( layout%fStencil%QQ, method%nMaxSources )
    ! interpolated moments
    real(kind=rk) :: targetMom( layout%fStencil%QQ )
    real(kind=rk) :: targetCoord(3) ! coordinates to interpolate to
    ! Coordinates in rotated reference cube
    real(kind=rk), parameter :: targetCoord_trafo(3) = [div1_4, div1_4, div1_4]
    ! shear stress scaling factor
    real(kind=rk) :: omfac
    real(kind=rk) :: tMom( layout%fStencil%QQ )  ! temp moment calculation
    integer :: iCase, iDir
    integer :: QQ
    real(kind=rk)  :: velocity(3)
    type(mus_fluid_type), pointer :: fluid
    integer :: nScalars
    real(kind=rk), allocatable :: momBuf(:,:)
    ! ---------------------------------------------------------------------------
    fluid => fieldProp(1)%fluid
    nScalars = varSys%nScalars

    QQ = layout%fStencil%QQ
    sourceLevel = level
    targetLevel = level + 1

    ! KM: temp fix: assuming omega is constant
    omfac = (1._rk - 0.5*fluid%viscKine%omLvl( sourceLevel )%val(1)) &
      &   / (1._rk - 0.5*fluid%viscKine%omLvl( targetLevel )%val(1))

    allocate( momBuf( QQ, tLevelDesc%sourceFromCoarser%nVals ) )

    ! First calculate all the required moments for all the source elements
    do iSourceElem = 1, tLevelDesc%sourceFromCoarser%nVals

      ! Get the source element's position in the state vector
      sourceElem = tLevelDesc%sourceFromCoarser%val( iSourceElem )

      do iDir = 1, QQ
        f(iDir)  = sState( ?IDX?( iDir, sourceElem, QQ, snSize ))
      enddo

      ! Calculate the raw moments form the pdf
      tMom( : ) = matmul( layout%moment%toMoments%A, f)

      velocity = [ tMom(2)/tMom(1), tMom(3)/tMom(1), tMom(4)/tMom(1) ]

      do iDir = 1, QQ
?? copy :: EqDir_comp( fEq(iDir), tMom(1), velocity, layout, iDir, cs2inv, rk )
      end do

      ! remove pi(0) parts from the second order moment
      ! that is ((rho u_\alpha * rho u_\beta) / rho)
      ! and additionally rho*cs2 from diagonal entries
      ! Pab = sum( c_ia*c_ib *f_i ) - kronecker_ab*(rho*cs2) - rho*u_a*u_b
      tMom( 5 ) = tMom( 5 ) - cs2*tMom( 1) & !Pxx
                   - tMom( 2)*tMom( 2)/tMom( 1)
      tMom( 6 ) = tMom( 6 ) - cs2*tMom( 1) & !Pyy
                   - tMom( 3)*tMom( 3)/tMom( 1)
      tMom( 7 ) = tMom( 7 ) - cs2*tMom( 1) & !Pzz
                   - tMom( 4)*tMom( 4)/tMom( 1)
      tMom( 8 ) = tMom( 8 )                & !Pxy
                   - tMom( 2)*tMom( 3)/tMom( 1)
      tMom( 9 ) = tMom( 9 )                & !Pyz
                   - tMom( 3)*tMom( 4)/tMom( 1)
      tMom(10 ) = tMom(10 )                & !Pzx
                   - tMom( 2)*tMom( 4)/tMom( 1)
      ! Transfer momentum to velocity
      tMom( 2:4 ) = tMom( 2:4 )/ tMom( 1 )

      ! Set the further higher moments to their equilibrium distribution
      tMom(11 ) = get_moment( QQ, layout%fStencil%cxDir, (/2,1,0/), fEq ) !mxxy
      tMom(12 ) = get_moment( QQ, layout%fStencil%cxDir, (/2,0,1/), fEq ) !mxxz
      tMom(13 ) = get_moment( QQ, layout%fStencil%cxDir, (/1,2,0/), fEq ) !myyx
      tMom(14 ) = get_moment( QQ, layout%fStencil%cxDir, (/0,2,1/), fEq ) !myyz
      tMom(15 ) = get_moment( QQ, layout%fStencil%cxDir, (/1,0,2/), fEq ) !mzzx
      tMom(16 ) = get_moment( QQ, layout%fStencil%cxDir, (/0,1,2/), fEq ) !mzzy
      tMom(17 ) = get_moment( QQ, layout%fStencil%cxDir, (/2,2,0/), fEq ) !mxxyy
      tMom(18 ) = get_moment( QQ, layout%fStencil%cxDir, (/0,2,2/), fEq ) !myyzz
      tMom(19 ) = get_moment( QQ, layout%fStencil%cxDir, (/2,0,2/), fEq ) !mzzxx

      momBuf( :, iSourceElem ) = tMom

    enddo ! nSourceElems


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

      iElem = targetList( indElem )

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

      ! get the target coordinates
      targetCoord = tLevelDesc%depFromCoarser( iElem )%coord

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

      !Compute rho and velocity for each field and interpolate pdf field wise
      ! Now loop over all fine source elements for this target:
      do iSourceElem = 1, nSourceElems
        ! Get the source element in the levelwise buffer
        sourceElem = tLevelDesc%depFromCoarser( iElem )%  &
          &                    elemBuffer%val( iSourceElem )
        ! and read its entries
        sourceMom(:, iSourceElem ) = momBuf(:, sourceElem)
      end do

      ! interpolate the density linearly
      targetMom(1)  = tem_intp_trilinearReduced(                             &
        &                    sourceMom(1,1:nSourceElems), targetCoord_trafo )

      ! get the rotation case from the child number of the target element
      iCase = tLevelDesc%depFromCoarser( iElem )%childNum

      ! targetMom(2:4) = mus_interpolate_quadraticCompact3d_leastSq(           &
      !   &             sourceMom(2:4, 1:nSourceElems),                        &
      !   &             targetCoord, sourceMom(5:10, 1:nSourceElems),          &
      !   &             nSourceElems,                                          &
      !   &             intp%fromCoarser( targetLevel )%matrixInv( iCase )%A )

      ! linear interpolation for higher moments
      do iVal = 5,19
        targetMom(iVal) = tem_intp_trilinearReduced( sourceMom(              &
          &                         iVal,1:nSourceElems), targetCoord_trafo )
      enddo

      ! Rescale the shear stress moments
      targetMom(5:10) = targetMom(5:10)*omfac

      ! and add back the pi(0) part to the tensor
      ! No need to rescale higher order moments as they were set
      ! to equilibrium
      targetMom( 5 ) = targetMom( 5 ) + cs2*targetMom( 1) & !Pxx
                   + targetMom( 2)*targetMom( 2)/targetMom( 1)
      targetMom( 6 ) = targetMom( 6 ) + cs2*targetMom( 1) & !Pyy
                   + targetMom( 3)*targetMom( 3)/targetMom( 1)
      targetMom( 7 ) = targetMom( 7 ) + cs2*targetMom( 1) & !Pzz
                   + targetMom( 4)*targetMom( 4)/targetMom( 1)
      targetMom( 8 ) = targetMom( 8 )                     & !Pxy
                   + targetMom( 2)*targetMom( 3)/targetMom( 1)
      targetMom( 9 ) = targetMom( 9 )                     & !Pyz
                   + targetMom( 3)*targetMom( 4)/targetMom( 1)
      targetMom(10 ) = targetMom(10 )                     & !Pzx
                   + targetMom( 2)*targetMom( 4)/targetMom( 1)

      ! Transfer back velocity field to momentum field
      targetMom(2:4) = targetMom(2:4)*targetMom(1)

      ! transform the moments back to the pdfs
      f = matmul( layout%moment%toPdf%A, targetMom )

      ! Now write the resulting pdf in the current direction to the target
      ! Element position
      do iVal = 1, QQ
        tState( ?IDX?( iVal, targetElem, QQ, tnSize )) = f(iVal)
      enddo
    enddo ! indElem loop

     deallocate( momBuf )

  end subroutine fillFinerGhostsFromMe_quadraticCompact3d
! ****************************************************************************** !


! ****************************************************************************** !
  !> No comment yet!
  !!
  !! @TODO add comment
  !!
  !! 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_quadraticCompact3d( 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 :: iVal           ! value counter
    integer :: iElem          ! element counter 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) :: f( layout%fStencil%QQ )   ! pdf to reconstruct from
    real(kind=rk) :: fEq( layout%fStencil%QQ ) ! source equilirbium part
    ! moments of the source elements' pdf
    real(kind=rk) :: sourceMom( layout%fStencil%QQ, method%nMaxSources )
    ! interpolated moments
    real(kind=rk) :: targetMom( layout%fStencil%QQ )
    ! shear stress scaling factor
    real(kind=rk) :: omfac
    integer :: iDir
    real(kind=rk) :: tMom( layout%fStencil%QQ )  ! temp moment calculation
    integer :: QQ
    type(mus_fluid_type), pointer :: fluid
    integer :: nScalars
    real(kind=rk), allocatable :: momBuf(:,:)
    ! ---------------------------------------------------------------------------
    fluid => fieldProp(1)%fluid
    nScalars = varSys%nScalars
    QQ = layout%fStencil%QQ
    targetLevel = level
    sourceLevel = level+1

    ! KM: temp fix: assuming omega is constant
    omfac = (1._rk - 0.5*fluid%viscKine%omLvl( sourceLevel )%val(1)) &
      &   / (1._rk - 0.5*fluid%viscKine%omLvl( targetLevel )%val(1))

    allocate( momBuf( layout%fStencil%QQ,                 &
       &               tLevelDesc%sourceFromFiner%nVals ) )

    ! First calculate all the required moments for all the source elements
    do iSourceElem = 1, tLevelDesc%sourceFromFiner%nVals
      ! Get the source element's position in the state vector
      sourceElem = tLevelDesc%sourceFromFiner%val( iSourceElem )
      do iDir = 1, QQ
        f(iDir)  = sState( ?IDX?( iDir, sourceElem, QQ, snSize ))
      enddo
      ! Calculate the raw moments form the pdf
      tMom( : ) = matmul( layout%moment%toMoments%A, f)

      ! remove pi(0) parts from the second order moment
      ! that is ((rho u_\alpha * rho u_\beta) / rho)
      ! and additionally rho*cs2 from diagonal entries
      ! Pab = sum( c_ia*c_ib *f_i ) - kronecker_ab*(rho*cs2) - rho*u_a*u_b
      tMom( 5 ) = tMom( 5 ) - cs2*tMom( 1) & !Pxx
                   - tMom( 2)*tMom( 2)/tMom( 1)
      tMom( 6 ) = tMom( 6 ) - cs2*tMom( 1) & !Pyy
                   - tMom( 3)*tMom( 3)/tMom( 1)
      tMom( 7 ) = tMom( 7 ) - cs2*tMom( 1) & !Pzz
                   - tMom( 4)*tMom( 4)/tMom( 1)
      tMom( 8 ) = tMom( 8 )                & !Pxy
                   - tMom( 2)*tMom( 3)/tMom( 1)
      tMom( 9 ) = tMom( 9 )                & !Pyz
                   - tMom( 3)*tMom( 4)/tMom( 1)
      tMom(10 ) = tMom(10 )                & !Pzx
                   - tMom( 2)*tMom( 4)/tMom( 1)
      ! Transfer momentum to velocity
      tMom( 2:4 ) = tMom( 2:4 )/ tMom( 1 )

      ! And calculate the equilibrium distribution from the macroscopic
      ! quantities in tMom
      fEq = getEquilibrium( tMom(1),    & ! density
        &             [ tMom(2), tMom(3), tMom(4) ],  & ! velocity
        &                layout )

      ! Set the further higher moments to their equilibrium distribution
      tMom(11) = get_moment( QQ, layout%fStencil%cxDir, [2,1,0], fEq ) !mxxy
      tMom(12) = get_moment( QQ, layout%fStencil%cxDir, [2,0,1], fEq ) !mxxz
      tMom(13) = get_moment( QQ, layout%fStencil%cxDir, [1,2,0], fEq ) !myyx
      tMom(14) = get_moment( QQ, layout%fStencil%cxDir, [0,2,1], fEq ) !myyz
      tMom(15) = get_moment( QQ, layout%fStencil%cxDir, [1,0,2], fEq ) !mzzx
      tMom(16) = get_moment( QQ, layout%fStencil%cxDir, [0,1,2], fEq ) !mzzy
      tMom(17) = get_moment( QQ, layout%fStencil%cxDir, [2,2,0], fEq ) !mxxyy
      tMom(18) = get_moment( QQ, layout%fStencil%cxDir, [0,2,2], fEq ) !myyzz
      tMom(19) = get_moment( QQ, layout%fStencil%cxDir, [2,0,2], fEq ) !mzzxx

      momBuf( :, iSourceElem ) = tMom

    end do ! nSourceElems

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

      iElem = targetList( indElem )

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

      nSourceElems = tLevelDesc%depFromFiner( iElem)%elem%nVals

      ! Now loop over all fine source elements for this target:
      do iSourceElem = 1, nSourceElems
        ! Get the source element in the levelwise buffer
        sourceElem = tLevelDesc%depFromFiner( iElem )%    &
          &                     elemBuffer%val( iSourceElem )
        ! and read its entries
        sourceMom(:, iSourceElem ) = momBuf(:, sourceElem )
      enddo

      ! interpolate the density linearly
      targetMom(1) = sum(sourceMom(1,1:nSourceElems)) / dble(nSourceElems)

      ! quadratic interpolation for velX, velY, velZ

      ! targetMom(2:4) = mus_interpolate_quadraticCompact3d_leastSq(           &
      !   &                 sourceMom(2:4, 1:nSourceElems),                    &
      !   &                 targetCoord, sourceMom(5:10, 1:nSourceElems),      &
      !   &                 nSourceElems,                                      &
      !   &                 intp%fromFiner( targetLevel )%matrixInv( 8 )%A )

      ! linear interpolation for higher moments
      do iVal = 5,19
        targetMom(iVal) =   sum(sourceMom(iVal,1:nSourceElems)) &
          &               / dble(nSourceElems)
      enddo

      ! shear stress part has to be scaled
      targetMom(5:10) = targetMom(5:10)*omfac

      ! and add back the pi(0) part to the tensor
      ! No need to rescale higher order moments as they were set to
      ! equilibrium
      targetMom( 5) =   targetMom( 5 ) + cs2*targetMom( 1) & !Pxx
                      + targetMom( 2)*targetMom( 2)/targetMom( 1)
      targetMom( 6) =   targetMom( 6 ) + cs2*targetMom( 1) & !Pyy
                      + targetMom( 3)*targetMom( 3)/targetMom( 1)
      targetMom( 7) =   targetMom( 7 ) + cs2*targetMom( 1) & !Pzz
                      + targetMom( 4)*targetMom( 4)/targetMom( 1)
      targetMom( 8) =   targetMom( 8 )                     & !Pxy
                      + targetMom( 2)*targetMom( 3)/targetMom( 1)
      targetMom( 9) =   targetMom( 9 )                     & !Pyz
                      + targetMom( 3)*targetMom( 4)/targetMom( 1)
      targetMom(10) =   targetMom(10 )                     & !Pzx
                      + targetMom( 2)*targetMom( 4)/targetMom( 1)

      ! Transfer back velocity field to momentum field
      targetMom(2:4) = targetMom(2:4)*targetMom(1)

      ! transform the moments back to the pdfs
      f = matmul( layout%moment%toPdf%A, targetMom )

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

    deallocate( momBuf )

  end subroutine fillMyGhostsFromFiner_quadraticCompact3d
! ****************************************************************************** !


! ****************************************************************************** !
  !> Fill parent from 4 children
  !!
  !! 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_quadraticCompact2d( 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 :: iVal           ! value counter
    integer :: iElem          ! element counter for outer loop
    integer :: iDir           ! link counter for stencil 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) :: targetPdf( layout%fStencil%QQ ) ! pdf to reconstruct
    real(kind=rk) :: sourcePdf( layout%fStencil%QQ ) ! pdf to reconstruct from
    ! moments of the source elements' pdf
    real(kind=rk) :: sourceMom( layout%fStencil%QQ, method%nMaxSources )
    ! interpolated moments
    real(kind=rk) :: targetMom( layout%fStencil%QQ )
    ! shear stress scaling factor
    real(kind=rk) :: omfac
    real(kind=rk) :: sourceEq( layout%fStencil%QQ )  ! source equilirbium part
    real(kind=rk) :: u_x, u_y, usq, ucx
    integer :: QQ
    type(mus_fluid_type), pointer :: fluid
    integer :: nScalars
    ! ---------------------------------------------------------------------------
    fluid => fieldProp(1)%fluid
    nScalars = varSys%nScalars
    QQ = layout%fStencil%QQ
    targetLevel = level
    sourceLevel = level+1

    ! KM: temp fix: assuming omega is constant
    omfac = (1._rk - 0.5*fluid%viscKine%omLvl( sourceLevel )%val(1)) &
      &   / (1._rk - 0.5*fluid%viscKine%omLvl( targetLevel )%val(1))

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

      iElem = targetList( indElem )

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

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

      ! Compute sourceRho and velocity for each field and interpolate pdf

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

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

          ! Read the pdf values from the state vector into sourcePdf
          do iDir = 1, QQ
            sourcePdf(iDir)  = sState( ?IDX?( iDir, sourceElem, QQ, snSize ))
          end do

          ! Calculate the raw moments form the pdf
          sourceMom( :,iSourceElem ) = matmul( layout%moment%toMoments%A, sourcePdf)

          ! Calculate fEq
          ! velX, velY
          u_x = sourceMom(2, iSourceElem) / sourceMom(1, iSourceElem)
          u_y = sourceMom(3, iSourceElem) / sourceMom(1, iSourceElem)

          ! square velocity and derived constants
          usq = u_x*u_x + u_y*u_y

          do iDir = 1, QQ
            ! velocity times lattice unit velocity
            ucx =   dble( layout%fStencil%cxDir( 1, iDir ) )*u_x &
              &   + dble( layout%fStencil%cxDir( 2, iDir ) )*u_y
            ! calculate equilibrium density
            sourceEq( iDir ) = layout%weight( iDir )                           &
              &              * sourceMom(1, iSourceElem)                       &
              &              * ( 1._rk + ucx * cs2inv + ucx * ucx * cs2inv     &
              &              * cs2inv * 0.5_rk - usq * cs2inv * 0.5_rk )
          end do

          ! remove pi(0) parts from the second order moment
          ! that is ((rho u_\alpha * rho u_\beta) / rho)
          ! and additionally rho*cs2 from diagonal entries
          ! Pab = sum( c_ia*c_ib *f_i ) - kronecker_ab*(rho*cs2) - rho*u_a*u_b
          sourceMom( 4, iSourceElem ) = sourceMom( 4, iSourceElem )            &
            &          - cs2*sourceMom( 1, iSourceElem)                        &
            &          - sourceMom( 2, iSourceElem)*sourceMom( 2, iSourceElem) &
            &          / sourceMom( 1, iSourceElem)
          sourceMom( 5, iSourceElem ) = sourceMom( 5, iSourceElem )            &
            &          - cs2*sourceMom( 1, iSourceElem)                        &
            &          - sourceMom( 3, iSourceElem)*sourceMom( 3, iSourceElem) &
            &          / sourceMom( 1, iSourceElem)
          sourceMom( 6, iSourceElem ) = sourceMom( 6, iSourceElem )            &
            &          - sourceMom( 2, iSourceElem)*sourceMom( 3, iSourceElem) &
            &          / sourceMom( 1, iSourceElem)

          ! Transfer momentum to velocity
          sourceMom( 2:3, iSourceElem ) = sourceMom( 2:3, iSourceElem )        &
            &                           / sourceMom( 1, iSourceElem )

          ! Set the further higher moments to their equilibrium distribution
          sourceMom( 7, iSourceElem ) = get_moment( QQ, layout%fStencil%cxDir, &
            &                                       (/2,1,0/), sourceEq        )
          sourceMom( 8, iSourceElem ) = get_moment( QQ, layout%fStencil%cxDir, &
            &                                       (/1,2,0/), sourceEq        )
          sourceMom( 9, iSourceElem ) = get_moment( QQ, layout%fStencil%cxDir, &
            &                                       (/2,2,0/), sourceEq        )
        end do

        ! interpolate the density by averaging
        targetMom(1) = sum(sourceMom(1,1:nSourceElems))                        &
          &          / real(nSourceElems,kind=rk)

        ! quadratic interpolation for momentumX, momentumY
        ! targetMom(2:3) = mus_interpolate_quadraticCompact2d_leastSq(           &
        !   &              sourceMom(2:3, 1:nSourceElems),                       &
        !   &              targetCoord(1:2), sourceMom(4:6, 1:nSourceElems),     &
        !   &              nSourceElems,                                         &
        !   &              intp%fromFiner( targetLevel )%matrixInv( 4 )%A )

        ! linear interpolation for higher moments
        do iVal = 4,9
          targetMom(iVal) = sum(sourceMom(iVal,1:nSourceElems)) &
            &             / real(nSourceElems,rk)
        enddo

        ! shear stress part has to be scaled with fac
        targetMom(4:6) = targetMom(4:6)*omfac
        ! and add back the pi(0) part to the tensor
        ! No need to rescale higher order moments as they were set to
        ! equilibrium
        targetMom(4) = targetMom(4) + cs2*targetMom(1)                         &
          &          + targetMom(2)*targetMom(2)/targetMom(1)
        targetMom(5) = targetMom(5) + cs2*targetMom(1)                         &
          &          + targetMom(3)*targetMom(3)/targetMom(1)
        targetMom(6) = targetMom(6)                                            &
          &          + targetMom(2)*targetMom(3)/targetMom(1)
        ! Transfer back velocity field to momentum field
        targetMom(2:3) = targetMom(2:3)*targetMom(1)

        ! transform the moments back to the pdfs
        targetPdf = matmul( layout%moment%toPdf%A, targetMom )

        ! Now write the resulting pdf in the current direction to the target
        ! Element position
        do iDir = 1, QQ
          tState( ?IDX?( iDir, targetElem, QQ, tnSize )) = targetPdf( iDir )
        end do ! iVal
    enddo

  end subroutine fillMyGhostsFromFiner_quadraticCompact2d
! ****************************************************************************** !


! ****************************************************************************** !
  !> Fill Finer Ghost elements from My Coarser element by 2D compact stencil.
  !!
  !! 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_quadraticCompact2d( 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 :: iElem          ! element counter for outer loop
    integer :: iDir           ! link counter for stencil 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) :: targetPdf( layout%fStencil%QQ ) ! pdf to reconstruct
    real(kind=rk) :: sourcePdf( layout%fStencil%QQ ) ! pdf to reconstruct from
    ! moments of the source elements' pdf
    real(kind=rk) :: sourceMom( layout%fStencil%QQ, method%nMaxSources )
    ! interpolated moments
    real(kind=rk) :: targetMom( layout%fStencil%QQ )
    real(kind=rk) :: targetCoord(3) ! coordinates to interpolate to
    ! Coordinates in rotated reference cube
    real(kind=rk) :: targetCoord_trafo(3)
    ! shear stress scaling factor
    real(kind=rk) :: omfac
    real(kind=rk) :: sourceEq( layout%fStencil%QQ ) ! source equilirbium part
    integer :: iCase ! child case for rotation of the unit cube
    real(kind=rk) :: u_x     ! local x-velocity
    real(kind=rk) :: u_y     ! local y-velocity
    real(kind=rk) :: usq, ucx  ! square velocity
    integer :: QQ
    type(mus_fluid_type), pointer :: fluid
    integer :: nScalars
    ! ---------------------------------------------------------------------------
    fluid => fieldProp(1)%fluid
    nScalars = varSys%nScalars

    QQ = layout%fStencil%QQ
    sourceLevel = level
    targetLevel = level + 1

    ! KM: temp fix: assuming omega is constant
    omfac = (1._rk - 0.5*fluid%viscKine%omLvl( sourceLevel )%val(1)) &
      &   / (1._rk - 0.5*fluid%viscKine%omLvl( targetLevel )%val(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)

      ! get the target coordinates
      targetCoord = tLevelDesc%depFromCoarser( iElem )%coord
      targetCoord_trafo = (/0.25_rk, 0.25_rk, 0._rk /)

      nSourceElems =4

      !Compute rho and velocity for each field and interpolate pdf field wise
      ! Now loop over all fine source elements for this target:
      do iSourceElem = 1, nSourceElems
        ! Get the source element's treeId
        sourceElem = tLevelDesc%depFromCoarser( iElem )%  &
          &                     elem%val( iSourceElem )
        ! get the pdf values for the source elements from the state vector
        ! and store in sourcePdf
        do iDir = 1, QQ
          sourcePdf(iDir)  = sState( ?IDX?( iDir, sourceElem, QQ, snSize ))
        enddo

        ! Calculate the raw moments from the pdfs
        sourceMom( :,iSourceElem ) = matmul( layout%moment%toMoments%A,      &
          &                                                       sourcePdf )

        ! velX, velY
        u_x = sourceMom( 2, iSourceElem)/sourceMom(1, iSourceElem)
        u_y = sourceMom( 3, iSourceElem)/sourceMom(1, iSourceElem)

        ! square velocity and derived constants
        usq  = u_x*u_x + u_y*u_y

        do iDir = 1, QQ
          ! velocity times lattice unit velocity
          ucx = dble( layout%fStencil%cxDir( 1, iDir ) )*u_x               &
            & + dble( layout%fStencil%cxDir( 2, iDir ) )*u_y
          ! calculate equilibrium density
          sourceEq( iDir ) = layout%weight( iDir )*sourceMom(1, iSourceElem) &
            &              * ( 1._rk + ucx*cs2inv + ucx * ucx * cs2inv       &
            &              * cs2inv*0.5_rk - usq * cs2inv * 0.5_rk )
        end do

        ! remove pi(0) parts from the second order moment
        sourceMom( 4, iSourceElem ) = sourceMom( 4, iSourceElem )            &
          &          - cs2*sourceMom( 1, iSourceElem)                        &
          &          - sourceMom( 2, iSourceElem)*sourceMom( 2, iSourceElem) &
          &          / sourceMom( 1, iSourceElem)
        sourceMom( 5, iSourceElem ) = sourceMom( 5, iSourceElem )            &
          &          - cs2*sourceMom( 1, iSourceElem)                        &
          &          - sourceMom( 3, iSourceElem)*sourceMom( 3, iSourceElem) &
          &          / sourceMom( 1, iSourceElem)
        sourceMom( 6, iSourceElem ) = sourceMom( 6, iSourceElem )            &
          &          - sourceMom( 2, iSourceElem)*sourceMom( 3, iSourceElem) &
          &          / sourceMom( 1, iSourceElem)
        ! Set the further higher moments to equilibrium
        sourceMom( 7, iSourceElem ) = get_moment( QQ, layout%fStencil%cxDir, &
          &                                       (/2,1,0/), sourceEq        )
        sourceMom( 8, iSourceElem ) = get_moment( QQ, layout%fStencil%cxDir, &
          &                                       (/1,2,0/), sourceEq        )
        sourceMom( 9, iSourceElem ) = get_moment( QQ, layout%fStencil%cxDir, &
          &                                       (/2,2,0/), sourceEq        )
        ! Transfer momentum field to velocity field
        sourceMom( 2:3, iSourceElem ) = sourceMom( 2:3, iSourceElem )        &
          &                           / sourceMom( 1, iSourceElem )
      end do
      ! interpolate the density linearly
      targetMom(1)  = tem_intp_bilinear( sourceMom(1,1:nSourceElems),        &
        &                                             targetCoord_trafo(1:2))

      ! get the rotation case from the child number of the target element
      iCase = tLevelDesc%depFromCoarser( iElem )%childNum

      ! targetMom(2:3) = mus_interpolate_quadraticCompact2d_leastSq(           &
      !   &              sourceMom(2:3, 1:nSourceElems),                       &
      !   &              targetCoord(1:2), sourceMom(4:6, 1:nSourceElems),     &
      !   &              nSourceElems,                                         &
      !   &              intp%fromCoarser( targetLevel )%matrixInv( iCase )%A )
      ! linear interpolation for higher moments
      targetMom(4:9) = tem_intp_bilinear( sourceMom(4:9,1:nSourceElems),     &
        &                                          targetCoord_trafo(1:2), 6)

      ! Rescale the shear stress moments
      targetMom(4:6) = targetMom(4:6)*omfac
      ! and add back the pi(0) part to the tensor
      targetMom(4) = targetMom(4) + cs2*targetMom(1)                         &
        &          + targetMom(2)*targetMom(2)/targetMom(1)
      targetMom(5) = targetMom(5) + cs2*targetMom(1)                         &
        &          + targetMom(3)*targetMom(3)/targetMom(1)
      targetMom(6) = targetMom(6)                                            &
        &          + targetMom(2)*targetMom(3)/targetMom(1)
      ! Transfer back velocity field to momentum field
      targetMom(2:3) = targetMom(2:3)*targetMom(1)

      ! Transform back into pdf space
      targetPdf = matmul( layout%moment%toPdf%A, targetMom )

      ! Now write the resulting pdf in the current direction to the target
      ! Element position
      do iDir = 1, QQ
        tState( ?IDX?( iDir, targetElem, QQ, tnSize )) = targetPdf( iDir )
      enddo
    end do

  end subroutine fillFinerGhostsFromMe_quadraticCompact2d
! ****************************************************************************** !


! ****************************************************************************** !
  !>
  !!
  subroutine unitTest_interpolate_quadraticCompact( iCase, passed )
    ! ---------------------------------------------------------------------------
    !>
    integer, intent(in)  :: iCase
    !>
    logical, intent(out) :: passed
    ! ---------------------------------------------------------------------------
    ! non-equilibrium part of the resulting pdf
    real(kind=rk) :: sourceMom( 9, 4 )
    real(kind=rk) :: targetMom( 9 ) ! non-equilibrium part of the resulting pdf
    real(kind=rk) :: refX(2)        ! reference values
    real(kind=rk) :: error(4)       ! error values
    integer :: nSourceElems, SE, SW, NE, NW
    ! ---------------------------------------------------------------------------
    passed = .false.
    error  = 0._rk
    nSourceElems = 4
    write(logUnit(1),*)' Performing unit test for interpolation moment '//     &
      &                'gradient routine'
    write(logUnit(1),*)' Current case for child ', iCase

    select case( iCase )
    case( 1, 5)
      SW = 4
      SE = 3
      NW = 2
      NE = 1
    case( 2, 6)
      SW = 2
      SE = 4
      NW = 1
      NE = 3
    case( 3, 7)
      SW = 3
      SE = 1
      NW = 4
      NE = 2
    case( 4, 8)
      SW = 1
      SE = 2
      NW = 3
      NE = 4
    end select

    ! sourceMom(9,4)
    sourceMom = 0._rk
    sourceMom( 1, : ) = 1._rk
    sourceMom( 3, : ) = 0.1_rk
?? if( PUSH ) then
      sourceMom( 6, SW ) = -0.074074074074074_rk  ! SW
      sourceMom( 6, SE ) = -0.037037037037037_rk  ! SE
      sourceMom( 6, NW ) =  0.074074074074074_rk  ! NW
      sourceMom( 6, NE ) =  0.037037037037037_rk  ! NE
?? else
      sourceMom( 6, SW ) =  0.059259259259259_rk  ! SW
      sourceMom( 6, SE ) =  0.029629629629630_rk  ! SE
      sourceMom( 6, NW ) = -0.059259259259259_rk  ! NW
      sourceMom( 6, NE ) = -0.029629629629630_rk  ! NE
?? endif
    refX(1) = 0.05625_rk
    refX(2) = 0.1_rk
    ! targetMom(1:2) = mus_interpolate_quadraticCompact2d_leastSq(               &
    !   &                                      sourceMom(2:3, 1:nSourceElems),   &
    !   &                (/0.25_rk, 0.25_rk/), sourceMom(4:6, 1:nSourceElems),   &
    !   &                nSourceElems, matrix )
    error(1:2) = (refX(1:2)-targetMom(1:2))/refX(1:2)*100._rk
?? if( PUSH ) then
?? else
    write(logUnit(6),*)'X  reference result ', refX(1), ' intp ', targetMom(1),&
      &                ' Error [%] ', error(1)
    write(logUnit(6),*)'Y  reference result ', refX(2), ' intp ', targetMom(2),&
      &                ' Error [%] ', error(2)
?? endif

    refX(1) = 0.075_rk
    refX(2) = 0.1_rk
    ! targetMom(1:2) = mus_interpolate_quadraticCompact2d_leastSq(               &
    !   &                                      sourceMom(2:3, 1:nSourceElems),   &
    !   &                  (/0.5_rk, 0.5_rk/), sourceMom(4:6, 1:nSourceElems),   &
    !   &                  nSourceElems, matrix )
    error(3:4) = (refX(1:2)-targetMom(1:2))/refX(1:2)*100._rk

?? if( PUSH ) then
?? else
    write(logUnit(6),*)'X  reference result ', refX(1), ' intp ', targetMom(1),&
      &                ' Error [%] ', error(3)
    write(logUnit(6),*)'Y  reference result ', refX(2), ' intp ', targetMom(2),&
      &                ' Error [%] ', error(4)
?? endif

    if( abs( maxval( error ) ) > 0.00001_rk ) then
      passed = .false.
    else
      passed = .true.
    endif

  end subroutine unitTest_interpolate_quadraticCompact
! ****************************************************************************** !


! ****************************************************************************** !
  !> Calculate the equilibrium distribution function in all directions
  !!
  !! The equilibrium distribution function is:\n
  !! \[ f^{eq}_i = w_i \rho ( 1 + \frac{\vec c_i \cdot \vec u}{c^2_s}
  !!                      + \frac{ {(\vec c_i \cdot \vec u)}^2}{2c^4_s}
  !!                      - \frac{\vec u \cdot \vec u}{2c^2_s}) \]\n
  !! where \(w_i\) is the weight in each direction,\n
  !! \(\rho\) is the macroscopic value of density,\n
  !! \(c_s\) is the speed of sound,\n
  !! \(\vec c_i\) is the lattice unit velocity in each direction,\n
  !! \(\vec u\) is the macroscopic value of velocity.
  !!
  pure function getEquilibrium( density, velocity, layout ) result( equil )
    ! ---------------------------------------------------------------------------
    !> scheme layout
    type(mus_scheme_layout_type), intent(in) :: layout
    !>
    real(kind=rk), intent(in) :: density
    !>
    real(kind=rk), intent(in) :: velocity(3)
    !> return value
    real(kind=rk)             :: equil(layout%fStencil%QQ)
    ! ---------------------------------------------------------------------------
    ! local variables
    real(kind=rk) :: ucx, usq
    integer :: iDir
    ! ---------------------------------------------------------------------------

    ! square of velocity
    usq = velocity(1) * velocity(1)                                            &
      & + velocity(2) * velocity(2)                                            &
      & + velocity(3) * velocity(3)

    do iDir = 1, layout%fStencil%QQ

      ! velocity times lattice unit velocity
      ucx = layout%fStencil%cxDirRK( 1, iDir )*velocity(1)     &
        & + layout%fStencil%cxDirRK( 2, iDir )*velocity(2)     &
        & + layout%fStencil%cxDirRK( 3, iDir )*velocity(3)

      ! calculate equilibrium density
      equil( iDir ) = layout%weight( iDir ) * density * ( 1.d0 + ucx*cs2inv    &
        &           + ucx*ucx*cs2inv*cs2inv*0.5d0                              &
        &           - usq*cs2inv*0.5d0 )

    enddo

  end function getEquilibrium
! ****************************************************************************** !


end module mus_interpolate_compact_module
! ****************************************************************************** !
