mus_interpolate_quadratic_module.f90 Source File


This file depends on

sourcefile~~mus_interpolate_quadratic_module.f90~~EfferentGraph sourcefile~mus_interpolate_quadratic_module.f90 mus_interpolate_quadratic_module.f90 sourcefile~mus_field_prop_module.f90 mus_field_prop_module.f90 sourcefile~mus_interpolate_quadratic_module.f90->sourcefile~mus_field_prop_module.f90 sourcefile~mus_scheme_layout_module.f90 mus_scheme_layout_module.f90 sourcefile~mus_interpolate_quadratic_module.f90->sourcefile~mus_scheme_layout_module.f90 sourcefile~mus_relaxationparam_module.f90 mus_relaxationParam_module.f90 sourcefile~mus_interpolate_quadratic_module.f90->sourcefile~mus_relaxationparam_module.f90 sourcefile~mus_fluid_module.f90 mus_fluid_module.f90 sourcefile~mus_interpolate_quadratic_module.f90->sourcefile~mus_fluid_module.f90 sourcefile~mus_interpolate_header_module.f90 mus_interpolate_header_module.f90 sourcefile~mus_interpolate_quadratic_module.f90->sourcefile~mus_interpolate_header_module.f90 sourcefile~mus_interpolate_tools_module.f90 mus_interpolate_tools_module.f90 sourcefile~mus_interpolate_quadratic_module.f90->sourcefile~mus_interpolate_tools_module.f90 sourcefile~mus_dervarpos_module.f90 mus_derVarPos_module.f90 sourcefile~mus_interpolate_quadratic_module.f90->sourcefile~mus_dervarpos_module.f90 sourcefile~mus_physics_module.f90 mus_physics_module.f90 sourcefile~mus_interpolate_quadratic_module.f90->sourcefile~mus_physics_module.f90 sourcefile~mus_pdf_module.f90 mus_pdf_module.f90 sourcefile~mus_interpolate_quadratic_module.f90->sourcefile~mus_pdf_module.f90 sourcefile~mus_field_prop_module.f90->sourcefile~mus_fluid_module.f90 sourcefile~mus_field_prop_module.f90->sourcefile~mus_physics_module.f90 sourcefile~mus_species_module.f90 mus_species_module.f90 sourcefile~mus_field_prop_module.f90->sourcefile~mus_species_module.f90 sourcefile~mus_scheme_header_module.f90 mus_scheme_header_module.f90 sourcefile~mus_field_prop_module.f90->sourcefile~mus_scheme_header_module.f90 sourcefile~mus_poisson_module.f90 mus_poisson_module.f90 sourcefile~mus_field_prop_module.f90->sourcefile~mus_poisson_module.f90 sourcefile~mus_moments_type_module.f90 mus_moments_type_module.f90 sourcefile~mus_scheme_layout_module.f90->sourcefile~mus_moments_type_module.f90 sourcefile~mus_relaxationparam_module.f90->sourcefile~mus_scheme_layout_module.f90 sourcefile~mus_relaxationparam_module.f90->sourcefile~mus_dervarpos_module.f90 sourcefile~mus_relaxationparam_module.f90->sourcefile~mus_physics_module.f90 sourcefile~mus_nonnewtonian_module.f90 mus_nonNewtonian_module.f90 sourcefile~mus_relaxationparam_module.f90->sourcefile~mus_nonnewtonian_module.f90 sourcefile~mus_turbulence_module.f90 mus_turbulence_module.f90 sourcefile~mus_relaxationparam_module.f90->sourcefile~mus_turbulence_module.f90 sourcefile~mus_relaxationparam_module.f90->sourcefile~mus_moments_type_module.f90 sourcefile~mus_relaxationparam_module.f90->sourcefile~mus_scheme_header_module.f90 sourcefile~mus_graddata_module.f90 mus_gradData_module.f90 sourcefile~mus_relaxationparam_module.f90->sourcefile~mus_graddata_module.f90 sourcefile~mus_fluid_module.f90->sourcefile~mus_relaxationparam_module.f90 sourcefile~mus_fluid_module.f90->sourcefile~mus_physics_module.f90 sourcefile~mus_fluid_module.f90->sourcefile~mus_pdf_module.f90 sourcefile~mus_fluid_module.f90->sourcefile~mus_nonnewtonian_module.f90 sourcefile~mus_fluid_module.f90->sourcefile~mus_turbulence_module.f90 sourcefile~mus_turb_viscosity_module.f90 mus_turb_viscosity_module.f90 sourcefile~mus_fluid_module.f90->sourcefile~mus_turb_viscosity_module.f90 sourcefile~mus_fluid_module.f90->sourcefile~mus_scheme_header_module.f90 sourcefile~mus_mrtrelaxation_module.f90 mus_mrtRelaxation_module.f90 sourcefile~mus_fluid_module.f90->sourcefile~mus_mrtrelaxation_module.f90 sourcefile~mus_interpolate_header_module.f90->sourcefile~mus_field_prop_module.f90 sourcefile~mus_interpolate_header_module.f90->sourcefile~mus_scheme_layout_module.f90 sourcefile~mus_interpolate_header_module.f90->sourcefile~mus_dervarpos_module.f90 sourcefile~mus_interpolate_header_module.f90->sourcefile~mus_physics_module.f90 sourcefile~mus_directions_module.f90 mus_directions_module.f90 sourcefile~mus_interpolate_header_module.f90->sourcefile~mus_directions_module.f90 sourcefile~mus_interpolate_tools_module.f90->sourcefile~mus_scheme_layout_module.f90 sourcefile~mus_interpolate_tools_module.f90->sourcefile~mus_interpolate_header_module.f90 sourcefile~mus_interpolate_tools_module.f90->sourcefile~mus_pdf_module.f90 sourcefile~mus_interpolate_tools_module.f90->sourcefile~mus_scheme_header_module.f90 sourcefile~mus_dervarpos_module.f90->sourcefile~mus_scheme_layout_module.f90 sourcefile~mus_nonnewtonian_module.f90->sourcefile~mus_scheme_layout_module.f90 sourcefile~mus_nonnewtonian_module.f90->sourcefile~mus_physics_module.f90 sourcefile~mus_nonnewtonian_module.f90->sourcefile~mus_scheme_header_module.f90 sourcefile~mus_turbulence_module.f90->sourcefile~mus_scheme_layout_module.f90 sourcefile~mus_turbulence_module.f90->sourcefile~mus_graddata_module.f90 sourcefile~mus_species_module.f90->sourcefile~mus_physics_module.f90 sourcefile~mus_turb_viscosity_module.f90->sourcefile~mus_scheme_layout_module.f90 sourcefile~mus_turb_viscosity_module.f90->sourcefile~mus_turbulence_module.f90 sourcefile~mus_turb_viscosity_module.f90->sourcefile~mus_scheme_header_module.f90 sourcefile~mus_turb_viscosity_module.f90->sourcefile~mus_graddata_module.f90 sourcefile~mus_vreman_module.f90 mus_Vreman_module.f90 sourcefile~mus_turb_viscosity_module.f90->sourcefile~mus_vreman_module.f90 sourcefile~mus_wale_module.f90 mus_WALE_module.f90 sourcefile~mus_turb_viscosity_module.f90->sourcefile~mus_wale_module.f90 sourcefile~mus_smagorinsky_module.f90 mus_Smagorinsky_module.f90 sourcefile~mus_turb_viscosity_module.f90->sourcefile~mus_smagorinsky_module.f90 sourcefile~mus_poisson_module.f90->sourcefile~mus_physics_module.f90 sourcefile~mus_mrtrelaxation_module.f90->sourcefile~mus_moments_type_module.f90 sourcefile~mus_mrtrelaxation_module.f90->sourcefile~mus_scheme_header_module.f90 sourcefile~mus_vreman_module.f90->sourcefile~mus_turbulence_module.f90 sourcefile~mus_vreman_module.f90->sourcefile~mus_graddata_module.f90 sourcefile~mus_derivedquantities_module.f90 mus_derivedQuantities_module.f90 sourcefile~mus_vreman_module.f90->sourcefile~mus_derivedquantities_module.f90 sourcefile~mus_wale_module.f90->sourcefile~mus_turbulence_module.f90 sourcefile~mus_wale_module.f90->sourcefile~mus_graddata_module.f90 sourcefile~mus_wale_module.f90->sourcefile~mus_derivedquantities_module.f90 sourcefile~mus_smagorinsky_module.f90->sourcefile~mus_scheme_layout_module.f90 sourcefile~mus_smagorinsky_module.f90->sourcefile~mus_turbulence_module.f90 sourcefile~mus_smagorinsky_module.f90->sourcefile~mus_graddata_module.f90 sourcefile~mus_smagorinsky_module.f90->sourcefile~mus_derivedquantities_module.f90

Files dependent on this one

sourcefile~~mus_interpolate_quadratic_module.f90~~AfferentGraph sourcefile~mus_interpolate_quadratic_module.f90 mus_interpolate_quadratic_module.f90 sourcefile~mus_interpolate_d2q9_module.f90 mus_interpolate_d2q9_module.f90 sourcefile~mus_interpolate_d2q9_module.f90->sourcefile~mus_interpolate_quadratic_module.f90 sourcefile~mus_interpolate_module.f90 mus_interpolate_module.f90 sourcefile~mus_interpolate_module.f90->sourcefile~mus_interpolate_quadratic_module.f90 sourcefile~mus_dynloadbal_module.f90 mus_dynLoadBal_module.f90 sourcefile~mus_dynloadbal_module.f90->sourcefile~mus_interpolate_module.f90 sourcefile~mus_construction_module.f90 mus_construction_module.f90 sourcefile~mus_dynloadbal_module.f90->sourcefile~mus_construction_module.f90 sourcefile~mus_hvs_aux_module.f90 mus_hvs_aux_module.f90 sourcefile~mus_hvs_aux_module.f90->sourcefile~mus_interpolate_module.f90 sourcefile~mus_aux_module.f90 mus_aux_module.f90 sourcefile~mus_aux_module.f90->sourcefile~mus_interpolate_module.f90 sourcefile~mus_construction_module.f90->sourcefile~mus_interpolate_module.f90 sourcefile~mus_harvesting.f90 mus_harvesting.f90 sourcefile~mus_harvesting.f90->sourcefile~mus_hvs_aux_module.f90 sourcefile~mus_harvesting.f90->sourcefile~mus_construction_module.f90 sourcefile~mus_hvs_construction_module.f90 mus_hvs_construction_module.f90 sourcefile~mus_harvesting.f90->sourcefile~mus_hvs_construction_module.f90 sourcefile~mus_control_module.f90 mus_control_module.f90 sourcefile~mus_control_module.f90->sourcefile~mus_aux_module.f90 sourcefile~musubi.f90 musubi.f90 sourcefile~musubi.f90->sourcefile~mus_aux_module.f90 sourcefile~musubi.f90->sourcefile~mus_control_module.f90 sourcefile~mus_program_module.f90 mus_program_module.f90 sourcefile~musubi.f90->sourcefile~mus_program_module.f90 sourcefile~mus_program_module.f90->sourcefile~mus_dynloadbal_module.f90 sourcefile~mus_program_module.f90->sourcefile~mus_aux_module.f90 sourcefile~mus_program_module.f90->sourcefile~mus_construction_module.f90 sourcefile~mus_program_module.f90->sourcefile~mus_control_module.f90 sourcefile~mus_hvs_construction_module.f90->sourcefile~mus_construction_module.f90

Contents


Source Code

! Copyright (c) 2011-2013 Manuel Hasert <m.hasert@grs-sim.de>
! Copyright (c) 2011 Jan Hueckelheim <j.hueckelheim@grs-sim.de>
! Copyright (c) 2011-2013 Harald Klimach <harald.klimach@uni-siegen.de>
! Copyright (c) 2011-2016, 2018-2020 Kannan Masilamani <kannan.masilamani@uni-siegen.de>
! Copyright (c) 2011-2014 Simon Zimny <s.zimny@grs-sim.de>
! Copyright (c) 2012-2016 Jiaxing Qi <jiaxing.qi@uni-siegen.de>
! Copyright (c) 2013-2014 Kartik Jain <kartik.jain@uni-siegen.de>
! Copyright (c) 2016 Tobias Schneider <tobias1.schneider@student.uni-siegen.de>
! Copyright (c) 2016, 2019-2020 Peter Vitt <peter.vitt2@uni-siegen.de>
! Copyright (c) 2016, 2018 Raphael Haupt <raphael.haupt@uni-siegen.de>
! Copyright (c) 2019 Jana Gericke <jana.gericke@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.
! ****************************************************************************** !
! Copyright (c) 2011-2013 Manuel Hasert <m.hasert@grs-sim.de>
! Copyright (c) 2011 Harald Klimach <harald.klimach@uni-siegen.de>
! Copyright (c) 2011 Konstantin Kleinheinz <k.kleinheinz@grs-sim.de>
! Copyright (c) 2011-2012 Simon Zimny <s.zimny@grs-sim.de>
! Copyright (c) 2012, 2014-2016 Jiaxing Qi <jiaxing.qi@uni-siegen.de>
! Copyright (c) 2012 Kartik Jain <kartik.jain@uni-siegen.de>
! Copyright (c) 2013-2015, 2019 Kannan Masilamani <kannan.masilamani@uni-siegen.de>
! Copyright (c) 2016 Tobias Schneider <tobias1.schneider@student.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.
! Copyright (c) 2014-2015, 2019-2020 Kannan Masilamani <kannan.masilamani@uni-siegen.de>
! Copyright (c) 2015-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.







! Copyright (c) 2015-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.



!> 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_quadratic_module
  use iso_c_binding, only: c_loc, c_ptr, c_f_pointer

  ! include treelm modules
  use env_module,              only: rk
  use tem_aux_module,          only: tem_abort
  use tem_param_module,        only: cs2inv, cs2, div1_3, div1_6, div1_9, &
    &                                div2_9, div5_9, div2_3, div1_4, &
    &                                div1_21, div4_21, div5_21, div1_7, &
    &                                div3_7, div1_42, div5_42, div1_2, &
    &                                div1_36, div3_4h, c_x, c_y, c_z, &
    &                                rho0, rho0Inv
  use tem_element_module,      only: eT_GhostFromCoarser
  use tem_construction_module, only: depSource_type, tem_levelDesc_type
  use tem_debug_module,        only: dbgUnit
  use tem_logging_module,      only: logUnit
  use tem_stencil_module,      only: tem_stencilHeader_type
  use tem_matrix_module,       only: tem_matrix_type, tem_matrix_dump
  use tem_varSys_module,       only: tem_varSys_type
  use tem_time_module,         only: tem_time_type

  ! include musubi modules
  use mus_pdf_module,                only: pdf_data_type
  use mus_scheme_layout_module,      only: mus_scheme_layout_type
  use mus_interpolate_header_module, only: mus_interpolation_method_type
  use mus_interpolate_tools_module,  only: mus_intp_convertMomToPDF3D,        &
    &                                      mus_intp_convertMomToPDF3D_incomp, &
    &                                      mus_intp_convertMomToPDF2D,        &
    &                                      mus_intp_convertMomToPDF2D_incomp, &
    &                                      mus_intp_getMoments
  use mus_physics_module,            only: mus_physics_type
  use mus_field_prop_module,         only: mus_field_prop_type
  use mus_fluid_module,              only: mus_fluid_type
  use mus_relaxationParam_module,    only: mus_calcOmegaFromVisc
  use mus_derVarPos_module,          only: mus_derVarPos_type

  implicit none

  private

  public :: fillArbiFinerGhostsFromMe_quad
  public :: fillArbiFinerGhostsFromMe_quad2D
  public :: fillFinerGhostsFromMe_quad
  public :: fillFinerGhostsFromMe_quadLES
  public :: fillFinerGhostsFromMe_quadIncomp
  public :: fillFinerGhostsFromMe_quadIncompLES
  public :: fillFinerGhostsFromMe_quad2D
  public :: fillFinerGhostsFromMe_quad2DIncomp
  public :: mus_interpolate_quad2D_leastSq
  public :: biquadratic_interpolate

contains

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

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

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

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

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

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

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

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

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

      iElem = targetList( indElem )

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

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

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

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

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

      end do  ! iSourceElem

      ! interpolate all auxiliary variables by quadratic interpolation
      tArbi(1:nScalars) = mus_interpolate_quad3D_leastSq(         &
        &   srcMom      = sArbi(1:nScalars, 1:nSourceElems),      &
        &   targetCoord = tLevelDesc%depFromCoarser( iElem )%coord, &
        &   LSFmat      = method%intpMat_forLSF%matArray            &
        &                       %val(posInIntpMatLSF),              &
        &   nSources    = nSourceElems,                             &
        &   nVals       = nScalars                               )

      ! write interpolated value
      tVal((targetElem-1)*nScalars+1 : targetElem*nScalars) &
          & = tArbi
    enddo

  end subroutine fillArbiFinerGhostsFromMe_quad
! **************************************************************************** !

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

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

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

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

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

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

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

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

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

      iElem = targetList( indElem )

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

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

      posInIntpMatLSF = tLevelDesc%depFromCoarser( iElem )%posInIntpMatLSF
      ! Now loop over all fine source elements for this target:
      do iSourceElem = 1, nSourceElems

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

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

      end do  ! iSourceElem

      ! interpolate all auxiliary variables by quadratic interpolation
      tArbi(1:nScalars) = mus_interpolate_quad2D_leastSq(         &
        &   srcMom      = sArbi(1:nScalars, 1:nSourceElems),      &
        &   targetCoord = tLevelDesc%depFromCoarser( iElem )%coord, &
        &   LSFmat      = method%intpMat_forLSF%matArray            &
        &                       %val(posInIntpMatLSF),              &
        &   nSources    = nSourceElems,                             &
        &   nVals       = nScalars                               )

      ! write interpolated value
      tVal((targetElem-1)*nScalars+1 : targetElem*nScalars) &
          & = tArbi
    enddo

  end subroutine fillArbiFinerGhostsFromMe_quad2D
! **************************************************************************** !

  ! ************************************************************************** !
  !> Fill fine ghost from coarse fluid by quadratic interpolation.
  !! 1. Compute moments for all source elements, save in momBuf
  !! 2. For each target, interpolate all moments using quadratic least square fit
  !! 3. Store target auxilary field
  !! 4. Compute viscosity on target element and compute source and target omega
  !! 5. Get nonEq scaling factor depeding on scheme layout and relaxation
  !! 6. Calculate Equilibrium and nonEquilibrium
  !! 7. calculate target: Eq + Scale * nonEquilibrium
  !! This routine is used by 3D 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_quad( 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) :: tPDF( layout%fStencil%QQ ) ! pdf to reconstruct from
    ! moments of the source elements' pdf
    real(kind=rk) :: srcMom( layout%fStencil%QQ, layout%fStencil%QQ )
    real(kind=rk) :: tMom(layout%fStencil%QQ)  ! target moment calculation
    integer :: QQ
    integer :: posInIntpMatLSF
    type(mus_fluid_type), pointer :: fluid
    real(kind=rk) :: nonEqScalingFacs(layout%fStencil%QQ)
    real(kind=rk) :: sOmegaKine, tOmegaKine, tVisc, invRho
    integer :: nScalars, tOffset
    integer :: dens_pos, vel_pos(3)
    ! ---------------------------------------------------------------------------
    fluid => fieldProp(1)%fluid
    nScalars = varSys%nScalars
    QQ = layout%fStencil%QQ
    dens_pos = varSys%method%val(derVarPos(1)%density)%auxField_varPos(1)
    vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:3)

    sourceLevel = level
    targetLevel = level + 1

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

      iElem = targetList( indElem )
      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 position
        sourceElem = tLevelDesc%depFromCoarser( iElem ) &
          &                    %elem%val( iSourceElem )

        ! Get macroscopic moments
        ! KM: \todo momBUF: precompute moments for all source elements before target
        ! loop because a source coarse element can be a source to multiple
        ! target fine elements, Use momBuf in mus_pdf_type in store the moments
        srcMom(:, iSourceElem) = mus_intp_getMoments(                         &
          &                               state     = sstate,                 &
          &                               elem      = sourceElem,             &
          &                               QQ        = QQ,                     &
          &                               nScalars  = nScalars,               &
          &                               nSize     = snSize,                 &
          &                               toMoments = layout%moment%toMoments )
      end do

      !interpolate moments to target element
      tMom(1:QQ) = mus_interpolate_quad3D_leastSq(                  &
        &   srcMom      = srcMom(1:QQ, 1:nSourceElems),             &
        &   targetCoord = tLevelDesc%depFromCoarser( iElem )%coord, &
        &   LSFmat      = method%intpMat_forLSF%matArray            &
        &                       %val(posInIntpMatLSF),              &
        &   nSources    = nSourceElems,                             &
        &   nVals       = QQ )

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

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

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

!      write(dbgUnit(1),*) 'sOmegaKine ', sOmegaKine
!      write(dbgUnit(1),*) 'tOmegaKine ', tOmegaKine

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

!      write(dbgUnit(1),*) 'Quad nonEqFacs', nonEqScalingFacs
      ! Convert moment to PDF
      tPDF = mus_intp_convertMomToPDF3D( moments          = tMom,             &
        &                                nonEqScalingFacs = nonEqScalingFacs, &
        &                                layout           = layout            )

      ! Now write the resulting pdf in the current direction to the target
      ! Element position
      do iVal = 1, QQ
        tState( ( targetelem-1)* nscalars+ ival) = tPDF(iVal)
      enddo

    end do ! indElem

  end subroutine fillFinerGhostsFromMe_quad
! ****************************************************************************** !

  ! ************************************************************************** !
  !> Fill fine ghost from coarse fluid by quadratic interpolation.
  !! 1. Compute moments for all source elements, save in momBuf
  !! 2. For each target, interpolate all moments using quadratic least square fit
  !! 3. Store target auxilary field
  !! 4. Compute viscosity on target element and compute source and target omega
  !! 5. Get nonEq scaling factor depeding on scheme layout and relaxation
  !! 6. Calculate Equilibrium and nonEquilibrium
  !! 7. calculate target: Eq + Scale * nonEquilibrium
  !! This routine is used by 3D 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_quadLES( 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) :: tPDF( layout%fStencil%QQ ) ! pdf to reconstruct from
    ! moments of the source elements' pdf
    real(kind=rk) :: srcMom( layout%fStencil%QQ, layout%fStencil%QQ )
    real(kind=rk) :: tMom(layout%fStencil%QQ)  ! target moment calculation
    integer :: QQ
    integer :: posInIntpMatLSF
    type(mus_fluid_type), pointer :: fluid
    real(kind=rk) :: nonEqScalingFacs(layout%fStencil%QQ)
    real(kind=rk) :: sOmegaKine, tOmegaKine, sVisc, tVisc, invRho
    real(kind=rk) :: sTurbVisc(layout%fStencil%QQ), tTurbVisc(1)
    integer :: nScalars, tOffset
    integer :: dens_pos, vel_pos(3)
    ! ---------------------------------------------------------------------------
    fluid => fieldProp(1)%fluid
    nScalars = varSys%nScalars
    QQ = layout%fStencil%QQ
    dens_pos = varSys%method%val(derVarPos(1)%density)%auxField_varPos(1)
    vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(:)

    sourceLevel = level
    targetLevel = level + 1

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

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

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

      !interpolate moments to target element
      tMom(1:QQ) = mus_interpolate_quad3D_leastSq(                  &
        &   srcMom      = srcMom(1:QQ, 1:nSourceElems),             &
        &   targetCoord = tLevelDesc%depFromCoarser( iElem )%coord, &
        &   LSFmat      = method%intpMat_forLSF%matArray            &
        &                       %val(posInIntpMatLSF),              &
        &   nSources    = nSourceElems,                             &
        &   nVals       = QQ )

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

      !interpolate turbulent viscosity to target element
      tTurbVisc = mus_interpolate_quad3D_leastSq(                   &
        &   srcMom      = sTurbVisc(1:nSourceElems),                &
        &   targetCoord = tLevelDesc%depFromCoarser( iElem )%coord, &
        &   LSFmat      = method%intpMat_forLSF%matArray            &
        &                       %val(posInIntpMatLSF),              &
        &   nSources    = nSourceElems,                             &
        &   nVals       = 1                                         )

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

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

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

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

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

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

      ! Now write the resulting pdf in the current direction to the target
      ! Element position
      do iVal = 1, QQ
        tState( ( targetelem-1)* nscalars+ ival) = tPDF(iVal)
      enddo

    end do ! indElem

  end subroutine fillFinerGhostsFromMe_quadLES
! ****************************************************************************** !

  ! ************************************************************************** !
  !> Fill fine ghost from coarse fluid by quadratic interpolation.
  !! 1. Compute moments for all source elements, save in momBuf
  !! 2. For each target, interpolate all moments using quadratic least square fit
  !! 3. Store target auxilary field
  !! 4. Get viscosity on target element and compute source and target omega
  !! 5. Get nonEq scaling factor depeding on scheme layout and relaxation
  !! 6. Convert moments to PDF
  !! 7. Calculate target: Eq + Scale * nonEquilibrium
  !! This routine is used by 3D acoustic quadratic interpolation for
  !! incompressible model.
  !!
  !! This subroutine's interface must match the abstract interface definition
  !! [[intpRoutine]] in intp/[[mus_interpolate_header_module]].f90 in order to
  !! be callable via [[mus_interpolation_method_type:do_intp]] function pointer.
  subroutine fillFinerGhostsFromMe_quadIncomp( 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) :: tPDF( layout%fStencil%QQ ) ! pdf to reconstruct from
    ! moments of the source elements' pdf
    real(kind=rk) :: srcMom( layout%fStencil%QQ, layout%fStencil%QQ )
    real(kind=rk) :: tMom(layout%fStencil%QQ)  ! target moment calculation
    integer :: QQ
    integer :: posInIntpMatLSF
    type(mus_fluid_type), pointer :: fluid
    real(kind=rk) :: nonEqScalingFacs(layout%fStencil%QQ)
    real(kind=rk) :: sOmegaKine, tOmegaKine, tVisc
    integer :: nScalars, tOffset
    integer :: dens_pos, vel_pos(3)
    ! ---------------------------------------------------------------------------
    fluid => fieldProp(1)%fluid
    QQ = layout%fStencil%QQ
    nScalars = varSys%nScalars
    dens_pos = varSys%method%val(derVarPos(1)%density)%auxField_varPos(1)
    vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:3)

    sourceLevel = level
    targetLevel = level + 1

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

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

      end do

      !interpolate moments to target element
      tMom(1:QQ) = mus_interpolate_quad3D_leastSq(                  &
        &   srcMom      = srcMom(1:QQ, 1:nSourceElems),             &
        &   targetCoord = tLevelDesc%depFromCoarser( iElem )%coord, &
        &   LSFmat      = method%intpMat_forLSF%matArray            &
        &                       %val(posInIntpMatLSF),              &
        &   nSources    = nSourceElems,                             &
        &   nVals       = QQ                                        )

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

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

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

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

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

      ! Now write the resulting pdf in the current direction to the target
      ! element position
      do iVal = 1, QQ
        tState( ( targetelem-1)* nscalars+ ival) = tPDF(iVal)
      enddo

    end do ! indElem

  end subroutine fillFinerGhostsFromMe_quadIncomp
! ****************************************************************************** !


  ! ************************************************************************** !
  !> Fill fine ghost from coarse fluid by quadratic interpolation.
  !! 1. Compute moments for all source elements, save in momBuf
  !! 2. For each target, interpolate all moments using quadratic least square fit
  !! 3. Store target auxilary field
  !! 4. Compute viscosity on target element and compute source and target omega
  !! 5. Get nonEq scaling factor depeding on scheme layout and relaxation
  !! 6. Calculate Equilibrium and nonEquilibrium
  !! 7. Calculate target: Eq + Scale * nonEquilibrium
  !! This routine is used by 3D acoustic quadratic interpolation for
  !! incompressible LES model.
  !!
  !! This subroutine's interface must match the abstract interface definition
  !! [[intpRoutine]] in intp/[[mus_interpolate_header_module]].f90 in order to
  !! be callable via [[mus_interpolation_method_type:do_intp]] function pointer.
  subroutine fillFinerGhostsFromMe_quadIncompLES( 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) :: tPDF( layout%fStencil%QQ ) ! pdf to reconstruct from
    ! moments of the source elements' pdf
    real(kind=rk) :: srcMom( layout%fStencil%QQ, layout%fStencil%QQ )
    real(kind=rk) :: tMom(layout%fStencil%QQ)  ! target moment calculation
    integer :: QQ
    integer :: posInIntpMatLSF
    type(mus_fluid_type), pointer :: fluid
    real(kind=rk) :: nonEqScalingFacs(layout%fStencil%QQ)
    real(kind=rk) :: sOmegaKine, tOmegaKine, sVisc, tVisc
    real(kind=rk) :: sTurbVisc(layout%fStencil%QQ), tTurbVisc(1)
    integer :: nScalars, tOffset
    integer :: dens_pos, vel_pos(3)
    ! ---------------------------------------------------------------------------
    fluid => fieldProp(1)%fluid
    nScalars = varSys%nScalars
    QQ = layout%fStencil%QQ
    dens_pos = varSys%method%val(derVarPos(1)%density)%auxField_varPos(1)
    vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:3)

    sourceLevel = level
    targetLevel = level + 1

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

      iElem = targetList( indElem )
      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 position
        sourceElem = tLevelDesc%depFromCoarser( iElem ) &
          &                    %elem%val( iSourceElem )

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

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

      !interpolate moments to target element
      tMom(1:QQ) = mus_interpolate_quad3D_leastSq(                  &
        &   srcMom      = srcMom(1:QQ, 1:nSourceElems),             &
        &   targetCoord = tLevelDesc%depFromCoarser( iElem )%coord, &
        &   LSFmat      = method%intpMat_forLSF%matArray            &
        &                       %val(posInIntpMatLSF),              &
        &   nSources    = nSourceElems,                             &
        &   nVals       = QQ )

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

      !interpolate turbulent viscosity to target element
      tTurbVisc = mus_interpolate_quad3D_leastSq(                   &
        &   srcMom      = sTurbVisc(1:nSourceElems),                &
        &   targetCoord = tLevelDesc%depFromCoarser( iElem )%coord, &
        &   LSFmat      = method%intpMat_forLSF%matArray            &
        &                       %val(posInIntpMatLSF),              &
        &   nSources    = nSourceElems,                             &
        &   nVals       = 1                                         )

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

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

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

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

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

!      write(dbgUnit(1),*) 'Quad nonEqFacs', nonEqScalingFacs
      ! Convert moment to PDF
      tPDF = mus_intp_convertMomToPDF3D_incomp(                              &
        &                               moments          = tMom,             &
        &                               nonEqScalingFacs = nonEqScalingFacs, &
        &                               layout           = layout            )

      ! Now write the resulting pdf in the current direction to the target
      ! Element position
      do iVal = 1, QQ
        tState( ( targetelem-1)* nscalars+ ival) = tPDF(iVal)
      enddo

    end do ! indElem

  end subroutine fillFinerGhostsFromMe_quadIncompLES
! ****************************************************************************** !


  ! ************************************************************************** !
  !> Fill fine ghost from coarse fluid by quadratic interpolation.
  !! 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 3D 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_quad2D( 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) :: tPDF( layout%fStencil%QQ ) ! pdf to reconstruct from
    ! moments of the source elements' pdf
    real(kind=rk) :: srcMom( layout%fStencil%QQ, layout%fStencil%QQ )
    real(kind=rk) :: tMom(layout%fStencil%QQ)  ! target moment calculation
    integer :: QQ
    integer :: posInIntpMatLSF
    type(mus_fluid_type), pointer :: fluid
    real(kind=rk) :: sOmegaKine, tOmegaKine, tVisc, invRho
    real(kind=rk) :: nonEqScalingFacs(layout%fStencil%QQ)
    integer :: nScalars, tOffset
    integer :: dens_pos, vel_pos(3)
    ! ---------------------------------------------------------------------------
    fluid => fieldProp(1)%fluid
    nScalars = varSys%nScalars
    QQ = layout%fStencil%QQ
    dens_pos = varSys%method%val(derVarPos(1)%density)%auxField_varPos(1)
    vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:3)

    sourceLevel = level
    targetLevel = level + 1

    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 position
        sourceElem = tLevelDesc%depFromCoarser( iElem ) &
          &                    %elem%val( iSourceElem )

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

      !interpolate moments to target element
      tMom(1:QQ) = mus_interpolate_quad2D_leastSq(                  &
        &   srcMom      = srcMom(1:QQ, 1:nSourceElems),             &
        &   targetCoord = tLevelDesc%depFromCoarser( iElem )%coord, &
        &   LSFmat      = method%intpMat_forLSF%matArray            &
        &                       %val(posInIntpMatLSF),              &
        &   nSources    = nSourceElems,                             &
        &   nVals       = QQ                                        )

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

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

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

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

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

      ! Now write the resulting pdf in the current direction to the target
      ! Element position
      do iVal = 1, QQ
        tState( ( targetelem-1)* nscalars+ ival) = tPDF(iVal)
      enddo

    end do ! indElem

  end subroutine fillFinerGhostsFromMe_quad2D
! ****************************************************************************** !


  ! ************************************************************************** !
  !> Fill fine ghost from coarse fluid by quadratic interpolation.
  !! 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 3D acoustic quadratic interpolation for
  !! incompressible model.
  !!
  !! This subroutine's interface must match the abstract interface definition
  !! [[intpRoutine]] in intp/[[mus_interpolate_header_module]].f90 in order to
  !! be callable via [[mus_interpolation_method_type:do_intp]] function pointer.
  subroutine fillFinerGhostsFromMe_quad2DIncomp( 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) :: tPDF( layout%fStencil%QQ ) ! pdf to reconstruct from
    ! moments of the source elements' pdf
    real(kind=rk) :: srcMom( layout%fStencil%QQ, layout%fStencil%QQ )
    real(kind=rk) :: tMom(layout%fStencil%QQ)  ! target moment calculation
    integer :: QQ
    integer :: posInIntpMatLSF
    type(mus_fluid_type), pointer :: fluid
    real(kind=rk) :: nonEqScalingFacs(layout%fStencil%QQ)
    real(kind=rk) :: sOmegaKine, tOmegaKine, tVisc
    integer :: nScalars, tOffset
    integer :: dens_pos, vel_pos(3)
    ! ---------------------------------------------------------------------------
    fluid => fieldProp(1)%fluid
    nScalars = varSys%nScalars
    QQ = layout%fStencil%QQ
    dens_pos = varSys%method%val(derVarPos(1)%density)%auxField_varPos(1)
    vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:3)

    sourceLevel = level
    targetLevel = level + 1

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

      iElem = targetList( indElem )
      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 position
        sourceElem = tLevelDesc%depFromCoarser( iElem ) &
          &                    %elem%val( iSourceElem )

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

      end do

      !interpolate moments to target element
      tMom(1:QQ) = mus_interpolate_quad2D_leastSq(                  &
        &   srcMom      = srcMom(1:QQ, 1:nSourceElems),             &
        &   targetCoord = tLevelDesc%depFromCoarser( iElem )%coord, &
        &   LSFmat      = method%intpMat_forLSF%matArray            &
        &                       %val(posInIntpMatLSF),              &
        &   nSources    = nSourceElems,                             &
        &   nVals       = QQ )

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

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

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

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

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

      ! Now write the resulting pdf in the current direction to the target
      ! Element position
      do iVal = 1, QQ
        tState( ( targetelem-1)* nscalars+ ival) = tPDF(iVal)
      enddo

    end do ! indElem

  end subroutine fillFinerGhostsFromMe_quad2DIncomp
! ****************************************************************************** !

! ****************************************************************************** !
  !> Biquadratic interpolation for a vector quantity phi
  !!
  ! real(kind=rk), dimension(6,9),parameter  :: intpMaxtrix =  &
  !           reshape((/ &
  !   &  div2_9,  -div1_6,    0.0d0,    div1_6,   -div1_3,   0.0d0, &
  !   &  div2_9,    0.0d0,  -div1_6,   -div1_3,    div1_6,   0.0d0, &
  !   &  div2_9,   div1_6,    0.0d0,    div1_6,   -div1_3,   0.0d0, &
  !   &  div2_9,    0.0d0,   div1_6,   -div1_3,    div1_6,   0.0d0, &
  !   & -div1_9,  -div1_6,  -div1_6,    div1_6,    div1_6,  0.25d0, &
  !   & -div1_9,  -div1_6,   div1_6,    div1_6,    div1_6, -0.25d0, &
  !   & -div1_9,   div1_6,  -div1_6,    div1_6,    div1_6, -0.25d0, &
  !   & -div1_9,   div1_6,   div1_6,    div1_6,    div1_6,  0.25d0, &
  !   &  div5_9,    0.0d0,    0.0d0,   -div1_3,   -div1_3,   0.0d0  &
  !  /),(/6,9/))
  pure function mus_interpolate_quad2D_leastSq( srcMom, targetCoord, LSFmat, &
    &                                           nSources, nVals ) result(phi)
    ! ---------------------------------------------------------------------------
    !> number of quantities to interpolation
    integer, intent(in) :: nVals
    !> Number of source elements
    integer, intent(in) :: nSources
    !> source values of the momentum on the square corners
    real(kind=rk), intent(in) :: srcMom(nVals, nSources)
    !> interpolation location within the square
    real(kind=rk), intent(in) :: targetCoord(3)
    !> matrix for least square fit
    type(tem_matrix_type), intent(in) :: LSFmat
    !> interpolated value
    real(kind=rk) :: phi( nVals )
    ! ---------------------------------------------------------------------------
    real(kind=rk) :: a(6) ! Coefficients
    integer :: iVal
    ! ---------------------------------------------------------------------------
    ! We extract momentum information completely on the view of the source
    ! coordinate system
    ! Set the right hand side of the equation system
    ! Solve the problem, where b = rhs, x = coefficients
    ! A*x = b
    ! overdetermined, solve the least Square fit problem
    ! (A^T)A*x = (A^T)b
    ! x = ((A^T)A)^-1*(A^T)b
    ! Solve linear system of equation with inverted matrix
    do iVal = 1, nVals
      a(:) = matmul( LSFmat%A, srcMom(iVal, :) )

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

  end function mus_interpolate_quad2D_leastSq
! ****************************************************************************** !


! ****************************************************************************** !
  !> Triquadratic interpolation for a vector quantity phi
  !! Each phi corresponds to each moment
  ! matrixInv%A(1:10,1:19) =
  !      1     2     3     4     5     6     7     8     9    10    11    12    13    14    15    16    17    18    19
  ! [ 4/21  4/21  4/21  4/21  4/21  4/21 -1/21 -1/21 -1/21 -1/21 -1/21 -1/21 -1/21 -1/21 -1/21 -1/21 -1/21 -1/21   3/7]
  ! [-1/10     0     0  1/10     0     0     0     0     0     0 -1/10  1/10 -1/10  1/10 -1/10 -1/10  1/10  1/10     0]
  ! [    0 -1/10     0     0  1/10     0 -1/10 -1/10  1/10  1/10     0     0     0     0 -1/10  1/10 -1/10  1/10     0]
  ! [    0     0 -1/10     0     0  1/10 -1/10  1/10 -1/10  1/10 -1/10 -1/10  1/10  1/10     0     0     0     0     0]
  ! [ 1/42  -1/7  -1/7  1/42  -1/7  -1/7 -1/21 -1/21 -1/21 -1/21  5/42  5/42  5/42  5/42  5/42  5/42  5/42  5/42 -5/21]
  ! [ -1/7  1/42  -1/7  -1/7  1/42  -1/7  5/42  5/42  5/42  5/42 -1/21 -1/21 -1/21 -1/21  5/42  5/42  5/42  5/42 -5/21]
  ! [ -1/7  -1/7  1/42  -1/7  -1/7  1/42  5/42  5/42  5/42  5/42  5/42  5/42  5/42  5/42 -1/21 -1/21 -1/21 -1/21 -5/21]
  ! [    0     0     0     0     0     0     0     0     0     0     0     0     0     0   1/4  -1/4  -1/4   1/4     0]
  ! [    0     0     0     0     0     0   1/4  -1/4  -1/4   1/4     0     0     0     0     0     0     0     0     0]
  ! [    0     0     0     0     0     0     0     0     0     0   1/4  -1/4  -1/4   1/4     0     0     0     0     0]
  !
! ****************************************************************************** !
  function mus_interpolate_quad3D_leastSq( srcMom, targetCoord, LSFmat, &
    &                                           nSources, nVals ) result(phi)
    ! ---------------------------------------------------------------------------
    !> number of quantities to interpolation
    integer, intent(in) :: nVals
    !> Number of source elements
    integer, intent(in) :: nSources
    !> source values of the momentum on the square corners
    real(kind=rk), intent(in) :: srcMom(nVals, nSources)
    !> interpolation location within the square
    real(kind=rk), intent(in) :: targetCoord(3)
    !> matrix for least square fit
    type(tem_matrix_type), intent(in) :: LSFmat
    !> interpolated value
    real(kind=rk) :: phi( nVals )
    ! ---------------------------------------------------------------------------
    real(kind=rk) :: a(10) ! Coefficients
    integer :: iVal
    ! ---------------------------------------------------------------------------
    ! We extract momentum information completely on the view of the source
    ! coordinate system
    ! Set the right hand side of the equation system
    ! Solve the problem, where b = rhs, x = coefficients
    ! A*x = b
    ! overdetermined, solve the least Square fit problem
    ! (A^T)A*x = (A^T)b
    ! x = ((A^T)A)^-1*(A^T)b
    ! Solve linear system of equation with inverted matrix
    do iVal = 1, nVals
      a(:) = matmul( LSFmat%A, srcMom(iVal, :) )

      ! Evaluate the bubble function with the above calculated coefficients
      ! m_pi(x) = a1+a2*xcoord+a3*ycoord+a4*xcoord*xcoord+a5*ycoord*ycoord+a6*xcoord*ycoord;
      phi( iVal ) =   a( 1) &
        &           + a( 2)*targetCoord(c_x)                  &
        &           + a( 3)*targetCoord(c_y)                  &
        &           + a( 4)*targetCoord(c_z)                  &
        &           + a( 5)*targetCoord(c_x)*targetCoord(c_x) &
        &           + a( 6)*targetCoord(c_y)*targetCoord(c_y) &
        &           + a( 7)*targetCoord(c_z)*targetCoord(c_z) &
        &           + a( 8)*targetCoord(c_x)*targetCoord(c_y) &
        &           + a( 9)*targetCoord(c_y)*targetCoord(c_z) &
        &           + a(10)*targetCoord(c_z)*targetCoord(c_x)
    enddo

  end function mus_interpolate_quad3D_leastSq
! ****************************************************************************** !


! ****************************************************************************** !
  !> 2D Biquadratic interpolation for fill fine from coarse.
  !! array f contains function values at 9 coarse elements.
  !! they are arranged as:
  !!   7  8  9
  !!   4  5  6
  !!   1  2  3
  !! position of element 1 is set as original, i.e. (0,0),
  !!             element 9 is set as (1,1)
  !! x,y are the coordinates of target position.
  !! fOut is the interpolated value at target position
  !!
  pure function biquadratic_interpolate( f, x, y ) result( fOut )
    ! ---------------------------------------------------------------------------
    real(kind=rk), intent(in) :: f(9)
    real(kind=rk), intent(in) :: x, y
    ! ---------------------------------------------------------------------------
    real(kind=rk) :: fOut
    real(kind=rk) :: a0, a1, a2, b0, b1, b2, c0, c1, c2, f1, f2, f3, m0, m1, m2
    ! ---------------------------------------------------------------------------

    a0 =          f(1)
    a1 = -3._rk * f(1) + 4._rk * f(2) -         f(3)
    a2 =  2._rk * f(1) - 4._rk * f(2) + 2._rk * f(3)

    b0 =          f(4)
    b1 = -3._rk * f(4) + 4._rk * f(5) -         f(6)
    b2 =  2._rk * f(4) - 4._rk * f(5) + 2._rk * f(6)

    c0 =          f(7)
    c1 = -3._rk * f(7) + 4._rk * f(8) -         f(9)
    c2 =  2._rk * f(7) - 4._rk * f(8) + 2._rk * f(9)

    f1 = a0 + a1 * x + a2 * x * x
    f2 = b0 + b1 * x + b2 * x * x
    f3 = c0 + c1 * x + c2 * x * x

    m0 =          f1
    m1 = -3._rk * f1 + 4._rk * f2 -         f3
    m2 =  2._rk * f1 - 4._rk * f2 + 2._rk * f3

    fOut = m0 + m1 * y + m2 * y * y

  end function biquadratic_interpolate
! ****************************************************************************** !

end module mus_interpolate_quadratic_module
! ****************************************************************************** !