! Copyright (c) 2013, 2016, 2019 Harald Klimach <harald.klimach@uni-siegen.de>
! Copyright (c) 2013 Manuel Hasert <m.hasert@grs-sim.de>
! Copyright (c) 2013-2020 Kannan Masilamani <kannan.masilamani@uni-siegen.de>
! Copyright (c) 2013-2014 Simon Zimny <s.zimny@grs-sim.de>
! Copyright (c) 2013-2016 Jiaxing Qi <jiaxing.qi@uni-siegen.de>
! Copyright (c) 2016 Verena Krupp <verena.krupp@uni-siegen.de>
! Copyright (c) 2016 Tobias Schneider <tobias1.schneider@student.uni-siegen.de>
! Copyright (c) 2016-2018 Raphael Haupt <raphael.haupt@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: Kannan Masilamani
!! author: Jiaxing Qi
!! This module provides the MUSUBI specific functions for calculating
!! macroscopic quantities from the state variables for incompressible LBM
!! models.\n
!! Notice that only those quantities that related to density should have a
!! formula which differs from normal LBM model.
!!
!! The depending common interface between MUSUBI and ATELES is defined in the
!! tem_derived_module. The functionality for accessing a variable from the state
!! and evaluating a lua function are also provided in the tem_derived module.
!!
!! Do not use get_Element or get_Point routines to update the state !
!!

?? include 'header/lbm_macros.inc'
?? include 'header/lbm_deriveMacros.inc'
?? include 'header/lbm_d3q19Macros.inc'
?? include 'treelm/source/deriveMacros.inc'
module mus_derQuanIncomp_module
  use iso_c_binding, only: c_loc, c_ptr, c_f_pointer

  ! include treelm modules
  use tem_param_module,         only: div1_2, div1_3, div1_54, div1_9, div3_4, &
    &                                 div1_36, div3_4h,                        &
    &                                 sqrt3, cs2inv, cs2, t2cs4inv, t2cs2inv,  &
    &                                 cs4inv, rho0, rho0Inv, q000
  use env_module,               only: rk, long_k, labelLen
  use tem_float_module,         only: operator(.feq.), operator(.fge.), &
    &                                 operator(.fle.)
  use tem_varSys_module,        only: tem_varSys_type, tem_varSys_op_type,     &
    &                                 tem_varSys_append_derVar,                &
    &                                 tem_varSys_proc_point,                   &
    &                                 tem_varSys_proc_element,                 &
    &                                 tem_varSys_proc_setParams,               &
    &                                 tem_varSys_proc_getParams,               &
    &                                 tem_varSys_proc_setupIndices,            &
    &                                 tem_varSys_proc_getValOfIndex,           &
    &                                 tem_varSys_getPoint_dummy,               &
    &                                 tem_varSys_getElement_dummy,             &
    &                                 tem_varSys_setupIndices_dummy,           &
    &                                 tem_varSys_getValOfIndex_dummy,          &
    &                                 tem_varSys_setParams_dummy,              &
    &                                 tem_varSys_getParams_dummy
  use tem_variable_module,      only: tem_variable_type
  use tem_stencil_module,       only: tem_stencilHeader_type
  use tem_topology_module,      only: tem_levelOf
  use tem_time_module,          only: tem_time_type
  use treelmesh_module,         only: treelmesh_type
  use tem_subTree_type_module,  only: tem_subTree_type, tem_treeIDfrom_subTree
  use tem_aux_module,           only: tem_abort
  use tem_logging_module,       only: logUnit
  use tem_operation_var_module, only: tem_evalMag_forElement,      &
    &                                 tem_evalMag_forPoint,        &
    &                                 tem_evalMag_fromIndex,       &
    &                                 tem_opVar_setupIndices,      &
    &                                 tem_get_new_varSys_data_ptr, &
    &                                 tem_opVar_setParams,         &
    &                                 tem_opVar_getParams
  use tem_debug_module,         only: dbgUnit
  use tem_grow_array_module,    only: grw_labelarray_type, append


  ! include musubi modules
  use mus_source_var_module,         only: mus_source_op_type
  use mus_pdf_module,                only: pdf_data_type
  use mus_varSys_module,             only: mus_varSys_data_type,               &
    &                                      mus_varSys_solverData_type,         &
    &                                      mus_get_new_solver_ptr,             &
    &                                      mus_deriveVar_ForPoint,             &
    &                                      mus_generic_varFromPDF_fromIndex,   &
    &                                      mus_generic_fromPDF_forElement,     &
    &                                      mus_derive_fromPDF
  use mus_stateVar_module,           only: mus_accessVar_setupIndices,         &
    &                                      mus_stateVar_Fetch_fromIndex,       &
    &                                      mus_stateVar_Fetch_now_fromIndex,   &
    &                                      mus_access_stateFetch_ForElement,   &
    &                                      mus_access_stateFetch_now_ForElement
  use mus_scheme_header_module,      only: mus_scheme_header_type
  use mus_scheme_layout_module,      only: mus_scheme_layout_type
  use mus_scheme_type_module,        only: mus_scheme_type
  use mus_derivedQuantities_module2, only: secondMom
  use mus_derQuan_module,            only: deriveDensity,                      &
    &                                      deriveDensity_fromIndex,            &
    &                                      derivePressure,                     &
    &                                      derivePressure_fromIndex,           &
    &                                      deriveShearStress,                  &
    &                                      deriveShearMag,                     &
    &                                      deriveWSS2D,                        &
    &                                      deriveWSS3D,                        &
    &                                      deriveTemp,                         &
    &                                      deriveShearRate,                    &
    &                                      deriveBndForce
  use mus_operation_var_module,      only: mus_opVar_setupIndices,         &
    &                                      mus_opVar_gradU_forElement,     &
    &                                      mus_opVar_vorticity_forElement, &
    &                                      mus_opVar_QCriterion_forElement
  use mus_derVarPos_module,          only: mus_derVarPos_type
  use mus_physics_module,            only: mus_convertFac_type

  implicit none

  private

  public :: mus_append_derVar_fluidIncomp
  ! equilbrium from macro uses different interface defined in
  ! mus_variable_module
  public :: deriveEquilIncomp_FromMacro
  public :: deriveEquilIncomp_FromMacro_d3q19
  public :: deriveVelIncomp_FromState
  public :: deriveVelIncomp_FromState_d3q19
  public :: deriveVelIncomp_FromPreColState
  public :: deriveEqIncomp_FromState
  public :: deriveAuxIncomp_fromState
  public :: deriveEquilIncomp_fromAux

  ! source variables
  public :: derive_absorbLayerIncomp
  public :: applySrc_absorbLayerIncomp

contains


  ! ************************************************************************* !
  !> subroutine to add derive variables for incompressible LBM
  !! (schemekind = 'fluid_incompressible') to the varsys.
  subroutine mus_append_derVar_fluidIncomp( varSys, solverData, schemeHeader, &
    &                                     stencil, fldLabel, derVarName )
    ! ---------------------------------------------------------------------- !
    !> global variable system
    type(tem_varSys_type), intent(inout)  :: varSys

    !> Contains pointer to solver data types
    type(mus_varSys_solverData_type), target, intent(in) :: solverData

    !> identifier of the scheme
    type(mus_scheme_header_type), intent(in)  :: schemeHeader

    !> compute stencil defintion
    type(tem_stencilHeader_type), intent(in)       :: stencil

    !> array of field label prefix. Size=nFields
    character(len=*), intent(in)              :: fldLabel

    !> array of derive physical variables
    type(grw_labelarray_type), intent(inout) :: derVarName
    ! ---------------------------------------------------------------------- !
    ! number of derive variables
    integer :: nDerVars, iVar, nComponents, addedPos, iIn
    logical :: wasAdded
    character(len=labelLen), allocatable ::  input_varname(:)
    character(len=labelLen)  ::  varName
    procedure(tem_varSys_proc_point), pointer :: get_point => NULL()
    procedure(tem_varSys_proc_element), pointer :: get_element => NULL()
    procedure(tem_varSys_proc_setParams), pointer :: set_params => null()
    procedure(tem_varSys_proc_getParams), pointer :: get_params => null()
    procedure(tem_varSys_proc_setupIndices), pointer :: &
      &                                      setup_indices => null()
    procedure(tem_varSys_proc_getValOfIndex), pointer :: &
      &                                       get_valOfIndex => null()
    type(c_ptr) :: method_data
    character(len=labelLen), allocatable :: derVarName_loc(:)
    ! ---------------------------------------------------------------------- !
    nullify(get_point, get_element, set_params, get_params, setup_indices, &
      &     get_valOfIndex)

    nDerVars = 19
    allocate(derVarName_loc(nDerVars))
    derVarName_loc = [ 'fetch_pdf      ', 'fetch_pdf_now  ', &
      &                'pressure       ', 'equilibrium    ', &
      &                'non_equilibrium', 'kinetic_energy ', &
      &                'shear_stress   ', 'strain_rate    ', &
      &                'shear_rate     ', 'wss            ', &
      &                'momentum       ', 'vel_mag        ', &
      &                'bnd_force      ', 'fetch_pdf      ', &
      &                'shear_mag      ', 'temperature    ', &
      &                'grad_velocity  ', 'vorticity      ', &
      &                'q_criterion    '                     ]

    do iVar = 1, nDerVars
      call append(derVarName, derVarName_loc(iVar))
      ! set default pointers, overwrite if neccessary
      get_element => tem_varSys_getElement_dummy
      get_point => mus_deriveVar_ForPoint
      setup_indices => mus_opVar_setupIndices
      get_valOfIndex => tem_varSys_getvalOfIndex_dummy
      set_params => tem_varSys_setParams_dummy
      get_params => tem_varSys_getParams_dummy
      method_data  = mus_get_new_solver_ptr(solverData)

      select case(trim(adjustl(derVarName_loc(iVar))))
      case ('fetch_pdf')
        get_element => mus_access_stateFetch_ForElement
        get_valOfIndex => mus_stateVar_Fetch_fromIndex
        setup_indices => mus_accessVar_setupIndices
        nComponents = stencil%QQ
        allocate(input_varname(1))
        input_varname(1) = 'pdf'

      case ('fetch_pdf_now')
        get_element => mus_access_stateFetch_now_ForElement
        get_valOfIndex => mus_stateVar_Fetch_now_fromIndex
        setup_indices => mus_accessVar_setupIndices
        nComponents = stencil%QQ
        allocate(input_varname(1))
        input_varname(1) = 'pdf'

      case ('pressure')
        get_element => derivePressure
        get_valOfIndex => derivePressure_fromIndex
        nComponents = 1
        allocate(input_varname(1))
        input_varname(1) = 'density'

      case ('bnd_force')
        get_element => deriveBndForce
        nComponents = 3
        allocate(input_varname(1))
        input_varname(1) = 'pdf'

      case ('equilibrium')
        get_element => deriveEquilIncomp
        get_valOfIndex => deriveEquilIncomp_fromIndex
        nComponents = stencil%QQ
        allocate(input_varname(1))
        input_varname(1) = 'pdf'

      case ('non_equilibrium')
        get_element => deriveNonEquilIncomp
        get_valOfIndex => deriveNonEquilIncomp_fromIndex
        nComponents = stencil%QQ
        allocate(input_varname(1))
        input_varname(1) = 'fetch_pdf_now'

      case ('kinetic_energy')
        get_element => deriveKeIncomp
        get_ValOfIndex => deriveKEIncomp_fromIndex
        nComponents = 1
        allocate(input_varname(1))
        input_varname(1) = 'velocity'

      case ('temperature')
        get_element => deriveTemp
        nComponents = 1
        allocate(input_varname(0))

      case ('shear_stress')
        nComponents = 6
        get_element => deriveShearStress
        allocate(input_varname(1))
        input_varname(1) = 'non_equilibrium'

      case ('strain_rate')
        nComponents = 6
        get_element => deriveStrainRateIncomp
        get_ValOfIndex => deriveStrainRateIncomp_fromIndex
        allocate(input_varname(1))
        input_varname(1) = 'fetch_pdf_now'

      case ('shear_rate')
        get_element => deriveShearRate
        nComponents = 1
        allocate(input_varname(1))
        input_varname(1) = 'strain_rate'

      case ('wss')
        nComponents = 1
        allocate(input_varname(1))
        input_varname(1) = 'shear_stress'
        if (stencil%nDims == 2) then
          get_element => deriveWSS2D
        else if (stencil%nDims == 3) then
          get_element => deriveWSS3D
        else
          write(logUnit(1),*) 'WARNING: WSS does not support 1D'
        end if

      case ('momentum')
        get_element => deriveMomentumIncomp
        get_valOfIndex => deriveMomentumIncomp_fromIndex
        nComponents = 3
        allocate(input_varname(1))
        input_varname(1) = 'velocity'

      case ('grad_velocity')
        get_element => mus_opVar_gradU_forElement
        nComponents = 9
        allocate(input_varname(1))
        input_varname(1) = 'velocity'

      case ('vorticity')
        get_element => mus_opVar_vorticity_forElement
        nComponents = 3
        allocate(input_varname(1))
        input_varname(1) = 'velocity'

      case ('q_criterion')
        get_element => mus_opVar_QCriterion_forElement
        nComponents = 1
        allocate(input_varname(1))
        input_varname(1) = 'velocity'

      case ('vel_mag')
        get_element => tem_evalMag_forElement
        get_point => tem_evalMag_forPoint
        get_valOfIndex => tem_evalMag_fromIndex
        setup_indices => tem_opVar_setupIndices
        set_params => tem_opVar_setParams
        get_params => tem_opVar_getParams
        method_data = tem_get_new_varSys_data_ptr(method_data)
        nComponents = 1
        allocate(input_varname(1))
        input_varname(1) = 'velocity'

      case ('shear_mag')
        get_element => deriveShearMag
        nComponents = 1
        allocate(input_varname(1))
        input_varname(1) = 'shear_stress'

      case default
        write(logUnit(1),*) 'WARNING: Unknown variable: '//&
          &                 trim(derVarName_loc(iVar))
        cycle !go to next variable
      end select

      ! update variable names with field label
      varname = trim(fldLabel)//trim(adjustl(derVarName_loc(iVar)))
      do iIn = 1, size(input_varname)
        input_varname(iIn) = trim(fldLabel)//trim(input_varname(iIn))
      end do

      ! append variable to varSys
      call tem_varSys_append_derVar( me             = varSys,         &
        &                            varName        = trim(varname),  &
        &                            nComponents    = nComponents,    &
        &                            input_varname  = input_varname,  &
        &                            method_data    = method_data,    &
        &                            get_point      = get_point,      &
        &                            get_element    = get_element,    &
        &                            set_params     = set_params,     &
        &                            get_params     = get_params,     &
        &                            setup_indices  = setup_indices,  &
        &                            get_valOfIndex = get_valOfIndex, &
        &                            pos            = addedPos,       &
        &                            wasAdded       = wasAdded        )

      if (wasAdded) then
        write(logUnit(10),*) ' Appended variable: '//trim(varname)
      else if (addedpos < 1) then
        write(logUnit(1),*) 'Error: variable '//trim(varname)// &
          &                 ' is not added to variable system'
      end if

      deallocate(input_varname)
    end do

  end subroutine mus_append_derVar_fluidIncomp
  ! **************************************************************************** !


! ************************************************************************** !
!       Subroutines with common interface for the function pointers            !
! ************************************************************************** !

! ************************************************************************** !
  !> Calculate the velocity on given elements for incompressible model
  !! \( \vec{u} = \frac{\sum(f_i c_i)}{\rho_0} \)
  !!
  !! The interface has to comply to the abstract interface
  !! [[tem_varSys_module:tem_varSys_proc_element]].
  !!
?? copy :: get_element_headtxt(deriveVelocityIncomp)
    ! ---------------------------------------------------------------------- !
    !> Function pointer to perform specific operation.
    procedure(mus_derive_fromPDF), pointer :: fnCalcPtr
    ! ---------------------------------------------------------------------------
    fnCalcPtr => mus_deriveVelocityIncomp

    call mus_generic_fromPDF_forElement( &
      &  fun       = fun,                &
      &  varSys    = varSys,             &
      &  elempos   = elempos,            &
      &  tree      = tree,               &
      &  time      = time,               &
      &  nVals     = nElems,             &
      &  fnCalcPtr = fnCalcPtr,          &
      &  nDofs     = nDofs,              &
      &  res       = res                 )

  end subroutine deriveVelocityIncomp
! ************************************************************************** !

! ****************************************************************************** !
  !> Calculate momentum from velocity stored in auxField
  !!
  !! The interface has to comply to the abstract interface
  !! [[tem_varSys_module:tem_varSys_proc_element]].
  !!
?? copy :: get_element_headtxt(deriveMomentumIncomp)
    ! ---------------------------------------------------------------------------
    !> Function pointer to perform specific operation.
    integer :: statePos, iElem, iComp, iLevel, elemOff
    type(mus_varSys_data_type), pointer :: fPtr
    integer :: dens_pos, vel_pos(3)
    real(kind=rk) :: rho
    ! ---------------------------------------------------------------------------
    call C_F_POINTER( fun%method_Data, fPtr )
    res = 0.0_rk
    associate( scheme => fPtr%solverData%scheme,                      &
      &        levelPointer => fPtr%solverData%geometry%levelPointer, &
      &        auxField => fPtr%solverData%scheme%auxField            )

      do iElem = 1, nElems
        ! if state array is defined level wise then use levelPointer(pos)
        ! to access state array
        statePos = levelPointer( elemPos(iElem) )
        iLevel = tem_levelOf( tree%treeID( elemPos(iElem) ) )

        ! element offset for auxField
        elemoff = (statePos-1)*varSys%nAuxScalars
        ! position of velocity in auxField array
        vel_pos = varSys%method%val( fun%input_varPos(1) ) &
          &                     %auxField_varPos(1:3)

        ! compute and store momentum
        res((iElem-1)*3 + 1) = rho0 * auxField(iLevel)%val(elemOff + vel_pos(1))
        res((iElem-1)*3 + 2) = rho0 * auxField(iLevel)%val(elemOff + vel_pos(2))
        res((iElem-1)*3 + 3) = rho0 * auxField(iLevel)%val(elemOff + vel_pos(3))

      end do ! iElem
    end associate

  end subroutine deriveMomentumIncomp
! ****************************************************************************** !



! ****************************************************************************** !
  !> Initiates the calculation of equlibrium
  !! This routine sets the function Pointer for equlibrium calcualtion and calls
  !! the generice get Element from PDF routine
  !!
  !! The interface has to comply to the abstract interface
  !! [[tem_varSys_module:tem_varSys_proc_element]].
  !!
?? copy :: get_element_headtxt(deriveEquilIncomp)
    ! ---------------------------------------------------------------------------
    !> Function pointer to perform specific operation.
    procedure(mus_derive_fromPDF), pointer :: fnCalcPtr
    ! ---------------------------------------------------------------------------
    fnCalcPtr => mus_deriveEquilIncomp

    call mus_generic_fromPDF_forElement( &
      &  fun       = fun,                &
      &  varSys    = varSys,             &
      &  elempos   = elempos,            &
      &  tree      = tree,               &
      &  time      = time,               &
      &  nVals     = nElems,             &
      &  fnCalcPtr = fnCalcPtr,          &
      &  nDofs     = nDofs,              &
      &  res       = res                 )

  end subroutine deriveEquilIncomp
! ************************************************************************** !


! ****************************************************************************** !
  !> Initiates the calculation of NonEquil
  !! This routine sets the function Pointer for NonEquil calcualtion and calls
  !! the generice get Element from PDF routine
  !!
  !! The interface has to comply to the abstract interface
  !! [[tem_varSys_module:tem_varSys_proc_element]].
  !!
?? copy :: get_element_headtxt(deriveNonEquilIncomp)
    ! ---------------------------------------------------------------------------
    !> Function pointer to perform specific operation.
    procedure(mus_derive_fromPDF), pointer :: fnCalcPtr
    ! ---------------------------------------------------------------------------
    fnCalcPtr => mus_deriveNonEquilIncomp

    call mus_generic_fromPDF_forElement( &
      &  fun       = fun,                &
      &  varSys    = varSys,             &
      &  elempos   = elempos,            &
      &  tree      = tree,               &
      &  time      = time,               &
      &  nVals     = nElems,             &
      &  fnCalcPtr = fnCalcPtr,          &
      &  nDofs     = nDofs,              &
      &  res       = res                 )

  end subroutine deriveNonEquilIncomp
! ************************************************************************** !


! ****************************************************************************** !
  !> Initiates the calculation of StrainRate
  !! This routine sets the function Pointer for StrainRate calcualtion and calls
  !! the generice get Element from PDF routine
  !!
  !! The interface has to comply to the abstract interface
  !! [[tem_varSys_module:tem_varSys_proc_element]].
  !!
?? copy :: get_element_headtxt(deriveStrainRateIncomp)
    ! ---------------------------------------------------------------------------
    !> Function pointer to perform specific operation.
    procedure(mus_derive_fromPDF), pointer :: fnCalcPtr
    ! ---------------------------------------------------------------------------
    fnCalcPtr => mus_deriveStrainRateIncomp

    call mus_generic_fromPDF_forElement( &
      &  fun       = fun,                &
      &  varSys    = varSys,             &
      &  elempos   = elempos,            &
      &  tree      = tree,               &
      &  time      = time,               &
      &  nVals     = nElems,             &
      &  fnCalcPtr = fnCalcPtr,          &
      &  nDofs     = nDofs,              &
      &  res       = res                 )


  end subroutine deriveStrainRateIncomp
! ****************************************************************************** !

! **************************************************************************** !
!       Subroutines with common interface for the function pointers            !
!                              getValOfIndex                                   !
! **************************************************************************** !

! ****************************************************************************** !
  !> Initiates the calculation of Velocity.
  !! This routine sets the function Pointer for velocity calcualtion and calls
  !! the generice get Value of Index routine
  !!
  !! The interface has to comply to the abstract interface
  !! [[tem_varSys_module:tem_varSys_proc_getValOfIndex]].
  !!
  recursive subroutine deriveVelocityIncomp_fromIndex( fun, varSys, time, &
     &                                    iLevel, idx, idxLen, nVals, res )
      !> Description of the method to obtain the variables, here some preset
      !! values might be stored, like the space time function to use or the
      !! required variables.
      class(tem_varSys_op_type), intent(in) :: fun

      !> The variable system to obtain the variable from.
      type(tem_varSys_type), intent(in)     :: varSys

      !> Point in time at which to evaluate the variable.
      type(tem_time_type), intent(in)       :: time

      !> Level on which values are requested
      integer, intent(in)                   :: iLevel

      !> Index of points in the growing array and variable val array to
      !! return.
      !! Size: most times nVals, if contiguous arrays are used it depends
      !! on the number of first indices
      integer, intent(in)                   :: idx(:)

      !> With idx as start index in contiguous memory,
      !! idxLength defines length of each contiguous memory
      !! Size: dependes on number of first index for contiguous array,
      !! but the sum of all idxLen is equal to nVals
      integer, optional, intent(in)         :: idxLen(:)

      !> Number of values to obtain for this variable (vectorized access).
      integer, intent(in)                   :: nVals

      !> Resulting values for the requested variable.
      !!
      !! Dimension: n requested entries x nComponents of this variable
      !! Access: (iElem-1)*fun%nComponents + iComp
      real(kind=rk), intent(out)            :: res(:)
    ! ---------------------------------------------------------------------------
    !> Function pointer to perform specific operation.
    procedure(mus_derive_fromPDF), pointer  :: fnCalcPtr
    ! ---------------------------------------------------------------------------

    fnCalcPtr => mus_derivevelocityIncomp

    call mus_generic_varFromPDF_fromIndex( &
      &  fun       = fun,                  &
      &  varSys    = varSys,               &
      &  time      = time,                 &
      &  iLevel    = iLevel,               &
      &  idx       = idx,                  &
      &  nVals     = nVals,                &
      &  fnCalcPtr = fnCalcPtr,            &
      &  res       = res                   )

  end subroutine deriveVelocityIncomp_fromIndex
! ****************************************************************************** !


! ****************************************************************************** !
  !> Calculate Momentum from density and velocity in auxField.
  !!
  !! The interface has to comply to the abstract interface
  !! [[tem_varSys_module:tem_varSys_proc_getValOfIndex]].
  !!
  recursive subroutine deriveMomentumIncomp_fromIndex( fun, varSys, time, &
    & iLevel, idx, idxLen, nVals, res )
    ! ---------------------------------------------------------------------------
    !> Description of the method to obtain the variables, here some preset
    !! values might be stored, like the space time function to use or the
    !! required variables.
    class(tem_varSys_op_type), intent(in) :: fun

    !> The variable system to obtain the variable from.
    type(tem_varSys_type), intent(in)     :: varSys

    !> Point in time at which to evaluate the variable.
    type(tem_time_type), intent(in)       :: time

    !> Level on which values are requested
    integer, intent(in)                   :: iLevel

    !> Index
    !! return.
    !! Size: most times nVals, if contiguous arrays are used it depends
    !! on the number of first indices
    integer, intent(in)                   :: idx(:)

    !> With idx as start index in contiguous memory,
    !! idxLength defines length of each contiguous memory
    !! Size: dependes on number of first index for contiguous array,
    !! but the sum of all idxLen is equal to nVals
    integer, optional, intent(in)         :: idxLen(:)

    !> Number of values to obtain for this variable (vectorized access).
    integer, intent(in)                   :: nVals

    !> Resulting values for the requested variable.
    !!
    !! Dimension: n requested entries x nComponents of this variable
    !! Access: (iElem-1)*fun%nComponents + iComp
    real(kind=rk), intent(out)            :: res(:)
    ! ---------------------------------------------------------------------------
    type(mus_varSys_data_type), pointer :: fPtr
    integer :: vel_pos, iVal
    ! ---------------------------------------------------------------------------

?? if( DEBUG ) then
    ! check if number of index are the same as number of values asked for
    if (size(idx) /= nVals) then
      call tem_abort(                                                         &
        & 'in deriveChargeDensity_fromIndex: not the same amount of index and'&
        & //'values are provided, stopping')
    end if
?? endif

    call C_F_POINTER( fun%method_Data, fPtr )
    res = 0.0_rk

    ! get velocity values for IDX
    vel_pos = fun%input_varPos(1)
    call varSys%method%val( vel_pos )%get_ValOfIndex( &
      &     varSys = varSys,                          &
      &     time   = time,                            &
      &     iLevel = iLevel,                          &
      &     idx    = fPtr%opData%input_pntIndex(1)    &
      &              %indexLvl(iLevel)%val( idx(:) ), &
      &     nVals  = nVals,                           &
      &     res    = res(:)                           )

    ! convert velocity to momentum
    do iVal = 1, nVals
      res( (iVal-1)*3 + 1 : iVal*3 ) = rho0 * res( (iVal-1)*3 + 1 : iVal*3 )
    end do

  end subroutine deriveMomentumIncomp_fromIndex
! ****************************************************************************** !


! ****************************************************************************** !
  !> Initiates the calculation of equilibrium.
  !! This routine sets the function Pointer for equilibrium calcualtion and calls
  !! the generice get Value of Index routine
  !!
  !! The interface has to comply to the abstract interface
  !! [[tem_varSys_module:tem_varSys_proc_getValOfIndex]].
  !!
  recursive subroutine deriveEquilIncomp_fromIndex( fun, varSys, time, iLevel, &
     &                                                 idx, idxLen, nVals, res )
      !> Description of the method to obtain the variables, here some preset
      !! values might be stored, like the space time function to use or the
      !! required variables.
      class(tem_varSys_op_type), intent(in) :: fun

      !> The variable system to obtain the variable from.
      type(tem_varSys_type), intent(in)     :: varSys

      !> Point in time at which to evaluate the variable.
      type(tem_time_type), intent(in)       :: time

      !> Level on which values are requested
      integer, intent(in)                   :: iLevel

      !> Index of points in the growing array and variable val array to
      !! return.
      !! Size: most times nVals, if contiguous arrays are used it depends
      !! on the number of first indices
      integer, intent(in)                   :: idx(:)

      !> With idx as start index in contiguous memory,
      !! idxLength defines length of each contiguous memory
      !! Size: dependes on number of first index for contiguous array,
      !! but the sum of all idxLen is equal to nVals
      integer, optional, intent(in)         :: idxLen(:)

      !> Number of values to obtain for this variable (vectorized access).
      integer, intent(in)                   :: nVals

      !> Resulting values for the requested variable.
      !!
      !! Dimension: n requested entries x nComponents of this variable
      !! Access: (iElem-1)*fun%nComponents + iComp
      real(kind=rk), intent(out)            :: res(:)
    ! ---------------------------------------------------------------------------
    !> Function pointer to perform specific operation.
    procedure(mus_derive_fromPDF), pointer  :: fnCalcPtr
    ! ---------------------------------------------------------------------------

    fnCalcPtr => mus_deriveEquilIncomp

    call mus_generic_varFromPDF_fromIndex( &
      &  fun       = fun,                  &
      &  varSys    = varSys,               &
      &  time      = time,                 &
      &  iLevel    = iLevel,               &
      &  idx       = idx,                  &
      &  nVals     = nVals,                &
      &  fnCalcPtr = fnCalcPtr,            &
      &  res       = res                   )

  end subroutine deriveEquilIncomp_fromIndex
! ****************************************************************************** !


! ****************************************************************************** !
  !> Initiates the calculation of non_equilibrium.
  !! This routine sets the function Pointer for non_equilibrium calcualtion and
  !! calls the generice get Value of Index routine
  !!
  !! The interface has to comply to the abstract interface
  !! [[tem_varSys_module:tem_varSys_proc_getValOfIndex]].
  !!
  recursive subroutine deriveNonEquilIncomp_fromIndex( fun, varSys, time, iLevel, idx, &
     &                                        idxLen, nVals, res )
      !> Description of the method to obtain the variables, here some preset
      !! values might be stored, like the space time function to use or the
      !! required variables.
      class(tem_varSys_op_type), intent(in) :: fun

      !> The variable system to obtain the variable from.
      type(tem_varSys_type), intent(in)     :: varSys

      !> Point in time at which to evaluate the variable.
      type(tem_time_type), intent(in)       :: time

      !> Level on which values are requested
      integer, intent(in)                   :: iLevel

      !> Index of points in the growing array and variable val array to
      !! return.
      !! Size: most times nVals, if contiguous arrays are used it depends
      !! on the number of first indices
      integer, intent(in)                   :: idx(:)

      !> With idx as start index in contiguous memory,
      !! idxLength defines length of each contiguous memory
      !! Size: dependes on number of first index for contiguous array,
      !! but the sum of all idxLen is equal to nVals
      integer, optional, intent(in)         :: idxLen(:)

      !> Number of values to obtain for this variable (vectorized access).
      integer, intent(in)                   :: nVals

      !> Resulting values for the requested variable.
      !!
      !! Dimension: n requested entries x nComponents of this variable
      !! Access: (iElem-1)*fun%nComponents + iComp
      real(kind=rk), intent(out)            :: res(:)
    ! ---------------------------------------------------------------------------
    !> Function pointer to perform specific operation.
    procedure(mus_derive_fromPDF), pointer  :: fnCalcPtr
    ! ---------------------------------------------------------------------------

    fnCalcPtr => mus_deriveNonEquilIncomp

    call mus_generic_varFromPDF_fromIndex( &
      &  fun       = fun,                  &
      &  varSys    = varSys,               &
      &  time      = time,                 &
      &  iLevel    = iLevel,               &
      &  idx       = idx,                  &
      &  nVals     = nVals,                &
      &  fnCalcPtr = fnCalcPtr,            &
      &  res       = res                   )

  end subroutine deriveNonEquilIncomp_fromIndex
! ****************************************************************************** !


! ****************************************************************************** !
  !> Initiates the calculation of kinetic_energy.
  !! This routine sets the function Pointer for kinetic_energy calcualtion and
  !! calls the generice get Value of Index routine
  !!
  !! The interface has to comply to the abstract interface
  !! [[tem_varSys_module:tem_varSys_proc_getValOfIndex]].
  !!
  recursive subroutine deriveKeIncomp_fromIndex( fun, varSys, time, iLevel,    &
    &                                            idx, idxLen, nVals, res )
    ! ---------------------------------------------------------------------------
    !> Description of the method to obtain the variables, here some preset
    !! values might be stored, like the space time function to use or the
    !! required variables.
    class(tem_varSys_op_type), intent(in) :: fun

    !> The variable system to obtain the variable from.
    type(tem_varSys_type), intent(in)     :: varSys

    !> Point in time at which to evaluate the variable.
    type(tem_time_type), intent(in)       :: time

    !> Level on which values are requested
    integer, intent(in)                   :: iLevel

    !> Index of points in the growing array and variable val array to
    !! return.
    !! Size: most times nVals, if contiguous arrays are used it depends
    !! on the number of first indices
    integer, intent(in)                   :: idx(:)

    !> With idx as start index in contiguous memory,
    !! idxLength defines length of each contiguous memory
    !! Size: dependes on number of first index for contiguous array,
    !! but the sum of all idxLen is equal to nVals
    integer, optional, intent(in)         :: idxLen(:)

    !> Number of values to obtain for this variable (vectorized access).
    integer, intent(in)                   :: nVals

    !> Resulting values for the requested variable.
    !!
    !! Dimension: n requested entries x nComponents of this variable
    !! Access: (iElem-1)*fun%nComponents + iComp
    real(kind=rk), intent(out)            :: res(:)
    ! ---------------------------------------------------------------------------
    type(mus_varSys_data_type), pointer :: fPtr
    integer :: vel_pos, iVal
    real(kind=rk) :: vel(3)
    ! ---------------------------------------------------------------------------

?? if( DEBUG ) then
    ! check if number of index are the same as number of values asked for
    if (size(idx) /= nVals) then
      call tem_abort(                                                         &
        & 'in deriveChargeDensity_fromIndex: not the same amount of index and'&
        & //'values are provided, stopping')
    end if
?? endif

    call C_F_POINTER( fun%method_Data, fPtr )
    res = 0.0_rk

    ! get velocity values for IDX
    vel_pos = fun%input_varPos(1)
    call varSys%method%val( vel_pos )%get_ValOfIndex( &
      &     varSys = varSys,                          &
      &     time   = time,                            &
      &     iLevel = iLevel,                          &
      &     idx    = fPtr%opData%input_pntIndex(1)    &
      &              %indexLvl(iLevel)%val( idx(:) ), &
      &     nVals  = nVals,                           &
      &     res    = res(:)                           )

    ! convert velocity to momentum
    do iVal = 1, nVals
      vel(1) = res((iVal-1)*3 + 1)
      vel(2) = res((iVal-1)*3 + 2)
      vel(3) = res((iVal-1)*3 + 3)
      res(iVal) = sum( vel(:)*vel(:) ) * 0.5_rk * rho0
    end do


  end subroutine deriveKeIncomp_fromIndex
! ****************************************************************************** !


! ****************************************************************************** !
  !> Initiates the calculation of StrainRate.
  !! This routine sets the function Pointer for StrainRate calcualtion and
  !! calls the generice get Value of Index routine
  !!
  !! The interface has to comply to the abstract interface
  !! [[tem_varSys_module:tem_varSys_proc_getValOfIndex]].
  !!
  recursive subroutine deriveStrainRateIncomp_fromIndex( fun, varSys, time,    &
     &                                         iLevel, idx, idxLen, nVals, res )
      !> Description of the method to obtain the variables, here some preset
      !! values might be stored, like the space time function to use or the
      !! required variables.
      class(tem_varSys_op_type), intent(in) :: fun

      !> The variable system to obtain the variable from.
      type(tem_varSys_type), intent(in)     :: varSys

      !> Point in time at which to evaluate the variable.
      type(tem_time_type), intent(in)       :: time

      !> Level on which values are requested
      integer, intent(in)                   :: iLevel

      !> Index of points in the growing array and variable val array to
      !! return.
      !! Size: most times nVals, if contiguous arrays are used it depends
      !! on the number of first indices
      integer, intent(in)                   :: idx(:)

      !> With idx as start index in contiguous memory,
      !! idxLength defines length of each contiguous memory
      !! Size: dependes on number of first index for contiguous array,
      !! but the sum of all idxLen is equal to nVals
      integer, optional, intent(in)         :: idxLen(:)

      !> Number of values to obtain for this variable (vectorized access).
      integer, intent(in)                   :: nVals

      !> Resulting values for the requested variable.
      !!
      !! Dimension: n requested entries x nComponents of this variable
      !! Access: (iElem-1)*fun%nComponents + iComp
      real(kind=rk), intent(out)            :: res(:)
    ! ---------------------------------------------------------------------------
    !> Function pointer to perform specific operation.
    procedure(mus_derive_fromPDF), pointer  :: fnCalcPtr
    ! ---------------------------------------------------------------------------

    fnCalcPtr => mus_deriveStrainRateIncomp

    call mus_generic_varFromPDF_fromIndex( &
      &  fun       = fun,                  &
      &  varSys    = varSys,               &
      &  time      = time,                 &
      &  iLevel    = iLevel,               &
      &  idx       = idx,                  &
      &  nVals     = nVals,                &
      &  fnCalcPtr = fnCalcPtr,            &
      &  res       = res                   )

  end subroutine deriveStrainRateIncomp_fromIndex
! ****************************************************************************** !

! **************************************************************************** !
!                           Calculation routines                               !
! **************************************************************************** !


! ****************************************************************************** !
  !> Calculate kinetic energy from velocity in auxField
  !!
  !! The interface has to comply to the abstract interface
  !! [[tem_varSys_module:tem_varSys_proc_element]].
  !!
?? copy :: get_element_headtxt(deriveKeIncomp)
    ! ---------------------------------------------------------------------------
    !> Function pointer to perform specific operation.
    integer :: statePos, iElem, iLevel, elemOff
    type(mus_varSys_data_type), pointer :: fPtr
    integer :: vel_pos(3)
    real(kind=rk) :: vel(3)
    ! ---------------------------------------------------------------------------
    call C_F_POINTER( fun%method_Data, fPtr )
    res = 0.0_rk
    associate( scheme => fPtr%solverData%scheme,                      &
      &        levelPointer => fPtr%solverData%geometry%levelPointer, &
      &        auxField => fPtr%solverData%scheme%auxField            )

      do iElem = 1, nElems
        ! if state array is defined level wise then use levelPointer(pos)
        ! to access state array
        statePos = levelPointer( elemPos(iElem) )
        iLevel = tem_levelOf( tree%treeID( elemPos(iElem) ) )

        ! element offset for auxField
        elemoff = (statePos-1)*varSys%nAuxScalars

        ! position of velocity in auxField array
        vel_pos = varSys%method%val( fun%input_varPos(1) ) &
          &                     %auxField_varPos(1:3)

        ! velocity
        vel(1) = auxField(iLevel)%val(elemOff + vel_pos(1))
        vel(2) = auxField(iLevel)%val(elemOff + vel_pos(2))
        vel(3) = auxField(iLevel)%val(elemOff + vel_pos(3))

        ! compute and store kinetic energy
        res(iElem) = sum( vel(:)*vel(:) ) * 0.5_rk * rho0
      end do ! iElem
    end associate

  end subroutine deriveKeIncomp
! ************************************************************************** !


! ****************************************************************************** !
   !> Derive absorb layer variable defined as a source term.
?? copy :: get_element_headtxt(derive_absorbLayerIncomp)
    ! ---------------------------------------------------------------------------
    call tem_abort('Not implemented')
    ! ---------------------------------------------------------------------------
  end subroutine derive_absorbLayerIncomp
! ****************************************************************************** !


! **************************************************************************** !
  !> Update state with source variable "absorb_layer".
  !! absorb_layer is used to absorb the flow and gradually reduce the flow
  !! quantities like pressure and velocity to a fixed value for incompressible
  !! model. It is based on:
  !! Xu, H., & Sagaut, P. (2013). Analysis of the absorbing layers for the
  !! weakly-compressible lattice Boltzmann methods. Journal of Computational
  !! Physics, 245(x), 14–42.
  !!
  !! This subroutine's interface must match the abstract interface definition
  !! [[proc_apply_source]] in derived/[[mus_source_var_module]].f90 in order to
  !! be callable via [[mus_source_op_type:applySrc]] function pointer.
  subroutine applySrc_absorbLayerIncomp( fun, inState, outState, neigh,      &
    &                                    auxField, nPdfSize, iLevel, varSys, &
    &                                    time, phyConvFac, derVarPos         )
    ! -------------------------------------------------------------------- !
    !> Description of method to apply source terms
    class(mus_source_op_type), intent(in) :: fun

    !> input  pdf vector
    real(kind=rk), intent(in) :: inState(:)

    !> output pdf vector
    real(kind=rk), intent(inout) :: outState(:)

    !> connectivity Array corresponding to state vector
    integer,intent(in) :: neigh(:)

    !> auxField array
    real(kind=rk), intent(in) :: auxField(:)

    !> number of elements in state Array
    integer, intent(in) :: nPdfSize

    !> current level
    integer, intent(in) :: iLevel

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

    !> Point in time at which to evaluate the variable.
    type(tem_time_type), intent(in)  :: time

    !> Physics conversion factor for current level
    type(mus_convertFac_type), intent(in) :: phyConvFac

    !> position of derived quantities in varsys
    type(mus_derVarPos_type), intent(in) :: derVarPos(:)
    ! -------------------------------------------------------------------- !
    call tem_abort('Error: Absorb layer is not yet implemented')
  end subroutine applySrc_absorbLayerIncomp
! **************************************************************************** !


! ************************************************************************** !
  !> This routine computes equilbrium from density and velocity
  !!
  !! This subroutine's interface must match the abstract interface definition
  !! [[derive_FromMacro]] in derived/[[mus_derVarPos_module]].f90 in order to be
  !! callable via [[mus_derVarPos_type:equilFromMacro]] function pointer.
  subroutine deriveEquilIncomp_FromMacro( density, velocity, iField, nElems, &
    &                                     varSys, layout, res                )
    ! -------------------------------------------------------------------- !
    !> Array of density.
    !! Single species: dens_1, dens_2 .. dens_n
    !! multi-species: dens_1_sp1, dens_1_sp2, dens_2_sp1, dens_2_sp2 ...
    !!                dens_n_sp1, dens_n_sp2
    real(kind=rk), intent(in) :: density(:)

    !> Array of velocity.
    !! Size: dimension 1: n*nFields. dimension 2: 3 (nComp)
    !! 1st dimension arrangement for multi-species is same as density
    real(kind=rk), intent(in) :: velocity(:, :)

    !> Current field
    integer, intent(in) :: iField

    !> number of elements
    integer, intent(in) :: nElems

    !> variable system which is required to access fieldProp
    !! information via variable method data c_ptr
    type(tem_varSys_type), intent(in) :: varSys

    !> scheme layout contains stencil definition and lattice weights
    type(mus_scheme_layout_type), intent(in) :: layout

    !> Output of this routine
    !! Dimension: n*nComponents of res
    real(kind=rk), intent(out) :: res(:)
    ! -------------------------------------------------------------------- !
    real(kind=rk) :: fEq(layout%fStencil%QQ), vel(3), usq, ucx
    integer :: QQ, iElem, iDir
    ! ---------------------------------------------------------------------- !

    QQ = layout%fStencil%QQ
    do iElem = 1, nElems
      vel = velocity(:,iElem)

      usq = dot_product(vel, vel)*t2cs2inv

      do iDir = 1, QQ
        ucx = dot_product( layout%fStencil%cxDirRK(:, iDir), vel )

        ! calculate equilibrium density
        fEq( iDir ) = layout%weight( iDir ) * ( density(iElem) &
          & + ( cs2inv * ucx + ucx * ucx * t2cs4inv - usq ) )
      enddo

      res( (iElem-1)*QQ+1: iElem*QQ ) = fEq
    end do
  end subroutine deriveEquilIncomp_FromMacro
! ************************************************************************** !


! ************************************************************************** !
  !> This routine computes equilbrium from density and velocity
  !!
  !! This subroutine's interface must match the abstract interface definition
  !! [[derive_FromMacro]] in derived/[[mus_derVarPos_module]].f90 in order to be
  !! callable via [[mus_derVarPos_type:equilFromMacro]] function pointer.
  subroutine deriveEquilIncomp_FromMacro_d3q19( density, velocity, iField,  &
    &                                           nElems, varSys, layout, res )
    ! -------------------------------------------------------------------- !
    !> Array of density.
    !! Single species: dens_1, dens_2 .. dens_n
    !! multi-species: dens_1_sp1, dens_1_sp2, dens_2_sp1, dens_2_sp2 ...
    !!                dens_n_sp1, dens_n_sp2
    real(kind=rk), intent(in) :: density(:)

    !> Array of velocity.
    !! Size: dimension 1: n*nFields. dimension 2: 3 (nComp)
    !! 1st dimension arrangement for multi-species is same as density
    real(kind=rk), intent(in) :: velocity(:, :)

    !> Current field
    integer, intent(in) :: iField

    !> number of elements
    integer, intent(in) :: nElems

    !> variable system which is required to access fieldProp
    !! information via variable method data c_ptr
    type(tem_varSys_type), intent(in) :: varSys

    !> scheme layout contains stencil definition and lattice weights
    type(mus_scheme_layout_type), intent(in) :: layout

    !> Output of this routine
    !! Dimension: n*nComponents of res
    real(kind=rk), intent(out) :: res(:)
    ! -------------------------------------------------------------------- !
    integer :: iElem
    real(kind=rk) :: fEq(19), u_x, u_y, u_z
?? copy :: fEq_d3q19_var
    ! ---------------------------------------------------------------------- !

    do iElem = 1, nElems

      u_x = velocity(1,iElem)
      u_y = velocity(2,iElem)
      u_z = velocity(3,iElem)

?? copy :: fEq_d3q19_a( fEq, u_x, u_y, u_z, density(iElem), rho0 )

      res( (iElem-1)*19 +  1) = fEq( 1)
      res( (iElem-1)*19 +  2) = fEq( 2)
      res( (iElem-1)*19 +  3) = fEq( 3)
      res( (iElem-1)*19 +  4) = fEq( 4)
      res( (iElem-1)*19 +  5) = fEq( 5)
      res( (iElem-1)*19 +  6) = fEq( 6)
      res( (iElem-1)*19 +  7) = fEq( 7)
      res( (iElem-1)*19 +  8) = fEq( 8)
      res( (iElem-1)*19 +  9) = fEq( 9)
      res( (iElem-1)*19 + 10) = fEq(10)
      res( (iElem-1)*19 + 11) = fEq(11)
      res( (iElem-1)*19 + 12) = fEq(12)
      res( (iElem-1)*19 + 13) = fEq(13)
      res( (iElem-1)*19 + 14) = fEq(14)
      res( (iElem-1)*19 + 15) = fEq(15)
      res( (iElem-1)*19 + 16) = fEq(16)
      res( (iElem-1)*19 + 17) = fEq(17)
      res( (iElem-1)*19 + 18) = fEq(18)
      res( (iElem-1)*19 + 19) = fEq(19)

    end do

  end subroutine deriveEquilIncomp_FromMacro_d3q19
! ************************************************************************** !

! ************************************************************************** !
  !> This routine computes equilbrium from auxField
  !!
  !! This subroutine's interface must match the abstract interface definition
  !! [[derive_equilFromAux]] in derived/[[mus_derVarPos_module]].f90 in order to
  !! be callable via [[mus_derVarPos_type:equilFromAux]] function pointer.
  subroutine deriveEquilIncomp_fromAux( derVarPos, auxField, iField, nElems, &
    &                                   varSys, layout, fEq                  )
    ! -------------------------------------------------------------------- !
    !> Position of derive variable in variable system
    class(mus_derVarPos_type), intent(in) :: derVarPos
    !> Array of auxField.
    !! Single species: dens_1, vel_1, dens_2, vel_2, .. dens_n, vel_n
    !! multi-species: dens_1_sp1, vel_1_spc1, dens_1_sp2, vel_1_spc2,
    !!                dens_2_sp1, vel_2_spc2, dens_2_sp2, vel_2_spc2 ...
    !!                dens_n_sp1, vel_n_sp1, dens_n_sp2, vel_n_spc2
    !! Access: (iElem-1)*nAuxScalars + auxField_varPos
    real(kind=rk), intent(in) :: auxField(:)

    !> Current field
    integer, intent(in) :: iField

    !> number of elements
    integer, intent(in) :: nElems

    !> variable system which is required to access fieldProp
    !! information via variable method data c_ptr
    type(tem_varSys_type), intent(in) :: varSys

    !> scheme layout contains stencil definition and lattice weights
    type(mus_scheme_layout_type), intent(in) :: layout

    !> Output of this routine
    !! Dimension: n*QQ of res
    real(kind=rk), intent(out) :: fEq(:)
    ! -------------------------------------------------------------------- !
    real(kind=rk) :: rho, vel(3), usq, ucx
    integer :: QQ, iElem, iDir, elemOff
    integer :: dens_pos, vel_pos(3)
    ! ---------------------------------------------------------------------- !
    dens_pos = varSys%method%val(derVarPos%density)%auxField_varPos(1)
    vel_pos = varSys%method%val(derVarPos%velocity)%auxField_varPos(1:3)

    QQ = layout%fStencil%QQ
    !NEC$ ivdep
    do iElem = 1, nElems
      ! element offset
      elemoff = (iElem-1)*varSys%nAuxScalars
      ! density
      rho = auxField(elemOff + dens_pos)
      ! velocity
      vel(1) = auxField(elemOff + vel_pos(1))
      vel(2) = auxField(elemOff + vel_pos(2))
      vel(3) = auxField(elemOff + vel_pos(3))

      usq = ( vel(1)*vel(1) + vel(2)*vel(2) + vel(3)*vel(3) )*t2cs2inv

      !NEC$ shortloop
      do iDir = 1, QQ
        ucx = layout%fStencil%cxDirRK(1, iDir) * vel(1) &
          & + layout%fStencil%cxDirRK(2, iDir) * vel(2) &
          & + layout%fStencil%cxDirRK(3, iDir) * vel(3)

        ! calculate equilibrium density
        fEq( (iElem-1)*QQ + iDir ) = layout%weight( iDir ) * ( rho  &
          & + rho0*( cs2inv * ucx + ucx * ucx * t2cs4inv - usq ) )
      enddo

    end do
  end subroutine deriveEquilIncomp_fromAux
! ************************************************************************** !


! ************************************************************************** !
  !> This routine computes velocity from state array
  !!
  !! This subroutine's interface must match the abstract interface definition
  !! [[derive_FromState]] in derived/[[mus_derVarPos_module]].f90 in order to be
  !! callable via [[mus_derVarPos_type:velFromState]],
  !! [[mus_derVarPos_type:equilFromState]], [[mus_derVarPos_type:momFromState]],
  !! [[mus_derVarPos_type:velocitiesFromState]], and
  !! [[mus_derVarPos_type:momentaFromState]] function pointers.
  subroutine deriveVelIncomp_FromState( state, iField, nElems, varSys, layout, &
    &                                   res                                    )
    ! -------------------------------------------------------------------- !
    !> Array of state
    !! n * layout%fStencil%QQ * nFields
    real(kind=rk), intent(in) :: state(:)

    !> Current field
    integer, intent(in) :: iField

    !> number of elements
    integer, intent(in) :: nElems

    !> variable system which is required to access fieldProp
    !! information via variable method data c_ptr
    type(tem_varSys_type), intent(in) :: varSys

    !> scheme layout contains stencil definition and lattice weights
    type(mus_scheme_layout_type), intent(in) :: layout

    !> Output of this routine
    !! Dimension: n * nComponents of res
    real(kind=rk), intent(out) :: res(:)
    ! -------------------------------------------------------------------- !
    integer :: iElem, iDir
    real(kind=rk) :: f( layout%fStencil%QQ )

    do iElem = 1, nElems
      do iDir = 1, layout%fStencil%QQ
        f(iDir) = state( iDir+(iElem-1)* varSys%nScalars )
      end do
      res( (iElem-1)*3+1 ) = sum( f * layout%fStencil%cxDirRK(1,:) )
      res( (iElem-1)*3+2 ) = sum( f * layout%fStencil%cxDirRK(2,:) )
      res( (iElem-1)*3+3 ) = sum( f * layout%fStencil%cxDirRK(3,:) )
    end do

  end subroutine deriveVelIncomp_FromState
! ************************************************************************** !


! ************************************************************************** !
  !> This routine computes velocity from state array
  !!
  !! This subroutine's interface must match the abstract interface definition
  !! [[derive_FromState]] in derived/[[mus_derVarPos_module]].f90 in order to be
  !! callable via [[mus_derVarPos_type:velFromState]],
  !! [[mus_derVarPos_type:equilFromState]], [[mus_derVarPos_type:momFromState]],
  !! [[mus_derVarPos_type:velocitiesFromState]], and
  !! [[mus_derVarPos_type:momentaFromState]] function pointers.
  subroutine deriveVelIncomp_FromState_d3q19( state, iField, nElems, varSys, &
    &                                         layout, res                    )
    ! -------------------------------------------------------------------- !
    !> Array of state
    !! n * layout%fStencil%QQ * nFields
    real(kind=rk), intent(in) :: state(:)

    !> Current field
    integer, intent(in) :: iField

    !> number of elements
    integer, intent(in) :: nElems

    !> variable system which is required to access fieldProp
    !! information via variable method data c_ptr
    type(tem_varSys_type), intent(in) :: varSys

    !> scheme layout contains stencil definition and lattice weights
    type(mus_scheme_layout_type), intent(in) :: layout

    !> Output of this routine
    !! Dimension: n * nComponents of res
    real(kind=rk), intent(out) :: res(:)
    ! -------------------------------------------------------------------- !
    integer :: iElem
    real(kind=rk) :: f(18), u_x, u_y, u_z

    !NEC$ ivdep
    do iElem = 1, nElems

      f( 1) = state(  1+(iElem-1)*varSys%nScalars )
      f( 2) = state(  2+(iElem-1)*varSys%nScalars )
      f( 3) = state(  3+(iElem-1)*varSys%nScalars )
      f( 4) = state(  4+(iElem-1)*varSys%nScalars )
      f( 5) = state(  5+(iElem-1)*varSys%nScalars )
      f( 6) = state(  6+(iElem-1)*varSys%nScalars )
      f( 7) = state(  7+(iElem-1)*varSys%nScalars )
      f( 8) = state(  8+(iElem-1)*varSys%nScalars )
      f( 9) = state(  9+(iElem-1)*varSys%nScalars )
      f(10) = state( 10+(iElem-1)*varSys%nScalars )
      f(11) = state( 11+(iElem-1)*varSys%nScalars )
      f(12) = state( 12+(iElem-1)*varSys%nScalars )
      f(13) = state( 13+(iElem-1)*varSys%nScalars )
      f(14) = state( 14+(iElem-1)*varSys%nScalars )
      f(15) = state( 15+(iElem-1)*varSys%nScalars )
      f(16) = state( 16+(iElem-1)*varSys%nScalars )
      f(17) = state( 17+(iElem-1)*varSys%nScalars )
      f(18) = state( 18+(iElem-1)*varSys%nScalars )

?? copy :: vel_d3q19( u_x, u_y, u_z, f, 1.0_rk )

      res( (iElem-1)*3+1 ) = u_x
      res( (iElem-1)*3+2 ) = u_y
      res( (iElem-1)*3+3 ) = u_z

    end do

  end subroutine deriveVelIncomp_FromState_d3q19
! ************************************************************************** !


! ************************************************************************** !
  !> This routine computes velocity from precollision state array using FETCH
  !! macro.
  !!
  !! This subroutine's interface must match the abstract interface definition
  !! [[derive_FromPreColState]] in derived/[[mus_derVarPos_module]].f90 in order
  !! to be callable via [[mus_derVarPos_type:velFromPreColState]] function
  !! pointer.
  subroutine deriveVelIncomp_FromPreColState( state, neigh, iField, nSize, &
    &                                         nElems, varSys, layout, res  )
    ! -------------------------------------------------------------------- !
    !> Array of state
    !! n * layout%fStencil%QQ * nFields
    real(kind=rk), intent(in) :: state(:)

    !> connectivity array
    integer, intent(in) :: neigh(:)

    !> Current field
    integer, intent(in) :: iField

    !> number of elements in state array
    integer, intent(in) :: nSize

    !> number of elements
    integer, intent(in) :: nElems

    !> variable system which is required to access fieldProp
    !! information via variable method data c_ptr
    type(tem_varSys_type), intent(in) :: varSys

    !> scheme layout contains stencil definition and lattice weights
    type(mus_scheme_layout_type), intent(in) :: layout

    !> Output of this routine
    !! Dimension: n * nComponents of res
    real(kind=rk), intent(out) :: res(:)
    ! -------------------------------------------------------------------- !
    integer :: iElem, iDir
    real(kind=rk) :: pdf( layout%fStencil%QQ ), vel(3)
    integer :: QQ, nScalars, nDims
    ! --------------------------------------------------------------------------
    QQ = layout%fStencil%QQ
    nScalars = varSys%nScalars
    nDims = layout%fStencil%nDims

    do iElem = 1, nElems
      do iDir = 1, QQ
        pdf(iDir) = state(                                          &
          & ?FETCH?(iDir, iField, iElem, QQ, nScalars, nSize, neigh))
      end do
      ! momentum
      vel( 1 ) = sum( pdf * layout%fStencil%cxDirRK(1,:) )
      vel( 2 ) = sum( pdf * layout%fStencil%cxDirRK(2,:) )
      vel( 3 ) = sum( pdf * layout%fStencil%cxDirRK(3,:) )
      ! return velocity field according on stencil dimensions
      res( (iElem-1)*nDims+1:iElem*nDims) = vel(1:nDims)
    end do
    ! convert to velocity
    res = res / rho0

  end subroutine deriveVelIncomp_FromPreColState
! ************************************************************************** !

! ************************************************************************** !
  !> This routine computes auxField from state array
  !!
  !! This subroutine's interface must match the abstract interface definition
  !! [[derive_auxFromState]] in derived/[[mus_derVarPos_module]].f90 in order to
  !! be callable via [[mus_derVarPos_type:auxFieldFromState]] function pointer.
  subroutine deriveAuxIncomp_fromState( derVarPos, state, iField, nElems, &
    &                                   nSize, iLevel, stencil, varSys,   &
    &                                   auxField                          )
    ! -------------------------------------------------------------------- !
    !> Position of derive variable in variable system
    class(mus_derVarPos_type), intent(in) :: derVarPos
    !> Array of state
    !! n * layout%stencil(1)%QQ * nFields
    real(kind=rk), intent(in) :: state(:)

    !> Current field
    integer, intent(in) :: iField

    !> number of elements
    integer, intent(in) :: nElems

    !> number of elements in state array
    integer, intent(in) :: nSize

    !> current level
    integer, intent(in) :: iLevel

    !> stencil header contains discrete velocity vectors
    type(tem_stencilHeader_type), intent(in) :: stencil

    !> variable system which is required to access fieldProp
    !! information via variable method data c_ptr
    type(tem_varSys_type), intent(in) :: varSys

    !> Output of this routine
    !! Size: nElems*nAuxScalars
    real(kind=rk), intent(inout) :: auxField(:)
    ! -------------------------------------------------------------------- !
    integer :: dens_pos, vel_pos(3), pdfPos
    integer :: iElem, iDir, elemOff
    real(kind=rk) :: pdf( stencil%QQ )
    ! -------------------------------------------------------------------------- !
    dens_pos = varSys%method%val(derVarPos%density)%auxField_varPos(1)
    vel_pos = varSys%method%val(derVarPos%velocity)%auxField_varPos(1:3)

    !NEC$ ivdep
    do iElem = 1, nElems
      !NEC$ shortloop
      do iDir = 1, stencil%QQ
        pdfPos = varSys%method%val(derVarPos%pdf)%state_varPos(iDir)
        pdf(iDir) = state( ?IDX?(pdfPos, iElem, varSys%nScalars, nSize) )
      end do

      ! element offset
      elemoff = (iElem-1)*varSys%nAuxScalars

      ! density
      auxField(elemOff+dens_pos) = sum( pdf )

      ! velocity
      auxField(elemOff+vel_pos(1)) = sum( pdf * stencil%cxDirRK(1,:) ) * rho0Inv
      auxField(elemOff+vel_pos(2)) = sum( pdf * stencil%cxDirRK(2,:) ) * rho0Inv
      auxField(elemOff+vel_pos(3)) = sum( pdf * stencil%cxDirRK(3,:) ) * rho0Inv
    end do

  end subroutine deriveAuxIncomp_fromState
! ************************************************************************** !

! ************************************************************************** !
  !> This routine computes equilibirium from state array
  !! Here it is assumed that rho0 = 1.0
  !!
  !! This subroutine's interface must match the abstract interface definition
  !! [[derive_FromState]] in derived/[[mus_derVarPos_module]].f90 in order to be
  !! callable via [[mus_derVarPos_type:velFromState]],
  !! [[mus_derVarPos_type:equilFromState]],
  !! [[mus_derVarPos_type:momFromState]],
  !! [[mus_derVarPos_type:velocitiesFromState]], and
  !! [[mus_derVarPos_type:momentaFromState]] function pointers.
  subroutine deriveEqIncomp_FromState( state, iField, nElems, varSys, layout, &
    &                                  res                                    )
    ! -------------------------------------------------------------------- !
    !> Array of state
    !! n * layout%fStencil%QQ * nFields
    real(kind=rk), intent(in) :: state(:)

    !> Current field
    integer, intent(in) :: iField

    !> number of elements
    integer, intent(in) :: nElems

    !> variable system which is required to access fieldProp
    !! information via variable method data c_ptr
    type(tem_varSys_type), intent(in) :: varSys

    !> scheme layout contains stencil definition and lattice weights
    type(mus_scheme_layout_type), intent(in) :: layout

    !> Output of this routine
    !! Dimension: n * nComponents of res
    real(kind=rk), intent(out) :: res(:)
    ! -------------------------------------------------------------------- !
    integer :: QQ, iElem, iDir
    real(kind=rk) :: f( layout%fStencil%QQ )
    real(kind=rk) :: rho, vel(3), usq, ucx

    QQ = layout%fStencil%QQ

    do iElem = 1, nElems

      do iDir = 1, QQ
        f(iDir) = state( iDir+(iElem-1)*varSys%nScalars )
      end do

      rho = sum( f )
      vel(1) = sum( f * layout%fStencil%cxDirRK(1,:) )
      vel(2) = sum( f * layout%fStencil%cxDirRK(2,:) )
      vel(3) = sum( f * layout%fStencil%cxDirRK(3,:) )

      usq = dot_product(vel, vel) * t2cs2inv

      do iDir = 1, QQ
        ucx = dot_product( layout%fStencil%cxDirRK(:, iDir), vel )

        ! calculate equilibrium density
        res( (iElem-1)*QQ+iDir ) = layout%weight( iDir ) * &
          &  ( rho + rho0 * ( cs2inv * ucx + ucx * ucx * t2cs4inv - usq ) )
      enddo

    end do

  end subroutine deriveEqIncomp_FromState
! ************************************************************************** !


! ************************************************************************** !
  !> Calculate the velocity on given elements for incompressible model
  !! \( \vec{u} = \frac{\sum(f_i c_i)}{\rho_0} \)
  !!
  !! The interface has to comply to the abstract interface
  !! [[tem_varSys_module:tem_varSys_proc_element]].
  !!
  recursive subroutine mus_deriveVelocityIncomp(fun, varsys, stencil, iLevel, &
    &                                             posInState, pdf, res, nVals )
    !> description of the method to obtain the variables, here some preset
    !! values might be stored, like the space time function to use or the
    !! required variables.
    class(tem_varsys_op_type), intent(in)     :: fun
    !> the variable system to obtain the variable from.
    type(tem_varsys_type), intent(in)         :: varsys
    !> fluid stencil defintion
    type(tem_stencilHeader_type), intent(in)  :: stencil
    !> current Level
    integer, intent(in)                       :: iLevel
    !> Position of element in levelwise state array
    integer, intent(in)                       :: posInState(:)
    !> pdf array
    real(kind=rk), intent(in)                 :: pdf(:)
    !> results
    real(kind=rk), intent(out)                :: res(:)
    !> nVals to get
    integer, intent(in)                       :: nVals
    ! ---------------------------------------------------------------------- !
    type(mus_varSys_data_type), pointer       :: fPtr
    type(mus_scheme_type), pointer            :: scheme
    real(kind=rk), allocatable                :: tmpPDF(:)
    integer                                   :: iComp, iVal, pdfPos, nCompsPDF
    ! ---------------------------------------------------------------------- !
    call C_F_POINTER( fun%method_Data, fPtr )
    scheme => fPtr%solverData%scheme

    pdfPos = fun%input_varPos(1)
    nCompsPDF = varSys%method%val( pdfPos )%nComponents
    allocate( tmpPDF( nCompsPDF ) )
    res = 0.0_rk

    do iVal = 1, nVals
      tmpPDF = pdf( (iVal-1)*nCompsPDF+1 : iVal*nCompsPDF )
      do iComp = 1, fun%nComponents
        res( iComp+ (iVal-1)*fun%nComponents) =                 &
          &  sum(tmpPDF * scheme%layout%fStencil%cxDirRK(iComp,:))
      end do ! iComp
    end do ! iVal
    deallocate( tmpPDF )

  end subroutine mus_deriveVelocityIncomp
! ************************************************************************** !


! ************************************************************************** !
  !> Calculate the equlibrium of a given element number with the given input
  !! state vector for incompressible model
  !! \( \rho0 = 1.0 \)
  !! \( f_eq = w_\alpha*(\rho + \rho0( (u . c_i) + (u . c_i)^2 - u^2 ) ) \)
  !!
  !! The interface has to comply to the abstract interface
  !! [[tem_varSys_module:tem_varSys_proc_element]].
  !!
  recursive subroutine mus_deriveEquilIncomp(fun, varsys, stencil, iLevel, &
    &                                          posInState, pdf, res, nVals )
    !> description of the method to obtain the variables, here some preset
    !! values might be stored, like the space time function to use or the
    !! required variables.
    class(tem_varsys_op_type), intent(in)     :: fun
    !> the variable system to obtain the variable from.
    type(tem_varsys_type), intent(in)         :: varsys
    !> fluid stencil defintion
    type(tem_stencilHeader_type), intent(in)  :: stencil
    !> current Level
    integer, intent(in)                       :: iLevel
    !> Position of element in levelwise state array
    integer, intent(in)                       :: posInState(:)
    !> pdf array
    real(kind=rk), intent(in)                 :: pdf(:)
    !> results
    real(kind=rk), intent(out)                :: res(:)
    !> nVals to get
    integer, intent(in)                       :: nVals
    ! ---------------------------------------------------------------------------
    type(mus_varSys_data_type), pointer       :: fPtr
    type(mus_scheme_type), pointer            :: scheme
    real(kind=rk), allocatable                :: tmpPDF(:)
    real(kind=rk), allocatable                :: fEq(:)
    real(kind=rk)                             :: dens, vel(3)
    integer                                   :: pdfPos, nCompsPDF, iVal, iDir
    ! ---------------------------------------------------------------------------
    call C_F_POINTER( fun%method_Data, fPtr )
    scheme => fPtr%solverData%scheme
    pdfPos = fun%input_varPos(1)
    nCompsPDF = varSys%method%val( pdfPos )%nComponents
    allocate( tmpPDF( nCompsPDF ) )
    allocate( fEq( fun%nComponents ) )
    res = 0.0_rk

    do iVal = 1, nVals
      tmpPDF = pdf( (iVal-1)*nCompsPDF+1 : iVal*nCompsPDF )
      dens   = sum(tmpPDF)
      vel(1) = sum(tmpPDF * scheme%layout%fStencil%cxDirRK(1,:))
      vel(2) = sum(tmpPDF * scheme%layout%fStencil%cxDirRK(2,:))
      vel(3) = sum(tmpPDF * scheme%layout%fStencil%cxDirRK(3,:))

      do iDir = 1, scheme%layout%fStencil%QQ
?? copy :: EqDir_incomp(fEq(iDir), dens, vel, scheme%layout, iDir, rho0, cs2inv, rk)
      end do
      res( (iVal-1)*fun%nComponents+1: iVal*fun%nComponents ) = fEq
    end do !iVal
    deallocate( tmpPDF )
    deallocate( fEq )

  end subroutine mus_deriveEquilIncomp
! ************************************************************************** !


! ************************************************************************** !
  !> Calculate the Non-Equlibrium
  !!
  recursive subroutine mus_deriveNonEquilIncomp(fun, varsys, stencil, iLevel, &
    &                                             posInState, pdf, res, nVals )
    !> description of the method to obtain the variables, here some preset
    !! values might be stored, like the space time function to use or the
    !! required variables.
    class(tem_varsys_op_type), intent(in)     :: fun
    !> the variable system to obtain the variable from.
    type(tem_varsys_type), intent(in)         :: varsys
    !> fluid stencil defintion
    type(tem_stencilHeader_type), intent(in)  :: stencil
    !> current Level
    integer, intent(in)                       :: iLevel
    !> Position of element in levelwise state array
    integer, intent(in)                       :: posInState(:)
    !> pdf array
    real(kind=rk), intent(in)                 :: pdf(:)
    !> results
    real(kind=rk), intent(out)                :: res(:)
    !> nVals to get
    integer, intent(in)                       :: nVals
    ! ---------------------------------------------------------------------------
    type(mus_varSys_data_type), pointer       :: fPtr
    type(mus_scheme_type), pointer            :: scheme
    real(kind=rk), allocatable                :: tmpPDF(:)
    real(kind=rk), allocatable                :: fEq(:)
    real(kind=rk)                             :: dens, vel(3)
    integer                                   :: iDir, iVal
    integer                                   :: pdfPos, nCompsPDF
    ! ---------------------------------------------------------------------------
    call C_F_POINTER( fun%method_Data, fPtr )
    scheme => fPtr%solverData%scheme

    pdfPos = fun%input_varPos(1)
    nCompsPDF = varSys%method%val( pdfPos )%nComponents
    allocate( fEq( fun%nComponents ) )
    allocate( tmpPDF( nCompsPDF ) )
    res = 0.0_rk

    do iVal = 1 , nVals
      tmpPDF = pdf( (iVal-1)*nCompsPDF+1 : iVal*nCompsPDF )
      ! computes density and velocity
      dens   = sum(tmpPDF)
      vel(1) = sum(tmpPDF * scheme%layout%fStencil%cxDirRK(1,:))
      vel(2) = sum(tmpPDF * scheme%layout%fStencil%cxDirRK(2,:))
      vel(3) = sum(tmpPDF * scheme%layout%fStencil%cxDirRK(3,:))

      ! computes equilibrium
      do iDir = 1, scheme%layout%fStencil%QQ
?? copy :: EqDir_incomp(fEq(iDir), dens, vel, scheme%layout, iDir, rho0, cs2inv, rk)
      end do

      res( (iVal-1)*fun%nComponents+1: iVal*fun%nComponents ) = tmpPDF - fEq
    end do ! iVal

    deallocate( tmpPDF )
    deallocate( fEq )


  end subroutine mus_deriveNonEquilIncomp
! ************************************************************************** !


! ************************************************************************** !
  !> Calculate the kinetic energy as
  !! \[
  !!  E = \frac{1}{2} \rho ( u_x^2 + u_y^2 + u_z^2 )
  !! \]
  !!
  !! The interface has to comply to the abstract interface
  !! [[tem_varSys_module:tem_varSys_proc_element]].
  !!
  recursive subroutine mus_deriveKeIncomp(fun, varsys, stencil, iLevel, &
    &                                       posInState, pdf, res, nVals )
    !> description of the method to obtain the variables, here some preset
    !! values might be stored, like the space time function to use or the
    !! required variables.
    class(tem_varsys_op_type), intent(in)     :: fun
    !> the variable system to obtain the variable from.
    type(tem_varsys_type), intent(in)         :: varsys
    !> fluid stencil defintion
    type(tem_stencilHeader_type), intent(in)  :: stencil
    !> current Level
    integer, intent(in)                       :: iLevel
    !> Position of element in levelwise state array
    integer, intent(in)                       :: posInState(:)
    !> pdf array
    real(kind=rk), intent(in)                 :: pdf(:)
    !> results
    real(kind=rk), intent(out)                :: res(:)
    !> nVals to get
    integer, intent(in)                       :: nVals
    ! ---------------------------------------------------------------------------
    type(mus_varSys_data_type), pointer       :: fPtr
    type(mus_scheme_type), pointer            :: scheme
    real(kind=rk), allocatable                :: tmpPDF(:)
    real(kind=rk)                             :: vel(3)
    integer                                   :: pdfPos, nCompsPDF, iVal
    ! ---------------------------------------------------------------------------
    call C_F_POINTER( fun%method_Data, fPtr )
    scheme => fPtr%solverData%scheme
    pdfPos = fun%input_varPos(1)
    nCompsPDF = varSys%method%val( pdfPos )%nComponents
    allocate( tmpPDF( nCompsPDF ) )
    res = 0.0_rk

    do iVal = 1, nVals
      tmpPDF = pdf( (iVal-1)*nCompsPDF+1 : iVal*nCompsPDF )
      vel(1) = sum(tmpPDF * scheme%layout%fStencil%cxDirRK(1,:))
      vel(2) = sum(tmpPDF * scheme%layout%fStencil%cxDirRK(2,:))
      vel(3) = sum(tmpPDF * scheme%layout%fStencil%cxDirRK(3,:))
      res( iVal ) = sum( vel(:)*vel(:) ) * 0.5_rk
    end do !iVal

    deallocate( tmpPDF )

  end subroutine mus_deriveKeIncomp
! ************************************************************************** !

! ************************************************************************** !
  !> author: Jiaxing Qi
  !! Calculate the strain rate ( or rate of strain, or rate of deformation)
  !!
  !! The formula is:
  !! \[
  !!  \tau_{\alpha \beta}=
  !!    -\frac{3\omega}{2{\rho}_0} \sum_{i} f^{neq}_{i} c_{i\alpha} c_{i\beta}
  !! \]
  !! where \( \tau_{\alpha \beta}\) is the stress
  !! in the \(\beta\)-direction on a face normal to the \(\alpha\)-axis,\n
  !! \( f^{neq}_i = f_i - f^{eq}_i\) is the non-equilibrium pdf.\n
  !! For more information, please refer to: equation 45 in\n
  !! Krueger T, Varnik F, Raabe D. Shear stress in lattice Boltzmann
  !! simulations. Physical Review E. 2009;79(4):1-14.\n
  !!
  !! For multi-level mesh, Omega on finer level needs to be adjusted in order to
  !! get the correct shearstress calculation.\n
  !! First, we defines c as the dx ratio between finer and coarse level.\n
  !! \( c={ \Delta dx }_{ c }/{ \Delta dx }_{ f } \)
  !! Then the viscosity on the different levels must satisfy:\n
  !! \( \frac { { \nu  }_{ f } }{ { \nu  }_{ c } } =c \)
  !! This constrain leads to a relationship of omega on different levels:\n
  !! \( {\omega}_f = \frac {1}{ {\lambda}(\frac{1}{{\omega}_c}-0.5)+0.5 } \)
  !! For more information, please refer to:\n
  !! Manuel H, Harald K, Joerg B, Sabine R. Aeroacoustic validation of the
  !! lattice boltzmann method on non-uniform grids. ECCOMAS 2012
  !!
  recursive subroutine mus_deriveStrainRateIncomp(fun, varsys, stencil, iLevel,&
      &                                            posInState, pdf, res, nVals )
    !> description of the method to obtain the variables, here some preset
    !! values might be stored, like the space time function to use or the
    !! required variables.
    class(tem_varsys_op_type), intent(in)     :: fun
    !> the variable system to obtain the variable from.
    type(tem_varsys_type), intent(in)         :: varsys
    !> fluid stencil defintion
    type(tem_stencilHeader_type), intent(in)  :: stencil
    !> current Level
    integer, intent(in)                       :: iLevel
    !> Position of element in levelwise state array
    integer, intent(in)                       :: posInState(:)
    !> pdf array
    real(kind=rk), intent(in)                 :: pdf(:)
    !> results
    real(kind=rk), intent(out)                :: res(:)
    !> nVals to get
    integer, intent(in)                       :: nVals
    ! ---------------------------------------------------------------------------
    type(mus_varSys_data_type), pointer :: fPtr
    type(mus_scheme_type), pointer :: scheme
    integer :: nCompsPDF, iVal, pdfPos, iDir
    real(kind=rk) :: dens, vel(3), omega
    real(kind=rk), allocatable :: nonEq(:)
    real(kind=rk), allocatable :: tmpPDF(:)
    real(kind=rk), allocatable :: fEq(:)
    real(kind=rk), allocatable :: tau(:)
    integer :: nScalars
    ! ---------------------------------------------------------------------- !
    call C_F_POINTER( fun%method_Data, fPtr )
    scheme => fPtr%solverData%scheme

    pdfPos = fun%input_varPos(1)
    nCompsPDF = varSys%method%val( pdfPos )%nComponents
    nScalars = varSys%nScalars

    allocate( tmpPDF( nCompsPDF ) )
    allocate( fEq( nCompsPDF ) )
    allocate( nonEq( nCompsPDF ) )
    allocate( tau( fun%nComponents ) )

    do iVal = 1, nVals
      tmpPDF = pdf( (iVal-1)*nCompsPDF+1 : iVal*nCompsPDF )

      ! computes density and velocity
      dens   = sum(tmpPDF)
      vel(1) = sum(tmpPDF * scheme%layout%fStencil%cxDirRK(1,:))
      vel(2) = sum(tmpPDF * scheme%layout%fStencil%cxDirRK(2,:))
      vel(3) = sum(tmpPDF * scheme%layout%fStencil%cxDirRK(3,:))

      ! computes equilibrium
      do iDir = 1, scheme%layout%fStencil%QQ
?? copy :: EqDir_incomp(fEq(iDir), dens, vel, scheme%layout, iDir, rho0, cs2inv, rk)
      end do

      ! Non-Eq
      nonEq = tmpPDF - fEq

      ! get the correct omega value
      omega = scheme%field(1)%fieldProp%fluid%viscKine%omLvl( iLevel ) &
        &                                    %val(posInState(iVal))

      ! compute shear stress
      tau(:) = secondMom( cxcx = scheme%layout%fStencil%cxcx(:,:), &
        &                 f    = nonEq(:),                         &
        &                 QQ   = scheme%layout%fStencil%QQ         )

      res( (iVal-1)*fun%nComponents+1: iVal*fun%nComponents ) = &
        &                   tau(:) * (-1.5_rk) * omega / rho0

    end do ! iVal

    deallocate( tmpPDF )
    deallocate( fEq )
    deallocate( nonEq )
    deallocate( tau )

  end subroutine mus_deriveStrainRateIncomp
! ****************************************************************************** !

end module mus_derQuanIncomp_module
! ************************************************************************** !
