! Copyright (c) 2012-2013 Manuel Hasert <m.hasert@grs-sim.de>
! Copyright (c) 2013-2014 Simon Zimny <s.zimny@grs-sim.de>
! Copyright (c) 2013-2015, 2018-2020 Kannan Masilamani <kannan.masilamani@uni-siegen.de>
! Copyright (c) 2014-2016 Jiaxing Qi <jiaxing.qi@uni-siegen.de>
! Copyright (c) 2016 Tobias Schneider <tobias1.schneider@student.uni-siegen.de>
! Copyright (c) 2020 Peter Vitt <peter.vitt2@uni-siegen.de>
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions are met:
!
! 1. Redistributions of source code must retain the above copyright notice,
! this list of conditions and the following disclaimer.
!
! 2. Redistributions in binary form must reproduce the above copyright notice,
! this list of conditions and the following disclaimer in the documentation
! and/or other materials provided with the distribution.
!
! THIS SOFTWARE IS PROVIDED BY THE UNIVERSITY OF SIEGEN “AS IS” AND ANY EXPRESS
! OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
! OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
! IN NO EVENT SHALL UNIVERSITY OF SIEGEN OR CONTRIBUTORS BE LIABLE FOR ANY
! DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
! (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
! ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
! ****************************************************************************** !
!> author: Manuel Hasert
!! 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
!! for sourceLevel = targetLevel+1
!! sourceLevel = targetLevel-1
!!
!! For an overview over implemented interpolation methods, see
!! [Interpolation methods](../page/features/intp_methods.html)
!!
?? include 'header/lbm_macros.inc'
?? include 'treelm/source/logMacros.inc'
module mus_interpolate_debug_module
  use iso_c_binding, only: c_f_pointer

  ! include treelm modules
  use env_module,            only: rk, long_k
  use tem_param_module,      only: cs2inv, PI
  use tem_topology_module,   only: tem_coordOfID
  use tem_element_module,    only: eT_GhostFromCoarser,eT_ghostFromFiner
  use tem_logging_module,    only: logUnit
  use tem_debug_module,      only: dbgUnit
  use tem_construction_module, only: tem_levelDesc_type
  use tem_varSys_module,       only: tem_varSys_type
  use tem_time_module,       only: tem_time_type
  use tem_stencil_module,    only: tem_stencilHeader_type
  use tem_construction_module, only: tem_levelDesc_type

  ! include musubi modules
  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_derivedQuantities_module2, only: set_pdfDiffusive
  use mus_interpolate_header_module, only: mus_interpolation_method_type
  use mus_field_prop_module,         only: mus_field_prop_type
  use mus_fluid_module,              only: mus_fluid_type

  implicit none

  private

  public :: do_nothing
  public :: do_nothing_arbi

?? if ( .not. SOA) then
  ! public :: fillMyGhostsFromFiner_average
  ! public :: fillFinerGhostsFromMe_average

  ! public :: fillMyGhostsFromFiner_debug
  ! public :: fillFinerGhostsFromMe_debug

  ! public :: fillMyGhostsFromFiner_copyFirst
  ! public :: fillFinerGhostsFromMe_copyFirst

  public :: fillMyGhostsFromFiner_TGV2D
  public :: fillFinerGhostsFromMe_TGV2D

  ! public :: fillMyGhostsFromFiner_TGV3D
  ! public :: fillFinerGhostsFromMe_TGV3D
?? endif

  public :: TGV_2D

  real(kind=rk), parameter :: debugValue = 1.0d0

 contains

! ****************************************************************************** !
  !> Fill GhostFromFiner elements on my level with debug value
  !!
  !! 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 do_nothing( 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(:)
    ! -------------------------------------------------------------------- !

  end subroutine do_nothing
! ****************************************************************************** !

! ****************************************************************************** !
  !> Fill GhostFromFiner elements on my level with debug value
  !!
  !! This subroutine's interface must match the abstract interface definition
  !! [[intpRoutine_arbitraryVal]] in intp/[[mus_interpolate_header_module]].f90
  !! in order to be callable via
  !! [[mus_interpolation_method_type:do_intpArbiVal]] function pointer.
  subroutine do_nothing_arbi( method, tLevelDesc, level, stencil, sVal,   &
    &                         snSize, tVal, tnSize, nTargets, targetList, &
    &                         nScalars                                    )
    ! ------------------------------------------------------------------ !
    class(mus_interpolation_method_type), intent(in) :: method

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

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

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

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

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

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

    !> number of scalars to interpolate
    integer, intent(in) :: nScalars
    ! ------------------------------------------------------------------ !

  end subroutine do_nothing_arbi
! ****************************************************************************** !


?? if( .not. SOA) then
! ****************************************************************************** !
  !> Interpolation routine that is based on a simple weighted average of source
  !! nodes
  !!
  !! This is the interpolation fine-> coarse. It does not need weights, as the
  !! distances source <-> target are equal for all source nodes. Thus, we can
  !! add all source pdf's and divide by the number of source nodes.
  !!
  !! This subroutine's interface must match the abstract interface definition
  !! [[intpRoutine]] in intp/[[mus_interpolate_header_module]].f90 in order to
  !! be callable via [[mus_interpolation_method_type:do_intp]] function pointer.
  subroutine fillMyGhostsFromFiner_Average( method, fieldProp, tLevelDesc, &
    &                                       level, sState, snSize, tState, &
    &                                       tnSize, tAuxField, layout,     &
    &                                       nTargets, targetList, physics, &
    &                                       time, varSys, derVarPos        )
    ! -------------------------------------------------------------------- !
    class(mus_interpolation_method_type), intent(in) :: method

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

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

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

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

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

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

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

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

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

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

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

    !> position of all derive variable in varSys for all fields
    type(mus_derVarPos_type), intent(in) :: derVarPos(:)
    ! -------------------------------------------------------------------- !
    integer :: sourceLevel    ! level of source elements
    integer :: sourceElem     ! treeId of current source element
    integer :: targetLevel    ! level of target elements
    integer :: targetElem     ! treeId of current source element
    integer :: iDir           ! current direction (discrete velocity) for loop
    integer :: iElem          ! current target element (for outer loop)
    integer :: indElem        ! element counter for indirection list
    integer :: iSourceElem    ! current source element (for inner loop)
    integer :: nSourceElems   ! number of source elements for the current target
    real(kind=rk) :: tVal     ! temporary value: accumulator for pdf density
    integer :: QQ
    ! --------------------------------------------------------------------------
    QQ = layout%fStencil%QQ

    targetLevel = level
    sourceLevel = level+1

    ! Treat every coarse target element:
    do indElem = 1, nTargets

      iElem = targetList( indElem )

      targetElem = iElem + tLevelDesc%offset( 1, eT_ghostFromFiner)

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

      ! Now that we know the target element, loop over all nVars and
      ! discrete velocity directions for each field variable
      do iDir = 1, QQ
        tVal = 0._rk
        !... and all the fine source elements...
        do iSourceElem = 1, nSourceElems
          ! Find out the source element's treeId
          sourceElem = tLevelDesc%depFromFiner( iElem )%elem%val( iSourceElem )
          ! Accumulate all source pdf's in tVal. We are just adding up values
          ! here...
          tVal = tVal+sstate(?IDX?(iDir, sourceElem, QQ, snSize))
        enddo ! iSourceElem

        ! Now we have all pdf's for this discrete velocity (iDir) added up.
        ! Divide this by the number of sources and write the result to the
        ! target element.
        tstate(?IDX?(iDir, targetElem, QQ, tnSize )) = tVal/dble( nSourceElems )
      enddo ! iDir
    enddo ! indElem
  end subroutine fillMyGhostsFromFiner_Average
! ****************************************************************************** !


! ****************************************************************************** !
  !> Interpolation routine that is based on a simple weighted average of
  !! source nodes
  !!
  !! This is the interpolation coarse-> fine. Weights are needed here, as
  !! the distances source <-> target are different for the source nodes.
  !!
  !! 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_Average( method, fieldProp, tLevelDesc, &
    &                                       level, sState, snSize, tState, &
    &                                       tnSize, tAuxField, layout,     &
    &                                       nTargets, targetList, physics, &
    &                                       time, varSys, derVarPos        )
    ! -------------------------------------------------------------------- !
    class(mus_interpolation_method_type), intent(in) :: method

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

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

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

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

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

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

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

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

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

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

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

    !> position of all derive variable in varSys for all fields
    type(mus_derVarPos_type), intent(in) :: derVarPos(:)
    ! -------------------------------------------------------------------- !
    integer :: sourceLevel    ! level of source elements
    integer :: sourceElem     ! treeId of current source element
    integer :: targetLevel    ! level of target elements
    integer :: targetElem     ! treeId of current source element
    integer :: iDir           ! current direction (discrete velocity) for loop
    integer :: iElem          ! current target element (for outer loop)
    integer :: indElem        ! element counter for indirection list
    integer :: iSourceElem    ! current source element (for inner loop)
    integer :: nSourceElems   ! number of source elements for the current target
    real(kind=rk) :: tVal     ! temporary value: accumulator for pdf density
    real(kind=rk) :: fTmp     ! temporary pdf value
    ! weight for source elements depending on distance to target element
    real(kind=rk) :: weight(method%nMaxSources)
    integer :: QQ
    logical :: error
    ! --------------------------------------------------------------------------
    QQ = layout%fStencil%QQ

    sourceLevel = level
    targetLevel = level + 1

    ! Treat every coarse target element:
    do indElem = 1, nTargets

      iElem = targetList( indElem )

      targetElem = iElem + tLevelDesc%offset( 1, eT_ghostFromCoarser)

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

      weight( 1:nSourceElems ) = 1._rk / dble(nSourceElems)

      error = .false.
      do iDir = 1, QQ

        tVal = 0._rk
        do iSourceElem = 1, nSourceElems
          sourceElem = tLevelDesc%depFromCoarser( iElem )%elem%val( iSourceElem)
          fTmp = sstate(?IDX?( iDir, sourceElem, QQ, snSize))
          if( fTmp < 0._rk ) then
            error = .true.
          end if

          ! Accumulate all source pdf's with weighting in tVal
          tVal = tVal + weight(iSourceElem) * fTmp
        end do ! iSourceElem

        ! Now we have all pdf's for this discrete velocity (iDir) weighted
        ! and added up. The weights add up to 1, so we don't have to
        ! divide here by anything.
        tstate(?IDX?( iDir, targetElem, QQ, tnSize)) = tVal
      end do ! iDir

      if( error ) then
        write(dbgUnit(1),*) ' error for tID ', tLevelDesc%total( targetElem )
      end if

    enddo

  end subroutine fillFinerGhostsFromMe_Average
! ****************************************************************************** !


! ****************************************************************************** !
  !> Fill GhostFromFiner elements on my level with debug value
  !!
  !! This subroutine's interface must match the abstract interface definition
  !! [[intpRoutine]] in intp/[[mus_interpolate_header_module]].f90 in order to
  !! be callable via [[mus_interpolation_method_type:do_intp function pointer.
  subroutine fillMyGhostsFromFiner_Debug( 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 :: targetLevel      ! level of target elements
    integer :: targetElem       ! treeId of current source element
    integer :: iDir             ! current direction (discrete velocity) for loop
    integer :: iElem            ! current target element (for outer loop)
    integer :: indElem          ! element counter for indirection list
    integer :: QQ
    ! --------------------------------------------------------------------------
    QQ = layout%fStencil%QQ

    targetLevel = level

    ! Treat every coarse target element:
    do indElem = 1, nTargets

      iElem = targetList( indElem )

      targetElem = iElem + tLevelDesc%offset( 1, eT_ghostFromFiner)

      ! Now that we know the target element, loop over all nVars and
      ! discrete velocity directions for each field variable
      do iDir = 1, QQ
        ! ... and simply copy the pdf from source to target.
        tstate(?IDX?(iDir, targetElem, QQ, tnSize))=debugValue
      enddo
    enddo

  end subroutine fillMyGhostsFromFiner_Debug
! ****************************************************************************** !


! ****************************************************************************** !
  !> Fill the ghost elements on the finer level with debug value
  !!
  !! 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_Debug( 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 :: targetLevel      ! level of target elements
    integer :: targetElem       ! treeId of current source element
    integer :: iDir             ! current direction (discrete velocity) for loop
    integer :: iElem            ! current target element (for outer loop)
    integer :: indElem          ! element counter for indirection list
    integer :: QQ
    ! --------------------------------------------------------------------------
    QQ = layout%fStencil%QQ
    targetLevel = level + 1

    ! Treat every fine target element:
    do indElem = 1, nTargets

      iElem = targetList( indElem )

      targetElem = iElem + tLevelDesc%offset( 1, eT_ghostFromCoarser)

      ! Now that we know the target element, loop over all nVars and
      ! discrete velocity directions for each field variable
      do iDir = 1, QQ
          ! ... and simply copy the pdf from source to target.
        tstate(?IDX?(iDir, targetElem, QQ, tnSize))=debugValue
      enddo
    enddo

  end subroutine fillFinerGhostsFromMe_Debug
! ****************************************************************************** !


! ****************************************************************************** !
  !> This simple method does not perform actual interpolation, but simply copies
  !! one of the source elements pdf values to the target element.
  !! It is for DEBUG purpose, the quality is not good enough for production use!
  !!
  !! This subroutine's interface must match the abstract interface definition
  !! [[intpRoutine]] in intp/[[mus_interpolate_header_module]].f90 in order to
  !! be callable via [[mus_interpolation_method_type:do_intp]] function pointer.
  subroutine fillMyGhostsFromFiner_CopyFirst( method, fieldProp, tLevelDesc, &
    &                                         level, sState, snSize, tState, &
    &                                         tnSize, tAuxField, layout,     &
    &                                         nTargets, targetList, physics, &
    &                                         time, varSys, derVarPos        )
    ! -------------------------------------------------------------------- !
    class(mus_interpolation_method_type), intent(in) :: method

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

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

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

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

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

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

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

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

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

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

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

    !> position of all derive variable in varSys for all fields
    type(mus_derVarPos_type), intent(in) :: derVarPos(:)
    ! -------------------------------------------------------------------- !
    integer :: sourceLevel      ! level of source elements
    integer :: sourceElem       ! treeId of current source element
    integer :: targetLevel      ! level of target elements
    integer :: targetElem       ! treeId of current source element
    integer :: iDir             ! current direction (discrete velocity) for loop
    integer :: iElem            ! current target element (for outer loop)
    integer :: indElem          ! element counter for indirection list
    integer :: QQ
    ! --------------------------------------------------------------------------
    QQ = layout%fStencil%QQ

    targetLevel = level
    sourceLevel = level+1

    ! Treat every coarse target element:
    do indElem = 1, nTargets

      iElem = targetList( indElem )

      targetElem = iElem + tLevelDesc%offset( 1, eT_ghostFromFiner)

      sourceElem = tLevelDesc%depFromFiner( iElem )%elem%val(1)

      ! Now that we know the target element, loop over all nVars and
      ! discrete velocity directions for each field variable
      do iDir = 1, QQ

        ! ... and simply copy the pdf from source to target.
        tstate( ?IDX?( iDir, targetElem, QQ, tnSize)) = &
          & sstate(?IDX?( iDir, sourceElem, QQ, snSize))

      enddo
    enddo
  end subroutine fillMyGhostsFromFiner_CopyFirst
! ****************************************************************************** !


! ****************************************************************************** !
  !> This simple method does not perform actual interpolation, but simply copies
  !! one of the source element's pdf values to the target element. It is mostly
  !! there for debugging purposes, the quality is not good enough for
  !! production use!
  !!
  !! 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_CopyFirst( method, fieldProp, tLevelDesc, &
    &                                         level, sState, snSize, tState, &
    &                                         tnSize, tAuxField, layout,     &
    &                                         nTargets, targetList, physics, &
    &                                         time, varSys, derVarPos        )
    ! -------------------------------------------------------------------- !
    class(mus_interpolation_method_type), intent(in) :: method

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

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

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

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

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

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

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

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

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

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

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

    !> position of all derive variable in varSys for all fields
    type(mus_derVarPos_type), intent(in) :: derVarPos(:)
    ! -------------------------------------------------------------------- !
    integer :: sourceLevel      ! level of source elements
    integer :: sourceElem       ! treeId of current source element
    integer :: targetLevel      ! level of target elements
    integer :: targetElem       ! treeId of current source element
    integer :: iDir             ! current direction (discrete velocity) for loop
    integer :: iElem            ! current target element (for outer loop)
    integer :: indElem          ! element counter for indirection list
    integer :: QQ
    ! --------------------------------------------------------------------------
    QQ = layout%fStencil%QQ
    sourceLevel = level
    targetLevel = level + 1

    ! Treat every fine target element:
    do indElem = 1, nTargets

      iElem = targetList( indElem )

      sourceElem = tLevelDesc%depFromCoarser( iElem )%elem%val( 1)

      targetElem = iElem + tLevelDesc%offset( 1, eT_ghostFromCoarser)

      ! Now that we know the target element, loop over all nVars and
      ! discrete velocity directions for each field variable
      do iDir = 1, QQ
          ! ... and simply copy the pdf from source to target.
        tstate(?IDX?( iDir, targetElem, QQ, tnSize)) = &
          &  sstate(?IDX?( iDir, sourceElem, QQ, snSize))
      enddo ! iDir
    enddo

  end subroutine fillFinerGhostsFromMe_CopyFirst
! ****************************************************************************** !


! ****************************************************************************** !
  !> This simple method does not perform actual interpolation,
  !! but get analytical solution and set pdf to target element.
  !! It is for debugging purposes, not for production use!
  !!
  !! 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_TGV2D( 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 :: targetLevel ! level of target elements
    integer :: targetElem  ! treeId of current source element
    integer :: iDir, iElem, indElem
    real(kind=rk) :: targetCoord(3), omega, dx
    real(kind=rk) :: mom(layout%fStencil%QQ)
    real(kind=rk) :: tPDF(layout%fStencil%QQ)
    integer :: coord(4)
    integer(kind=long_k) :: targetID
    integer :: QQ
    type(mus_fluid_type), pointer :: fluid
    ! --------------------------------------------------------------------------
    fluid => fieldProp(1)%fluid
    QQ = layout%fStencil%QQ

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

    ! Treat every fine target element:
    do indElem = 1, nTargets

      iElem = targetList( indElem )

      ! target position in state array
      targetElem = iElem + tLevelDesc%offset( 1, eT_ghostFromCoarser)
      ! target treeID
      targetID   = tLevelDesc%total( targetElem )

      ! get bary
      coord       = tem_CoordOfId(targetID)
      dx          = 2.0_rk * PI / dble( 2**coord(4) )
      targetCoord = (dble(coord(:3)) + 0.5_rk) * dx

      ! get pressure, velocity, strainRate
      mom(1:6) = TGV_2D( targetCoord, time%sim )

! -----------------------------------------------------------------------------
! Debug output
?? IF (loglvl >= 1) THEN
write(dbgUnit(1),"(A, 3F12.8)") 'target pos: ', targetCoord(1:3)
write(dbgUnit(1),"(A, F12.8, A, 2F12.8, A, 3F12.8)") &
  &            'P: ', mom(1), ', V: ', mom(2:3), ', S: ', mom(4:6)
write(dbgUnit(1),"(A)") ' '
?? ENDIF
! -----------------------------------------------------------------------------

      ! convert PHY to LB
      mom(1)   = mom(1) * cs2inv / physics%fac( targetLevel )%press + 1.0_rk
      mom(2:3) = mom(2:3) / physics%fac( targetLevel )%vel
      mom(4:6) = mom(4:6) / physics%fac( targetLevel )%strainRate

      tPDF = set_pdfDiffusive( layout, omega, 1.0_rk, mom )

      ! fill state vector
      do iDir = 1, QQ
        tstate(?IDX?(iDir, targetElem, QQ, tnSize)) = tPDF(iDir)
      enddo

    enddo

  end subroutine fillFinerGhostsFromMe_TGV2D
! ****************************************************************************** !

! ****************************************************************************** !
  !> This simple method does not perform actual interpolation,
  !! but get analytical solution and set pdf to target element
  !! It is mostly
  !! there for debugging purposes, the quality is not good enough for
  !! production use!
  !!
  !! This subroutine's interface must match the abstract interface definition
  !! [[intpRoutine]] in intp/[[mus_interpolate_header_module]].f90 in order to
  !! be callable via [[mus_interpolation_method_type:do_intp]] function pointer.
  subroutine fillMyGhostsFromFiner_TGV2D( 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 :: targetLevel ! level of target elements
    integer :: targetElem  ! treeId of current source element
    integer :: iDir, iElem, indElem
    real(kind=rk) :: targetCoord(3), omega, dx
    real(kind=rk) :: mom(layout%fStencil%QQ)
    real(kind=rk) :: tPDF(layout%fStencil%QQ)
    integer :: coord(4)
    integer(kind=long_k) :: targetID
    integer :: QQ
    type(mus_fluid_type), pointer :: fluid
    ! --------------------------------------------------------------------------
    fluid => fieldProp(1)%fluid
    QQ = layout%fStencil%QQ

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

    ! Treat every fine target element:
    do indElem = 1, nTargets

      iElem = targetList( indElem )

      ! target position in state array
      targetElem = iElem + tLevelDesc%offset( 1, eT_ghostFromFiner)
      targetID   = tLevelDesc%total( targetElem )

      ! get bary
      coord       = tem_CoordOfId(targetID)
      dx          = 2.0_rk * PI / real(2**coord(4), rk)
      targetCoord = (real(coord(:3), rk) + 0.5_rk) * dx

      ! get pressure, velocity, strainRate
      mom(1:6) = TGV_2D( targetCoord, time%sim )

! -----------------------------------------------------------------------------
! Debug output
?? IF (loglvl >= 1) THEN
write(dbgUnit(1),"(A, 3F12.8)") 'target pos: ', targetCoord(1:3)
write(dbgUnit(1),"(A, F12.8, A, 2F12.8, A, 3F12.8)") &
  &            'P: ', mom(1), ', V: ', mom(2:3), ', S: ', mom(4:6)
write(dbgUnit(1),"(A)") ' '
?? ENDIF
! -----------------------------------------------------------------------------

      ! convert PHY to LB
      mom(1)   = mom(1) * cs2inv / physics%fac( targetLevel )%press + 1.0_rk
      mom(2:3) = mom(2:3) / physics%fac( targetLevel )%vel
      mom(4:6) = mom(4:6) / physics%fac( targetLevel )%strainRate

      tPDF = set_pdfDiffusive( layout, omega, 1.0_rk, mom )

      ! fill state vector
      do iDir = 1, QQ
        tstate(?IDX?(iDir, targetElem, QQ, tnSize)) = tPDF(iDir)
      enddo

    enddo

  end subroutine fillMyGhostsFromFiner_TGV2D
! ****************************************************************************** !


! ****************************************************************************** !
  !> No comment yet!
  !!
  !! @TODO add comment
  !!
  !! This subroutine's interface must match the abstract interface definition
  !! [[intpRoutine]] in intp/[[mus_interpolate_header_module]].f90 in order to
  !! be callable via [[mus_interpolation_method_type:do_intp]] function pointer.
  subroutine fillFinerGhostsFromMe_TGV3D( 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 :: targetLevel ! level of target elements
    integer :: targetElem  ! treeId of current source element
    integer :: iDir, iElem, indElem
    real(kind=rk) :: targetCoord(3), omega, dx
    real(kind=rk) :: mom(layout%fStencil%QQ)
    real(kind=rk) :: tPDF(layout%fStencil%QQ)
    real(kind=rk) :: m(6)
    integer :: coord(4)
    integer(kind=long_k) :: targetID
    integer :: QQ
    type(mus_fluid_type), pointer :: fluid
    ! --------------------------------------------------------------------------
    fluid => fieldProp(1)%fluid
    QQ = layout%fStencil%QQ

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

    ! Treat every fine target element:
    do indElem = 1, nTargets

      iElem = targetList( indElem )

      ! target position in state array
      targetElem = iElem + tLevelDesc%offset( 1, eT_ghostFromCoarser)
      ! target treeID
      targetID   = tLevelDesc%total( targetElem )

      ! get bary
      coord       = tem_CoordOfId(targetID)
      dx          = 2.0_rk * PI / dble( 2**coord(4) )
      targetCoord = (dble(coord(:3)) + 0.5_rk) * dx

      ! get pressure, velocity, strainRate
      m = TGV_2D( targetCoord, time%sim )

      ! convert PHY to LB
      mom(1) = m(1) * cs2inv / physics%fac( targetLevel )%press + 1.0_rk
      mom(2) = m(2) / physics%fac( targetLevel )%vel
      mom(3) = m(3) / physics%fac( targetLevel )%vel
      mom(4) = 0._rk
      mom(5) = m(4) / physics%fac( targetLevel )%strainRate
      mom(6) = m(5) / physics%fac( targetLevel )%strainRate
      mom(7) = 0._rk
      mom(8) = m(6) / physics%fac( targetLevel )%strainRate
      mom(9:) = 0._rk

      tPDF = set_pdfDiffusive( layout, omega, 1.0_rk, mom )

      ! fill state vector
      do iDir = 1, QQ
        tstate(?IDX?(iDir, targetElem, QQ, tnSize)) = tPDF(iDir)
      enddo

    enddo

  end subroutine fillFinerGhostsFromMe_TGV3D
! ****************************************************************************** !

! ****************************************************************************** !
  !> No comment yet!
  !!
  !! @TODO add comment
  !!
  !! This subroutine's interface must match the abstract interface definition
  !! [[intpRoutine]] in intp/[[mus_interpolate_header_module]].f90 in order to
  !! be callable via [[mus_interpolation_method_type:do_intp]] function pointer.
  subroutine fillMyGhostsFromFiner_TGV3D( 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 :: targetLevel ! level of target elements
    integer :: targetElem  ! treeId of current source element
    integer :: iDir, iElem, indElem
    real(kind=rk) :: targetCoord(3), omega, dx
    real(kind=rk) :: mom(layout%fStencil%QQ)
    real(kind=rk) :: tPDF(layout%fStencil%QQ), m(6)
    integer :: coord(4)
    integer(kind=long_k) :: targetID
    integer :: QQ
    type(mus_fluid_type), pointer :: fluid
    ! --------------------------------------------------------------------------
    fluid => fieldProp(1)%fluid
    QQ = layout%fStencil%QQ

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

    ! Treat every fine target element:
    do indElem = 1, nTargets

      iElem = targetList( indElem )

      ! target position in state array
      targetElem = iElem + tLevelDesc%offset( 1, eT_ghostFromFiner)
      targetID   = tLevelDesc%total( targetElem )

      ! get bary
      coord       = tem_CoordOfId(targetID)
      dx          = 2.0_rk * PI / real(2**coord(4), rk)
      targetCoord = (real(coord(:3), rk) + 0.5_rk) * dx

      ! get pressure, velocity, strainRate
      m = TGV_2D( targetCoord, time%sim )

      ! convert PHY to LB
      mom(1)  = m(1) * cs2inv / physics%fac( targetLevel )%press + 1.0_rk
      mom(2)  = m(2) / physics%fac( targetLevel )%vel
      mom(3)  = m(3) / physics%fac( targetLevel )%vel
      mom(4)  = 0._rk
      mom(5)  = m(4) / physics%fac( targetLevel )%strainRate
      mom(6)  = m(5) / physics%fac( targetLevel )%strainRate
      mom(7)  = 0._rk
      mom(8)  = m(6) / physics%fac( targetLevel )%strainRate
      mom(9:) = 0._rk

      tPDF = set_pdfDiffusive( layout, omega, 1.0_rk, mom )

      ! fill state vector
      do iDir = 1, QQ
        tstate(?IDX?(iDir, targetElem, QQ, tnSize)) = tPDF(iDir)
      enddo

    enddo

  end subroutine fillMyGhostsFromFiner_TGV3D
! ****************************************************************************** !
?? ENDIF

! ****************************************************************************** !
  !> This routine returns the analytical solution of TGV 2D testcase for a given
  !! position and time (coord, t)
  !!
  pure function TGV_2D( coord, t ) result ( res )
    ! ---------------------------------------------------------------------------
    !> position and time
    real(kind=rk), intent(in) :: coord(3), t
    !> pressure, velX, velY, Sxx, Syy, Sxy
    real(kind=rk)             :: res(6)
    ! ---------------------------------------------------------------------------
    real(kind=rk) :: nuPhy, tD, x, y
    ! ---------------------------------------------------------------------------

    x = coord(1)
    y = coord(2)
    ! viscosity = 1 / Re
    nuPhy = 1.0 / 25.0_rk
    tD = 1.0 / nuPhy / 2

    ! pressure
    res(1) = -(cos(2.0*x)+cos(2.0*y)) * exp(-2.0*t/tD) / 4.0
    ! velocity
    res(2) = -cos(x) * sin(y) * exp(-t/tD)
    res(3) =  sin(x) * cos(y) * exp(-t/tD)
    ! strain rate S
    res(4) =  sin(x) * sin(y) * exp(-t/tD)
    res(5) = -sin(x) * sin(y) * exp(-t/tD)
    res(6) = 0.0_rk

  end function TGV_2D
! ****************************************************************************** !
end module mus_interpolate_debug_module
! ****************************************************************************** !
