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

  ! include treelm modules
  use env_module,              only: rk, long_k, newUnit, pathLen
  use treelmesh_module,        only: treelmesh_type
  use tem_grow_array_module,   only: grw_intArray_type
  use tem_construction_module, only: tem_levelDesc_type
  use tem_element_module,      only: eT_GhostFromCoarser, &
    &                                eT_ghostFromFiner
  use tem_geometry_module,     only: tem_baryOfId, &
    &                                tem_elemSize
  use tem_param_module,        only: cs2inv, cs2, PI, div1_2, div1_9, div4_9,  &
    &                                div1_36, div2_3, div1_18, div2_9, div1_18,&
    &                                rho0
  use tem_matrix_module,       only: tem_matrix_type
  ! use tem_varSys_module,       only: tem_varSys_type

  ! include musubi modules
  use mus_interpolate_header_module, only: mus_interpolation_type
  use mus_scheme_layout_module,      only: mus_scheme_layout_type
  use mus_pdf_module,                only: pdf_data_type
  use mus_scheme_header_module,      only: mus_scheme_header_type

  implicit none

  private

  public :: dump_intpLists
  public :: debug_dependencies
  public :: fillMomBuf
  public :: get_fNeqFac_f2c
  public :: get_fNeqFac_c2f

  public :: mus_intp_getMoments
  public :: mus_intp_getSrcMoments
  public :: mus_intp_convertMomToPDF3D
  public :: mus_intp_convertMomToPDF3D_incomp
  public :: mus_intp_convertMomToPDF2D
  public :: mus_intp_convertMomToPDF2D_incomp
  public :: mus_intp_scaleMoments2D
  public :: mus_intp_scaleMoments2D_incomp
  public :: mus_intp_scaleMoments3D
  public :: mus_intp_scaleMoments3D_incomp

  contains
! ****************************************************************************** !
  !> check the dependencies from Finer
  !!
  subroutine debug_dependencies( intp, leveldesc, tree, rank )
    ! ---------------------------------------------------------------------------
    !> state properties
    type( tem_leveldesc_type ), intent(in) :: levelDesc(:)
    !> interpolation method info
    type( mus_interpolation_type ), intent(in) :: intp
    !> global tree information
    type( treelmesh_type ), intent(in) :: tree
    !> musubi mpi communicator environment
    integer, intent(in) :: rank
    ! ---------------------------------------------------------------------------
    integer :: nUnit
    character(len=17)  :: buffer
    ! ---------------------------------------------------------------------------
    nUnit = newunit()
    write(buffer,'(a,i7.7)') 'debugDep. rank: ', rank
    open(file = trim(buffer), unit= nUnit, recl=1024)

    write(nUnit,*)
    write(nUnit,*) 'GhostFromFiner'
    write(nUnit,*)
    write(nUnit,*) '-----------------------------------------------------------'
    call dump_MyGhostsFromFiner( intp, levelDesc, nUnit, tree )
    write(nUnit,*) '-----------------------------------------------------------'
    write(nUnit,*) '-----------------------------------------------------------'
    write(nUnit,*)
    write(nUnit,*) 'GhostFromCoarser'
    write(nUnit,*)
    write(nUnit,*) '-----------------------------------------------------------'
    call dump_FinerGhostsFromMe( intp, levelDesc, nUnit, tree )
    write(nUnit,*) '-----------------------------------------------------------'
    write(nUnit,*) '-----------------------------------------------------------'
    write(nUnit,*)
    write(nUnit,*) 'GhostFromCoarserBuffer'
    write(nUnit,*)
    write(nUnit,*) '-----------------------------------------------------------'
    call dump_FinerGhostsFromMeBuffer( intp, leveldesc, nUnit, tree )
    write(nUnit,*) '-----------------------------------------------------------'

    close( nUnit )

  end subroutine debug_dependencies
! ****************************************************************************** !


! ****************************************************************************** !
  !> check the dependencies from Finer
  !!
  subroutine dump_intpLists( minLevel, maxLevel, order, levelDesc, rank )
    ! ---------------------------------------------------------------------------
    !> global pdf information
    integer, intent(in) :: minLevel, maxLevel
    !> state properties
    type( tem_leveldesc_type ), intent(in) :: levelDesc(minLevel:maxLevel)
    integer, intent(in) :: order
    !> musubi mpi communicator environment
    integer, intent(in) :: rank
    ! ---------------------------------------------------------------------------
    integer :: nUnit, iLevel, iOrder
    character(len=17)  :: buffer
    ! ---------------------------------------------------------------------------
    nUnit = newunit()
    write(buffer,'(a4,i7.7)') 'dumpIntpLists.', rank
    open(file = trim(buffer), unit= nUnit, recl=1024)

    do iLevel = minLevel, maxLevel
      write(nUnit,*)
      write(nUnit,*) 'GhostFromFiner      '
      write(nUnit,*)
      write(nUnit,*) '---------------------------------------------------------'
      call dump_intpList( eType = eT_ghostFromFiner, &
        &    levelDesc = levelDesc(iLevel), &
        &    ind = levelDesc( iLevel)%intpFromFiner, nUnit = nUnit )
      write(nUnit,*) '---------------------------------------------------------'
      write(nUnit,*) '---------------------------------------------------------'
      write(nUnit,*)
      write(nUnit,*) 'GhostFromCoarser    '
      write(nUnit,*)
      do iOrder = 0, order
        write(nUnit,*) '--------------------------------------------------------'
        write(nUnit,*) 'order: ',iOrder, '--------------------------------------'
        write(nUnit,*) '--------------------------------------------------------'
        call dump_intpList( eType = eT_ghostFromCoarser, &
          &    levelDesc = levelDesc(iLevel), &
          &    ind = levelDesc(iLevel)%intpFromCoarser(iOrder), nUnit = nUnit )
        write(nUnit,*) '--------------------------------------------------------'
      end do
    end do
    close( nUnit )

  end subroutine dump_intpLists
! ****************************************************************************** !


! ****************************************************************************** !
  !> check the dependencies from Finer and write them out so we can compare
  !!
  subroutine dump_intpList( eType, levelDesc, ind, nUnit )
    ! ---------------------------------------------------------------------------
    !> state properties
    type( tem_levelDesc_type ), intent(in) :: levelDesc
    !> indirectio list
    type( grw_intArray_type ), intent(in)  :: ind
    !>
    integer, intent(in) :: nUnit
    !>
    integer, intent(in) :: eType
    ! ---------------------------------------------------------------------------
    integer :: iElem          ! current target element (for outer loop)
    integer :: indElem        ! element counter for indirection list
    integer :: targetElem     ! element counter for indirection list
    integer(kind=long_k) :: tID
    character(len=pathLen) :: buffer
    ! ---------------------------------------------------------------------------

    buffer = ''
    do indElem = 1, ind%nVals
      iElem = ind%val( indElem )
      targetElem = iElem + levelDesc%offset( 1, eType )
      tID = levelDesc%total( targetElem )
      write(buffer, '(2a, i9)') trim(buffer), ' ', tID
      if( mod( indElem, 32 ) == 0 .or. indElem == ind%nVals ) then
        write(nUnit, *) trim(buffer)
        buffer = ''
      end if
    enddo

  end subroutine dump_intpList
! ****************************************************************************** !


! ****************************************************************************** !
  !> check the dependencies from Finer and write them out so we can compare
  !!
  subroutine dump_MyGhostsFromFiner( intp, levelDesc, nUnit, tree)
    ! ---------------------------------------------------------------------------
    !> state properties
    type( tem_leveldesc_type ), intent(in) :: leveldesc(:)
    !> interpolation method info
    type( mus_interpolation_type ), intent(in) :: intp
    !> global tree information
    type( treelmesh_type ), intent(in) :: tree
    !> unit to write to
    integer, intent(in) :: nUnit
    ! --------------------------------------------------------------------------
    integer :: ilevel         ! grid refinement level
    integer :: sourceLevel    ! level of source elements
    integer :: sourceElem     ! treeId of current source element
    integer :: targetLevel    ! level of target elements
    integer :: targetElem     ! treeId of current source element
    integer :: iElem          ! current target element (for outer loop)
    integer :: iSourceElem    ! current source element (for inner loop)
    integer :: nSourceElems   ! number of source elements for the current target
    integer(kind=long_k) :: tID(0:intp%fillMineFromFiner%nMaxSources )
    ! --------------------------------------------------------------------------

    do iLevel = tree%global%minLevel, tree%global%maxLevel - 1

      sourceLevel = iLevel + 1
      targetLevel = iLevel

      ! Treat all coarse target elements
      do iElem = 1, levelDesc( targetLevel )%elem%nElems( eT_ghostFromFiner )
        ! Read the target element treeId
        targetElem = iElem + levelDesc( targetlevel )%offset( 1, eT_ghostFromFiner )
        ! Find out how many fine source elements we have for interpolation.
        ! Usually 8, but can differ at corners, obstacles, boundaries...
        nSourceELems = levelDesc( ilevel )%depFromFiner( iElem )%elem%nVals
        tID = 0
        tID( 0 ) = levelDesc( targetLevel )%total( targetElem )

        write(nUnit, '(a,i9,3f10.5)') ' targetGhost  ', tID( 0 ),              &
          &                           tem_baryOfId( tree, tID(0) )

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

          ! Get the source element's treeId
          sourceElem = levelDesc( targetLevel )%depFromFiner( iElem )     &
            &                                       %elem%val( iSourceElem )
          tID( iSourceElem ) = levelDesc( sourceLevel )%total( sourceElem )

          call dump_elemDep( targetElem = tID(0),                              &
            &                sourceElem = tID( iSourceElem ),                  &
            &                nUnit = nUnit,                                    &
            &                tree = tree)
        end do  ! iSourceElem
      enddo
    enddo

  end subroutine dump_MyGhostsFromFiner
! ****************************************************************************** !


! ****************************************************************************** !
  !> check the dependencies from Coarser
  !!
  subroutine dump_FinerGhostsFromMe( intp, leveldesc, nUnit, tree)
    ! ---------------------------------------------------------------------------
    !> state properties
    type( tem_levelDesc_type ), intent(in) :: levelDesc(:)
    !> interpolation method info
    type( mus_interpolation_type ), intent(in) :: intp
    !> global tree information
    type( treelmesh_type ), intent(in) :: tree
    !> unit to write to
    integer, intent(in) :: nUnit
    ! ---------------------------------------------------------------------------
    integer :: ilevel         ! grid refinement level
    integer :: sourceLevel    ! level of source elements
    integer :: sourceElem     ! treeId of current source element
    integer :: targetLevel    ! level of target elements
    integer :: targetElem     ! treeId of current source element
    integer :: iElem          ! current target element (for outer loop)
    integer :: iSourceElem    ! current source element (for inner loop)
    integer :: nSourceElems   ! number of source elements for the current target
    integer(kind=long_k), allocatable :: tID(:)
    logical :: weights
    integer :: nMaxSources
    ! ---------------------------------------------------------------------------
    nMaxSources = maxval(intp%fillFinerFromMe(:)%nMaxSources)
    allocate(tID(0:nMaxSources))

    do iLevel = tree%global%minLevel, tree%global%maxLevel - 1

      sourceLevel = ilevel
      targetLevel = ilevel + 1

      ! Treat all coarse target elements
      do iElem = 1, levelDesc( targetLevel )%elem%nElems( eT_ghostFromCoarser )
        ! Read the target element treeId
        targetElem = iElem + levelDesc( targetlevel )%offset( 1, eT_ghostFromCoarser)
        ! Find out how many fine source elements we have for interpolation.
        ! Usually 8, but can differ at corners, obstacles, boundaries...
        nSourceELems = levelDesc( targetlevel )%depFromCoarser( iElem )%  &
          &                                                          elem%nVals
        tID = 0
        tID( 0 ) = levelDesc( targetLevel )%total( targetElem )
        write(nUnit, '(a,i9,3f10.5)') ' targetGhost  ', tID( 0 ),              &
          &                                        tem_baryOfId( tree, tID(0) )

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

          ! Get the source element's treeId
          sourceElem = levelDesc( targetLevel )%depFromCoarser( iElem )   &
            &                                       %elem%val( iSourceElem )
          tID( iSourceElem ) = levelDesc( sourceLevel )%total( sourceElem )

          if( allocated( levelDesc( targetLevel )                         &
              &       %depFromCoarser( iElem )%weight)) then
            weights = .true.
          else
            weights = .false.
          end if
          if( weights ) then
            call dump_elemDep( targetElem = tID(0),                            &
              &                sourceElem = tID( iSourceElem ),                &
              &                nUnit = nUnit,                                  &
              &                tree = tree,                                    &
              &                weight = levelDesc( targetLevel )          &
              &                   %depFromCoarser( iElem )%weight( iSourceElem))
          else
            call dump_elemDep( targetElem = tID(0),                            &
              &                sourceElem = tID( iSourceElem ),                &
              &                nUnit = nUnit,                                  &
              &                tree = tree )
          end if

        end do ! iSourceElem
      end do ! iElem
    end do ! iLevel

  end subroutine dump_FinerGhostsFromMe
! ****************************************************************************** !


! ****************************************************************************** !
  !> check the dependencies from Coarser
  !!
  subroutine dump_FinerGhostsFromMeBuffer( intp, levelDesc, nUnit, tree )
    ! ---------------------------------------------------------------------------
    !> state properties
    type( tem_leveldesc_type ), intent(in) :: leveldesc(:)
    !> interpolation method info
    type( mus_interpolation_type ), intent(in) :: intp
    !> global tree information
    type( treelmesh_type ), intent(in) :: tree
    !> unit to write to
    integer, intent(in) :: nUnit
    ! ---------------------------------------------------------------------------
    integer :: ilevel         ! grid refinement level
    integer :: sourceLevel    ! level of source elements
    integer :: sourceElem     ! treeId of current source element
    integer :: targetLevel    ! level of target elements
    integer :: targetElem     ! treeId of current source element
    integer :: iElem          ! current target element (for outer loop)
    integer :: iSourceElem    ! current source element (for inner loop)
    integer :: nSourceElems   ! number of source elements for the current target
    integer :: totalPos ! Position of the element in the total list
    integer(kind=long_k), allocatable :: tID(:)
    logical :: weights
    integer :: nMaxSources
    ! ---------------------------------------------------------------------------
    nMaxSources = maxval(intp%fillFinerFromMe(:)%nMaxSources)
    allocate(tID(0:nMaxSources))

    do iLevel = tree%global%minLevel, tree%global%maxLevel - 1
      sourceLevel = ilevel
      targetLevel = ilevel + 1

      ! Treat all coarse target elements
      do iElem = 1, levelDesc( targetLevel )%elem%nElems( eT_ghostFromCoarser )
        ! Read the target element treeId
        targetElem = iElem + levelDesc( targetlevel )                     &
          &                      %offset( 1, eT_ghostFromCoarser)
        ! Find out how many fine source elements we have for interpolation.
        ! Usually 8, but can differ at corners, obstacles, boundaries...
        nSourceELems = levelDesc( targetlevel )%depFromCoarser( iElem )%  &
          &                                                          elem%nVals
        tID = 0
        tID( 0 ) = levelDesc( targetLevel )%total( targetElem )
        write(nUnit, '(a,i9,3f10.5)') ' targetGhost  ', tID( 0 ),              &
          &                           tem_baryOfId( tree, tID(0) )

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

          ! Get the source element's treeId
          sourceElem = levelDesc( targetLevel )%depFromCoarser( iElem )   &
            &                                    %elemBuffer%val( iSourceElem )
          totalPos = levelDesc( targetLevel )%sourceFromCoarser%          &
            &                                                 val( sourceElem )
          tID( iSourceElem ) = levelDesc( sourceLevel )%total( totalPos )

          if( allocated( levelDesc( targetLevel )                         &
              &       %depFromCoarser( iElem )%weight)) then
              weights = .true.
          else
            weights = .false.
          end if
          if( weights ) then
            call dump_elemDep( targetElem = tID(0),                            &
              &                sourceElem = tID( iSourceElem ),                &
              &                nUnit = nUnit,                                  &
              &                tree = tree,                                    &
              &                weight = levelDesc( targetLevel )          &
              &                   %depFromCoarser( iElem )%weight( iSourceElem))
          else
            call dump_elemDep( targetElem = tID(0),                            &
              &                sourceElem = tID( iSourceElem ),                &
              &                nUnit = nUnit,                                  &
              &                tree = tree )
          end if

        end do  ! iSourceElem
      enddo
    enddo

  end subroutine dump_FinerGhostsFromMeBuffer
! ****************************************************************************** !


! ****************************************************************************** !
  !> dump dependencies for one element
  !!
  subroutine dump_elemDep( targetElem, sourceElem, nUnit, tree, weight )
    ! ---------------------------------------------------------------------------
    !>
    integer(kind=long_k), intent(in) :: targetElem
    !>
    integer(kind=long_k), intent(in) :: sourceElem
    !>
    integer, intent(in) :: nUnit
    !>
    type( treelmesh_type ), intent(in) :: tree
    !>
    real(kind=rk), optional :: weight
    ! ---------------------------------------------------------------------------
    real(kind=rk) :: xTarget(3), xSource(3)
    real(kind=rk) :: elemSize
    character(len=pathLen) :: buffer
    ! ---------------------------------------------------------------------------
    buffer = ''
    elemSize = tem_elemSize( tree, targetElem )
    xTarget = tem_baryOfId( tree, targetElem )
    xSource = tem_baryOfId( tree, sourceElem )
    write(buffer,'(a,i9, 6f10.5)') '  ', sourceElem,                           &
      &          ((xSource(1:3) - xTarget(1:3))/elemSize), xSource(1:3)
    if( present( weight )) then
      write(buffer, '(a,a,f8.5)') trim( buffer), '    weight: ', weight
    end if
    write(nUnit, '(a)') trim(buffer)

  end subroutine dump_elemDep
! ****************************************************************************** !


! ****************************************************************************** !
  !> Convert PDF of a list of elements to moment, save to buf
  !! @todo: use IDX instead of SAVE macro?
  pure function fillMomBuf( elemList, nElems, state, QQ,  nSize, A) result ( buf )
    ! ---------------------------------------------------------------------------
    integer,       intent(in) :: nElems !< number of source elements
    integer,       intent(in) :: elemList( nElems ) !< elements list
    real(kind=rk), intent(in) :: state(:) !< state array
    integer,       intent(in) :: QQ !< nComponents
    integer,       intent(in) :: nSize
    real(kind=rk), intent(in) :: A(QQ,QQ) !< moments matrix
    real(kind=rk) :: buf(QQ,nElems)
    ! ---------------------------------------------------------------------------
    integer :: iDir, iElem, elemPos
    real(kind=rk) :: PDF(QQ)
    ! ---------------------------------------------------------------------------

    do iElem = 1, nElems
      ! Get the source element's position in the state vector
      elemPos = elemList( iElem )

      ! transfer PDF
      do iDir = 1, QQ
        PDF( iDir ) = state( ?idx?( iDir,elemPos,QQ,nSize))
      end do ! iDir

      ! convert to moment
      buf(:,iElem) = matmul( A, PDF )

    end do ! iElem

  end function fillMomBuf
! ****************************************************************************** !

! ****************************************************************************** !
  !> Calculate fNeq scale factor from fine level to(2) coarse level.
  !! Such scale make sure continuuity of strain rate and shear stress over
  !! levels.
  !! To achive equal strain rate, we want:
  !!   fNeq_c * omega_c / dt_c = fNeq_f * omega_f / dt_f
  !! rearrange this, we get:
  !!   fNeq_c = fNeq_f * ( omega_f / omega_c ) * ( dt_c / dt_f )
  !! dt_c/dt_f is strain rate scale factor (i.e. sFac)
  !!
  !! This function can be used converting fNeq_f to fNeq_c, such as fill coarse
  !! from finer interpolation routines.
  !! How to use this function:
  !!  sFac   = physics%sfac( fineLevel, coarseLevel )
  !!  fNeq_c = fNeq_f * fNeqFac( omega_f, omega_c, sFac )
  !!
  pure function get_fNeqFac_f2c( omegaF, omegaC, sFac ) result ( fac )
    ! ---------------------------------------------------------------------------
    !> omega on fine and coarse level
    real(kind=rk), intent(in) :: omegaF, omegaC
    !> strain rate scaling factor that from fine to coarse
    real(kind=rk), intent(in) :: sFac
    ! ---------------------------------------------------------------------------
    real(kind=rk) :: fac
    ! ---------------------------------------------------------------------------

    fac = omegaF / omegaC * sFac

  end function get_fNeqFac_f2c
! ****************************************************************************** !

! ****************************************************************************** !
  !> Calculate fNeq scale factor from coarse level to(2) fine level.
  !! Such scale make sure continuuity of strain rate and shear stress over
  !! levels.
  !! To achive equal strain rate, we want:
  !!   fNeq_c * omega_c / dt_c = fNeq_f * omega_f / dt_f
  !! rearrange this, we get:
  !!   fNeq_f = fNeq_c * ( omega_c / omega_f ) * ( dt_f / dt_c )
  !! dt_f/dt_c is strain rate from coarse to fine scale factor (i.e. sFac)
  !!
  !! This function can be used converting fNeq_c to fNeq_f,
  !! such as fill fine from coarse interpolation routines.
  !! How to use this function:
  !!  sFac   = physics%sfac( coarseLevel, fineLevel )
  !!  fNeq_f = fNeq_c * fNeqFac( omega_c, omega_f, sFac )
  !!
  pure function get_fNeqFac_c2f( omegaC, omegaF, sFac ) result ( fac )
    ! ---------------------------------------------------------------------------
    !> omega on fine and coarse level
    real(kind=rk), intent(in) :: omegaF, omegaC
    !> strain rate scaling factor that from fine to coarse
    real(kind=rk), intent(in) :: sFac
    ! ---------------------------------------------------------------------------
    real(kind=rk) :: fac
    ! ---------------------------------------------------------------------------

    fac = omegaC / omegaF * sFac

  end function get_fNeqFac_c2f
! ****************************************************************************** !

! **************************************************************************** !
  !> This routine computes moments for all sources elements in momBuf
  !! which can be accessed using depFromCoarser(iElem)%elemBuffer.
  !! In FillFiner, fine elements can have same sources so to avoid calculating
  !! moments for same source elements multiple times, momBuf is used.
  subroutine mus_intp_getSrcMoments(momBuf, state, nSourceElems, srcElemList, &
    &                               QQ, nScalars, nSize, momTransMat, nMoms)
    ! --------------------------------------------------------------------------
    !> Source moments
    real(kind=rk), intent(out) :: momBuf(:,:)
    !> State vector of source elemnents
    real(kind=rk), intent(in) :: state(:)
    !> Total number of source elements
    integer, intent(in) :: nSourceElems
    !> element position in source state array
    integer, intent(in) :: srcElemList(1:nSourceElems)
    !> number of direction in stencil
    integer, intent(in) :: QQ
    !> number of scalars in state vector
    integer, intent(in) :: nScalars
    !> Size of state vector
    integer, intent(in) :: nSize
    !> Moments transformation matrix
    type(tem_matrix_type), intent(in) :: momTransMat
    !> Number of moments to return
    integer, intent(in) :: nMoms
    ! --------------------------------------------------------------------------
    integer :: iSourceElem, iDir, sourceElem
    real(kind=rk) :: fTmp(QQ)
    ! --------------------------------------------------------------------------
    do iSourceElem = 1, nSourceElems
      ! Get the source element's position in the state vector
      sourceElem = srcElemList(iSourceElem)

      do iDir = 1, QQ
        fTmp(iDir) = state(?IDX?(iDir, sourceElem, nScalars, nSize))
      end do

      ! Calculate the raw moments form the pdf
      momBuf(1:nMoms,iSourceElem) = matmul( momTransMat%A(1:nMoms,:), fTmp )

      ! Remove reference density before interpolation and add it after
      ! interpolation and scaling to target element
      momBuf(1,iSourceElem) = momBuf(1,iSourceElem) - rho0
    end do

  end subroutine mus_intp_getSrcMoments
! **************************************************************************** !


! **************************************************************************** !
  !> This function returns macroscopic moments from state
  pure function mus_intp_getMoments(state, elem, QQ, nScalars, nSize, &
    &                               toMoments) result(moments)
    ! --------------------------------------------------------------------------
    !> State vector
    real(kind=rk), intent(in) :: state(:)
    !> element position in state array
    integer, intent(in) :: elem
    !> number of direction in stencil
    integer, intent(in) :: QQ
    !> number of scalars in state vector
    integer, intent(in) :: nScalars
    !> Size of state vector
    integer, intent(in) :: nSize
    !> Moments transformation matrix
    type(tem_matrix_type), intent(in) :: toMoments
    !> moments
    real(kind=rk) :: moments(QQ)
    ! --------------------------------------------------------------------------
    integer :: iDir
    real(kind=rk) :: fTmp(QQ)
    ! --------------------------------------------------------------------------
    do iDir = 1, QQ
      fTmp(iDir) = state(?IDX?(iDir, elem, nScalars, nSize))
    end do

    ! Calculate the raw moments form the pdf
    moments = matmul( toMoments%A, fTmp )

  end function mus_intp_getMoments
! **************************************************************************** !


! **************************************************************************** !
  !> This function scales interpolated moments like density, momentum and
  !! shear stress to target level for 3D compressible model
  pure function mus_intp_scaleMoments3D(momIn, pFac, vFac, nEqFac) &
    & result(momOut)
    ! --------------------------------------------------------------------------
    !> Interpolated moments, dens, m_1, m_2, m_3, SS_XX, SS_YY, S_ZZ,
    !! SS_XY, SS_YZ, SS_ZX
    real(kind=rk), intent(in) :: momIn(10)
    !> Scaling factor for Pressure
    real(kind=rk), intent(in) :: pFac
    !> Scaling factor for velocity
    real(kind=rk), intent(in) :: vFac
    !> Scaling factor for nonEq Shear stress momIn
    real(kind=rk), intent(in) :: nEqFac
    !> Scaled moments, dens, m_1, m_2, m_3, SS_XX, SS_YY, S_ZZ,
    !! SS_XY, SS_YZ, SS_ZX
    real(kind=rk) :: momOut(10)
    ! --------------------------------------------------------------------------
    real(kind=rk) :: tDens, inv_tDens
    ! --------------------------------------------------------------------------
    ! total density
    tDens = momIn(1) + rho0
    inv_tDens = 1.0_rk / tDens
    ! remove pi(0) parts from the second order momInent
    ! that is ((rho u_\alpha * rho u_\beta) / rho)
    ! and additionally rho*cs2 from diagonal entries
    ! Pab = sum( c_ia*c_ib *f_i ) - kronecker_ab*(rho*cs2) - rho*u_a*u_b
    ! KM: its multiplied by inv_tDens because momIn(2,3) are momentum not
    ! velocity
    momOut( 5) = momIn( 5) - momIn(2)*momIn(2)*inv_tDens - cs2*tDens !Pxx
    momOut( 6) = momIn( 6) - momIn(3)*momIn(3)*inv_tDens - cs2*tDens !Pyy
    momOut( 7) = momIn( 7) - momIn(4)*momIn(4)*inv_tDens - cs2*tDens !Pzz
    momOut( 8) = momIn( 8) - momIn(2)*momIn(3)*inv_tDens             !Pxy
    momOut( 9) = momIn( 9) - momIn(3)*momIn(4)*inv_tDens             !Pyz
    momOut(10) = momIn(10) - momIn(4)*momIn(2)*inv_tDens             !Pxz

    ! Rescale density and momentum
    momOut(1) = momIn(1) * pFac + rho0
    momOut(2:4) = momIn(2:4) * vFac

    inv_tDens = 1.0_rk / momOut(1)
    ! Rescale the nonEq moments
    ! and add back the scaled pi(0) part to the tensor
    momOut( 5) = momOut( 5)*nEqFac + momOut(2)*momOut(2)*inv_tDens &
      &         + cs2*momOut(1)
    momOut( 6) = momOut( 6)*nEqFac + momOut(3)*momOut(3)*inv_tDens &
      &         + cs2*momOut(1)
    momOut( 7) = momOut( 7)*nEqFac + momOut(4)*momOut(4)*inv_tDens &
      &         + cs2*momOut(1)
    momOut( 8) = momOut( 8)*nEqFac + momOut(2)*momOut(3)*inv_tDens
    momOut( 9) = momOut( 9)*nEqFac + momOut(3)*momOut(4)*inv_tDens
    momOut(10) = momOut(10)*nEqFac + momOut(4)*momOut(2)*inv_tDens

  end function mus_intp_scaleMoments3D
! **************************************************************************** !

! **************************************************************************** !
  !> This function scales interpolated moments like density, velocity and
  !! shear stress to target level for 3D incompressible model
  pure function mus_intp_scaleMoments3D_incomp(momIn, pFac, vFac, nEqFac) &
    & result(momOut)
    ! --------------------------------------------------------------------------
    !> Interpolated moments, dens, m_1, m_2, m_3, SS_XX, SS_YY, S_ZZ,
    !! SS_XY, SS_YZ, SS_ZX
    real(kind=rk), intent(in) :: momIn(10)
    !> Scaling factor for Pressure
    real(kind=rk), intent(in) :: pFac
    !> Scaling factor for velocity
    real(kind=rk), intent(in) :: vFac
    !> Scaling factor for nonEq Shear stress moments
    real(kind=rk), intent(in) :: nEqFac
    !> Scaled moments, dens, m_1, m_2, m_3, SS_XX, SS_YY, S_ZZ,
    !! SS_XY, SS_YZ, SS_ZX
    real(kind=rk) :: momOut(10)
    ! --------------------------------------------------------------------------
    real(kind=rk) :: tDens, inv_rho
    ! --------------------------------------------------------------------------
    ! total density
    tDens = momIn(1) + rho0
    inv_rho = 1.0_rk / rho0

    ! remove pi(0) parts from the second order momentsent
    ! that is ((rho u_\alpha * rho u_\beta) / rho)
    ! and additionally rho*cs2 from diagonal entries
    ! Pab = sum( c_ia*c_ib *f_i ) - kronecker_ab*(rho*cs2) - rho*u_a*u_b
    momOut( 5) = momIn( 5) - momIn(2)*momIn(2)*inv_rho - cs2*tDens !Pxx
    momOut( 6) = momIn( 6) - momIn(3)*momIn(3)*inv_rho - cs2*tDens !Pyy
    momOut( 7) = momIn( 7) - momIn(4)*momIn(4)*inv_rho - cs2*tDens !Pzz
    momOut( 8) = momIn( 8) - momIn(2)*momIn(3)*inv_rho             !Pxy
    momOut( 9) = momIn( 9) - momIn(3)*momIn(4)*inv_rho             !Pyz
    momOut(10) = momIn(10) - momIn(4)*momIn(2)*inv_rho             !Pxz

    ! Rescale density and momentum
    momOut(1) = momIn(1) * pFac + rho0
    momOut(2:4) = momIn(2:4) * vFac

    ! Rescale the nonEq momOut
    ! and add back the scaled pi(0) part to the tensor
    momOut( 5) = momOut( 5)*nEqFac + momOut(2)*momOut(2)*inv_rho &
      &         + cs2*momOut(1)
    momOut( 6) = momOut( 6)*nEqFac + momOut(3)*momOut(3)*inv_rho &
      &         + cs2*momOut(1)
    momOut( 7) = momOut( 7)*nEqFac + momOut(4)*momOut(4)*inv_rho &
      &         + cs2*momOut(1)
    momOut( 8) = momOut( 8)*nEqFac + momOut(2)*momOut(3)*inv_rho
    momOut( 9) = momOut( 9)*nEqFac + momOut(3)*momOut(4)*inv_rho
    momOut(10) = momOut(10)*nEqFac + momOut(4)*momOut(2)*inv_rho

  end function mus_intp_scaleMoments3D_incomp
! **************************************************************************** !

! **************************************************************************** !
  !> This function scales interpolated moments like density, momentum and
  !! shear stress to target level for 2D compressible model
  pure function mus_intp_scaleMoments2D(momIn, pFac, vFac, nEqFac) &
    & result(momOut)
    ! --------------------------------------------------------------------------
    !> Interpolated moments, dens, m_1, m_2, SS_XX, SS_YY, SS_XY
    real(kind=rk), intent(in) :: momIn(6)
    !> Scaling factor for Pressure
    real(kind=rk), intent(in) :: pFac
    !> Scaling factor for velocity
    real(kind=rk), intent(in) :: vFac
    !> Scaling factor for nonEq Shear stress moments
    real(kind=rk), intent(in) :: nEqFac
    !> Scaled moments, dens, m_1, m_2, SS_XX, SS_YY, SS_XY
    real(kind=rk) :: momOut(6)
    ! --------------------------------------------------------------------------
    real(kind=rk) :: tDens, inv_tDens
    ! --------------------------------------------------------------------------
    ! total density
    tDens = momIn(1) + rho0
    inv_tDens = 1.0_rk / tDens
    ! remove pi(0) parts from the second order moment
    ! that is ((rho u_\alpha * rho u_\beta) / rho)
    ! and additionally rho*cs2 from diagonal entries
    ! Pab = sum( c_ia*c_ib *f_i ) - kronecker_ab*(rho*cs2) - rho*u_a*u_b
    ! KM: its multiplied by inv_tDens because moments(2,3) are momentum not
    ! velocity
    momOut( 4) = momIn( 4) - momIn(2)*momIn(2)*inv_tDens - cs2*tDens
    momOut( 5) = momIn( 5) - momIn(3)*momIn(3)*inv_tDens - cs2*tDens
    momOut( 6) = momIn( 6) - momIn(2)*momIn(3)*inv_tDens

    ! Rescale density and momentum
    momOut(1) = momIn(1) * pFac + rho0
    momOut(2:3) = momIn(2:3) * vFac

    inv_tDens = 1.0_rk / momOut(1)
    ! Rescale the nonEq momOut
    ! and add back the scaled pi(0) part to the tensor
    momOut( 4) = momOut( 4)*nEqFac + momOut(2)*momOut(2)*inv_tDens &
      &         + cs2*momOut(1)
    momOut( 5) = momOut( 5)*nEqFac + momOut(3)*momOut(3)*inv_tDens &
      &         + cs2*momOut(1)
    momOut( 6) = momOut( 6)*nEqFac + momOut(2)*momOut(3)*inv_tDens

  end function mus_intp_scaleMoments2D
! **************************************************************************** !

! **************************************************************************** !
  !> This function scales interpolated moments like density, velocity and
  !! shear stress to target level for 2D incompressible model
  pure function mus_intp_scaleMoments2D_incomp(momIn, pFac, vFac, nEqFac) &
    & result(momOut)
    ! --------------------------------------------------------------------------
    !> Interpolated momentsents, dens, m_1, m_2, SS_XX, SS_YY, SS_XY
    real(kind=rk), intent(in) :: momIn(6)
    !> Scaling factor for Pressure
    real(kind=rk), intent(in) :: pFac
    !> Scaling factor for velocity
    real(kind=rk), intent(in) :: vFac
    !> Scaling factor for nonEq Shear stress momentsents
    real(kind=rk), intent(in) :: nEqFac
    !> Scaled moments, dens, m_1, m_2, SS_XX, SS_YY, SS_XY
    real(kind=rk) :: momOut(6)
    ! --------------------------------------------------------------------------
    real(kind=rk) :: tDens, inv_rho
    ! --------------------------------------------------------------------------
    ! total density
    tDens = momIn(1) + rho0
    inv_rho = 1.0_rk / rho0

    ! remove pi(0) parts from the second order moments
    ! that is ((rho u_\alpha * rho u_\beta) / rho)
    ! and additionally rho*cs2 from diagonal entries
    ! Pab = sum( c_ia*c_ib *f_i ) - kronecker_ab*(rho*cs2) - rho*u_a*u_b
    momOut( 4) = momIn( 4) - momIn(2)*momIn(2)*inv_rho - cs2*tDens
    momOut( 5) = momIn( 5) - momIn(3)*momIn(3)*inv_rho - cs2*tDens
    momOut( 6) = momIn( 6) - momIn(2)*momIn(3)*inv_rho

    ! Rescale density and momentum
    momOut(1) = momIn(1) * pFac + rho0
    momOut(2:3) = momIn(2:3) * vFac

    ! Rescale the nonEq moments
    ! and add back the scaled pi(0) part to the tensor
    momOut( 4) = momOut( 4)*nEqFac + momOut(2)*momOut(2)*inv_rho &
      &         + cs2*momOut(1)
    momOut( 5) = momOut( 5)*nEqFac + momOut(3)*momOut(3)*inv_rho &
      &         + cs2*momOut(1)
    momOut( 6) = momOut( 6)*nEqFac + momOut(2)*momOut(3)*inv_rho

  end function mus_intp_scaleMoments2D_incomp
! **************************************************************************** !


! **************************************************************************** !
  !> This function computes converts moments to pdf by computing
  !! equilibrium from dens and vel and use it to compute higher moments
  !! and transform moments to PDF
  pure function mus_intp_convertMomToPDF3D(moments, nonEqScalingFacs, layout) &
    & result(PDF)
    ! --------------------------------------------------------------------------
    !> Scheme layout
    type(mus_scheme_layout_type), intent(in) :: layout
    !> nonEquilibrium scaling factor
    real(kind=rk), intent(in) :: nonEqScalingFacs(layout%fStencil%QQ)
    !> All moments
    real(kind=rk), intent(in) :: moments(layout%fStencil%QQ)
    !> output pdf
    real(kind=rk) :: PDF(layout%fStencil%QQ)
    ! --------------------------------------------------------------------------
    integer :: QQ, iDir
    real(kind=rk) :: rho, inv_rho, vel(3)
    real(kind=rk) :: tmpMom(layout%fStencil%QQ)
    real(kind=rk) :: mEq(layout%fStencil%QQ)
    ! --------------------------------------------------------------------------
    QQ = layout%fStencil%QQ
    rho = moments(1)
    inv_rho = 1.0_rk/rho
    vel(1) = moments(layout%moment%first_moments(1))*inv_rho
    vel(2) = moments(layout%moment%first_moments(2))*inv_rho
    vel(3) = moments(layout%moment%first_moments(3))*inv_rho
    ! calculate Eq
    do iDir = 1, QQ
?? copy :: EqDir_comp( PDF(iDir), rho, vel, layout, iDir, cs2inv, rk )
    end do

    ! equilibrium moments
    mEq = matmul( layout%moment%toMoments%A, PDF )

    ! Scaled moments
    tmpMom = mEq + nonEqScalingFacs*(moments - mEq)

    ! transform the moments back to the pdfs
    PDF = matmul( layout%moment%toPdf%A, tmpMom )

  end function mus_intp_convertMomToPDF3D
! **************************************************************************** !

! **************************************************************************** !
  !> This function computes converts moments to pdf by computing
  !! equilibrium from dens and vel and use it to compute higher moments
  !! and transform moments to PDF
  pure function mus_intp_convertMomToPDF3D_incomp(moments, nonEqScalingFacs, &
    &                                             layout) result(PDF)
    ! --------------------------------------------------------------------------
    !> Scheme layout
    type(mus_scheme_layout_type), intent(in) :: layout
    !> nonEquilibrium scaling factor
    real(kind=rk), intent(in) :: nonEqScalingFacs(layout%fStencil%QQ)
    !>  all moments
    real(kind=rk), intent(in) :: moments(layout%fStencil%QQ)
    !> output pdf
    real(kind=rk) :: PDF(layout%fStencil%QQ)
    ! --------------------------------------------------------------------------
    integer :: QQ, iDir
    real(kind=rk) :: rho, vel(3)
    real(kind=rk) :: tmpMom(layout%fStencil%QQ)
    real(kind=rk) :: mEq(layout%fStencil%QQ)
    ! --------------------------------------------------------------------------
    QQ = layout%fStencil%QQ
    rho = moments(1)
    vel(1) = moments(layout%moment%first_moments(1))
    vel(2) = moments(layout%moment%first_moments(2))
    vel(3) = moments(layout%moment%first_moments(3))
    ! calculate Eq
    do iDir = 1, layout%fStencil%QQ
?? copy :: EqDir_incomp( PDF(iDir), rho, vel, layout, iDir, rho0, cs2inv, rk )
    end do

    ! equilibrium moments
    mEq = matmul( layout%moment%toMoments%A, PDF )

    ! Scaled moments
    tmpMom = mEq + nonEqScalingFacs*(moments - mEq)

    ! transform the moments back to the pdfs
    PDF = matmul( layout%moment%toPdf%A, tmpMom )

  end function mus_intp_convertMomToPDF3D_incomp
! **************************************************************************** !

! **************************************************************************** !
  !> This function computes converts moments to pdf by computing
  !! equilibrium from dens and vel and use it to compute higher moments
  !! and transform moments to PDF
  pure function mus_intp_convertMomToPDF2D(moments, nonEqScalingFacs, layout) &
    &                                      result(PDF)
    ! --------------------------------------------------------------------------
    !> Scheme layout
    type(mus_scheme_layout_type), intent(in) :: layout
    !> nonEquilibrium scaling factor
    real(kind=rk), intent(in) :: nonEqScalingFacs(layout%fStencil%QQ)
    !> All moments
    real(kind=rk), intent(in) :: moments(layout%fStencil%QQ)
    !> output pdf
    real(kind=rk) :: PDF(layout%fStencil%QQ)
    ! --------------------------------------------------------------------------
    integer :: QQ
    real(kind=rk) :: rho, vel(2)
    real(kind=rk) :: tmpMom(layout%fStencil%QQ)
    real(kind=rk) :: mEq(layout%fStencil%QQ)
?? copy :: var_fEq_d2q9
    ! --------------------------------------------------------------------------
    QQ = layout%fStencil%QQ
    rho = moments(1)
    vel(1) = moments(layout%moment%first_moments(1))/rho
    vel(2) = moments(layout%moment%first_moments(2))/rho
    ! calculate Eq
?? copy :: fEq_d2q9_a( PDF, vel(1), vel(2), rho, rho )

    ! equilibrium moments
    mEq = matmul( layout%moment%toMoments%A, PDF )

    ! Scaled moments
    tmpMom = mEq + nonEqScalingFacs*(moments - mEq)

    ! transform the moments back to the pdfs
    PDF = matmul( layout%moment%toPdf%A, tmpMom )

  end function mus_intp_convertMomToPDF2D
! **************************************************************************** !

! **************************************************************************** !
  !> This function computes converts moments to pdf by computing
  !! equilibrium from dens and vel and use it to compute higher moments
  !! and transform moments to PDF
  pure function mus_intp_convertMomToPDF2D_incomp(moments, nonEqScalingFacs,   &
    &                                             layout) result(PDF)
    ! --------------------------------------------------------------------------
    !> Scheme layout
    type(mus_scheme_layout_type), intent(in) :: layout
    !> nonEquilibrium scaling factor
    real(kind=rk), intent(in) :: nonEqScalingFacs(layout%fStencil%QQ)
    !> All moments
    real(kind=rk), intent(in) :: moments(layout%fStencil%QQ)
    !> output pdf
    real(kind=rk) :: PDF(layout%fStencil%QQ)
    ! --------------------------------------------------------------------------
    real(kind=rk) :: tmpMom(layout%fStencil%QQ)
    real(kind=rk) :: mEq(layout%fStencil%QQ)
    integer :: QQ
    real(kind=rk) :: vel(2)
?? copy :: var_fEq_d2q9
    ! --------------------------------------------------------------------------
    QQ = layout%fStencil%QQ
    vel(1) = moments(layout%moment%first_moments(1))
    vel(2) = moments(layout%moment%first_moments(2))
    ! calculate Eq
?? copy :: fEq_d2q9_a( PDF, vel(1), vel(2), moments(1), rho0 )

    ! equilibrium moments
    mEq = matmul( layout%moment%toMoments%A, PDF )

    ! Scaled moments
    tmpMom = mEq + nonEqScalingFacs*(moments - mEq)

    ! transform the moments back to the pdfs
    PDF = matmul( layout%moment%toPdf%A, tmpMom )

  end function mus_intp_convertMomToPDF2D_incomp
! **************************************************************************** !

end module mus_interpolate_tools_module
! ****************************************************************************** !
