! Copyright (c) 2011-2013 Manuel Hasert <m.hasert@grs-sim.de>
! Copyright (c) 2011 Jan Hueckelheim <j.hueckelheim@grs-sim.de>
! Copyright (c) 2011-2012 Harald Klimach <harald.klimach@uni-siegen.de>
! Copyright (c) 2011-2013, 2015, 2018-2020 Kannan Masilamani <kannan.masilamani@uni-siegen.de>
! Copyright (c) 2011-2014 Simon Zimny <s.zimny@grs-sim.de>
! Copyright (c) 2012, 2014-2016 Jiaxing Qi <jiaxing.qi@uni-siegen.de>
! Copyright (c) 2014 Kartik Jain <kartik.jain@uni-siegen.de>
! Copyright (c) 2016-2017 Tobias Schneider <tobias1.schneider@student.uni-siegen.de>
! Copyright (c) 2018 Raphael Haupt <raphael.haupt@uni-siegen.de>
! Copyright (c) 2019-2020 Peter Vitt <peter.vitt2@uni-siegen.de>
! Copyright (c) 2019 Jana Gericke <jana.gericke@uni-siegen.de>
! Copyright (c) 2021-2022 Gregorio Gerardo Spinelli <gregoriogerardo.spinelli@dlr.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 <a id="interpolation-of"></a>
!!
!! 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)
!!
! 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.
module mus_interpolate_module

  ! include treelm modules
  use mpi
  use env_module,            only: rk
  use tem_aux_module,        only: tem_abort
  use tem_debug_module,      only: main_debug, dbgUnit
  use tem_logging_module,    only: logUnit
  use tem_construction_module, only: tem_levelDesc_type, depSource_append
  use tem_dyn_array_module,  only: dyn_intArray_type, init, append, truncate
  use tem_grow_array_module, only: grw_intArray_type, init, append, truncate
  use tem_element_module,    only: eT_ghostFromCoarser, eT_ghostFromFiner
  use tem_stencil_module,    only: tem_stencilHeader_type
  use tem_matrix_module,     only: init, append
  use tem_tools_module,      only: tem_horizontalSpacer
  use tem_property_module,   only: prp_fluid, prp_fineGhostClosestToFluid

  ! include musubi modules
  use mus_field_prop_module,            only: mus_field_prop_type
  use mus_interpolate_header_module,    only: mus_interpolation_type,         &
    &                                         mus_interpolation_config_type,  &
    &                                         mus_interpolation_stencil_type, &
    &                                         weighted_average, linear,       &
    &                                         quadratic
  use mus_interpolate_debug_module,     only: do_nothing, do_nothing_arbi
  use mus_interpolate_average_module,   only:           &
    & fillArbiMyGhostsFromFiner_avg,                    &
    & fillMyGhostsFromFiner_avgIncomp_feq_fneq,         &
    & fillMyGhostsFromFiner_avgIncompLES_feq_fneq,      &
    & fillMyGhostsFromFiner_avg_feq_fneq,               &
    & fillMyGhostsFromFiner_avgLES_feq_fneq,            &
    & fillMyGhostsFromFiner_avg2DIncomp_feq_fneq,       &
    & fillMyGhostsFromFiner_avg2D_feq_fneq,             &
    & fillArbiFinerGhostsFromMe_weighAvg,               &
    & fillFinerGhostsFromMe_weighAvgIncomp_feq_fneq,    &
    & fillFinerGhostsFromMe_weighAvgIncompLES_feq_fneq, &
    & fillFinerGhostsFromMe_weighAvg_feq_fneq,          &
    & fillFinerGhostsFromMe_weighAvgLES_feq_fneq,       &
    & fillFinerGhostsFromMe_weighAvg2DIncomp_feq_fneq,  &
    & fillFinerGhostsFromMe_weighAvg2D_feq_fneq
  use mus_interpolate_linear_module,    only:         &
    & fillArbiFinerGhostsFromMe_linear,               &
    & fillArbiFinerGhostsFromMe_linear2D,             &
    & fillFinerGhostsFromMe_linear_feq_fneq,          &
    & fillFinerGhostsFromMe_linearLES_feq_fneq,       &
    & fillFinerGhostsFromMe_linear2D_feq_fneq,        &
    & fillFinerGhostsFromMe_linearIncomp_feq_fneq,    &
    & fillFinerGhostsFromMe_linearIncompLES_feq_fneq, &
    & fillFinerGhostsFromMe_linear2DIncomp_feq_fneq
  use mus_interpolate_quadratic_module, only:       &
    & fillArbiFinerGhostsFromMe_quad,               &
    & fillArbiFinerGhostsFromMe_quad2D,             &
    & fillFinerGhostsFromMe_quad_feq_fneq,          &
    & fillFinerGhostsFromMe_quadLES_feq_fneq,       &
    & fillFinerGhostsFromMe_quadIncomp_feq_fneq,    &
    & fillFinerGhostsFromMe_quadIncompLES_feq_fneq, &
    & fillFinerGhostsFromMe_quad2D_feq_fneq,        &
    & fillFinerGhostsFromMe_quad2DIncomp_feq_fneq
  use mus_scheme_header_module,     only: mus_scheme_header_type
  use mus_scheme_layout_module,     only: mus_scheme_layout_type

  implicit none

  private

  public :: mus_init_interpolate
  public :: mus_init_levelDescIntpArrays
  public :: mus_dump_levelDescIntp_nElems
  public :: mus_intp_update_depFromCoarser
  public :: mus_intp_update_depFromFiner

  contains

  ! ************************************************************************** !
  !> This subroutine initialzes the interpolation
  !!
  !! - setting the fillMineFromFiner and fillFinerFromMe function pointers
  !!   to the interpolation function chosen by the user and
  !! - pre-calculating the weights based on the distances between source
  !!   and target nodes
  !! - allocate buffers if required
  !!
  subroutine mus_init_interpolate( intp, levelDesc, schemeHeader, stencil, &
    &                              minLevel, maxLevel, fieldProp)
    ! ---------------------------------------------------------------------------
    !> interpolation type
    type(mus_interpolation_type), intent(inout) :: intp
    integer,                        intent(in) :: minLevel, maxLevel
    !> level descriptor is actually used
    type(tem_levelDesc_type), intent(inout) :: levelDesc(minLevel:maxLevel)
    !> the stencil header
    type(tem_stencilHeader_type), intent(in) :: stencil
    !> scheme header
    type(mus_scheme_header_type), intent(in) :: schemeHeader
    !> field properties
    type(mus_field_prop_type), intent(in) :: fieldProp(:)
    ! ---------------------------------------------------------------------------
    integer :: iOrder
    ! ---------------------------------------------------------------------------

    write(logUnit(1),"(A)") 'Initialize interpolation: '//trim(intp%config%method)
    write(logUnit(1),"(A)") '  with scheme: '//trim(schemeHeader%relaxation)
    flush(logUnit(1))

    select case (trim(schemeHeader%kind))
    case ('fluid')
      if (fieldProp(1)%fluid%turbulence%active) then
        call assign_intp_fluidLES( intp  = intp,          &
          &                        nDims = stencil%nDims, &
          &                        QQ    = stencil%QQ     )
      else
        call assign_intp_fluid( intp  = intp,          &
          &                     nDims = stencil%nDims, &
          &                     QQ    = stencil%QQ     )
      end if
    case ('fluid_incompressible')
      if (fieldProp(1)%fluid%turbulence%active) then
        call assign_intp_fluidIncompLES( intp  = intp,          &
          &                              nDims = stencil%nDims, &
          &                              QQ    = stencil%QQ     )
      else
        call assign_intp_fluidIncomp( intp  = intp,          &
          &                           nDims = stencil%nDims, &
          &                           QQ    = stencil%QQ     )
      end if
    case default
        write(logUnit(1),"(A)") ' Multilevel interpolation is not supported' &
          &                     // 'for scheme: '//trim(schemeHeader%kind)
        call tem_abort()
    end select

    ! Used for debugging. Configurable variable, Default is False
    if (intp%config%noIntpFromFiner) then
      intp%fillMineFromFiner%do_intp => do_nothing
      intp%fillMineFromFiner%do_intpArbiVal => do_nothing_arbi
    end if
    if (intp%config%noIntpFromCoarser) then
      do iOrder = 0, intp%config%order
        intp%fillFinerFromMe(iOrder)%do_intp => do_nothing
        intp%fillFinerFromMe(iOrder)%do_intpArbiVal => do_nothing_arbi
      end do
    end if

    write(logUnit(1),*) 'Done with interpolation method.'
  end subroutine mus_init_interpolate
  ! ************************************************************************** !

! ****************************************************************************** !
  !> Set up interpolation routines for fluid (weakly compressible) scheme
  subroutine assign_intp_fluid( intp, nDims, QQ )
    ! ---------------------------------------------------------------------------
    !> interpolation type
    type(mus_interpolation_type), intent(inout) :: intp
    !> Number of dimensions
    integer,                      intent(in) :: nDims
    !> number of stencil directions
    integer,                      intent(in) :: QQ
    ! ---------------------------------------------------------------------------
    integer :: ii
    ! ---------------------------------------------------------------------------
    write(logUnit(3),*)   ' Assign intp function pointer for fluid:'

    ! assign pointers to fillCoarser
    intp%fillMineFromFiner%do_intpArbiVal => fillArbiMyGhostsFromFiner_avg
    select case(nDims)
    case(2)
      write(logUnit(3),*) 'fillCoarser average 2D'
      intp%fillMineFromFiner%do_intp => fillMyGhostsFromFiner_avg2D_feq_fneq

    case(3)
      write(logUnit(3),*) 'fillCoarser average 3D'
      intp%fillMineFromFiner%do_intp => fillMyGhostsFromFiner_avg_feq_fneq

    case default
      call tem_abort('Unsupported nDims')
    end select

    ! assign pointers to fillFiner
    do ii = 0, intp%config%order
      select case(ii)
      case(weighted_average)
        intp%fillFinerFromMe(ii)%do_intpArbiVal &
          & => fillArbiFinerGhostsFromMe_weighAvg
        select case(nDims)
        case(2)
          write(logUnit(3),*) 'fillFiner weighted average 2D'
          intp%fillFinerFromMe(ii)                           &
            & %do_intp => fillFinerGhostsFromMe_weighAvg2D_feq_fneq

        case(3)
          write(logUnit(3),*) 'fillFiner weighted average 3D'
          intp%fillFinerFromMe(ii)%do_intp &
            & => fillFinerGhostsFromMe_weighAvg_feq_fneq
        end select

      case(linear)
        select case(nDims)
        case(2)
          write(logUnit(3),*) 'fillFiner linear 2D'
          intp%fillFinerFromMe(ii)%do_intpArbiVal    &
            & => fillArbiFinerGhostsFromMe_linear2D
          intp%fillFinerFromMe(ii)                         &
            & %do_intp => fillFinerGhostsFromMe_linear2D_feq_fneq

        case(3)
          write(logUnit(3),*) 'fillFiner linear 3D'
          intp%fillFinerFromMe(ii)%do_intpArbiVal &
            & => fillArbiFinerGhostsFromMe_linear
          intp%fillFinerFromMe(ii)%do_intp => &
            & fillFinerGhostsFromMe_linear_feq_fneq
        end select

      case(quadratic)
        select case(nDims)
        case(2)
          write(logUnit(3),*) 'fillFiner quadratic 2D'
          intp%fillFinerFromMe(ii)%do_intpArbiVal &
            & => fillArbiFinerGhostsFromMe_quad2D
          intp%fillFinerFromMe(ii)                       &
            & %do_intp => fillFinerGhostsFromMe_quad2D_feq_fneq

        case(3)
          write(logUnit(3),*) 'fillFiner quadratic 3D'
          intp%fillFinerFromMe(ii)%do_intpArbiVal &
            & => fillArbiFinerGhostsFromMe_quad
          intp%fillFinerFromMe(ii)%do_intp => fillFinerGhostsFromMe_quad_feq_fneq

        end select

      case default
        call tem_abort('Unknown interpolation method')
      end select

    end do

  end subroutine assign_intp_fluid
! ****************************************************************************** !

! ****************************************************************************** !
  !> Set up interpolation routines for fluid_incomp scheme
  subroutine assign_intp_fluidIncomp( intp, nDims, QQ )
    ! ---------------------------------------------------------------------------
    !> interpolation type
    type(mus_interpolation_type), intent(inout) :: intp
    !> Number of dimensions
    integer,                      intent(in) :: nDims
    !> number of stencil directions
    integer,                      intent(in) :: QQ
    ! ---------------------------------------------------------------------------
    integer :: ii
    ! ---------------------------------------------------------------------------
    write(logUnit(3),*)   ' Assign intp function pointer for fluid_incomp:'

    ! assign pointers to fillCoarser
    intp%fillMineFromFiner%do_intpArbiVal => fillArbiMyGhostsFromFiner_avg
    select case(nDims)
    case(2)
      write(logUnit(3),*) 'fillCoarser average 2D'
      intp%fillMineFromFiner%do_intp => fillMyGhostsFromFiner_avg2DIncomp_feq_fneq

    case(3)
      write(logUnit(3),*) 'fillCoarser average 3D'
      intp%fillMineFromFiner%do_intp => fillMyGhostsFromFiner_avgIncomp_feq_fneq

    end select

    ! assign pointers to fillFiner
    do ii = 0, intp%config%order
      select case(ii)
      case(weighted_average)
        intp%fillFinerFromMe(ii)%do_intpArbiVal &
          & => fillArbiFinerGhostsFromMe_weighAvg
        select case(nDims)
        case(2)
          write(logUnit(3),*) 'fillFiner weighted average 2D'
          intp%fillFinerFromMe(ii)                             &
            & %do_intp => fillFinerGhostsFromMe_weighAvg2DIncomp_feq_fneq

        case(3)
          write(logUnit(3),*) 'fillFiner weighted average 3D'
          intp%fillFinerFromMe(ii)                           &
            & %do_intp => fillFinerGhostsFromMe_weighAvgIncomp_feq_fneq

        end select

      case(linear)
        select case(nDims)
        case(2)
          write(logUnit(3),*) 'fillFiner linear 2D'
          intp%fillFinerFromMe(ii)%do_intpArbiVal    &
            & => fillArbiFinerGhostsFromMe_linear2D
          intp%fillFinerFromMe(ii)                           &
            & %do_intp => fillFinerGhostsFromMe_linear2DIncomp_feq_fneq
        case(3)
          write(logUnit(3),*) 'fillFiner linear 3D'
          intp%fillFinerFromMe(ii)%do_intpArbiVal &
            & => fillArbiFinerGhostsFromMe_linear
          intp%fillFinerFromMe(ii)%do_intp => fillFinerGhostsFromMe_linearIncomp_feq_fneq

        end select

      case(quadratic)
        select case(nDims)
        case(2)
          write(logUnit(3),*) 'fillFiner quadratic 2D'
          intp%fillFinerFromMe(ii)%do_intpArbiVal &
            & => fillArbiFinerGhostsFromMe_quad2D
          intp%fillFinerFromMe(ii)                         &
            & %do_intp => fillFinerGhostsFromMe_quad2DIncomp_feq_fneq

        case(3)
          write(logUnit(3),*) 'fillFiner quadratic 3D'
          intp%fillFinerFromMe(ii)%do_intpArbiVal => fillArbiFinerGhostsFromMe_quad
          intp%fillFinerFromMe(ii)%do_intp => fillFinerGhostsFromMe_quadIncomp_feq_fneq

        end select

      case default
        call tem_abort('Unknown interpolation method')
      end select

    end do

  end subroutine assign_intp_fluidIncomp
! ****************************************************************************** !

  !> Set up interpolation routines for fluid (weakly compressible) scheme and
  !! turbulence active
  subroutine assign_intp_fluidLES( intp, nDims, QQ )
    ! ---------------------------------------------------------------------------
    !> interpolation type
    type(mus_interpolation_type), intent(inout) :: intp
    !> Number of dimensions
    integer,                      intent(in) :: nDims
    !> number of stencil directions
    integer,                      intent(in) :: QQ
    ! ---------------------------------------------------------------------------
    integer :: ii
    ! ---------------------------------------------------------------------------
    write(logUnit(3),*)   ' Assign intp function pointer for fluid LES:'

    ! assign pointers to fillCoarser
    intp%fillMineFromFiner%do_intpArbiVal => fillArbiMyGhostsFromFiner_avg
    select case(nDims)
    case(3)
      write(logUnit(3),*) 'fillCoarser average 3D LES'
      intp%fillMineFromFiner%do_intp => fillMyGhostsFromFiner_avgLES_feq_fneq

    case default
      call tem_abort('Unsupported nDims')
    end select

    ! assign pointers to fillFiner
    do ii = 0, intp%config%order
      select case(ii)
      case(weighted_average)
        intp%fillFinerFromMe(ii)%do_intpArbiVal &
          & => fillArbiFinerGhostsFromMe_weighAvg
        select case(nDims)
        case(3)
          write(logUnit(3),*) 'fillFiner weighted average 3D LES'
          intp%fillFinerFromMe(ii)%do_intp &
            & => fillFinerGhostsFromMe_weighAvgLES_feq_fneq

        end select

      case(linear)
        select case(nDims)
        case(3)
          write(logUnit(3),*) 'fillFiner linear 3D LES'
          intp%fillFinerFromMe(ii)%do_intpArbiVal &
            & => fillArbiFinerGhostsFromMe_linear
          intp%fillFinerFromMe(ii)%do_intp => fillFinerGhostsFromMe_linearLES_feq_fneq

        end select

      case(quadratic)
        select case(nDims)
        case(3)
          intp%fillFinerFromMe(ii)%do_intpArbiVal &
            & => fillArbiFinerGhostsFromMe_quad
          write(logUnit(3),*) 'fillFiner quadratic 3D LES'
          intp%fillFinerFromMe(ii)%do_intp => fillFinerGhostsFromMe_quadLES_feq_fneq

        end select
      case default
        call tem_abort('Unknown interpolation method')
      end select

    end do

  end subroutine assign_intp_fluidLES
! ****************************************************************************** !

  !> Set up interpolation routines for fluid (incompressible) scheme and
  !! turbulence active
  subroutine assign_intp_fluidIncompLES( intp, nDims, QQ )
    ! ---------------------------------------------------------------------------
    !> interpolation type
    type(mus_interpolation_type), intent(inout) :: intp
    !> Number of dimensions
    integer,                      intent(in) :: nDims
    !> number of stencil directions
    integer,                      intent(in) :: QQ
    ! ---------------------------------------------------------------------------
    integer :: ii
    ! ---------------------------------------------------------------------------
    write(logUnit(3),*)   ' Assign intp function pointer for fluid LES:'

    ! assign pointers to fillCoarser
    intp%fillMineFromFiner%do_intpArbiVal => fillArbiMyGhostsFromFiner_avg
    select case(nDims)
    case(3)
      write(logUnit(3),*) 'fillCoarser average 3D LES'
      intp%fillMineFromFiner%do_intp => fillMyGhostsFromFiner_avgIncompLES_feq_fneq

    case default
      call tem_abort('Unsupported nDims')
    end select

    ! assign pointers to fillFiner
    do ii = 0, intp%config%order
      select case(ii)
      case(weighted_average)
        intp%fillFinerFromMe(ii)%do_intpArbiVal &
          & => fillArbiFinerGhostsFromMe_weighAvg
        select case(nDims)
        case(3)
          write(logUnit(3),*) 'fillFiner weighted average 3D LES'
          intp%fillFinerFromMe(ii)%do_intp &
            & => fillFinerGhostsFromMe_weighAvgIncompLES_feq_fneq

        end select

      case(linear)
        select case(nDims)
        case(3)
          write(logUnit(3),*) 'fillFiner linear 3D LES'
          intp%fillFinerFromMe(ii)%do_intpArbiVal &
            & => fillArbiFinerGhostsFromMe_linear
          intp%fillFinerFromMe(ii)%do_intp => fillFinerGhostsFromMe_linearIncompLES_feq_fneq

        end select

      case(quadratic)
        select case(nDims)
        case(3)
          intp%fillFinerFromMe(ii)%do_intpArbiVal &
            & => fillArbiFinerGhostsFromMe_quad
          write(logUnit(3),*) 'fillFiner quadratic 3D LES'
          intp%fillFinerFromMe(ii)%do_intp => fillFinerGhostsFromMe_quadIncompLES_feq_fneq

        end select
      case default
        call tem_abort('Unknown interpolation method')
      end select

    end do

  end subroutine assign_intp_fluidIncompLES
! ****************************************************************************** !

! ****************************************************************************** !
  !> Initialize levelwise ghost and source list for interpolation
  subroutine mus_init_levelDescIntpArrays( intp, levelDesc, minLevel, maxLevel )
    ! --------------------------------------------------------------------------
    !> interpolation type
    type(mus_interpolation_type), intent(inout)  :: intp
    integer, intent(in) :: minLevel
    integer, intent(in) :: maxLevel
    !> levelDesc
    type(tem_levelDesc_type), intent(inout) :: levelDesc(minLevel:maxLevel)
    ! --------------------------------------------------------------------------
    integer :: iLevel
    ! --------------------------------------------------------------------------
    write(logUnit(10),*) 'Initiale levelDesc intp arrays'
    do iLevel = minLevel, maxLevel
      call init_intpArraysPerLevel(                                   &
        & sourceFromCoarser    = levelDesc(iLevel)%sourceFromCoarser, &
        & intpFromCoarser      = levelDesc(iLevel)%intpFromCoarser,   &
        & sourceFromFiner      = levelDesc(iLevel)%sourceFromFiner,   &
        & intpFromFiner        = levelDesc(iLevel)%intpFromFiner,     &
        & nGhostFromCoarser    = levelDesc(iLevel)%elem               &
        &                        %nElems( eT_ghostFromCoarser ),      &
        & nGhostFromFiner      = levelDesc(iLevel)%elem               &
        &                        %nElems( eT_ghostFromFiner ),        &
        & intp_order           = intp%config%order,                   &
        & nMaxSourcesFromFiner = intp%fillMineFromFiner%nMaxSources,  &
        & iLevel               = iLevel                               )
    end do !iLevel
  end subroutine mus_init_levelDescIntpArrays
! ****************************************************************************** !


! ****************************************************************************** !
  subroutine init_intpArraysPerLevel( sourceFromCoarser, intpFromCoarser,      &
    &                                 sourceFromFiner, intpFromFiner,          &
    &                                 nGhostFromCoarser, nGhostFromFiner,      &
    &                                 intp_order, nMaxSourcesFromFiner, iLevel )
    ! --------------------------------------------------------------------------
    !> dynamic array of source elements from coarser
    type(dyn_intArray_type), intent(out) :: sourceFromCoarser
    !> growing array of ghost elements from coarser for different intp order
    type(grw_intArray_type), allocatable, intent(out) :: intpFromCoarser(:)
    !> dynamic array of source elements from finer
    type(dyn_intArray_type), intent(out) :: sourceFromFiner
    !> growing array of ghost elements from finer
    type(grw_intArray_type), intent(out) :: intpFromFiner
    !> Number of ghost from coarser on this level
    integer, intent(in) :: nGhostFromCoarser
    !> Number of ghost from finer on this level
    integer, intent(in) :: nGhostFromFiner
    !> current level
    integer, intent(in) :: iLevel
    !> interpolation order defined by user
    integer, intent(in) :: intp_order
    !> nMaxSources for fillMineFromFiner
    integer, intent(in) :: nMaxSourcesFromFiner
    ! --------------------------------------------------------------------------
    integer :: iOrder
    ! --------------------------------------------------------------------------
    write(logUnit(10),*) 'Initiale levelDesc intp array on level:', iLevel

    ! fromCoarser
    call init( me = sourceFromCoarser, length = nGhostFromCoarser )

    allocate(intpFromCoarser(0:intp_order))
    ! initialize only user defined order with nGhostFromCoarser since
    ! order below user defined order are fall back which will be small in size
    call init( me     = intpFromCoarser(intp_order), &
      &        length = nGhostFromCoarser            )
    do iOrder = 0, intp_order-1
      call init( me = intpFromCoarser(iOrder) )
    end do

    ! fromFiner
    call init( me     = sourceFromFiner,                       &
      &        length = nGhostFromFiner * nMaxSourcesFromFiner )

    call init( me  = intpFromFiner, length = nGhostFromFiner )

  end subroutine init_intpArraysPerLevel
! ****************************************************************************** !


! ****************************************************************************** !
  !> This routine dumps global nElems in intpFromCoarser, intpFromFiner,
  !! sourcesFromCoarser ans sourcesFromFiner
  subroutine mus_dump_levelDescIntp_nElems(intp, levelDesc, minLevel, maxLevel,&
    &                                      root, comm                          )
    ! --------------------------------------------------------------------------
    !> interpolation type
    type(mus_interpolation_type), intent(in)  :: intp
    integer, intent(in) :: minLevel
    integer, intent(in) :: maxLevel
    !> levelDesc
    type(tem_levelDesc_type), intent(inout) :: levelDesc(minLevel:maxLevel)
    integer, intent(in) :: root !< root process
    integer, intent(in) :: comm !< mpi communicator
    ! --------------------------------------------------------------------------
    integer :: iLevel, iOrder
    ! amount of High and Low Oder Ghost elements form Fine and Coarser
    integer :: nGFC(0:intp%config%order), tSFC, nGFF, tSFF, iErr
    ! ---------------------------------------------------------------------------
    do iLevel = minLevel, maxLevel
      call truncate( me = levelDesc(iLevel)%sourceFromCoarser    )
      do iOrder = 0, intp%config%order
        call truncate( me = levelDesc(iLevel)%intpFromCoarser(iOrder) )
      end do

      call truncate( me = levelDesc(iLevel)%sourceFromFiner    )
      call truncate( me = levelDesc(iLevel)%intpFromFiner )

      ! get ghostFromCoarser for every order ( nGFC )
      do iOrder = 0, intp%config%order
        call MPI_REDUCE( levelDesc(iLevel)%intpFromCoarser(iOrder)%nVals, &
          &              nGFC(iOrder), 1, MPI_INTEGER, mpi_sum,           &
          &              root, comm, iErr)
      end do

      ! get total sourceFromCoarser ( tSFC )
      call MPI_REDUCE( levelDesc( iLevel )%sourceFromCoarser%nVals, &
        &              tSFC, 1, MPI_INTEGER, mpi_sum,               &
        &              root, comm, iErr)

      ! get ghostFromFiner ( nGFF )
      call MPI_REDUCE( levelDesc( iLevel )%intpFromFiner%nVals, &
        &              nGFF, 1, MPI_INTEGER, mpi_sum,           &
        &              root, comm, iErr)

      ! get total sourceFromFiner ( tSFF )
      call MPI_REDUCE( levelDesc( iLevel )%sourceFromFiner%nVals, &
        &              tSFF, 1, MPI_INTEGER, mpi_sum,             &
        &              root, comm, iErr)

      write(logUnit(3),"(A,I0)") '                       level: ', iLevel
      do iOrder = 0, intp%config%order
        write(logUnit(3),"(A,I0)") '         Interpolation order: ', iOrder
        write(logUnit(3),"(A,I0)") '              ghostFromCoarser: ', &
          &                        nGFC(iOrder)
      end do
      write(logUnit(3),"(A,I0)") '     total sourceFromCoarser: ', tSFC
      write(logUnit(3),"(A,I0)") '              ghostFromFiner: ', nGFF
      write(logUnit(3),"(A,I0)") '     total   sourceFromFiner: ', tSFF
    end do ! iLevel = minLevel, maxLevel
    ! ------------------------------------------------------------------------

  end subroutine mus_dump_levelDescIntp_nElems
! ****************************************************************************** !


! ****************************************************************************** !
  !> The required source elements for ghost from coarser elements are
  !! identified in this routine. Moreover, the weights for each sources based
  !! on distance are calculated. \n
  !!
  !! The depdendencies for ghostFromCoarser elements initially only include its
  !! parent element.
  !! Additional dependencies Y for the element x are found which were not
  !! included in the construction routine before hand. However, they are
  !! defined to be inherently part of the stencil and hence exist locally and
  !! can simply be added to the dependencies. Two different ways are
  !! implemented.
  !! Full interpolation stencil
  !!   the complete compute stencil is included into the dependency list
  !!```
  !!  Y       Y       Y
  !!
  !!        x   o
  !!  Y       O       Y
  !!        o   o
  !!
  !!  Y       Y       Y
  !!```
  !!
  !! Note: The view here is according to how the interpolation will be done.
  !!       The coarse level does the interpolation for the finer level.
  !!       Thus, we access the information from the levelDesc with the
  !!       targetLevel and sourceLevel perspective.
  !!
  !! Steps of finding sources:
  !! 1.  loop over all neighbors: iNeigh = 1, QQ
  !! 2.  get neighbor posInTotal: elemPos = neigh(1)%nghElems
  !! 3.  if elemPos > 0,
  !!       add this source into depFromCoarser%elem
  !!       add this source into sourceFromCoarser
  !!       save its position in sourceFromCoarser into depFromCoarser%elemBuffer
  !!       calculate its weight based on relative coordinates in stencil
  !!       weight is inverse proportional to the distance from source to target
  !!       weight is normalized at last, so that sum of weights equals 1
  !! 4.  save its offset on stencil
  !!
  !! After finding sources, append target to low or high interpolation list
  !! based on if all required sources are found.
  !!
  subroutine mus_intp_update_depFromCoarser( intp, levelDesc, stencil, &
    &                                        minLevel, maxLevel )
    ! ---------------------------------------------------------------------------
    !> interpolation type
    type(mus_interpolation_type), intent(inout)  :: intp
    integer,                        intent(in) :: minLevel
    integer,                        intent(in) :: maxLevel
    !> Level Descriptor
    type( tem_levelDesc_type )                 :: levelDesc(minLevel:maxLevel)
    type( tem_stencilHeader_type ), intent(in) :: stencil
    ! ---------------------------------------------------------------------------
    integer :: sourceLevel  ! level of source elements
    integer :: targetLevel  ! level of target elements
    integer :: targetElem   ! position of target element in total list
    integer :: iElem, iNeigh
    integer :: nGhostElems
    integer :: nFoundSources
    integer :: posInTotal, parentPosInTotal
    real(kind=rk) :: targetBary(3) ! barycenter of current target element
    integer :: mySources(stencil%QQ), myNeighDir(stencil%QQ)
    ! ---------------------------------------------------------------------------
    integer :: posInMatArray, intpOrder, childNum
    integer :: iOrder
    logical :: success
    integer :: sourceElem, iSourceElem
    ! ---------------------------------------------------------------------------

    write(logUnit(1),"(A)") 'Updating dependences of ghost from coarser'

! DEBUG output  ----------------------------------------------------------------
call tem_horizontalSpacer( fUnit = dbgUnit(1), before = 1 )
write(dbgUnit(3),"(A)") "Inside mus_intp_update_depFromCoarser routine"
write(dbgUnit(3),"(A)") "Up to now, for each ghostFromCoarser we only have its parent."
write(dbgUnit(3),"(A)") "Here we are trying to find all neighbors of this parent."
write(dbgUnit(3),"(A,I0)") ''
! DEBUG output  ----------------------------------------------------------------

    do sourceLevel = minLevel, maxLevel - 1

      targetLevel = sourceLevel + 1

      nGhostElems = levelDesc( targetLevel )%elem%nElems( eT_ghostFromCoarser )
write(dbgUnit(4),"(A,I0)") " source level: ", sourceLevel
write(dbgUnit(4),"(A,I0)") " target level: ", targetLevel
write(dbgUnit(4),"(A,I0)") ' nGhostFromCoarser: ', nGhostElems

      !!suppChild = 1
      !do iElem = 1, nGhostElems
      !  !childNum = levelDesc(targetLevel)%depFromCoarser(iElem)%childNum
      !  !allocate(levelDesc( targetLevel )%depFromCoarser( iElem )%bitmask(stencil%QQ))
      !  !levelDesc( targetLevel )%depFromCoarser( iElem )%bitmask = .true.
      !
      !  !if (suppChild == 1) supp_bitmask = .true.
      !
      !  targetElem = iElem + levelDesc( targetLevel )%             &
      !    &                                offset( 1, eT_ghostFromCoarser)

!write(dbgUnit(5),"(A, I0)") '      iElem: ', iElem
!!write(dbgUnit(5),"(A, I0)") '      suppChild: ', suppChild
!!write(dbgUnit(5),"(A, I0)") '      childNum: ', childNum
!flush(dbgUnit(5))

      !  do iNeigh = 1, stencil%QQ-1 ! no rest position
      !
      !    !if ( iNeigh == stencil%restPosition ) then
      !      !! for rest position
      !      !supp_bitmask(iNeigh) = .true.
      !    !else
      !      ! Find neighor using neigh array because for boundary, neigh array
      !      ! points to current element so no need to check for existence of the
      !      ! neighbor element
      !      ! check for both first and second layer of ghost cells
      !      ! the second condition is true only if the secon neigh is fluid
      !      nghElem = levelDesc( targetLevel )%neigh(1)%nghElems( iNeigh, targetElem )
      !      if (nghElem > 0) then ! nghElem exists
      !        if ( btest(levelDesc( targetLevel )%property( nghElem ), prp_fluid) ) then
      !          ! nghElem is
      !          !invDir = stencil%cxDirInv( iNeigh )
      !          !supp_bitmask(invDir) = .false.
      !          ! then this is the closest fine ghost cell to the fluid cell
      !          levelDesc( targetLevel )%property( targetElem ) = &
      !            & ibset( levelDesc( targetLevel )%property( targetElem ), prp_fineGhostClosestToFluid )
      !         else
      !          ! find the second layer of closest fine ghosts to fluid
      !          nghElem_2 = levelDesc( targetLevel )%neigh(1)%nghElems( iNeigh, nghElem )
      !          if ( nghElem_2 > 0 ) then
      !            if ( btest(levelDesc( targetLevel )%property( nghElem_2 ), prp_fluid) ) then
      !      !        supp_bitmask = .true.
      !              levelDesc( targetLevel )%property( targetElem ) = &
      !            & ibset( levelDesc( targetLevel )%property( targetElem ), prp_fineGhostClosestToFluid )
      !            end if
      !          end if
      !        end if
      !      end if
      !    !end if
      !  end do ! iNeigh
      !
      !
      !  !if (childNum == 8) then
      !  !  do iChild = 1, suppChild
      !  !    levelDesc( targetLevel )%depFromCoarser( iElem - suppChild + iChild )%bitmask = supp_bitmask
      !  !  end do
      !  !  suppChild = 1
      !  !else
      !  !  suppChild = suppChild + 1
      !  !end if
      !
      !end do

      ! Treat all the fromCoarser elements
      do iElem = 1, nGhostElems

        parentPosInTotal = &
          & levelDesc( targetLevel )%depFromCoarser( iElem )%elem%val(1)
        targetBary = levelDesc( targetLevel )%depFromCoarser( iElem )%coord(:)
        childNum = levelDesc(targetLevel)%depFromCoarser(iElem)%childNum

        ! target element treeID
write(dbgUnit(4),"(A,I0)") ''
        ! target element position in total list
        targetElem = iElem + levelDesc( targetLevel )%             &
          &                                offset( 1, eT_ghostFromCoarser)

write(dbgUnit(4),"(A,I0)") ' target treeID: ', leveldesc( targetlevel )%total( targetElem )
write(dbgUnit(4),"(A,I0)") ' child number: ', childNum
write(dbgUnit(4),"(A,3F8.3)") ' coord relative to parent: ', targetBary(1:3)
write(dbgUnit(4),"(A,I0)") ' number of sources (parent): ', &
  &             levelDesc( targetLevel )%depFromCoarser( iElem )%elem%nVals
write(dbgUnit(4),"(A,I0)") ' parent treeID: ', levelDesc( sourceLevel )%total( parentPosInTotal )
write(dbgUnit(4),"(A   )") ''

        ! --------------- treat parent element (first source) --------------
        ! start with treating the parent element, which is located at the
        ! first position in the dependency list
        ! Later on, proceed to siblings of parent in the target element's
        ! directions from its parent element
        ! add the parent to the list of source elements on this level
        nFoundSources = 0
        mySources = 0
        myNeighDir = -1
        ! Now loop over all neighbors and find maximum number of source elements
        do iNeigh = 1, stencil%QQ

write(dbgUnit(5),"(A, I0)") '      iNeigh: ', iNeigh

          if ( iNeigh == stencil%restPosition ) then
            ! normally parent is the last source due to stencil definition
            posInTotal = parentPosInTotal
          else
            ! this is a neihgbor of the parent element ->
            ! access neighbor information.
            posInTotal = levelDesc( sourceLevel )%neigh(1)%       &
              &                   nghElems( iNeigh, parentPosInTotal )
          end if

write(dbgUnit(5),"(A, I0)") '  posInTotal: ', posInTotal
          if ( posInTotal > 0 ) then
write(dbgUnit(5),"(A, I0)") '  sourceID: ', levelDesc(sourceLevel)%total(posInTotal)
            nFoundSources = nFoundSources + 1
            mySources( nFoundSources ) = posInTotal
            myNeighDir( nFoundSources ) = iNeigh
          end if ! posIntotal > 0, i.e. this source is a valid one
        end do ! iNeigh

!write(dbgUnit(3),*) '  bitmask = ', levelDesc( targetLevel )%depFromCoarser( iElem )%bitmask(:)

        ! if sources are found then determine which interpolation to use
        ! depending on available sources
        if ( nFoundSources /= 0 ) then
write(dbgUnit(4),"(A, I0)")   '      nFound: ', nFoundSources
write(dbgUnit(4),*) '   myNeighDir: ', myNeighDir(1:nFoundSources)
write(dbgUnit(4),*) '    mySources: ', mySources(1:nFoundSources)

          ! Find maximum possible interpolation order for current nFoundSources
          ! and update mySources, neighDir and nFoundSources if all sources
          ! in stencil are found.
          call find_possIntpOrderAndUpdateMySources(          &
            &    intp          = intp,                        &
            &    mySources     = mySources(1:nFoundSources),  &
            &    neighDir      = myNeighDir(1:nFoundSources), &
            &    nFoundSources = nFoundSources,               &
            &    childNum      = childNum,                    &
            &    intpOrder     = intpOrder                    )
write(dbgUnit(4),"(A,I0)") '   intpOrder: ', intpOrder
write(dbgUnit(4),"(A, I0)")   '      nFound_new: ', nFoundSources


          ! Compute interpolation matrix for least square fit using stencil
          ! direction of available sources
          ! The parent of target childs coord is 0,0,0 so we could
          ! just use of stencil%cxDir to build up this matrix entries
          ! Every row in matrix is evaluated with coord of source element
          if (intpOrder == quadratic) then
            call append( me       = intp%fillFinerFromMe(intpOrder) &
              &                          %intpMat_forLSF,           &
              &          order    = intpOrder,                      &
              &          QQ       = stencil%QQ,                     &
              &          nDims    = stencil%nDims,                  &
              &          nSources = nFoundSources,                  &
              &          cxDirRK  = stencil%cxDirRK,                &
              &          neighDir = myNeighDir(1:nFoundSources),    &
              &          pos      = posInMatArray,                  &
              &          success  = success                         )
             if (success) then
               levelDesc( targetLevel )%depFromCoarser(iElem)         &
                 &                     %posInIntpMatLSF = posInMatArray
             else
               ! matrix is singular then fall back to linear
               intpOrder = linear
             end if
          end if

          if (intpOrder == linear) then
            call append( me       = intp%fillFinerFromMe(intpOrder) &
              &                          %intpMat_forLSF,           &
              &          order    = intpOrder,                      &
              &          QQ       = stencil%QQ,                     &
              &          nDims    = stencil%nDims,                  &
              &          nSources = nFoundSources,                  &
              &          cxDirRK  = stencil%cxDirRK,                &
              &          neighDir = myNeighDir(1:nFoundSources),    &
              &          pos      = posInMatArray,                  &
              &          success  = success                         )
             if (success) then
               levelDesc( targetLevel )%depFromCoarser(iElem)         &
                 &                     %posInIntpMatLSF = posInMatArray
             else
               ! matrix is singular then fall back to weighted average
               intpOrder = weighted_average
             end if
          end if

          if (intpOrder == weighted_average) then
            ! Compute weights for weighted average interperation.
            ! Weights are computed according to available sources.
            call compute_weight( &
              &    weights       = levelDesc(targetLevel)         &
              &                    %depFromCoarser(iElem)%weight, &
              &    nFoundSources = nFoundSources,                 &
              &    neighDir      = myNeighDir(1:nFoundSources),   &
              &    targetBary    = targetBary,                    &
              &    intp_config   = intp%config,                   &
              &    cxDirRK       = stencil%cxDirRK                )
          end if

write(dbgUnit(4),"(A,I0)") '   Final intpOrder: ', intpOrder
          ! assign current element to correct interpolation scheme according
          ! to available sources
          call append( me  = levelDesc( targetLevel )     &
            &                %intpFromCoarser(intpOrder), &
            &          val = iElem                        )

        else ! no Sources
          write(logUnit(1),*) 'Can NOT find any source for this ghostFromCoarser.'
          write(logUnit(1),*) 'Thus can NOT do interpolation.'
          call tem_abort()
        end if

        call depSource_append(                                                &
          &    me         = levelDesc( targetLevel )%depFromCoarser( iElem ), &
          &    sourceList = levelDesc( targetLevel )%sourceFromCoarser,       &
          &    mySources  = mySources(1:nFoundSources),                       &
          &    n          = nFoundSources                                     )

      end do ! iElem = 1, nGhostElems

      write(dbgUnit(4),"(A,I0)") ' level: ', targetLevel
      write(dbgUnit(4),"(A,I0)") ' number of ghost elements: ', nGhostElems
      do iOrder = 0, intp%config%order
        write(dbgUnit(4),"(A,I0)") ' ghostFromCoarser: ', &
          &      levelDesc( targetLevel )%intpFromCoarser(iOrder)%nVals
      end do
      write(dbgUnit(4),"(A,I0)") ' total sources: ',      &
        &      levelDesc( targetLevel )%sourceFromCoarser%nVals

    end do ! level loop

    do sourceLevel = minLevel, maxLevel-1
      targetLevel = sourceLevel + 1
      write(dbgUnit(4),*) 'sourceIDs for target level: ', targetLevel
      do iSourceElem = 1, levelDesc(targetLevel)%sourceFromCoarser%nVals
        sourceElem = levelDesc(targetLevel)%sourceFromCoarser%val(iSourceElem)
        write(dbgUnit(4),*) iSourceElem, sourceElem, &
          & levelDesc(sourceLevel)%total(sourceElem)
      end do
    end do

    write(dbgUnit(4),"(A)") "Outside of mus_intp_update_depFromCoarser routine"
    call tem_horizontalSpacer( fUnit = dbgUnit(4), after = 1 )
  end subroutine mus_intp_update_depFromCoarser
! ****************************************************************************** !


! **************************************************************************** !
  !> Find maximum possible interpolation order which can be used for
  !! fillFinerFromMe by comparing nFoundSources with nMaxSources of different
  !! interpolation order starting from interpolation order defined by user
  subroutine find_possIntpOrderAndUpdateMySources( intp, mySources, neighDir, &
    &                                              nFoundSources, childNum,   &
    &                                              intpOrder                  )
    ! ---------------------------------------------------------------------------
    !> interpolation type
    type(mus_interpolation_type), intent(in)  :: intp
    !> Number of source elements found
    integer, intent(inout) :: nFoundSources
    !> position of found source elements in total list
    !! Update this list if intpStencil%isActive
    integer, intent(inout) :: mySources(:)
    !> cxDir for found sounce elements
    integer, intent(inout) :: neighDir(:)
    !> Curent finer ghost child number
    integer, intent(in) :: childNum
    !> interpolation order
    integer, intent(out) :: intpOrder
    ! --------------------------------------------------------------------------
    integer :: iOrder, intpOrder_tmp
    logical :: allSrcFound
    ! --------------------------------------------------------------------------
    ! Start with interpolation order defined by user and if requirement
    ! for that interpolation is not satisfied then fall back to lower
    ! order
    intpOrder = intp%config%order

    ! if interpolation order <= quadratic then
    ! Find which order to use for this ghost element depending on
    ! order defined in config and number of sources found
    ! Set to weighted_average if nSources does not satisfy
    ! for any higher order
    intpOrder_tmp = weighted_average
    do iOrder = intpOrder, linear, -1
      if ( nFoundSources >= intp%fillFinerFromMe(iOrder)%nMinSources ) then
        intpOrder_tmp = iOrder
        exit
      end if
    end do
    intpOrder = intpOrder_tmp

    ! if weightedAvgStencil is defined then use only neighDir from
    ! stencilNeighDir as found sources for weighted_average.
    ! If all sources required by the stencil are not found then use found
    ! sources
    if ( (intpOrder == weighted_average  .or. intpOrder == linear ) &
      & .and. intp%weightedAvgStencil%isActive ) then
      call updateMySources(                                         &
        &    mySources     = mySources,                             &
        &    neighDir      = neighDir,                              &
        &    nFoundSources = nFoundSources,                         &
        &    childNum      = childNum,                              &
        &    nMaxSources   = intp%fillFinerFromMe(weighted_average) &
        &                        %nMaxSources,                      &
        &    intpStencil   = intp%weightedAvgStencil,               &
        &    allSrcFound   = allSrcFound                            )
    end if

  contains

    ! ************************************************************************ !
    !> Update mySources if all sources in required by the stencil are found
    subroutine updateMySources(mySources, neighDir, nFoundSources, childNum, &
      &                        nMaxSources, intpStencil, allSrcFound)
      ! ------------------------------------------------------------------------
      !> Number of source elements found
      integer, intent(inout) :: nFoundSources
      !> position of found source elements in total list
      !! Update this list if intpStencil%isActive
      integer, intent(inout) :: mySources(:)
      !> cxDir for found sounce elements
      integer, intent(inout) :: neighDir(:)
      !> ChildNum of current target fine element
      integer, intent(in) :: childNum
      !> Maximum number of sources required for this intpStencil
      integer, intent(in) :: nMaxSources
      !> Interpolation stencil
      type(mus_interpolation_stencil_type), intent(in) :: intpStencil
      !> is true if all sources required by the stencil are found
      logical, intent(out) :: allSrcFound
      ! ------------------------------------------------------------------------
      integer :: iSrc, nFoundSources_new
      integer :: mySources_tmp(nFoundSources), neighDir_tmp(nFoundSources)
      ! ------------------------------------------------------------------------
      allSrcFound = .false.
      nFoundSources_new = 0

      do iSrc = 1, nFoundSources
        if ( any(intpStencil%neighDir(1:nMaxSources, childNum) &
          & == neighDir(iSrc)) ) then
          nFoundSources_new = nFoundSources_new + 1
          neighDir_tmp(nFoundSources_new) = neighDir(iSrc)
          mySources_tmp(nFoundSources_new) = mySources(iSrc)
        end if
      end do
      ! all sources required by this stencil are found. update sources, neighdir
      ! and nFoundSources
      if (nFoundSources_new == nMaxSources) then
        allSrcFound = .true.
        nFoundSources = nFoundSources_new
        mySources(1:nFoundSources) = mySources_tmp(1:nFoundSources)
        neighDir(1:nFoundSources) = neighDir_tmp(1:nFoundSources)
      end if

    end subroutine updateMySources
    ! ************************************************************************ !

  end subroutine find_possIntpOrderAndUpdateMySources
! **************************************************************************** !

! ****************************************************************************** !
  !> This routine computes weights for weighted_average interpolation.
  subroutine compute_weight( weights, nFoundSources, neighDir, targetBary, &
    &                        intp_config, cxDirRK )
    ! ---------------------------------------------------------------------------
    !> computed weight
    real(kind=rk), allocatable, intent(out) :: weights(:)
    !> Number of source elements found
    integer, intent(in) :: nFoundSources
    !> cxDir for found sounce elements
    integer, intent(in) :: neighDir(nFoundSources)
    !> child bary relative to parent
    real(kind=rk), intent(in) :: targetBary(3)
    !> Interpolation config info
    type(mus_interpolation_config_type), intent(in) :: intp_config
    !> cxDir of current source
    real(kind=rk), intent(in) :: cxDirRK(:,:)
    ! ---------------------------------------------------------------------------
    real(kind=rk) :: scalar_dist
    !> absolute distances
    real(kind=rk) :: distance(3)
    real(kind=rk)  :: sumWeight
    integer :: iSrc, iNeigh
    ! ---------------------------------------------------------------------------
write(dbgUnit(4),"(A)") "Inside compute weight"
    allocate( weights( nFoundSources ) )
    do iSrc = 1, nFoundSources
      iNeigh = neighDir(iSrc)
      ! distance from this source to target
      ! coordinates of source is indeed its offset to parent of target
      ! For 2D: use only X and Y
      distance(1:3) = abs( cxDirRK(1:3, iNeigh) - targetBary )

      scalar_dist = sqrt(dot_product(distance, distance))

      select case ( trim(intp_config%weights_method) )
        case( 'inverse_distance' )
          weights(iSrc) = 1.0_rk / scalar_dist**intp_config%IDW_powerfac
        case( 'modified_shepards' )
          weights(iSrc) = ( max( 0.0_rk, 1.0_rk - scalar_dist ) &
            &           / scalar_dist )**2
        case( 'linear_distance')
          weights(iSrc) = (1.0_rk - distance(1)) * (1.0_rk-distance(2)) &
            &           * (1.0_rk - distance(3))
        case default
          write(*,*) 'Unknown weighting method, check inputs. Aborting ...'
          write(logUnit(0),*) 'weighting_method in lua file should be:'
          write(logUnit(0),*) 'inverse_distance'
          write(logUnit(0),*) 'modified_shepards'
          write(logUnit(0),*) 'Aborting ...'
          call tem_abort()
      end select
write(dbgUnit(4),"(A,3F7.3)") '    distance: ', distance(1:3)
write(dbgUnit(4),"(A,3F7.3)") '      weight: ', weights(iSrc)
    end do !iSrc

    ! Normalize weight, so that sum(weight) = 1
    sumWeight = sum( weights(1:nFoundSources) )
    weights(:) = weights( 1:nFoundSources ) / sumWeight
write(dbgUnit(7),*) '   Weights: ', weights

  end subroutine compute_weight
! ****************************************************************************** !

! ****************************************************************************** !
  !> All sources (children) have been found in treelm,
  !! updated number of sources needed, based on nDims
  !! collect all source elements into sourceFromFiner
  !! assign ghost intp list. Currently only average interpolation is implemented
  !! for fillMineFromFiner
  !!
  subroutine mus_intp_update_depFromFiner( intp, levelDesc, minLevel, maxLevel, &
    & stencil )
    ! ---------------------------------------------------------------------------
    !> interpolation type
    type(mus_interpolation_type), intent(inout)  :: intp
    integer,                        intent(in) :: minLevel
    integer,                        intent(in) :: maxLevel
    !> Level Descriptor
    type( tem_levelDesc_type )                 :: levelDesc(minLevel:maxLevel)
    !> Stencil
    type( tem_stencilHeader_type ), intent(in) :: stencil
    ! ---------------------------------------------------------------------------
    integer :: iLevel, iElem
    integer :: nSources   ! number of source elements actually appended
    integer :: mySources( intp%fillMineFromFiner%nMaxSources )
    integer :: iSourceElem, sourceElem
    !integer :: nghElem, targetElem, iNeigh, invDir
    ! ---------------------------------------------------------------------------

    write(logUnit(1),"(A)") 'Updating dependencies of ghost from finer.'
call tem_horizontalSpacer( fUnit = dbgUnit(4), before = 1 )
    write(dbgUnit(4),"(A)") 'Updating dependencies of ghost from finer.'

    ! all levels except highest level
    do iLevel = minLevel, maxLevel-1

      do iElem = 1, levelDesc( iLevel )%elem%nElems( eT_ghostFromFiner )

        ! Depending of dimension elements are added to source list
        ! nMaxSources in fillMineFromFiner is set according to nDims
        ! for 3D, all 8 elements are considered
        ! for 2D, only first 4 elements are considered
        ! for 1D, only first 2 elements are considered.
        nSources = min( intp%fillMineFromFiner%nMaxSources,                  &
          &             levelDesc( iLevel )%depFromFiner( iElem )%elem%nVals )
        mySources(1:nSources) = levelDesc( iLevel )%depFromFiner( iElem ) &
          &                                        %elem%val(1:nSources)

        call depSource_append(                                       &
          &  me         = levelDesc( iLevel )%depFromFiner( iElem ), &
          &  sourceList = levelDesc( iLevel )%sourceFromFiner,       &
          &  mySources  = mySources(1:nSources),                     &
          &  n          = nSources                                   )

        ! add into interpolation list
        call append( me  = levelDesc( iLevel )%intpFromFiner, &
          &          val = iElem )

        !targetElem = iElem + levelDesc( iLevel )%             &
        !  &                                offset( 1, eT_ghostFromFiner)

        !write(logUnit(1),"(A)") 'Create mask for coarse ghosts.'
        !flush(logunit(1))
        !allocate(levelDesc( iLevel )%depFromFiner( iElem )%bitmask(stencil%QQ))
        !levelDesc( iLevel )%depFromFiner( iElem )%bitmask = .false.

        ! Now loop over all neighbors
        !do iNeigh = 1, stencil%QQ-1 !no rest position
        !  !if ( iNeigh == stencil%restPosition ) then
        !    ! for rest position coarse ghost has the correct value
        !    !levelDesc( iLevel )%depFromFiner( iElem )%bitmask(iNeigh) = .true.
        !  !else if ( .not. iNeigh == stencil%restPosition ) then
        !    ! Find neighor using neigh array because for boundary, neigh array
        !    ! points to current element so no need to check for existence of the
        !    ! neighbor element
        !    ! check for both first and second layer of ghost cells
        !    ! the second condition is true only if the secon neigh is fluid
        !    nghElem = levelDesc( iLevel )%neigh(1)%nghElems( iNeigh, targetElem )
        !    if (nghElem > 0) then
        !      if ( btest(levelDesc( iLevel )%property( nghElem ), prp_fluid) ) then
        !        levelDesc( iLevel )%depFromFiner( iElem )%bitmask(iNeigh) = .true.
        !      end if
        !    end if
        !  !end if
        !end do

!write(dbgUnit(3),*) '  bitmask = ', levelDesc( iLevel )%depFromFiner( iElem )%bitmask(:)


      end do ! elem%nElems( eT_ghostFromFiner )

      write(dbgUnit(4),*) 'SourceIDs on targetlevel:', iLevel
      do iSourceElem = 1, levelDesc(iLevel)%sourceFromFiner%nVals
        sourceElem = levelDesc(iLevel)%sourceFromFiner%val(iSourceElem)
        write(dbgUnit(4),*) iSourceElem, sourceElem, levelDesc(iLevel+1)%total(sourceElem)
      end do

    end do ! iLevel = minLevel, maxLevel-1

    write(dbgUnit(4),"(A)") "Outside of mus_intp_update_depFromFiner routine"
    call tem_horizontalSpacer( fUnit = dbgUnit(4), after = 1 )

  end subroutine mus_intp_update_depFromFiner
! ****************************************************************************** !


end module mus_interpolate_module
! ****************************************************************************** !
!> \page intp_methods Interpolation methods
!!
!! \section intp_average Averaging interpolation
!!\verbatim interpolation_method = 'linear'
!!\endverbatim
!! Simple averaging of all quantities with weighted summation
!! \section intp_linear Linear interpolation
!! Linear interpolation of the first two moments (density and momentum) and
!! separation of equilibrium and non-equilibrium part to account for correct
!! shear stress calculation
!! This fine->coarse interpolation method is based on Dazhi Yu et. al. (2002):
!! "A multi-block lattice Boltzmann method for viscous fluid flows".
!! The basic idea is to split up the source pdfs into equilibrium- and
!! non-equilibrum- parts, then interpolate both parts with average interpolation,
!! and apply a corrective scaling factor to the non-equilibrium part before
!! adding both parts up and writing it to the target element.
!! Values are interpolated in a linear fashion by multiplying the quantities
!! of the source elements by weighting factors
