! Copyright (c) 2015-2016 Jiaxing Qi <jiaxing.qi@uni-siegen.de>
! Copyright (c) 2016 Tobias Schneider <tobias1.schneider@student.uni-siegen.de>
! Copyright (c) 2018-2020 Kannan Masilamani <kannan.masilamani@uni-siegen.de>
! Copyright (c) 2018 Raphael Haupt <raphael.haupt@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_d3q19Macros.inc'
?? include 'header/lbm_d3q27Macros.inc'
?? include 'treelm/source/logMacros.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_d3q27_module
  use iso_c_binding, only: c_f_pointer

  ! include treelm modules
  use env_module,              only: rk
  use tem_param_module,        only: cs2inv, cs2, div1_2, &
    &                                div1_36, div1_8, div3_4h, div1_4, div2_3, &
    &                                rho0
  use tem_element_module,      only: eT_GhostFromCoarser, &
    &                                eT_ghostFromFiner
  use tem_construction_module, only: tem_levelDesc_type
  use tem_debug_module,        only: dbgUnit
  use tem_logging_module,      only: logUnit
  use tem_varSys_module,       only: tem_varSys_type
  use tem_time_module,         only: tem_time_type

  ! include musubi modules
  ! use mus_field_prop_module,         only: mus_field_prop_type
  use mus_scheme_layout_module,      only: mus_scheme_layout_type
  use mus_derVarPos_module,          only: mus_derVarPos_type
  use mus_physics_module,            only: mus_physics_type
  use mus_interpolate_header_module, only: mus_interpolation_method_type
  use mus_interpolate_tools_module,  only: get_fNeqFac_f2c, &
    &                                      get_fNeqFac_c2f
  use mus_derivedQuantities_module2, only: secondMom
  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 ::  get_polynomial_3D
  public :: eval_polynomial_3D
  public :: fillCoarser_compactD3Q27
  public :: fillFiner_compactD3Q27

  integer :: QQ = 27
?? copy :: q27

contains

! ****************************************************************************** !
  pure function eval_polynomial_3D( a, b, c, coord ) result( res )
    ! ---------------------------------------------------------------------------
    real(kind=rk), intent(in) :: a(10), b(10), c(10)
    ! coordinates(x,y,z) of target position
    real(kind=rk), intent(in) :: coord(3)
    real(kind=rk)             :: res(3)
    real(kind=rk)             :: xx, yy, zz, xy, yz, xz
    ! ---------------------------------------------------------------------------

    xx = coord(1)*coord(1)
    yy = coord(2)*coord(2)
    zz = coord(3)*coord(3)
    xy = coord(1)*coord(2)
    yz = coord(2)*coord(3)
    xz = coord(1)*coord(3)

    ! Evaluate polynomial
    res(1) = a(1) + a(2)*coord(1) + a(3)*coord(2) + a(4)*coord(3) &
      &           + a(5)*xx       + a(6)*yy       + a(7)*zz &
      &           + a(8)*xy       + a(9)*yz       + a(10)*xz
    res(2) = b(1) + b(2)*coord(1) + b(3)*coord(2) + b(4)*coord(3) &
      &           + b(5)*xx       + b(6)*yy       + b(7)*zz &
      &           + b(8)*xy       + b(9)*yz       + b(10)*xz
    res(3) = c(1) + c(2)*coord(1) + c(3)*coord(2) + c(4)*coord(3) &
      &           + c(5)*xx       + c(6)*yy       + c(7)*zz &
      &           + c(8)*xy       + c(9)*yz       + c(10)*xz

  end function eval_polynomial_3D
! ****************************************************************************** !

! ****************************************************************************** !
  !> Construct polynomial and interpolate at target position.
  !! To construct polynomial, velocity and strain rate of 4 source locations
  !! are needed:
  !! H(0,0,0), K(1,1,0), M(1,0,1), N(0,1,1)
  subroutine get_polynomial_3D( u, s, a, b, c )
    ! ---------------------------------------------------------------------------
    ! Velocity (Ux, Uy, Uz)
    real(kind=rk), intent(in) :: u(3,4)
    ! Strain rate (Sxx,Syy,Szz,Sxy,Syz,Sxz)
    real(kind=rk), intent(in) :: s(6,4)
    real(kind=rk), intent(out) :: a(10), b(10), c(10)
    ! ---------------------------------------------------------------------------
    real(kind=rk) :: a0,ax,ay,az,axx,ayy,azz,axy,ayz,axz
    real(kind=rk) :: b0,bx,by,bz,bxx,byy,bzz,bxy,byz,bxz
    real(kind=rk) :: c0,cx,cy,cz,cxx,cyy,czz,cxy,cyz,cxz
    real(kind=rk) :: P, Q, R
    real(kind=rk) :: Pa, Qa, Ra
    real(kind=rk) :: Pb, Qb, Rb
    real(kind=rk) :: Pc, Qc, Rc
    ! ---------------------------------------------------------------------------

    a0 = u(1,1)
    b0 = u(2,1)
    c0 = u(3,1)

    axx = ( -s(1,1) + s(1,2) + s(1,3) - s(1,4) ) * div1_4
    axy = ( -s(1,1) + s(1,2) - s(1,3) + s(1,4) ) * div1_2
    axz = ( -s(1,1) - s(1,2) + s(1,3) + s(1,4) ) * div1_2

    bxy = ( -s(2,1) + s(2,2) + s(2,3) - s(2,4) ) * div1_2
    byy = ( -s(2,1) + s(2,2) - s(2,3) + s(2,4) ) * div1_4
    byz = ( -s(2,1) - s(2,2) + s(2,3) + s(2,4) ) * div1_2

    cxz = ( -s(3,1) + s(3,2) + s(3,3) - s(3,4) ) * div1_2
    cyz = ( -s(3,1) + s(3,2) - s(3,3) + s(3,4) ) * div1_2
    czz = ( -s(3,1) - s(3,2) + s(3,3) + s(3,4) ) * div1_4

    !  axy + 2bxx = ( -s(4,1) + s(4,2) + s(4,3) - s(4,4) )
              bxx = ( -s(4,1) + s(4,2) + s(4,3) - s(4,4) - axy ) * div1_2
    ! 2ayy +  bxy = ( -s(4,1) + s(4,2) - s(4,3) + s(4,4) )
       ayy        = ( -s(4,1) + s(4,2) - s(4,3) + s(4,4) - bxy ) * div1_2
                P =   -s(4,1) - s(4,2) + s(4,3) + s(4,4)

    !  byz + 2cyy = ( -s(5,1) + s(5,2) - s(5,3) + s(5,4) )
              cyy = ( -s(5,1) + s(5,2) - s(5,3) + s(5,4) - byz ) * div1_2
    ! 2bzz +  cyz = ( -s(5,1) - s(5,2) + s(5,3) + s(5,4) )
       bzz        = ( -s(5,1) - s(5,2) + s(5,3) + s(5,4) - cyz ) * div1_2
                Q =   -s(5,1) + s(5,2) + s(5,3) - s(5,4)

    !  axz + 2cxx = ( -s(6,1) + s(6,2) + s(6,3) - s(6,4) )
              cxx = ( -s(6,1) + s(6,2) + s(6,3) - s(6,4) - axz ) * div1_2
    ! 2azz +  cxz = ( -s(6,1) - s(6,2) + s(6,3) + s(6,4) )
       azz        = ( -s(6,1) - s(6,2) + s(6,3) + s(6,4) - cxz ) * div1_2
                R =   -s(6,1) + s(6,2) - s(6,3) + s(6,4)

    ! P == ayz + bzx
    ! Q == bzx + cxy
    ! R == ayz + cxy
    ayz = (   P - Q + R ) * div1_2
    bxz = (   P + Q - R ) * div1_2
    cxy = ( - P + Q + R ) * div1_2

    Pa = u(1,2) - a0 - axx - ayy - axy
    Qa = u(1,3) - a0 - axx - azz - axz
    Ra = u(1,4) - a0 - ayy - azz - ayz
    ax = (  Pa + Qa - Ra ) * div1_2
    ay = (  Pa - Qa + Ra ) * div1_2
    az = ( -Pa + Qa + Ra ) * div1_2

    Pb = u(2,2) - b0 - bxx - byy - bxy
    Qb = u(2,3) - b0 - bxx - bzz - bxz
    Rb = u(2,4) - b0 - byy - bzz - byz
    bx = (  Pb + Qb - Rb ) * div1_2
    by = (  Pb - Qb + Rb ) * div1_2
    bz = ( -Pb + Qb + Rb ) * div1_2

    Pc = u(3,2) - c0 - cxx - cyy - cxy
    Qc = u(3,3) - c0 - cxx - czz - cxz
    Rc = u(3,4) - c0 - cyy - czz - cyz
    cx = (  Pc + Qc - Rc ) * div1_2
    cy = (  Pc - Qc + Rc ) * div1_2
    cz = ( -Pc + Qc + Rc ) * div1_2

    a = [a0,ax,ay,az,axx,ayy,azz,axy,ayz,axz]
    b = [b0,bx,by,bz,bxx,byy,bzz,bxy,byz,bxz]
    c = [c0,cx,cy,cz,cxx,cyy,czz,cxy,cyz,cxz]

  end subroutine get_polynomial_3D
! ****************************************************************************** !

! ****************************************************************************** !
  !> 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_compactd3q27( 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, indElem, iDir, iSourceElem, ii
    ! macroscopic velocities for all source elements
    real(kind=rk) :: sourceV(3,4)
    real(kind=rk) :: sourceS(6,4)
    ! macroscopic velocities for all source elements
    real(kind=rk) :: f(QQ), fEq(QQ), fNeq(QQ), tempNeq(QQ)
    real(kind=rk) :: rho, vel(3), p, inv_rho, XYZ(-1:1,1:3)
    real(kind=rk) :: omegaSource, omegaTarget
    real(kind=rk) :: pFac, vFac, sFac, fNeqFac
    real(kind=rk) :: sNeqPost2Pre, tNeqPre2Post
    real(kind=rk) :: a(10), b(10), c(10)
    !> 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
    type(mus_fluid_type), pointer :: fluid
    integer :: nScalars
    ! ---------------------------------------------------------------------------
    fluid => fieldProp(1)%fluid
    nScalars = varSys%nScalars

    sourceLevel = level + 1
    targetLevel = level
    ! KM: temp fix: assuming omega is constant
    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 )
    fNeqFac = get_fNeqFac_f2c( omegaSource, omegaTarget, sFac )
    sNeqPost2Pre = 1.0d0 / ?post2pre?( omegaSource )
    tNeqPre2Post = ?pre2post?( omegaTarget )

?? IF (loglvl >= 10) THEN
write(dbgUnit(10), *) ''
write(dbgUnit(10), *) '---- Enter fillMyGhostsFromFiner_Compact2Dincomp --------'
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 )

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

        posInTotal = tLevelDesc%depFromFiner( iElem )%elem%val( iSourceElem )
        do iDir = 1, QQ
          f(iDir) = sState( ?IDX?( iDir, posInTotal, nScalars, snSize ))
        end do

        ! calculate rho and take average
        rho = sum( f )
        inv_rho = 1.0_rk / rho
        p = p + ( rho - rho0 )

        ! calculate vel
        vel(1) = sum( layout%fStencil%cxDirRK(1,:) * f(:) ) * inv_rho
        vel(2) = sum( layout%fStencil%cxDirRK(2,:) * f(:) ) * inv_rho
        vel(3) = sum( layout%fStencil%cxDirRK(3,:) * f(:) ) * inv_rho
        sourceV(1:3,iSourceElem) = vel(1:3)

        ! calculate feq and fNeq
        do ii = 1,3
          XYZ( 0,ii) = -div2_3 + vel(ii) * vel(ii)
          XYZ( 1,ii) = - ( XYZ(0,ii) + 1.0d0 + vel(ii) ) * 0.5d0
          XYZ(-1,ii) = XYZ(1,ii) + vel(ii)
        end do
        do iDir = 1, QQ
          fEq(iDir) = - rho * XYZ(layout%fStencil%cxDir(1,iDir),1) &
            &               * XYZ(layout%fStencil%cxDir(2,iDir),2) &
            &               * XYZ(layout%fStencil%cxDir(3,iDir),3)
        enddo
        tempNeq = (f - fEq) * sNeqPost2Pre
        fNeq = fNeq + tempNeq

        ! calculate S
        sourceS(1:6,iSourceElem) = -1.5_rk * omegaSource * inv_rho &
          &                        * secondMom( layout%fStencil%cxcx, tempNeq, QQ )

      end do ! iSourceElem
      ! ------------------------------------------------------------------------

      ! ------------------ Compact INTERPOLATION -------------------------------
      ! pressure and fNeq by average
      ! velocity: quadratic compact interpolate
      p = p * div1_4
      call get_polynomial_3D( sourceV, sourceS, a, b, c)
      vel(:) = eval_polynomial_3D( a, b, c, [0.5_rk, 0.5_rk, 0.5_rk] )
      fNeq = fNeq * div1_4
      ! ------------------ Compact INTERPOLATION -------------------------------

! DEBUG --------------------------------------------------------------
?? IF (loglvl >= 10) THEN
write(dbgUnit(10),*) 'iElem: ', iElem !, ', nSources: ', nSourceElems
write(dbgUnit(10),"(A,3F12.8)") '1: ', sourceV(:,1)
write(dbgUnit(10),"(A,3F12.8)") '2: ', sourceV(:,2)
write(dbgUnit(10),"(A,3F12.8)") '3: ', sourceV(:,3)
write(dbgUnit(10),"(A,3F12.8)") '4: ', sourceV(:,4)
write(dbgUnit(10),"(A,3F12.8)") '   ', vel*vFac
write(dbgUnit(10),"(A)") ''
?? ENDIF
! DEBUG --------------------------------------------------------------

      ! Scale quantities on source level to target level
      rho  = p    * pFac + rho0
      vel  = vel  * vFac
      fNeq = fNeq * fNeqFac * tNeqPre2Post

      ! Calculate fEq
      do ii = 1,3
        XYZ( 0,ii) = -div2_3 + vel(ii) * vel(ii)
        XYZ( 1,ii) = - ( XYZ(0,ii) + 1.0d0 + vel(ii) ) * 0.5d0
        XYZ(-1,ii) = XYZ(1,ii) + vel(ii)
      end do

      ! Ghost pos in total
      posInTotal = iElem + tLevelDesc%offset( 1, eT_ghostFromFiner)
      do iDir = 1, QQ
        fEq(iDir) = - rho * XYZ(layout%fStencil%cxDir(1,iDir),1) &
          &               * XYZ(layout%fStencil%cxDir(2,iDir),2) &
          &               * XYZ(layout%fStencil%cxDir(3,iDir),3)
        tState(?IDX?( iDir, posInTotal, nScalars, tnSize ))=fEq(iDir)+fNeq(iDir)
      end do

    end do ! indElem

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

  end subroutine fillCoarser_compactd3q27
! ****************************************************************************** !

! ****************************************************************************** !
  !> 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_compactD3Q27( 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, iSource, ii, iDir
    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 :: sourceS(:,:)
    real(kind=rk), allocatable :: weight(:)
    ! temporary array
    real(kind=rk) :: f(QQ), fEq(QQ), fNeq(QQ)
    real(kind=rk) :: rho, vel(3), p, inv_rho, XYZ(-1:1,1:3)
    real(kind=rk) :: omegaSource, omegaTarget
    real(kind=rk) :: pFac, vFac, sFac, fNeqFac, sNeqPost2Pre, tNeqPre2Post
    real(kind=rk) :: a(10), b(10), c(10)
    integer :: nMaxSources
    ! real(kind=rk) :: sum1, sum2, sum3, sum4, sum5, sum6, sum7, sum8, sum9, sum0
    !> debug variables
?? IF (loglvl >= 1) THEN
    integer(kind=long_k) :: treeID
?? ENDIF
    type(mus_fluid_type), pointer :: fluid
    integer :: nScalars
    ! --------------------------------------------------------------------------
    fluid => fieldProp(1)%fluid
    nScalars = varSys%nScalars
    nMaxSources = method%nMaxSources

    sourceLevel = level
    targetLevel = level + 1
    ! KM: temp fix: assuming omega is constant
    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 fillFiner_compactD3Q27   -----'
write(dbgUnit(10), *) ''
write(dbgUnit(10), *) ' nMaxSources: ', 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(nMaxSources)      )
    allocate( sourceV(3,nMaxSources)    )
    allocate( sourcefNeq(QQ,nMaxSources))
    allocate( weight(nMaxSources)       )
    allocate( sourceS(6,nMaxSources)    )

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

      iElem = targetList( indElem )

      ! Find out how many fine source elements we have for interpolation.
      targetPos = iElem + tLevelDesc%offset( 1, eT_ghostFromCoarser)
      nSourceElems = tLevelDesc%depFromCoarser( iElem )%elem%nVals

      weight(1:nSourceElems) = &
        &    tLevelDesc%depFromCoarser( iElem )%weight(1:nSourceElems)

      do iSource = 1, nSourceElems

        sourcePos = tLevelDesc%depFromCoarser( iElem )%elem%val( iSource )
        do iDir = 1, QQ
          f(iDir) = sState( ?IDX?( iDir, sourcePos, nScalars, snSize ))
        end do

        ! calculate rho and vel
        rho = sum( f )
        inv_rho = 1.0_rk / rho
        vel(1) = sum( layout%fStencil%cxDirRK(1,:) * f(:) ) * inv_rho
        vel(2) = sum( layout%fStencil%cxDirRK(2,:) * f(:) ) * inv_rho
        vel(3) = sum( layout%fStencil%cxDirRK(3,:) * f(:) ) * inv_rho

        ! calculate feq and fNeq
        ! XYZ( -1:1, 1:3 )
        do ii = 1,3
          XYZ( 0,ii) = -div2_3 + vel(ii) * vel(ii)
          XYZ( 1,ii) = - ( XYZ(0,ii) + 1.0d0 + vel(ii) ) * 0.5d0
          XYZ(-1,ii) = XYZ(1,ii) + vel(ii)
        end do
        do iDir = 1, QQ
          fEq(iDir) = - rho * XYZ(layout%fStencil%cxDir(1,iDir),1) &
            &               * XYZ(layout%fStencil%cxDir(2,iDir),2) &
            &               * XYZ(layout%fStencil%cxDir(3,iDir),3)
          fNeq(iDir) = ( f(iDir) - fEq(iDir) ) * sNeqPost2Pre
        enddo

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

        ! calculate S
sourceS(1,iSource) = -1.5_rk*omegaSource*inv_rho*sum(layout%fStencil%cxcx(1,:) * fNeq(:))
sourceS(2,iSource) = -1.5_rk*omegaSource*inv_rho*sum(layout%fStencil%cxcx(2,:) * fNeq(:))
sourceS(3,iSource) = -1.5_rk*omegaSource*inv_rho*sum(layout%fStencil%cxcx(3,:) * fNeq(:))
sourceS(4,iSource) = -1.5_rk*omegaSource*inv_rho*sum(layout%fStencil%cxcx(4,:) * fNeq(:))
sourceS(5,iSource) = -1.5_rk*omegaSource*inv_rho*sum(layout%fStencil%cxcx(5,:) * fNeq(:))
sourceS(6,iSource) = -1.5_rk*omegaSource*inv_rho*sum(layout%fStencil%cxcx(6,:) * fNeq(:))
      end do ! iSource

      ! ------------------ LINEAR INTERPOLATION -------------------
      ! linearly interpolate Pressure and fNeq
      p = sum( weight(1:nSourceElems)*sourceP(1:nSourceElems) )
      do iDir = 1, QQ
        fNeq(iDir) = sum(   weight(         1:nSourceElems) &
          &               * sourcefNeq(iDir,1:nSourceElems) )
      end do

      ! interpolate vel by compact
      call get_polynomial_3D( sourceV, sourceS, a, b, c)
      vel(:) = eval_polynomial_3D( a, b, c, &
        &            tLevelDesc%depFromCoarser( iElem )%coord(1:3)  )
      ! ------------------ LINEAR INTERPOLATION -------------------

?? IF (loglvl >= 10) THEN
write(dbgUnit(10),"(2(A,I0))")' iElem: ', iElem, ', nSources: ', nSourceElems
write(dbgUnit(10),"(A,4F10.7)") 'weight: ', weight(:)
write(dbgUnit(10),"(A,5F10.7)") ' p: ', sourceP(:), p * pFac
write(dbgUnit(10),"(A,5F10.7)") 'Ux: ', sourceV(1,:), vel(1) * vFac
write(dbgUnit(10),"(A,5F10.7)") 'Uy: ', sourceV(2,:), vel(2) * vFac
write(dbgUnit(10),"(A,5F10.7)") 'Uz: ', sourceV(3,:), vel(3) * vFac
do iDir = 1, QQ
  write(dbgUnit(10),"(A,I2,A,5F10.7)") '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
      do ii = 1,3
        XYZ( 0,ii) = -div2_3 + vel(ii) * vel(ii)
        XYZ( 1,ii) = - ( XYZ(0,ii) + 1.0d0 + vel(ii) ) * 0.5d0
        XYZ(-1,ii) = XYZ(1,ii) + vel(ii)
      end do
      do iDir = 1, QQ
        fEq(iDir) = - rho * XYZ(layout%fStencil%cxDir(1,iDir),1) &
          &               * XYZ(layout%fStencil%cxDir(2,iDir),2) &
          &               * XYZ(layout%fStencil%cxDir(3,iDir),3)
        tState(?IDX?( iDir, targetPos, nScalars, tnSize ))=fEq(iDir)+fNeq(iDir)
      end do

    end do ! indElem

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

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

  end subroutine fillFiner_compactD3Q27
! ****************************************************************************** !

end module mus_interpolate_d3q27_module
