! Copyright (c) 2015-2016 Jiaxing Qi <jiaxing.qi@uni-siegen.de>
! Copyright (c) 2015-2016, 2018-2020 Kannan Masilamani <kannan.masilamani@uni-siegen.de>
! Copyright (c) 2016 Tobias Schneider <tobias1.schneider@student.uni-siegen.de>
! Copyright (c) 2017 Harald Klimach <harald.klimach@uni-siegen.de>
! Copyright (c) 2018 Raphael Haupt <raphael.haupt@uni-siegen.de>
! Copyright (c) 2019-2020 Peter Vitt <peter.vitt2@uni-siegen.de>
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions are met:
!
! 1. Redistributions of source code must retain the above copyright notice,
! this list of conditions and the following disclaimer.
!
! 2. Redistributions in binary form must reproduce the above copyright notice,
! this list of conditions and the following disclaimer in the documentation
! and/or other materials provided with the distribution.
!
! THIS SOFTWARE IS PROVIDED BY THE UNIVERSITY OF SIEGEN “AS IS” AND ANY EXPRESS
! OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
! OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
! IN NO EVENT SHALL UNIVERSITY OF SIEGEN OR CONTRIBUTORS BE LIABLE FOR ANY
! DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
! (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
! ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
! ****************************************************************************** !
!> author: Jiaxing Qi
!! Linear Interpolation of flow quantities between different grid levels
!! ---
!!
!! # Interpolation
!! The routines defined here, fill up the ghost elements with valid data.
!! 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
!!
!! # Workflow
!!
!! Each interpolation routine acts on a list of ghost elements.
!! This list contains pointers to the position in the total list.
!! For each of these ghost elements, the source elements are identified.
!! Before that, the sourceLevel is identified. However, the code is restricted
!! to work with a level jump of only one level, so the sourceLevel is:\n
!! for fillMyGhostsFromFiner: sourceLevel = targetLevel + 1
!! for fillFinerGhostsFromMe: sourceLevel = targetLevel - 1
!!
!! For an overview over implemented interpolation methods, see
!! [Interpolation methods](../page/features/intp_methods.html)
!!
?? include 'header/lbm_d2q9Macros.inc'
?? include 'header/lbm_macros.inc'
?? include 'treelm/source/logMacros.inc'
module mus_interpolate_d2q9_module
  use iso_c_binding, only: 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, div1_2, div1_9, div4_9, &
    &                              div2_9, div1_18, div1_36, div1_4, rho0
  use tem_construction_module, only: tem_levelDesc_type

  use tem_debug_module,      only: dbgUnit
  use tem_varSys_module,     only: tem_varSys_type
  use tem_time_module,       only: tem_time_type
?? IF (loglvl >= 1) THEN
  use tem_debug_module,      only: dbgUnit
  use env_module,            only: long_k
  use tem_topology_module,   only: tem_coordOfID
  use tem_param_module,      only: PI
?? ENDIF

  ! include musubi modules
  use mus_interpolate_tools_module,  only: get_fNeqFac_f2c, &
    &                                      get_fNeqFac_c2f
  use mus_interpolate_header_module, only: mus_interpolation_method_type
  use mus_interpolate_quadratic_module, only: mus_interpolate_quad2D_leastSq
  use mus_scheme_layout_module,      only: mus_scheme_layout_type
  use mus_physics_module,            only: mus_physics_type
  use mus_moments_module,             only: get_moment
  use mus_derivedQuantities_module2,  only: getNonEqFac
  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 :: mus_interpolate_compact2D

  ! used by linear interpolation
  public :: fillFinerGhostsFromMe_LinearIncomp_2D
  ! used by quadratic interpolation
  public :: fillFiner_Quadratic2D_diffusive
  public :: fillFinerGhostsFromMe_quadraticD2Q9
  public :: fillFinerGhostsFromMe_quadraticD2Q9Incomp

  public :: fillCoarser_compact2DIncomp
  public :: fillFiner_compact2DIncomp
  public :: fillFiner_compact2D

  integer, parameter :: QQ = 9 !< number of pdf directions
  real(kind=rk) :: sourceP(   QQ )
  real(kind=rk) :: sourceV(2, QQ )
  real(kind=rk) :: sourceS(3, QQ )
  real(kind=rk) ::  weight(   QQ )

  contains

! ****************************************************************************** !
  !> Fill Finer ghost from Coarser fluid by 2D linear 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_linearIncomp_2D( 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 :: sourcePos      ! source element position in total list
    integer :: targetLevel    ! level of target elements
    integer :: targetPos      ! target element position in total list
    integer :: iElem, indElem, iSourceElem
    integer :: nSourceElems   ! number of source elements for the current target
    ! macroscopic velocities for all source elements
    real(kind=rk), allocatable :: sourceP(:)
    real(kind=rk), allocatable :: sourceV(:,:)
    real(kind=rk), allocatable :: sourcefNeq(:,:)
    real(kind=rk), allocatable :: weight(:)
    ! temporary array
    real(kind=rk) :: f(9), fEq(9), fNeq(9)
    real(kind=rk) :: rho, vel(2), p
?? copy :: var_fEq_d2q9
    real(kind=rk) :: omegaSource, omegaTarget
    real(kind=rk) :: pFac, vFac, sFac, fNeqFac, sNeqPost2Pre, tNeqPre2Post
!    real(kind=rk), parameter :: div1_16 = 1.d0 / 16.d0
    !> debug variables
?? IF (loglvl >= 1) THEN
    integer(kind=long_k) :: treeID
?? ENDIF
?? IF (loglvl >= 10) THEN
    integer :: iDir
?? ENDIF
    type(mus_fluid_type), pointer :: fluid
    integer :: nScalars
    ! --------------------------------------------------------------------------
    fluid => fieldProp(1)%fluid
    nScalars = varSys%nScalars

    sourceLevel = level
    targetLevel = level + 1
    ! KM: Temp fix: Assuming omega is contant
    omegaSource = fluid%viscKine%omLvl( sourceLevel )%val(1)
    omegaTarget = fluid%viscKine%omLvl( targetLevel )%val(1)

    ! macroscopic quantity scale factor from source level to target level
    pFac = physics%pFac( sourceLevel, targetLevel )
    vFac = physics%vFac( sourceLevel, targetLevel )
    sFac = physics%sFac( sourceLevel, targetLevel )
    ! coarse = source, fine = target
    fNeqFac = get_fNeqFac_c2f( omegaSource, omegaTarget, sFac )
    sNeqPost2Pre = 1.0d0 / ?post2pre?( omegaSource )
    tNeqPre2Post = ?pre2post?( omegaTarget )

?? IF (loglvl >= 10) THEN
write(dbgUnit(10), *) ''
write(dbgUnit(10), *) '---- Enter fillFinerGhostsFromMe_linearIncomp_2D -----'
write(dbgUnit(10), *) ''
write(dbgUnit(10), *) '         nMaxSources: ', method%nMaxSources
write(dbgUnit(10), *) '        source level: ', sourceLevel
write(dbgUnit(10), *) '        source omega: ', omegaSource
write(dbgUnit(10), *) 'source fNeq post2pre: ', sNeqPost2Pre
write(dbgUnit(10), *) '        target level: ', targetLevel
write(dbgUnit(10), *) '        target omega: ', omegaTarget
write(dbgUnit(10), *) 'target fNeq pre2post: ', tNeqPre2Post
write(dbgUnit(10), *) '                pFac: ', pFac
write(dbgUnit(10), *) '                vFac: ', vFac
write(dbgUnit(10), *) '                sFac: ', sFac
write(dbgUnit(10), *) '             fNeqFac: ', fNeqFac
write(dbgUnit(10), *) ''
?? ENDIF

    allocate( sourceP(method%nMaxSources)      )
    allocate( sourceV(2,method%nMaxSources)    )
    allocate( sourcefNeq(9,method%nMaxSources) )
    allocate( weight(method%nMaxSources)       )

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

      iElem = targetList( indElem )

      ! 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
        sourcePos = tLevelDesc%depFromCoarser( iElem )%elem%val( iSourceElem )

        ! load pdf
        f(1) = sState( ?IDX?( 1, sourcePos, nScalars, snSize))
        f(2) = sState( ?IDX?( 2, sourcePos, nScalars, snSize))
        f(3) = sState( ?IDX?( 3, sourcePos, nScalars, snSize))
        f(4) = sState( ?IDX?( 4, sourcePos, nScalars, snSize))
        f(5) = sState( ?IDX?( 5, sourcePos, nScalars, snSize))
        f(6) = sState( ?IDX?( 6, sourcePos, nScalars, snSize))
        f(7) = sState( ?IDX?( 7, sourcePos, nScalars, snSize))
        f(8) = sState( ?IDX?( 8, sourcePos, nScalars, snSize))
        f(9) = sState( ?IDX?( 9, sourcePos, nScalars, snSize))

        ! calculate rho and vel
?? copy :: rho_d2q9_a( rho, f )
?? copy :: vel_d2q9_a( vel(1), vel(2), f, rho0 )

        ! calculate feq and fNeq
?? copy :: fEq_d2q9_a( fEq, vel(1), vel(2), rho, rho0 )
        fNeq(:) = ( f(:) - fEq(:) ) * sNeqPost2Pre

             sourceP(iSourceElem) = rho - rho0
           sourceV(:,iSourceElem) = vel(:)
        sourcefNeq(:,iSourceElem) = fNeq(:)

      end do ! iSourceElem

      ! ------------------ LINEAR INTERPOLATION -------------------
      ! linearly interpolate Pressure, Vel, and fNeq
      p       = sum( weight(1:nSourceElems)*sourceP(     1:nSourceElems) )
      vel(1)  = sum( weight(1:nSourceElems)*sourceV(   1,1:nSourceElems) )
      vel(2)  = sum( weight(1:nSourceElems)*sourceV(   2,1:nSourceElems) )
      fNeq(1) = sum( weight(1:nSourceElems)*sourcefNeq(1,1:nSourceElems) )
      fNeq(2) = sum( weight(1:nSourceElems)*sourcefNeq(2,1:nSourceElems) )
      fNeq(3) = sum( weight(1:nSourceElems)*sourcefNeq(3,1:nSourceElems) )
      fNeq(4) = sum( weight(1:nSourceElems)*sourcefNeq(4,1:nSourceElems) )
      fNeq(5) = sum( weight(1:nSourceElems)*sourcefNeq(5,1:nSourceElems) )
      fNeq(6) = sum( weight(1:nSourceElems)*sourcefNeq(6,1:nSourceElems) )
      fNeq(7) = sum( weight(1:nSourceElems)*sourcefNeq(7,1:nSourceElems) )
      fNeq(8) = sum( weight(1:nSourceElems)*sourcefNeq(8,1:nSourceElems) )
      fNeq(9) = sum( weight(1:nSourceElems)*sourcefNeq(9,1:nSourceElems) )

      ! Bilinear interpolation
      ! for each fine ghost, it has 4 sources
      ! C     D
      !
      !   x
      ! A     B
      ! x = ( 9A + 3B + 3C + D ) / 16
      !
      ! p       = (      sourceP(1) * 9.0d0 +      sourceP(2) * 3.0d0 +      sourceP(3) * 3.0d0 +      sourceP(4) ) * div1_16
      ! vel(1)  = (    sourceV(1,1) * 9.0d0 +    sourceV(1,2) * 3.0d0 +    sourceV(1,3) * 3.0d0 +    sourceV(1,4) ) * div1_16
      ! vel(2)  = (    sourceV(2,1) * 9.0d0 +    sourceV(2,2) * 3.0d0 +    sourceV(2,3) * 3.0d0 +    sourceV(2,4) ) * div1_16
      ! fNeq(1) = ( sourcefNeq(1,1) * 9.0d0 + sourcefNeq(1,2) * 3.0d0 + sourcefNeq(1,3) * 3.0d0 + sourcefNeq(1,4) ) * div1_16
      ! fNeq(2) = ( sourcefNeq(2,1) * 9.0d0 + sourcefNeq(2,2) * 3.0d0 + sourcefNeq(2,3) * 3.0d0 + sourcefNeq(2,4) ) * div1_16
      ! fNeq(3) = ( sourcefNeq(3,1) * 9.0d0 + sourcefNeq(3,2) * 3.0d0 + sourcefNeq(3,3) * 3.0d0 + sourcefNeq(3,4) ) * div1_16
      ! fNeq(4) = ( sourcefNeq(4,1) * 9.0d0 + sourcefNeq(4,2) * 3.0d0 + sourcefNeq(4,3) * 3.0d0 + sourcefNeq(4,4) ) * div1_16
      ! fNeq(5) = ( sourcefNeq(5,1) * 9.0d0 + sourcefNeq(5,2) * 3.0d0 + sourcefNeq(5,3) * 3.0d0 + sourcefNeq(5,4) ) * div1_16
      ! fNeq(6) = ( sourcefNeq(6,1) * 9.0d0 + sourcefNeq(6,2) * 3.0d0 + sourcefNeq(6,3) * 3.0d0 + sourcefNeq(6,4) ) * div1_16
      ! fNeq(7) = ( sourcefNeq(7,1) * 9.0d0 + sourcefNeq(7,2) * 3.0d0 + sourcefNeq(7,3) * 3.0d0 + sourcefNeq(7,4) ) * div1_16
      ! fNeq(8) = ( sourcefNeq(8,1) * 9.0d0 + sourcefNeq(8,2) * 3.0d0 + sourcefNeq(8,3) * 3.0d0 + sourcefNeq(8,4) ) * div1_16
      ! fNeq(9) = ( sourcefNeq(9,1) * 9.0d0 + sourcefNeq(9,2) * 3.0d0 + sourcefNeq(9,3) * 3.0d0 + sourcefNeq(9,4) ) * div1_16
      ! ------------------ LINEAR INTERPOLATION -------------------

?? IF (loglvl >= 10) THEN
write(dbgUnit(10),*) '   iElem: ', iElem, ', nSources: ', nSourceElems
write(dbgUnit(10),*) '  weight: ', weight(:)
write(dbgUnit(10),*) 'pressure: ', sourceP(:), p * pFac
write(dbgUnit(10),*) '   vel x: ', sourceV(1,:), vel(1) * vFac
write(dbgUnit(10),*) '   vel y: ', sourceV(2,:), vel(2) * vFac
do iDir = 1, QQ
  write(dbgUnit(10),*) 'fNeq', iDir, ': ', sourcefNeq(iDir,:), fNeq(iDir)*fNeqFac
end do
write(dbgUnit(10),*) ''
?? ENDIF

      ! scale macroscopic quantities (LB unit) from source level to target level
      rho  = p    * pFac + rho0
      vel  = vel  * vFac
      fNeq = fNeq * fNeqFac * tNeqPre2Post

      ! Calculate fEq
?? copy :: fEq_d2q9_a( fEq, vel(1), vel(2), rho, rho0 )

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

      ! Now write pdf (fEq+fNeq) to the target element
      tState(?IDX?( 1, targetPos, nScalars, tnSize))=fEq(1)+fNeq(1)
      tState(?IDX?( 2, targetPos, nScalars, tnSize))=fEq(2)+fNeq(2)
      tState(?IDX?( 3, targetPos, nScalars, tnSize))=fEq(3)+fNeq(3)
      tState(?IDX?( 4, targetPos, nScalars, tnSize))=fEq(4)+fNeq(4)
      tState(?IDX?( 5, targetPos, nScalars, tnSize))=fEq(5)+fNeq(5)
      tState(?IDX?( 6, targetPos, nScalars, tnSize))=fEq(6)+fNeq(6)
      tState(?IDX?( 7, targetPos, nScalars, tnSize))=fEq(7)+fNeq(7)
      tState(?IDX?( 8, targetPos, nScalars, tnSize))=fEq(8)+fNeq(8)
      tState(?IDX?( 9, targetPos, nScalars, tnSize))=fEq(9)+fNeq(9)

    end do ! indElem

?? IF (loglvl >= 1) THEN
write(dbgUnit(1),*) '----- Leave fillFinerGhostsFromMe_linearIncomp_2D -----'
write(dbgUnit(1),*) ''
?? ENDIF

    deallocate( sourceP   )
    deallocate( sourceV   )
    deallocate( sourcefNeq)

  end subroutine fillFinerGhostsFromMe_linearIncomp_2D
! ****************************************************************************** !

! ****************************************************************************** !
  !> Fill Finer ghost from Coarser fluid by 2D quadratic interpolation
  !! This routine serves as a HIGH order interpolation, which means every ghost
  !! should have 9 source elements. It only works for D2Q9.
  !!
  !! 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 fillFiner_quadratic2D_diffusive( 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 :: targetLevel    ! level of target elements
    integer :: targetPos      ! target element position in total list
    integer :: iElem, indElem, ii
    real(kind=rk) :: targetCoord(2)
    real(kind=rk) :: fEq(9), fNeq(9), f(9)
    real(kind=rk) :: rho, vel(2), S(2,2), sourceVals(6,9), targetVals(6)
?? copy :: var_fEq_d2q9
    real(kind=rk) :: omegaSource, omegaTarget
    real(kind=rk) :: pFac, vFac, sFac
    real(kind=rk) :: sBuf(6,tLevelDesc%sourceFromCoarser%nVals)
    !> debug variables
?? IF (loglvl >= 1) THEN
    integer(kind=long_k) :: treeID
?? ENDIF
?? IF (loglvl >= 10) THEN
    integer :: iDir
?? ENDIF
    type(mus_fluid_type), pointer :: fluid
    integer :: nScalars
    ! --------------------------------------------------------------------------
    fluid => fieldProp(1)%fluid
    nScalars = varSys%nScalars

    sourceLevel = level
    targetLevel = level + 1
    ! KM: Temp fix: Assuming omega is contant
    omegaSource = fluid%viscKine%omLvl( sourceLevel )%val(1)
    omegaTarget = fluid%viscKine%omLvl( targetLevel )%val(1)

    ! macroscopic quantity scale factor from source level to target level
    pFac = physics%pFac( sourceLevel, targetLevel )
    vFac = physics%vFac( sourceLevel, targetLevel )
    sFac = physics%sFac( sourceLevel, targetLevel )

?? IF (loglvl >= 10) THEN
write(dbgUnit(10), *) ''
write(dbgUnit(10), *) '---- Enter fillFiner_quadratic2D_diffusive -----'
write(dbgUnit(10), *) ''
write(dbgUnit(10), *) '                pFac: ', pFac
write(dbgUnit(10), *) '                vFac: ', vFac
write(dbgUnit(10), *) '                sFac: ', sFac
write(dbgUnit(10), *) '             fNeqFac: ', fNeqFac
write(dbgUnit(10), *) '        source level: ', sourceLevel
write(dbgUnit(10), *) '               omega: ', omegaSource
! write(dbgUnit(10), *) '       fNeq post2pre: ', sNeqPost2Pre
write(dbgUnit(10), *) '        target level: ', targetLevel
write(dbgUnit(10), *) '               omega: ', omegaTarget
write(dbgUnit(10), *) '       fNeq pre2post: ', tNeqPre2Post
write(dbgUnit(10), *) ''
?? ENDIF

    ! fill buffer for all sources  ----------------------------------
    do iElem = 1, tLevelDesc%sourceFromCoarser%nVals
      do ii = 1, QQ
        f(ii) = sState( ?IDX?( ii, tLevelDesc%sourceFromCoarser%val(iElem), nScalars, snSize ))
      end do

      rho = sum( f )
?? copy :: vel_d2q9_a( vel(1), vel(2), f, rho0 )
      S = getStrain( f, rho, vel, omegaSource, layout%fStencil%cxDir )

      sBuf(1,iElem) = rho - rho0
      sBuf(2,iElem) = vel(1)
      sBuf(3,iElem) = vel(2)
      sBuf(4,iElem) = S(1,1)
      sBuf(5,iElem) = S(2,2)
      sBuf(6,iElem) = S(1,2)

    end do ! nVals
    ! fill buffer for all sources  ----------------------------------

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

      iElem = targetList( indElem )
      do ii = 1, QQ
        sourceVals(:,ii) = sBuf(:,tLevelDesc%depFromCoarser(iElem)%elemBuffer%val(ii))
      end do

      ! ------------------ QUADRATIC INTERPOLATION -------------------
      ! get the target coordinates
      targetCoord = tLevelDesc%depFromCoarser( iElem )%coord(:2)
      call tem_abort('Deactivated: quadratic 2D for diffusive')
!KM!      targetVals = mus_interpolate_quadratic2d_leastSq(          &
!KM!        &            sourceVals(:,:), targetCoord(:), (1+2+3), 9 )
      ! ------------------ QUADRATIC INTERPOLATION -------------------

?? IF (loglvl >= 10) THEN
write(dbgUnit(10),*) '   iElem: ', iElem
write(dbgUnit(10),*) 'pressure: ', sourceVals(1,:), targetVals(1) * pFac
write(dbgUnit(10),*) '   vel x: ', sourceVals(2,:), targetVals(2) * vFac
write(dbgUnit(10),*) '   vel y: ', sourceVals(3,:), targetVals(3) * vFac
do iDir = 1, QQ
  write(dbgUnit(10),*) 'fNeq', iDir, ': ', sourceVals(iDir+3,:), &
    &                  targetVals(iDir+3)*fNeqFac
end do
write(dbgUnit(10),*) ''
?? ENDIF

      ! scale macroscopic quantities (LB unit) from source level to target level
      rho  = targetVals(1)   * pFac + rho0
      vel  = targetVals(2:3) * vFac
      S(1,1)  = targetVals(4) * sFac
      S(2,2)  = targetVals(5) * sFac
      S(1,2)  = targetVals(6) * sFac
      S(2,1)  = S(1,2)
      fNeq = getNEq( S, omegaTarget, layout%fStencil%cxDirRK, layout%weight )

      ! Calculate fEq
?? copy :: fEq_d2q9_a( fEq, vel(1), vel(2), rho, rho0 )

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

      ! Now write pdf (fEq+fNeq) to the target element
      tState(?IDX?( 1, targetPos, nScalars, tnSize))=fEq(1)+fNeq(1)
      tState(?IDX?( 2, targetPos, nScalars, tnSize))=fEq(2)+fNeq(2)
      tState(?IDX?( 3, targetPos, nScalars, tnSize))=fEq(3)+fNeq(3)
      tState(?IDX?( 4, targetPos, nScalars, tnSize))=fEq(4)+fNeq(4)
      tState(?IDX?( 5, targetPos, nScalars, tnSize))=fEq(5)+fNeq(5)
      tState(?IDX?( 6, targetPos, nScalars, tnSize))=fEq(6)+fNeq(6)
      tState(?IDX?( 7, targetPos, nScalars, tnSize))=fEq(7)+fNeq(7)
      tState(?IDX?( 8, targetPos, nScalars, tnSize))=fEq(8)+fNeq(8)
      tState(?IDX?( 9, targetPos, nScalars, tnSize))=fEq(9)+fNeq(9)

    end do ! indElem

?? IF (loglvl >= 1) THEN
write(dbgUnit(1),*) '----- Leave fillFiner_quadratic2D_diffusive -----------'
write(dbgUnit(1),*) ''
?? ENDIF

  end subroutine fillFiner_quadratic2D_diffusive
! ****************************************************************************** !


  ! ************************************************************************** !
  !> Fill fine ghost from coarse fluid by quadratic interpolation for D2Q9
  !! stencil.
  !! 1. Compute moments for all source elements, save in momBuf
  !! 2. For each target, interpolate moments (den, vel, tau)
  !!    (10 moments for 3D and 6 moments for 2D)
  !! 3. calculate fEq and use it to calculate high order moments
  !! 4. convert moments to PDF
  !! This routine is used by acoustic quadratic 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_quadraticD2Q9( 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, iElem, indElem
    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
    ! moments of the source elements' pdf
    real(kind=rk) :: sourceMom( 6, method%nMaxSources )
    ! shear stress scaling factor
    real(kind=rk) :: fac, inv_rho, rho, vel(2)
    real(kind=rk) :: tMom( layout%fStencil%QQ )  ! temp moment calculation
    integer :: QQ
    integer :: posInIntpMatLSF
?? copy :: var_fEq_d2q9
    type(mus_fluid_type), pointer :: fluid
    integer :: nScalars
    real(kind=rk), allocatable :: momBuf(:,:)
    ! ---------------------------------------------------------------------------
    fluid => fieldProp(1)%fluid
    QQ = layout%fStencil%QQ
    nScalars = varSys%nScalars

    sourceLevel = level
    targetLevel = level + 1
    ! KM: Temp fix: Assuming omega is constant
    fac = getNoneqFac( fluid%viscKine%omLvl( sourceLevel )%val(1), &
      &                fluid%viscKine%omLvl( targetLevel )%val(1)  )

    allocate( momBuf( 6, 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 iVal = 1, QQ
        f(iVal) = sState(?IDX?( iVal, sourceElem, nScalars, snSize ))
      enddo

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

      ! Transfer momentum to velocity
      tMom(2:3) = tMom(2:3) * inv_rho
      ! 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( 4) = tMom( 4) - tMom(2)*tMom(2)*tMom(1) - cs2*tMom(1) !Pxx
      tMom( 5) = tMom( 5) - tMom(3)*tMom(3)*tMom(1) - cs2*tMom(1) !Pyy
      tMom( 6) = tMom( 6) - tMom(2)*tMom(3)*tMom(1)               !Pxy

      ! Store all the moments to the current source element's position in the
      ! buffer
      momBuf( 1:6, iSourceElem ) = tMom(1:6)
    enddo

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

      iElem = targetList( indElem )
      targetElem = iElem + tLevelDesc%offset( 1, eT_ghostFromCoarser)
      nSourceElems = tLevelDesc%depFromCoarser( iElem )%elem%nVals
      posInIntpMatLSF = tLevelDesc%depFromCoarser( iElem )%posInIntpMatLSF

      do iSourceElem = 1, nSourceElems
        ! Get the source element in the levelwise buffer
        sourceElem = tLevelDesc%depFromCoarser( iElem )%elemBuffer%val( iSourceElem )
        ! and read its entries
        sourceMom(1:6, iSourceElem) = momBuf(1:6, sourceElem)
      end do

      tMom(1:6) = mus_interpolate_quad2D_leastSq(                   &
        &   srcMom      = sourceMom(1:6, 1:nSourceElems),           &
        &   targetCoord = tLevelDesc%depFromCoarser( iElem )%coord, &
        &   LSFmat      = method%intpMat_forLSF%matArray            &
        &                       %val(posInIntpMatLSF),              &
        &   nSources    = nSourceElems,                             &
        &   nVals       = 6 )

      inv_rho = 1.0_rk / tMom(1)
      ! Transfer velocity to momentum
      tMom(2:3) = tMom(2:3)*tMom(1)

      ! Rescale the shear stress moments
      ! and add back the pi(0) part to the tensor
      tMom( 4) = tMom( 4)*fac + tMom(2)*tMom(2)*inv_rho + cs2*tMom(1)
      tMom( 5) = tMom( 5)*fac + tMom(3)*tMom(3)*inv_rho + cs2*tMom(1)
      tMom( 6) = tMom( 6)*fac + tMom(2)*tMom(3)*inv_rho

      rho = tMom(1)
      vel(1:2) = tMom(2:3) * inv_rho

?? copy :: fEq_d2q9_a( f, vel(1), vel(2), rho, rho )

      tMom(7:QQ) = matmul( layout%moment%toMoments%A(7:QQ,:), f )

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

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

    end do ! indElem

  end subroutine fillFinerGhostsFromMe_quadraticD2Q9
! ****************************************************************************** !

  ! ************************************************************************** !
  !> Fill fine ghost from coarse fluid by quadratic interpolation for D2Q9
  !! stencil.
  !! 1. Compute moments for all source elements, save in momBuf
  !! 2. For each target, interpolate moments (den, vel, tau)
  !!    (10 moments for 3D and 6 moments for 2D)
  !! 3. calculate fEq and use it to calculate high order moments
  !! 4. convert moments to PDF
  !! This routine is used by acoustic quadratic interpolation for
  !! incompressible.
  !!
  !! 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_quadraticD2Q9Incomp( 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, iElem, indElem
    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
    ! moments of the source elements' pdf
    real(kind=rk) :: sourceMom( 6, method%nMaxSources )
    ! shear stress scaling factor
    real(kind=rk) :: fac
    real(kind=rk) :: tMom( layout%fStencil%QQ )  ! temp moment calculation
    integer :: QQ
    integer :: posInIntpMatLSF
?? copy :: var_fEq_d2q9
    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
    fac = getNoneqFac( fluid%viscKine%omLvl( sourceLevel )%val(1), &
      &                fluid%viscKine%omLvl( targetLevel )%val(1)  )

    allocate( momBuf( 6, 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 iVal = 1, QQ
        f(iVal) = sState(?IDX?( iVal, sourceElem, nScalars, 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( 4) = tMom( 4) - tMom(2)*tMom(2) - cs2*tMom(1) !Pxx
      tMom( 5) = tMom( 5) - tMom(3)*tMom(3) - cs2*tMom(1) !Pyy
      tMom( 6) = tMom( 6) - tMom(2)*tMom(3)               !Pxy

      ! Store all the moments to the current source element's position in the
      ! buffer
      momBuf( 1:6, iSourceElem ) = tMom(1:6)
    enddo

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

      iElem = targetList( indElem )
      targetElem = iElem + tLevelDesc%offset( 1, eT_ghostFromCoarser)
      nSourceElems = tLevelDesc%depFromCoarser( iElem )%elem%nVals
      posInIntpMatLSF = tLevelDesc%depFromCoarser( iElem )%posInIntpMatLSF

      do iSourceElem = 1, nSourceElems
        ! Get the source element in the levelwise buffer
        sourceElem = tLevelDesc%depFromCoarser( iElem )%elemBuffer%val( iSourceElem )
        ! and read its entries
        sourceMom(1:6, iSourceElem) = momBuf(1:6, sourceElem)
      end do

      tMom(1:6) = mus_interpolate_quad2D_leastSq(                   &
        &   srcMom      = sourceMom(1:6, 1:nSourceElems),           &
        &   targetCoord = tLevelDesc%depFromCoarser( iElem )%coord, &
        &   LSFmat      = method%intpMat_forLSF%matArray            &
        &                       %val(posInIntpMatLSF),              &
        &   nSources    = nSourceElems,                             &
        &   nVals       = 6 )

      ! Rescale the shear stress moments
      ! and add back the pi(0) part to the tensor
      tMom( 4) = tMom( 4)*fac + tMom(2)*tMom(2) + cs2*tMom(1)
      tMom( 5) = tMom( 5)*fac + tMom(3)*tMom(3) + cs2*tMom(1)
      tMom( 6) = tMom( 6)*fac + tMom(2)*tMom(3)

?? copy :: fEq_d2q9_a( f, tMom(2), tMom(3), tMom(1), rho0 )

      tMom(7:QQ) = matmul( layout%moment%toMoments%A(7:QQ,:), f )

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

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

    end do ! indElem

  end subroutine fillFinerGhostsFromMe_quadraticD2Q9Incomp
! ****************************************************************************** !

! ****************************************************************************** !


! ****************************************************************************** !
  !> Given the velocity and strain rate (S) at four corner,
  !! interpolate the velocity at a given position (cx, cy)
  !! vel and S contain values at the following position:
  !! [0,0], [1,0], [0,1], [1,1]
  !! Checked!
  pure function mus_interpolate_Compact2d( vel, S, cx, cy ) result( phi )
    ! ---------------------------------------------------------------------------
    !> velocity (Ux, Uy) of 4 source elements
    real(kind=rk), intent(in) :: vel(2,4)
    !> Strain rate (Sxx,Syy,Sxy) of 4 source elements
    real(kind=rk), intent(in) :: S(3,4)
    !> interpolation location within the square
    !! it should be in the range of [0,1]
    real(kind=rk), intent(in) :: cx, cy
    !> interpolated velocity
    real(kind=rk) :: phi(2)
    ! ---------------------------------------------------------------------------
    real(kind=rk) :: a0,ax,ay,axx,ayy,axy
    real(kind=rk) :: b0,bx,by,bxx,byy,bxy
    real(kind=rk) :: Ux(4)
    real(kind=rk) :: Uy(4)
    real(kind=rk) :: Sxxyy(4) ! Sxx - Syy
    real(kind=rk) :: Sxy(4)
    ! ---------------------------------------------------------------------------

    ! FILL RIGHT HAND SIDE
    ux(:)    = vel(1,:)
    uy(:)    = vel(2,:)
    Sxxyy(:) = S(1,:) - S(2,:)
    Sxy(:)   = S(3,:)

    ! ---------------------------------------------------------------------------
    a0  = ux(1)
    axy = ux(1) + ux(4) - ux(2) - ux(3)
    b0  = uy(1)
    bxy = uy(1) + uy(4) - uy(2) - uy(3)

    ayy = ( Sxy(4) + Sxy(3) - Sxy(2) - Sxy(1) - bxy ) * 0.5_rk
    bxx = ( Sxy(4) + Sxy(2) - Sxy(3) - Sxy(1) - axy ) * 0.5_rk

    axx = ( Sxxyy(4) + Sxxyy(2) - Sxxyy(3) - Sxxyy(1) ) * 0.25_rk + bxy * 0.5_rk
    byy = ( Sxxyy(1) + Sxxyy(2) - Sxxyy(3) - Sxxyy(4) ) * 0.25_rk + axy * 0.5_rk

    ax = ux(2) - ux(1) - axx
    ay = ux(3) - ux(1) - ayy
    bx = uy(2) - uy(1) - bxx
    by = uy(3) - uy(1) - byy
    ! ---------------------------------------------------------------------------

    ! Evaluate the bubble function with the above calculated coefficients
    phi(1) = ux(1) + ax*cx + ay*cy + axx*cx*cx + ayy*cy*cy + axy*cx*cy
    phi(2) = uy(1) + bx*cx + by*cy + bxx*cx*cx + byy*cy*cy + bxy*cx*cy

  end function mus_interpolate_compact2D
! ****************************************************************************** !

! ****************************************************************************** !
  !> 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 fillCoarser_compact2DIncomp( 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 :: targetLevel    ! level of target elements
    integer :: posInTotal     ! source element position in total list
    integer :: iElem          ! current target element (for outer loop)
    integer :: indElem        ! element counter for indirection list
    integer :: iSource        ! current source element (for inner loop)
    ! macroscopic velocities for all source elements
    real(kind=rk) :: f(9), fEq(9), fNeq(9) !, tempNeq(9)
    real(kind=rk) :: rho, vel(2), p, S(2,2) !,Sxx, Syy, Sxy,
?? copy :: var_fEq_d2q9
    real(kind=rk) :: omegaSource, omegaTarget
    real(kind=rk) :: pFac, vFac, sFac
    !> debug variables
?? IF (loglvl >= 1) THEN
    integer(kind=long_k) :: targetID
    integer :: coord(4)
    real(kind=rk) :: targetCoord(3), dx
    real(kind=rk) :: mom(QQ)
?? ENDIF
    integer :: nScalars
    type(mus_fluid_type), pointer :: fluid
    ! ---------------------------------------------------------------------------
    fluid => fieldProp(1)%fluid
    nScalars = varSys%nScalars

    sourceLevel = level + 1
    targetLevel = level
    ! KM: Temp fix: Assuming omega is contant
    omegaSource = fluid%viscKine%omLvl( sourceLevel )%val(1)
    omegaTarget = fluid%viscKine%omLvl( targetLevel )%val(1)
    ! pressure and velocity conversion factor between levels
    pFac = physics%pfac( sourceLevel, targetLevel )
    vFac = physics%vfac( sourceLevel, targetLevel )
    sFac = physics%sfac( sourceLevel, targetLevel )

?? IF (loglvl >= 10) THEN
write(dbgUnit(10), *) ''
write(dbgUnit(10), *) '---- Enter fillCoarser_Compact2D --------'
write(dbgUnit(10), *) ''
write(dbgUnit(10), *) '                pFac: ', pFac
write(dbgUnit(10), *) '                vFac: ', vFac
write(dbgUnit(10), *) '                sFac: ', sFac
write(dbgUnit(10), *) '             fNeqFac: ', fNeqFac
write(dbgUnit(10), *) '        source level: ', sourceLevel
write(dbgUnit(10), *) '               omega: ', omegaSource
write(dbgUnit(10), *) '       fNeq post2pre: ', sNeqPost2Pre
write(dbgUnit(10), *) '        target level: ', targetLevel
write(dbgUnit(10), *) '               omega: ', omegaTarget
write(dbgUnit(10), *) '       fNeq pre2post: ', tNeqPre2Post
write(dbgUnit(10), *) ''
?? ENDIF

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

      iElem = targetList( indElem )

      ! Find out how many source elements we have for interpolation.
      ! nSourceElems = levelDesc( level )%depFromFiner( iElem )%elem%nVals

      ! ------------------------------------------------------------------------
      ! Now loop over all fine source elements for this target:
      ! Every ghost should have exactly 4 source elements.
      p = 0._rk
      do iSource = 1, 4

        posInTotal = tLevelDesc%depFromFiner( iElem )%elem%val( iSource )

        ! load pdf
        f(1) = sState( ?IDX?( 1, posInTotal, nScalars, snSize))
        f(2) = sState( ?IDX?( 2, posInTotal, nScalars, snSize))
        f(3) = sState( ?IDX?( 3, posInTotal, nScalars, snSize))
        f(4) = sState( ?IDX?( 4, posInTotal, nScalars, snSize))
        f(5) = sState( ?IDX?( 5, posInTotal, nScalars, snSize))
        f(6) = sState( ?IDX?( 6, posInTotal, nScalars, snSize))
        f(7) = sState( ?IDX?( 7, posInTotal, nScalars, snSize))
        f(8) = sState( ?IDX?( 8, posInTotal, nScalars, snSize))
        f(9) = sState( ?IDX?( 9, posInTotal, nScalars, snSize))

        ! calculate rho and take average
        rho = sum( f )
        p = p + ( rho - rho0 )

        ! calculate vel
?? copy :: vel_d2q9_a( vel(1), vel(2), f, rho0 )
        sourceV(1,iSource) = vel(1)
        sourceV(2,iSource) = vel(2)

        S = getStrain( f, rho, vel, omegaSource, layout%fStencil%cxDir )
        sourceS(1,iSource) = S(1,1)
        sourceS(2,iSource) = S(2,2)
        sourceS(3,iSource) = S(1,2)

      end do ! iSource
      ! ------------------------------------------------------------------------

      ! ------------------ Compact INTERPOLATION -------------------------------
      ! pressure and fNeq by average
      ! velocity: quadratic compact interpolate
      p = p * div1_4
      vel(:) = mus_interpolate_compact2d( sourceV(:,1:4),  &
        &                                 sourceS(:,1:4),  &
        &                                 0.5_rk, 0.5_rk )
      S(1,1) = sum(sourceS(1,1:4)) * div1_4
      S(2,2) = sum(sourceS(2,1:4)) * div1_4
      S(1,2) = sum(sourceS(3,1:4)) * div1_4
      S(2,1) = S(1,2)
      ! ------------------ Compact INTERPOLATION -------------------------------

      ! Scale quantities on source level to target level
      rho = p   * pFac + rho0
      vel = vel * vFac
      S   = s   * sFac

      ! Calculate fEq
?? copy :: fEq_d2q9_a( fEq, vel(1), vel(2), rho, rho0 )
      fNeq = getNEq( S, omegaTarget, layout%fStencil%cxDirRK, layout%weight )

      ! Get its position in state array
      posInTotal = iElem + tLevelDesc%offset( 1, eT_ghostFromFiner)

      ! Now write the resulting pdf in the current direction to the target
      ! Element position resulting pdf = f_eq + f_neq
      tState(?IDX?( 1, posInTotal, nScalars, tnSize))=fEq(1)+fNeq(1)
      tState(?IDX?( 2, posInTotal, nScalars, tnSize))=fEq(2)+fNeq(2)
      tState(?IDX?( 3, posInTotal, nScalars, tnSize))=fEq(3)+fNeq(3)
      tState(?IDX?( 4, posInTotal, nScalars, tnSize))=fEq(4)+fNeq(4)
      tState(?IDX?( 5, posInTotal, nScalars, tnSize))=fEq(5)+fNeq(5)
      tState(?IDX?( 6, posInTotal, nScalars, tnSize))=fEq(6)+fNeq(6)
      tState(?IDX?( 7, posInTotal, nScalars, tnSize))=fEq(7)+fNeq(7)
      tState(?IDX?( 8, posInTotal, nScalars, tnSize))=fEq(8)+fNeq(8)
      tState(?IDX?( 9, posInTotal, nScalars, tnSize))=fEq(9)+fNeq(9)

    end do ! indElem

?? IF (loglvl >= 1) THEN
write(dbgUnit(1), *) '---- Leave fillCoarser_Compact2D -----------'
write(dbgUnit(1), *) ''
?? ENDIF

  end subroutine fillCoarser_compact2DIncomp
! ****************************************************************************** !

! ****************************************************************************** !
  !> Fill Finer ghost from Coarser fluid by 2D linear 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 fillFiner_compact2DIncomp( 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 :: sourcePos      ! source element position in total list
    integer :: targetLevel    ! level of target elements
    integer :: targetPos      ! target element position in total list
    integer :: iElem, indElem, iSourceElem
    integer :: nSourceElems   ! number of source elements for the current target
    ! temporary array
    real(kind=rk) :: f(QQ), fEq(QQ), fNeq(QQ)
    real(kind=rk) :: rho, vel(2), p, S(2,2) !, Sxx, Syy, Sxy
?? copy :: var_fEq_d2q9
    real(kind=rk) :: omegaSource, omegaTarget
    real(kind=rk) :: pFac, vFac, sFac !, fNeqFac, sNeqPost2Pre, tNeqPre2Post
    !> debug variables
?? IF (loglvl >= 1) THEN
    integer(kind=long_k) :: treeID
?? ENDIF
?? IF (loglvl >= 10) THEN
    integer :: iDir
?? ENDIF
    type(mus_fluid_type), pointer :: fluid
    integer :: nScalars
    ! --------------------------------------------------------------------------
    fluid => fieldProp(1)%fluid
    nScalars = varSys%nScalars

    sourceLevel = level
    targetLevel = level + 1
    ! KM: Temp fix: Assuming omega is contant
    omegaSource = fluid%viscKine%omLvl( sourceLevel )%val(1)
    omegaTarget = fluid%viscKine%omLvl( targetLevel )%val(1)

    ! macroscopic quantity scale factor from source level to target level
    pFac = physics%pFac( sourceLevel, targetLevel )
    vFac = physics%vFac( sourceLevel, targetLevel )
    sFac = physics%sFac( sourceLevel, targetLevel )

?? IF (loglvl >= 10) THEN
write(dbgUnit(10), *) ''
write(dbgUnit(10), *) '---- Enter fillFiner_compact2D -----'
write(dbgUnit(10), *) ''
write(dbgUnit(10), *) ' nMaxSources: ', method%nMaxSources
write(dbgUnit(10), *) 'source level: ', sourceLevel
write(dbgUnit(10), *) 'source omega: ', omegaSource
write(dbgUnit(10), *) 'source fNeq post2pre: ', sNeqPost2Pre
write(dbgUnit(10), *) 'target level: ', targetLevel
write(dbgUnit(10), *) 'target omega: ', omegaTarget
write(dbgUnit(10), *) 'target fNeq pre2post: ', tNeqPre2Post
write(dbgUnit(10), *) 'pFac: ', pFac
write(dbgUnit(10), *) 'vFac: ', vFac
write(dbgUnit(10), *) 'sFac: ', sFac
write(dbgUnit(10), *) 'fNeqFac: ', fNeqFac
write(dbgUnit(10), *) ''
?? ENDIF

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

      iElem = targetList( indElem )

      ! 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
        sourcePos = tLevelDesc%depFromCoarser( iElem )%elem%val( iSourceElem )

        ! load pdf
        f(1) = sState( ?IDX?( 1, sourcePos, nScalars, snSize))
        f(2) = sState( ?IDX?( 2, sourcePos, nScalars, snSize))
        f(3) = sState( ?IDX?( 3, sourcePos, nScalars, snSize))
        f(4) = sState( ?IDX?( 4, sourcePos, nScalars, snSize))
        f(5) = sState( ?IDX?( 5, sourcePos, nScalars, snSize))
        f(6) = sState( ?IDX?( 6, sourcePos, nScalars, snSize))
        f(7) = sState( ?IDX?( 7, sourcePos, nScalars, snSize))
        f(8) = sState( ?IDX?( 8, sourcePos, nScalars, snSize))
        f(9) = sState( ?IDX?( 9, sourcePos, nScalars, snSize))

        ! calculate rho and vel
        rho = sum( f )
        sourceP(iSourceElem) = rho - rho0

?? copy :: vel_d2q9_a( vel(1), vel(2), f, rho0 )
        sourceV(:,iSourceElem) = vel(:)

        S = getStrain( f, rho, vel, omegaSource, layout%fStencil%cxDir )
        sourceS(1,iSourceElem) = S(1,1)
        sourceS(2,iSourceElem) = S(2,2)
        sourceS(3,iSourceElem) = S(1,2)

      end do ! iSourceElem

      ! ------------------ LINEAR INTERPOLATION -------------------
      ! JQ: taking pressure from parent gave better results than linear
      ! linearly interpolate Pressure, Vel, and fNeq
      p      = sum( weight(1:nSourceElems)*sourceP(  1:nSourceElems) )
      S(1,1) = sum( weight(1:nSourceElems)*sourceS(1,1:nSourceElems) )
      S(2,2) = sum( weight(1:nSourceElems)*sourceS(2,1:nSourceElems) )
      S(1,2) = sum( weight(1:nSourceElems)*sourceS(3,1:nSourceElems) )
      S(2,1) = S(1,2)

      ! interpolate vel by compact
      vel(:) = mus_interpolate_compact2d( sourceV(1:2,1:nSourceElems),  &
        &                                 sourceS(1:3,1:nSourceElems),  &
        &                  tLevelDesc%depFromCoarser( iElem )%coord(1), &
        &                  tLevelDesc%depFromCoarser( iElem )%coord(2)  )
      ! ------------------ LINEAR INTERPOLATION -------------------

?? IF (loglvl >= 10) THEN
write(dbgUnit(10),*) '   iElem: ', iElem, ', nSources: ', nSourceElems
write(dbgUnit(10),*) '  weight: ', weight(:)
write(dbgUnit(10),*) 'pressure: ', sourceP(:), p * pFac
write(dbgUnit(10),*) '   vel x: ', sourceV(1,:), vel(1) * vFac
write(dbgUnit(10),*) '   vel y: ', sourceV(2,:), vel(2) * vFac
do iDir = 1, QQ
  write(dbgUnit(10),*) 'fNeq', iDir, ': ', sourcefNeq(iDir,:), fNeq(iDir)*fNeqFac
end do
write(dbgUnit(10),*) ''
?? ENDIF

      ! scale macroscopic quantities (LB unit) from source level to target level
      rho = p   * pFac + rho0
      vel = vel * vFac
      S   = s   * sFac

      ! Calculate fEq
?? copy :: fEq_d2q9_a( fEq, vel(1), vel(2), rho, rho0 )
      fNeq = getNEq( S, omegaTarget, layout%fStencil%cxDirRK, layout%weight )

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

      ! Now write pdf (fEq+fNeq) to the target element
      tState(?IDX?( 1, targetPos, nScalars, tnSize))=fEq(1)+fNeq(1)
      tState(?IDX?( 2, targetPos, nScalars, tnSize))=fEq(2)+fNeq(2)
      tState(?IDX?( 3, targetPos, nScalars, tnSize))=fEq(3)+fNeq(3)
      tState(?IDX?( 4, targetPos, nScalars, tnSize))=fEq(4)+fNeq(4)
      tState(?IDX?( 5, targetPos, nScalars, tnSize))=fEq(5)+fNeq(5)
      tState(?IDX?( 6, targetPos, nScalars, tnSize))=fEq(6)+fNeq(6)
      tState(?IDX?( 7, targetPos, nScalars, tnSize))=fEq(7)+fNeq(7)
      tState(?IDX?( 8, targetPos, nScalars, tnSize))=fEq(8)+fNeq(8)
      tState(?IDX?( 9, targetPos, nScalars, tnSize))=fEq(9)+fNeq(9)

    end do ! indElem

write(dbgUnit(1),*) ''
?? IF (loglvl >= 1) THEN
write(dbgUnit(1),*) '----- Leave fillFiner_compact2D -----------'
write(dbgUnit(1),*) ''
?? ENDIF

  end subroutine fillFiner_compact2DIncomp
! ****************************************************************************** !

! ****************************************************************************** !
  !> Fill Finer ghost from Coarser fluid by 2D linear 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 fillFiner_compact2D( 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 :: sourcePos      ! source element position in total list
    integer :: targetLevel    ! level of target elements
    integer :: targetPos      ! target element position in total list
    integer :: iElem, indElem, iSourceElem
    integer :: nSourceElems   ! number of source elements for the current target
    ! temporary array
    real(kind=rk) :: f(QQ), fEq(QQ), fNeq(QQ)
    real(kind=rk) :: rho, inv_rho, vel(2), p, S(2,2) !, Sxx, Syy, Sxy
?? copy :: var_fEq_d2q9
    real(kind=rk) :: omegaSource, omegaTarget
    real(kind=rk) :: pFac, vFac, sFac !, fNeqFac, sNeqPost2Pre, tNeqPre2Post
    !> debug variables
?? IF (loglvl >= 1) THEN
    integer(kind=long_k) :: treeID
?? ENDIF
?? IF (loglvl >= 10) THEN
    integer :: iDir
?? ENDIF
    type(mus_fluid_type), pointer :: fluid
    integer :: nScalars
    ! --------------------------------------------------------------------------
    fluid => fieldProp(1)%fluid
    nScalars = varSys%nScalars

    sourceLevel = level
    targetLevel = level + 1
    ! KM: Temp fix: Assuming omega is contant
    omegaSource = fluid%viscKine%omLvl( sourceLevel )%val(1)
    omegaTarget = fluid%viscKine%omLvl( targetLevel )%val(1)

    ! macroscopic quantity scale factor from source level to target level
    pFac = physics%pFac( sourceLevel, targetLevel )
    vFac = physics%vFac( sourceLevel, targetLevel )
    sFac = physics%sFac( sourceLevel, targetLevel )

?? IF (loglvl >= 10) THEN
write(dbgUnit(10), *) ''
write(dbgUnit(10), *) '---- Enter fillFiner_compact2D -----'
write(dbgUnit(10), *) ''
write(dbgUnit(10), *) ' nMaxSources: ', method%nMaxSources
write(dbgUnit(10), *) 'source level: ', sourceLevel
write(dbgUnit(10), *) 'source omega: ', omegaSource
write(dbgUnit(10), *) 'source fNeq post2pre: ', sNeqPost2Pre
write(dbgUnit(10), *) 'target level: ', targetLevel
write(dbgUnit(10), *) 'target omega: ', omegaTarget
write(dbgUnit(10), *) 'target fNeq pre2post: ', tNeqPre2Post
write(dbgUnit(10), *) 'pFac: ', pFac
write(dbgUnit(10), *) 'vFac: ', vFac
write(dbgUnit(10), *) 'sFac: ', sFac
write(dbgUnit(10), *) 'fNeqFac: ', fNeqFac
write(dbgUnit(10), *) ''
?? ENDIF

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

      iElem = targetList( indElem )

      ! 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
        sourcePos = tLevelDesc%depFromCoarser( iElem )%elem%val( iSourceElem )

        ! load pdf
        f(1) = sState( ?IDX?( 1, sourcePos, nScalars, snSize))
        f(2) = sState( ?IDX?( 2, sourcePos, nScalars, snSize))
        f(3) = sState( ?IDX?( 3, sourcePos, nScalars, snSize))
        f(4) = sState( ?IDX?( 4, sourcePos, nScalars, snSize))
        f(5) = sState( ?IDX?( 5, sourcePos, nScalars, snSize))
        f(6) = sState( ?IDX?( 6, sourcePos, nScalars, snSize))
        f(7) = sState( ?IDX?( 7, sourcePos, nScalars, snSize))
        f(8) = sState( ?IDX?( 8, sourcePos, nScalars, snSize))
        f(9) = sState( ?IDX?( 9, sourcePos, nScalars, snSize))

        ! calculate rho and vel
        rho = sum( f )
        inv_rho = 1.0_rk/rho
        sourceP(iSourceElem) = rho - rho0

?? copy :: vel_d2q9_a( vel(1), vel(2), f, inv_rho )
        sourceV(:,iSourceElem) = vel(:)

        S = getStrain( f, rho, vel, omegaSource, layout%fStencil%cxDir )
        sourceS(1,iSourceElem) = S(1,1)
        sourceS(2,iSourceElem) = S(2,2)
        sourceS(3,iSourceElem) = S(1,2)

      end do ! iSourceElem

      ! ------------------ LINEAR INTERPOLATION -------------------
      ! JQ: taking pressure from parent gave better results than linear
      ! p = sourceP( tLevelDesc%depFromCoarser( iElem )%parentPosInSource )
      ! linearly interpolate Pressure, Vel, and fNeq
      p      = sum( weight(1:nSourceElems)*sourceP(  1:nSourceElems) )
      S(1,1) = sum( weight(1:nSourceElems)*sourceS(1,1:nSourceElems) )
      S(2,2) = sum( weight(1:nSourceElems)*sourceS(2,1:nSourceElems) )
      S(1,2) = sum( weight(1:nSourceElems)*sourceS(3,1:nSourceElems) )
      S(2,1) = S(1,2)

      ! interpolate vel by compact
      vel(:) = mus_interpolate_compact2d( sourceV(1:2,1:nSourceElems),  &
        &                                 sourceS(1:3,1:nSourceElems),  &
        &                  tLevelDesc%depFromCoarser( iElem )%coord(1), &
        &                  tLevelDesc%depFromCoarser( iElem )%coord(2)  )
      ! ------------------ LINEAR INTERPOLATION -------------------

?? IF (loglvl >= 10) THEN
write(dbgUnit(10),*) '   iElem: ', iElem, ', nSources: ', nSourceElems
write(dbgUnit(10),*) '  weight: ', weight(:)
write(dbgUnit(10),*) 'pressure: ', sourceP(:), p * pFac
write(dbgUnit(10),*) '   vel x: ', sourceV(1,:), vel(1) * vFac
write(dbgUnit(10),*) '   vel y: ', sourceV(2,:), vel(2) * vFac
do iDir = 1, QQ
  write(dbgUnit(10),*) 'fNeq', iDir, ': ', sourcefNeq(iDir,:), fNeq(iDir)*fNeqFac
end do
write(dbgUnit(10),*) ''
?? ENDIF

      ! scale macroscopic quantities (LB unit) from source level to target level
      rho = p   * pFac + rho0
      vel = vel * vFac
      S   = s   * sFac

      ! Calculate fEq
?? copy :: fEq_d2q9_a( fEq, vel(1), vel(2), rho, rho )
      fNeq = getNEq( S, omegaTarget, layout%fStencil%cxDirRK, layout%weight )

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

      ! Now write pdf (fEq+fNeq) to the target element
      tState(?IDX?( 1, targetPos, nScalars, tnSize))=fEq(1)+fNeq(1)
      tState(?IDX?( 2, targetPos, nScalars, tnSize))=fEq(2)+fNeq(2)
      tState(?IDX?( 3, targetPos, nScalars, tnSize))=fEq(3)+fNeq(3)
      tState(?IDX?( 4, targetPos, nScalars, tnSize))=fEq(4)+fNeq(4)
      tState(?IDX?( 5, targetPos, nScalars, tnSize))=fEq(5)+fNeq(5)
      tState(?IDX?( 6, targetPos, nScalars, tnSize))=fEq(6)+fNeq(6)
      tState(?IDX?( 7, targetPos, nScalars, tnSize))=fEq(7)+fNeq(7)
      tState(?IDX?( 8, targetPos, nScalars, tnSize))=fEq(8)+fNeq(8)
      tState(?IDX?( 9, targetPos, nScalars, tnSize))=fEq(9)+fNeq(9)

    end do ! indElem

write(dbgUnit(1),*) ''
?? IF (loglvl >= 1) THEN
write(dbgUnit(1),*) '----- Leave fillFiner_compact2D -----------'
write(dbgUnit(1),*) ''
?? ENDIF

  end subroutine fillFiner_compact2D
! ****************************************************************************** !

! ****************************************************************************** !
  pure function getStrain( f, rho, vel, omega, cxDir ) result( S )
    ! ---------------------------------------------------------------------------
    !> stencil layout
    integer, intent(in) :: cxDir(3,QQ)
    !> pdf array ( post-collision value )
    real(kind=rk), intent(in) :: f(QQ), rho, vel(2)
    !> relaxation parameter
    real(kind=rk), intent(in) :: omega
    !> output array: strain rate tensor
    real(kind=rk) :: S(2,2)
    ! ---------------------------------------------------------------------------
    real(kind=rk) :: fac

    ! Get the second velocity moments of the source element's pdf
    ! do jVal = 1, 2
    !   do iVal = 1, 2

    !     order = 0
    !     order( iVal ) = order( iVal ) + 1
    !     order( jVal ) = order( jVal ) + 1

    !     ! Calculate second velocity moment
    !     S(iVal,jVal) = - get_moment( QQ, layout%fStencil%cxDir, order, f) &
    !       &            + vel( iVal ) * vel( jVal ) + diagVal( iVal, jVal )

    !   end do
    ! end do

    S(1,1) = - get_moment( QQ, cxDir, [2,0,0], f) + vel(1) * vel(1) + cs2 * rho
    S(2,2) = - get_moment( QQ, cxDir, [0,2,0], f) + vel(2) * vel(2) + cs2 * rho

    S(2,1) = - get_moment( QQ, cxDir, [1,1,0], f) + vel(2) * vel(1)
    S(1,2) = S(2,1)

    ! Multiply factor and convert to pre-collision value
    fac = omega * cs2inv / ?post2pre?( omega ) * div1_2
    s = s * fac

  end function getStrain
! ****************************************************************************** !

! ****************************************************************************** !
  pure function getNEq( S, omega, cxDirRK, weight ) result( nEq )
    ! ---------------------------------------------------------------------------
    !> Strain rate tensor. It is a symmetric 2x2 matrix
    real(kind=rk), intent(in) :: S(2,2)
    real(kind=rk), intent(in) :: cxDirRK(3,QQ)
    real(kind=rk), intent(in) :: weight(QQ)
    real(kind=rk), intent(in) :: omega
    real(kind=rk) :: nEq(QQ)
    ! ---------------------------------------------------------------------------
    integer :: iDir
    real(kind=rk) :: fac
    ! ---------------------------------------------------------------------------

    do iDir = 1, QQ
      nEq(iDir) =   s(1,1) * ( cxDirRK(1,iDir) * cxDirRK(1,iDir) - cs2 ) &
        &         + s(2,2) * ( cxDirRK(2,iDir) * cxDirRK(2,iDir) - cs2 ) &
        &         + s(2,1) * cxDirRK(2,iDir) * cxDirRK(1,iDir) * 2.0_rk
    end do
    fac = cs2inv / omega * ?pre2post?( omega )
    nEq(:) = - weight(:) * nEq(:) * fac

  end function getNEq
! ****************************************************************************** !

end module mus_interpolate_d2q9_module
! ****************************************************************************** !
