! 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-2015, 2017-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 Tobias Schneider <tobias1.schneider@student.uni-siegen.de>
! Copyright (c) 2017-2018, 2020 Raphael Haupt <raphael.haupt@uni-siegen.de>
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions are met:
!
! 1. Redistributions of source code must retain the above copyright notice,
! this list of conditions and the following disclaimer.
!
! 2. Redistributions in binary form must reproduce the above copyright notice,
! this list of conditions and the following disclaimer in the documentation
! and/or other materials provided with the distribution.
!
! THIS SOFTWARE IS PROVIDED BY THE UNIVERSITY OF SIEGEN “AS IS” AND ANY EXPRESS
! OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
! OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
! IN NO EVENT SHALL UNIVERSITY OF SIEGEN OR CONTRIBUTORS BE LIABLE FOR ANY
! DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
! (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON
! ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
! SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
! ****************************************************************************** !
!> author: Manuel Hasert
!! author: Jiaxing Qi
!! Interpolation header to load confiugration and type definition
module mus_interpolate_header_module

  ! include aotus modules
  use aotus_module,     only: flu_State, aot_get_val
  use aot_table_module, only: aot_table_open, aot_table_close, aot_table_length
  use aot_out_module,   only: aot_out_type, aot_out_val

  ! include treelm modules
  use env_module,            only: rk, labelLen
  use tem_tools_module,      only: tem_horizontalSpacer
  use tem_aux_module,        only: tem_abort
  use tem_logging_module,    only: logUnit
  use tem_debug_module,      only: dbgUnit
  use tem_varSys_module,     only: tem_varSys_type
  use tem_construction_module, only: tem_levelDesc_type
  use tem_param_module,      only: qN00, q0N0, q00N, q100, q010, q001, &
    &                              q0NN, q0N1, q01N, q011, qN0N, q10N, &
    &                              qN01, q101, qNN0, qN10, q1N0, q110, &
    &                              qNNN, qNN1, qN1N, qN11, q1NN, q1N1, &
    &                              q11N, q111, q000
  use tem_matrix_module,     only: tem_intpMatrixLSF_type, init
  use tem_time_module,       only: tem_time_type
  use tem_stencil_module,    only: tem_stencilHeader_type

  ! include musubi modules
  use mus_field_prop_module,     only: mus_field_prop_type
  use mus_scheme_layout_module,  only: mus_scheme_layout_type
  use mus_physics_module,        only: mus_physics_type
  use mus_derVarPos_module,      only: mus_derVarPos_type
  use mus_directions_module,     only: qN0, q0N, q10, q01, qNN, qN1, q1N, q11

  implicit none

  private

  public :: mus_interpolation_type
  public :: mus_interpolation_config_type
  public :: mus_interpolation_method_type
  public :: mus_interpolation_stencil_type
  public :: mus_load_interpolate
  public :: mus_interpolate_out
  public :: mus_set_nSources

  !> Interpolation parameter to choose fillFinerFromMe
  integer, parameter, public ::          no_intp = -1
  integer, parameter, public :: weighted_average = 0
  integer, parameter, public ::           linear = 1
  integer, parameter, public ::        quadratic = 2
  integer, parameter, public ::          compact = 3

  ! D2Q9 directions are imported from mus_directions_module, we only need to
  ! specify q00
  integer,parameter :: q00 = 9  !< rest

  !> This data types contains intpRoutine function pointer for FillFiner
  !! and FillCoarser.
  !! For fillFiner, it build least square fit matrix for linear
  !! quadratic interpolations
  !! For fillCoarser: currently we do simple average
  !!
  !! Why do we need different intpRoutine for fillFinerFromMe?
  !! The order of interpolation to finer depends on available number
  !! of coarser source elements so for every order we use different
  !! interpolation routines.
  !! We start with user defined interpolation order and
  !! If nMinSources for that order is not found then we fall back to lower order
  !! Weighted average will be the lowest level for which nMinSources = 1
  type mus_interpolation_method_type

    !> Routine to interpolate coarse to fine for ghostFromCoarser elements
    !! and interpolate fine to coarse for ghostFromFiner elements.
    !! Sets pdf for ghost elements by f_eq + f_neq
    !! The moments required to compute equilibrium function are obtained
    !! from auxField array and the auxField of ghost elements are interpolate
    !! seperately using do_intpArbitraryField
    procedure(intpRoutine), pointer :: do_intp => null()

    !> Routine to interpolate coarse to fine and fine to coarse for
    !! arbitrary variables
    procedure(intpRoutine_arbitraryVal), pointer :: do_intpArbiVal => null()


    !> Matrix entries for linear/Quadratic interpolation least square fit
    !! ((A^T)A)^-1*(A^T)
    !! Size: (6,9) for D2Q9 stencil
    !! Size: (10,QQ)  for D3Q19 and D3Q27
    type(tem_intpMatrixLSF_type) :: intpMat_forLSF

    !> how many source elements are required by this interpolation order
    integer :: nMinSources

    !> Max number of sources amoung target ghosts
    !! Computed in mus_contruction::mus_intp_complete_coarseDep
    integer :: nMaxSources
  end type mus_interpolation_method_type


  !> Contains stencil for interpolation
  type mus_interpolation_stencil_type
    !> Is active only for specific layouts like d2q9, d3q19, d3q27
    logical :: isActive = .false.
    !> cxDir for interpolation stencil for depFromCoarser
    integer, allocatable :: neighDir(:,:)
  end type mus_interpolation_stencil_type

  !> Contains information loaded from config file
  type mus_interpolation_config_type
    !> name of the interpolation method for fillFinerFromMe
    character(len=labelLen) :: label

    !> Order of the interpolation for fillFinerFromMe
    integer :: order

    !> name of used weighting method
    character(len=labelLen) :: weights_method = 'linear_distance'

    !> Power factor for inverse distance weighting
    integer :: IDW_powerfac = 6

    !> Stencil for linear interpolation.
    !! By default use stencil from weighted average
    logical :: useComputeStencil = .false.

    !> Interpolation test by comparing against the initial condition
    logical :: testInterpolation = .false.
    logical :: testEachElement   = .false.
    logical :: testFluids        = .false.
    logical :: noIntpFromFiner   = .false.
    logical :: noIntpFromCoarser = .false.

  end type mus_interpolation_config_type

  !> definition of the used interpolation method
  type mus_interpolation_type

    !> Information loaded from config file
    type(mus_interpolation_config_type) :: config

    !> Interpolation routines to fillFiner
    !! Size: interpolation order
    type(mus_interpolation_method_type), allocatable :: fillFinerFromME(:)

    !> Interpolation routines to fillCoarser
    type(mus_interpolation_method_type) :: fillMineFromFiner

    !> stencil for compact interpolation
    type(mus_interpolation_stencil_type) :: compactStencil

    !> stencil for weighted average interpolation
    type(mus_interpolation_stencil_type) :: weightedAvgStencil

  end type mus_interpolation_type

! ****************************************************************************** !
  abstract interface
    !> This is the interface for all interpolation methods that
    !  can be called from outside to set state variable for ghost elements
    subroutine intpRoutine( method,  fieldProp, tLevelDesc, level, sState, &
      & snSize, tState, tnSize, tAuxField, layout, nTargets, targetList,   &
      & physics, time, varSys, derVarPos )
      import :: mus_interpolation_method_type, mus_scheme_layout_type,         &
        &       mus_field_prop_type, tem_levelDesc_type, mus_physics_type, rk, &
        &       tem_varSys_type, tem_time_type, mus_derVarPos_type

      class(mus_interpolation_method_type), intent(in) :: method

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

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

      !> State vector of SOURCE 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(:)

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

      !> 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
      type(mus_physics_type), intent(in) :: physics

      !> time required to compute analytical solution for TGV case
      type(tem_time_type), intent(in) :: time

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

      !> position of all derive variable in varSys
      type(mus_derVarPos_type), intent(in) :: derVarPos(:)
    end subroutine intpRoutine

    !> This is the interface for all interpolation methods that
    !  can be called from outside to interpolate arbitrary variable
    subroutine intpRoutine_arbitraryVal( method, tLevelDesc, level, stencil,   &
      &                                 sVal, snSize, tVal, tnSize, nTargets, &
      &                                 targetList, nScalars    )
      import :: mus_interpolation_method_type, rk, tem_levelDesc_type, &
        &       tem_stencilHeader_type

      class(mus_interpolation_method_type), intent(in) :: method

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

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

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

      !> array of SOURCE elements
      real(kind=rk),    intent(in) :: sVal(:)

      !> size of sVal
      integer,          intent(in) :: snSize

      !> array of TARGET GHOST elements
      real(kind=rk), intent(inout) :: tVal(:)

      !> size of tVal
      integer,          intent(in) :: tnSize

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

      !> Number of scalars to interpolate
      integer, intent(in) :: nScalars
    end subroutine intpRoutine_arbitraryVal

  end interface
! ****************************************************************************** !

  contains

! ****************************************************************************** !
  !> Read in the type of interpolation scheme
  !!
  !!```lua
  !! interpolation_method = 'linear'               -- simple definition
  !! interpolation_method = {method ='debug', value = 1.} -- definition in a table
  !!```
  !!
  subroutine mus_load_interpolate( me, conf, parent )
    ! ---------------------------------------------------------------------------
    !> interpolation type to load info to
    type(mus_interpolation_config_type), intent(out)  :: me
    !> lua state to load from
    type(flu_state)                        :: conf
    !> optional parent table to load from
    integer, optional, intent(in)          :: parent
    ! ---------------------------------------------------------------------------
    character(len=32) :: localKey
    integer :: iError, nEntries, intp_handle
    ! ---------------------------------------------------------------------------

    localKey = 'interpolation_method'

    ! if more than one variable is defined in a table then
    ! ref_value should also be a table with reference value for each variable
    ! defined in the variable table.
    if( present( parent )) then
      call aot_table_open( L=conf, thandle=intp_handle, key=localKey,          &
        &                  parent= parent )
    else
      call aot_table_open( L=conf, thandle=intp_handle, key=localKey)
    endif

    ! not a table, just the name of the interpolation method
    if (intp_handle == 0) then
      call aot_get_val(L=conf, key=localKey,                                   &
        &              val=me%label,                                           &
        &              ErrCode=iError, default='quadratic')
    else
      ! interpolation is defined as a table with optional additional information
      ! Get size of table
      nEntries = aot_table_length(L=conf, thandle=intp_handle)
      call aot_get_val(L=conf, thandle=intp_handle,                            &
        &              val=me%label, ErrCode=iError,                           &
        &              key = 'method', default='quadratic')

      call aot_get_val(L = conf, thandle = intp_handle, &
        &              val     = me%useComputeStencil,  &
        &              ErrCode = iError,                &
        &              key     = 'use_compute_stencil', &
        &              default = .false.                )


      call aot_get_val(L = conf, thandle = intp_handle, &
        &              val     = me%weights_method,     &
        &              ErrCode = iError,                &
        &              key     = 'weights_method',      &
        &              default = 'linear_distance'      )

      call aot_get_val(L = conf, thandle = intp_handle,       &
        &              val     = me%IDW_powerfac,             &
        &              ErrCode = iError,                      &
        &              key     = 'inverse_distance_powerfac', &
        &              default = 6                            )

      call aot_get_val(L=conf, thandle=intp_handle,                            &
        &              val=me%testInterpolation, ErrCode=iError,               &
        &              key = 'test', default=.false.)
      if( me%testInterpolation ) then
        write(logUnit(1),*)' Activated the Interpolation test.'
        call aot_get_val(L=conf, thandle=intp_handle,                          &
          &              val=me%testEachElement, ErrCode=iError,               &
          &              key = 'testEach', default=.false.)
        if( me%testEachElement   ) then
          write(logUnit(1),*)' Testing each element... '
        end if
        call aot_get_val(L=conf, thandle=intp_handle,                          &
          &              val=me%testFluids, ErrCode=iError,                    &
          &              key = 'testFluids', default=.false.)
        if( me%testFluids ) then
          write(logUnit(1),*)' Testing fluid element... '
        end if
      end if
      call aot_get_val(L=conf, thandle=intp_handle,                            &
        &              val=me%noIntpFromFiner, ErrCode=iError,                 &
        &              key = 'noIntpFromFiner', default=.false.)

      call aot_get_val(L=conf, thandle=intp_handle,                            &
        &              val=me%noIntpFromCoarser, ErrCode=iError,               &
        &              key = 'noIntpFromCoarser', default=.false.)

    endif
    call aot_table_close(L=conf, thandle=intp_handle)

    ! set interpolation order for fillFinerFromMe
    select case( trim(me%label) )
    case( 'weighted_average' )
      me%order = weighted_average
    case( 'linear' )
      me%order = linear
    case( 'quadratic' )
      me%order = quadratic
    case( 'compact' )
      me%order = compact
      call tem_abort('Error: Compact interpolation is not supported!')
    case( 'none' )
      me%order = no_intp
    case default
      call tem_abort('From mus_load_interpolate: Unknown interpolation method')
    end select

    call interpolate_dump( me, logUnit(1) )

  end subroutine mus_load_interpolate
! ****************************************************************************** !

! ****************************************************************************** !
  subroutine mus_set_nSources( me, nDims, QQ, layout)
    ! ---------------------------------------------------------------------------
    type(mus_interpolation_type), intent(inout) :: me
    integer,     intent(in)  :: nDims, QQ
    character(len=*), intent(in) :: layout
    ! ---------------------------------------------------------------------------
    integer :: iOrder
    ! ---------------------------------------------------------------------------

    ! Set nSources for fillCoarse
    me%fillMineFromFiner%nMinSources = 1
    select case(nDims)
    case(1)
      me%fillMineFromFiner%nMaxSources = 2
    case(2)
      me%fillMineFromFiner%nMaxSources = 4
    case(3)
      me%fillMineFromFiner%nMaxSources = 8
    end select

    ! Set nSources for fillFiner
    if (allocated(me%fillFinerFromMe)) deallocate(me%fillFinerFromMe)
    allocate(me%fillFinerFromMe(0:me%config%order))
    do iOrder = 0, me%config%order
      select case(iOrder)
      case(no_intp)
        ! do nothing
      case(weighted_average)
        me%fillFinerFromMe(iOrder)%nMinSources = 1
        select case(nDims)
        case(1)
          me%fillFinerFromMe(iOrder)%nMaxSources = 2
        case(2)
          me%fillFinerFromMe(iOrder)%nMaxSources = 4
          if (trim(layout) == 'd2q9') me%weightedAvgStencil%isActive = .true.
        case(3)
          select case(trim(layout))
          case('d3q19')
            me%fillFinerFromMe(iOrder)%nMaxSources = 7
            me%weightedAvgStencil%isActive = .true.
          case('d3q27')
            me%fillFinerFromMe(iOrder)%nMaxSources = 8
            me%weightedAvgStencil%isActive = .true.
          case default
            me%fillFinerFromMe(iOrder)%nMaxSources = QQ
          end select
        end select

        ! Use weighted average stencil only if QQ stencil is D2Q9, D3Q19
        ! or D3Q27
        ! else compute weight from all available sources
        if (me%weightedAvgStencil%isActive) then
          if ( allocated(me%weightedAvgStencil%neighDir) ) &
            & deallocate(me%weightedAvgStencil%neighDir)
          allocate(me%weightedAvgStencil%neighDir(me%fillFinerFromMe(iOrder) &
            &                                        %nMaxSources, 8))
          me%weightedAvgStencil%neighDir                              &
            & = init_cxDirWeightedAvg( QQ, me%fillFinerFromMe(iOrder) &
            &                                %nMaxSources             )
        end if

      case(linear)
        ! initialize least square matrix for linear interpolation
        call init( me     = me%fillFinerFromMe(iOrder)%intpMat_forLSF, &
          &        length = 1,                                         &
          &        nDims  = nDims,                                     &
          &        order  = linear                                     )

        me%fillFinerFromMe(iOrder)%nMinSources = me%fillFinerFromMe(iOrder)  &
          &                                         %intpMat_forLSF%nCoeffs

        ! Number of sources in weighted average stencil is enough for linear
        ! interpolation so use this stencil. Using weighted average stencil
        ! showed better result then using compute stencil
        if (me%config%useComputeStencil) then
          me%fillFinerFromMe(iOrder)%nMaxSources = QQ
          me%weightedAvgStencil%isActive = .false.
        else
          select case(nDims)
          case(1)
            me%fillFinerFromMe(iOrder)%nMaxSources = 2
          case(2)
            me%fillFinerFromMe(iOrder)%nMaxSources = 4
            if (trim(layout) == 'd2q9') me%weightedAvgStencil%isActive = .true.
          case(3)
            select case(trim(layout))
            case('d3q19')
              me%fillFinerFromMe(iOrder)%nMaxSources = 7
              me%weightedAvgStencil%isActive = .true.
            case('d3q27')
              me%fillFinerFromMe(iOrder)%nMaxSources = 8
              me%weightedAvgStencil%isActive = .true.
            case default
              me%fillFinerFromMe(iOrder)%nMaxSources = QQ
            end select
          end select
        end if

      case(quadratic)
        ! initialize least square matrix for quadratic interpolation
        call init( me     = me%fillFinerFromMe(iOrder)%intpMat_forLSF, &
          &        length = 1,                                         &
          &        nDims  = nDims,                                     &
          &        order  = quadratic                                  )

        me%fillFinerFromMe(iOrder)%nMinSources = me%fillFinerFromMe(iOrder)  &
          &                                        %intpMat_forLSF%nCoeffs
        me%fillFinerFromMe(iOrder)%nMaxSources = QQ

      case(compact)
        me%fillFinerFromMe(iOrder)%nMinSources = 4
        me%fillFinerFromMe(iOrder)%nMaxSources = 4
        me%compactStencil%isActive = .true.
        allocate(me%compactStencil%neighDir(me%fillFinerFromMe(iOrder) &
          &                                   %nMaxSources, 8))
        me%compactStencil%neighDir = init_cxDirCompact( QQ )
      case default
        call tem_abort('From mus_set_nSources: Unknown interpolation order')
      end select
    end do

    write(logUnit(1),"(A)") 'Setting for interpolation scheme: ' &
      &                      //trim(me%config%label)
    write(logUnit(3),"(A)") '  Number of sources from coarser: '
    write(logUnit(3),"(A,I0)") '     nMinSources: ', &
      &                        me%fillMineFromFiner%nMinSources
    write(logUnit(3),"(A,I0)") '     nMaxSources: ', &
      &                        me%fillMineFromFiner%nMaxSources
    write(logUnit(3),"(A)") '  Number of sources from finer: '
    do iOrder = 0, me%config%order
      write(logUnit(3),"(A,I0)") '    order: ', iOrder
      write(logUnit(3),"(A,I0)") '      nMinSources: ', &
      &                          me%fillFinerFromMe(iOrder)%nMinSources
      write(logUnit(3),"(A,I0)") '      nMaxSources: ', &
      &                          me%fillFinerFromMe(iOrder)%nMaxSources
    end do
    if (me%weightedAvgStencil%isActive) then
      write(logUnit(3),"(A)") '  use WeightedAvg stencil: T'
    end if
    if (me%compactStencil%isActive) then
      write(logUnit(3),"(A)") '  use compact stencil: T'
    end if

  end subroutine mus_set_nSources
! ****************************************************************************** !


! ****************************************************************************** !
  !> Dump interpolation method to lua
  subroutine mus_interpolate_out( me, conf )
    ! ---------------------------------------------------------------------------
    !> interpolation type to dump info to
    type(mus_interpolation_type), intent(in)  :: me
    !> aotus type handling the output to the file in lua format
    type(aot_out_type), optional, intent(inout) :: conf
    ! ---------------------------------------------------------------------------
    call aot_out_val( put_conf = conf, vname = 'interpolation_method', &
      &               val = trim(me%config%label) )

  end subroutine mus_interpolate_out
! ****************************************************************************** !

! ****************************************************************************** !
  !> Dump interpolation method to logUnit
  subroutine interpolate_dump( me, outUnit )
    ! ---------------------------------------------------------------------------
    !> interpolation type to dump info to
    type(mus_interpolation_config_type), intent(in)  :: me
    !> File unit to write to
    integer, intent(in) :: outUnit
    ! ---------------------------------------------------------------------------
    call tem_horizontalSpacer(fUnit = outUnit)
    write(outUnit,"(A)") 'Interpolation:'
    write(outUnit,"(2A)") '  method: ',trim(me%label)
    write(outUnit,"(A,i2)")  '   order: ', me%order
    write(outUnit,"(2A)") '  weights_method: ',trim(me%weights_method)
    write(outUnit,"(A,i2)") '  inverse distance powefac: ', &
      &                           me%IDW_powerfac
    call tem_horizontalSpacer(fUnit = outUnit)

  end subroutine interpolate_dump
! ****************************************************************************** !

! ****************************************************************************** !
  !> Initialize compact setencil
  function init_cxDirCompact( QQ ) result( me )
    ! ---------------------------------------------------------------------------
    integer, intent(in)  :: QQ
    ! integer, intent(in)  :: cxDir(3,QQ)
    integer              :: me(4,8)
    ! ---------------------------------------------------------------------------
    integer :: iChild
    integer,parameter :: q000_19 = 19  !< rest density is last for d3q19 stencil
    ! ---------------------------------------------------------------------------

    select case (QQ)
    case (9)
      me(:,:) = reshape( [ qNN, q0N, qN0, q00, & ! A
        &                   q0N, q1N, q00, q10, & ! B
        &                   qN0, q00, qN1, q01, & ! C
        &                   q00, q10, q01, q11, & ! D
        &                   qNN, q0N, qN0, q00, & ! A
        &                   q0N, q1N, q00, q10, & ! B
        &                   qN0, q00, qN1, q01, & ! C
        &                   q00, q10, q01, q11  ], [4,8] ) ! D
    case (27)
      me(:,:) = reshape( [ qNNN, q00N, q0N0, qN00, & ! child 1
        &                   q0NN, q10N, q1N0, q000, & ! child 2
        &                   qN0N, q01N, q000, qN10, & ! child 3
        &                   q00N, q11N, q100, q010, & ! child 4
        &                   qNN0, q000, q0N1, qN01, & ! child 5
        &                   q0N0, q100, q1N1, q001, & ! child 6
        &                   qN00, q010, q001, qN11, & ! child 7
        &                   q000, q110, q101, q011 ], [4,8] )
!KM 20180503    case (19)
!KM 20180503      me(:,:) = reshape( [ qNN0, q00N, q0N0, qN00, & ! child 1
!KM 20180503        &                  q0NN, q10N, q1N0, q000_19, & ! child 2
!KM 20180503        &                  qN0N, q01N, q000_19, qN10, & ! child 3
!KM 20180503        &                  q00N, q000_19, q100, q010, & ! child 4
!KM 20180503        &                  qNN0, q000_19, q0N1, qN01, & ! child 5
!KM 20180503        &                  q0N0, q100, q000_19, q001, & ! child 6
!KM 20180503        &                  qN00, q010, q001, q000_19, & ! child 7
!KM 20180503        &                  q000_19, q110, q101, q011 ], [4,8])
    case default
      write(logUnit(1),"(A)") 'Compact interpolation requires D2Q9 or D3Q27!'
      call tem_abort()
    end select

    write(dbgUnit(1),"(A)") ''
    write(dbgUnit(1),"(A)") ' compact stencil dir '
    do iChild = 1, 8
      ! me(1:3,iNeigh,iChild) = cxDir(1:3, dir(iNeigh, iChild) )
      write(dbgUnit(1), "(A,I0,A,4I3)")  'childNum: ', iChild, &
        &                                ', dir: ', me(1:4,iChild)
    end do
    write(dbgUnit(1),"(A)") ''
  end function init_cxDirCompact
! ****************************************************************************** !


! ****************************************************************************** !
  !> Initialize stencil for weighted average interpolation
  function init_cxDirWeightedAvg( QQ, nSources ) result( me )
    ! ---------------------------------------------------------------------------
    integer, intent(in)  :: QQ
    integer, intent(in)  :: nSources
    integer              :: me(nSources,8)
    ! ---------------------------------------------------------------------------
    integer :: iChild
    integer,parameter :: q000_19 = 19  !< rest density is last for d3q19 stencil
    ! ---------------------------------------------------------------------------

    select case (QQ)
    case (9)
      me(:,:) = reshape( [ qNN, q0N, qN0, q00, & ! A
        &                   q0N, q1N, q00, q10, & ! B
        &                   qN0, q00, qN1, q01, & ! C
        &                   q00, q10, q01, q11, & ! D
        &                   qNN, q0N, qN0, q00, & ! A
        &                   q0N, q1N, q00, q10, & ! B
        &                   qN0, q00, qN1, q01, & ! C
        &                   q00, q10, q01, q11  ], [nSources,8] ) ! D
    case (19)
      me(:,:) = reshape( [ q000_19, qN00, q0N0, q00N, qNN0, qN0N, q0NN, & ! child 1
        &                  q000_19, q100, q0N0, q00N, q1N0, q10N, q0NN, & ! child 2
        &                  q000_19, qN00, q010, q00N, qN10, qN0N, q01N, & ! child 3
        &                  q000_19, q100, q010, q00N, q110, q10N, q01N, & ! child 4
        &                  q000_19, qN00, q0N0, q001, qNN0, qN01, q0N1, & ! child 5
        &                  q000_19, q100, q0N0, q001, q1N0, q101, q0N1, & ! child 6
        &                  q000_19, qN00, q010, q001, qN10, qN01, q011, & ! child 7
        &                  q000_19, q100, q010, q001, q110, q101, q011 ], & ! child 8
        &                                                      [nSources,8])
    case (27)
      me(:,:) = reshape( [ q000, qN00, q0N0, q00N, qNN0, qN0N, q0NN, qNNN, & ! child 1
        &                  q000, q100, q0N0, q00N, q1N0, q10N, q0NN, q1NN, & ! child 2
        &                  q000, qN00, q010, q00N, qN10, qN0N, q01N, qN1N, & ! child 3
        &                  q000, q100, q010, q00N, q110, q10N, q01N, q11N, & ! child 4
        &                  q000, qN00, q0N0, q001, qNN0, qN01, q0N1, qNN1, & ! child 5
        &                  q000, q100, q0N0, q001, q1N0, q101, q0N1, q1N1, & ! child 6
        &                  q000, qN00, q010, q001, qN10, qN01, q011, qN11, & ! child 7
        &                  q000, q100, q010, q001, q110, q101, q011, q111 ],& ! child 8
        &                                                      [nSources,8])
    case default
      write(logUnit(1),"(A)") 'Weighted average interpolation requires D2Q9, D3Q19 '&
        &                    // 'or D3Q27!'
      call tem_abort()
    end select

    write(dbgUnit(1),"(A)") ''
    write(dbgUnit(1),"(A)") ' linear stencil dir '
    do iChild = 1, 8
      ! me(1:3,iNeigh,iChild) = cxDir(1:3, dir(iNeigh, iChild) )
      write(dbgUnit(1), "(A,I0,A,4I3)")  'childNum: ', iChild, &
        &                                ', dir: ', me(1:nSources,iChild)
    end do
    write(dbgUnit(1),"(A)") ''
  end function init_cxDirWeightedAvg
! ****************************************************************************** !


end module mus_interpolate_header_module
! ****************************************************************************** !
