! Copyright (c) 2013, 2016, 2020 Harald Klimach <harald.klimach@uni-siegen.de>
! Copyright (c) 2013-2014 Kartik Jain <kartik.jain@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) 2015-2016 Tobias Schneider <tobias1.schneider@student.uni-siegen.de>
! Copyright (c) 2016 Verena Krupp <verena.krupp@uni-siegen.de>
! Copyright (c) 2016 Philipp Otte <otte@mathcces.rwth-aachen.de>
! Copyright (c) 2016-2017 Raphael Haupt <raphael.haupt@uni-siegen.de>
! Copyright (c) 2019-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.
!!
!! 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 'treelm/source/deriveMacros.inc'
module mus_derQuan_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, &
    &                                 sqrt3, cs2inv, cs2, t2cs2inv, t2cs4inv,  &
    &                                 cs4inv, q000, rho0
  use env_module,               only: rk, long_k, labelLen
  use tem_float_module,         only: operator(.feq.), operator(.fge.), &
    &                                 operator(.fle.)
  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_logging_module,       only: logUnit
  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_aux_module,           only: tem_abort
  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_property_module,      only: prp_hasBnd, prp_hasQval
  use tem_tools_module,         only: tem_PositionInSorted
  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_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, getShearRate
  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_operation_var_module,      only: mus_opVar_setupIndices,         &
    &                                      mus_opVar_gradU_forElement,     &
    &                                      mus_opVar_vorticity_forElement, &
    &                                      mus_opVar_QCriterion_forElement
  use mus_stateVar_module,           only: mus_accessVar_setupIndices
  use mus_derVarPos_module,          only: mus_derVarPos_type
  use mus_physics_module,            only: mus_convertFac_type

  implicit none

  private

  public :: mus_append_derVar_fluid

  public :: derivePressure, derivePressure_fromIndex
  public :: deriveDensity
  public :: deriveDensity_fromIndex
  public :: deriveShearStress
  public :: deriveShearMag
  public :: deriveWSS3D
  public :: deriveWSS2D
  public :: deriveTemp
  public :: deriveShearRate
  public :: deriveBndForce
  ! equilbrium from macro uses different interface defined in
  ! mus_variable_module
  public :: deriveEquil_FromMacro
  public :: deriveRho_fromState
  public :: deriveVel_fromState
  public :: deriveVel_FromPreColState
  public :: deriveEq_fromState
  public :: deriveEquil_fromAux
  public :: deriveAux_fromState

  ! source variable
  public :: derive_absorbLayer
  public :: derive_force
  public :: derive_force1stOrd

  ! source update
  public :: applySrc_absorbLayer
  public :: applySrc_force
  public :: applySrc_force_MRT
  public :: applySrc_force_MRT_d3q19
  public :: applySrc_force1stOrd

contains

  ! **************************************************************************** !
  !> subroutine to add derive variables for weakly compressible LBM
  !! (schemekind = 'fluid') to the varsys.
  subroutine mus_append_derVar_fluid( 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      ', 'moment         ', &
      &                   '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
      method_data  = mus_get_new_solver_ptr(solverData)
      set_params => tem_varSys_setParams_dummy
      get_params => tem_varSys_getParams_dummy

      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 => deriveEquil
        get_valOfIndex => deriveEquil_fromIndex
        nComponents = stencil%QQ
        allocate(input_varname(1))
        input_varname(1) = 'pdf'

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

      case ('kinetic_energy')
        get_element => deriveKE
        get_valOfIndex => deriveKe_fromIndex
        nComponents = 1
        allocate(input_varname(2))
        input_varname(1) = 'density'
        input_varname(2) = '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 => deriveStrainRate
        get_valOfIndex => deriveStrainRate_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 => deriveMomentum
        get_valOfIndex => deriveMomentum_fromIndex
        nComponents = 3
        allocate(input_varname(2))
        input_varname(1) = 'density'
        input_varname(2) = '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 ('moment') !ONLY FOR 2D
        get_element => deriveMoment
        nComponents = 9
        allocate(input_varname(1))
        input_varname(1) = 'pdf'

      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_fluid
  ! **************************************************************************** !


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

! ****************************************************************************** !
  !> Initiates the calculation of density
  !! This routine sets the function Pointer for density 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(deriveDensity)
    ! ---------------------------------------------------------------------------
    !> Function pointer to perform specific operation.
    procedure(mus_derive_fromPDF), pointer :: fnCalcPtr
    ! ---------------------------------------------------------------------------

    fnCalcPtr => mus_derivedensity

    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 deriveDensity
! ****************************************************************************** !


! ****************************************************************************** !
  !> Initiates the calculation of velocity
  !! This routine sets the function Pointer for velocity 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(deriveVelocity)
    ! ---------------------------------------------------------------------------
    !> Function pointer to perform specific operation.
    procedure(mus_derive_fromPDF), pointer :: fnCalcPtr
    ! ---------------------------------------------------------------------------
    fnCalcPtr => mus_deriveVelocity

    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 deriveVelocity
! ****************************************************************************** !


! ****************************************************************************** !
  !> Calculate the temperature of a given set of elements (sum up all links).
  !!
  !! The interface has to comply to the abstract interface
  !! [[tem_varSys_module:tem_varSys_proc_element]].
  !!
?? copy :: get_element_headtxt(deriveTemp)
    ! ---------------------------------------------------------------------------

    res( 1:nElems ) = div1_3

  end subroutine deriveTemp
! ****************************************************************************** !


! ****************************************************************************** !
  !> Calculate the pressure of a given set of elements (sum up all links).
  !!
  !! Pressure calculation according to the isentropic equation of state for
  !! the LBM \( p = \rho c_s^2 \)
  !! with the calculation of density as in deriveDensity
  !!
  !! The interface has to comply to the abstract interface
  !! [[tem_varSys_module:tem_varSys_proc_element]].
  !!
?? copy :: get_element_headtxt(derivePressure)
    ! ---------------------------------------------------------------------------
    integer :: dens_pos
    ! ---------------------------------------------------------------------------
    ! position of density in glob system
    dens_pos = fun%input_varPos(1)
    ! get density variable from auxField
    call varSys%method%val(dens_pos)%get_element( varSys  = varSys,  &
      &                                           elemPos = elemPos, &
      &                                           time    = time,    &
      &                                           tree    = tree,    &
      &                                           nElems  = nElems,  &
      &                                           nDofs   = nDofs,   &
      &                                           res     = res      )

    ! convert density to pressure
    res(1:nElems) = ( res(1:nElems) - rho0 ) * cs2

  end subroutine derivePressure
! ****************************************************************************** !


! ****************************************************************************** !
  !> Calculate the force on the boundary of a given set of elements
  !!
  !! Force calculation in LBM:
  !! \[ F = \sum_\alpha{\alpha \ne 0} \vec{e}_{\alphabar} (2f_\alphabar
  !!       + 6\omega_\alpha \vec{e}_\alpha \rho \vec{u_b})  \]
  !! \n where, \n
  !! \alpha = direction towards boundary
  !! \alphabar = direction opposite to boundary
  !! \vec{u_b} = velocity of the obstacle cell.
  !! For stationary object, this should be zero.
  !!
  !! The interface has to comply to the abstract interface
  !! [[tem_varSys_module:tem_varSys_proc_element]].
  !!
?? copy :: get_element_headtxt(deriveBndForce)
    ! ---------------------------------------------------------------------------
    integer :: statePos, iElem, iDir, iLevel, QQN, QQ, bcID
    type(mus_varSys_data_type), pointer :: fPtr
    type(mus_scheme_type), pointer :: scheme
    integer :: pdfPos, invDir, posInBndID, posInQVal
    real(kind=rk) :: force(3), fIn, fOut, fOutNeigh, fIntp, qVal
    real(kind=rk), allocatable :: f(:)
    integer, allocatable :: iMeshDir(:), varPos(:)
    integer :: nSize
    ! ---------------------------------------------------------------------------
    call C_F_POINTER( fun%method_Data, fPtr )
    scheme => fPtr%solverData%scheme

    pdfPos = fun%input_varPos(1)
    QQ = scheme%layout%fStencil%QQ
    QQN    = scheme%layout%fStencil%QQN
    res = 0.0_rk

    allocate(        f( QQN ) )
    allocate( iMeshDir( QQN ) )
    allocate(   varPos( QQN ) )

    do iDir = 1, QQN
      iMeshDir(iDir) = scheme%layout%fStencil%map(iDir)
      varPos(iDir)   = varSys%method%val(pdfPos)%state_varPos(iDir)
    end do

    ! compute bnd_force only if shape is defined as boundary kind
    do iElem = 1, nElems

! write(*,*) ''
! write(*,*) 'treeID: ', tree%treeID( elemPos(iElem) )
      force = 0.0_rk

      ! if state array is defined level wise then use levelPointer(pos)
      ! to access state array
      statePos = fPtr%solverData%geometry%levelPointer( elemPos(iElem) )
      iLevel   = tem_levelOf( tree%treeID( elemPos(iElem) ) )
      nSize = scheme%pdf( iLevel )%nSize
      ! get pdf
      do iDir = 1, QQN
        f( iDir ) =  &
          & scheme%state( iLevel )%val(                              &
          & ?IDX?( varPos(iDir), statePos, varSys%nScalars, nSize ), &
          & scheme%pdf( iLevel )%nNext )
      end do !iComp

      ! compute force only if current element has boundary property
      if ( btest( tree%elemPropertyBits( elemPos(iElem) ), &
        &         prp_hasBnd )  ) then

        ! position of current element in boundary_ID list
        posInBndID = fPtr%solverData%geometry%posInBndID( elempos(iElem) )

        if (posInBndID > 0) then
          do iDir = 1, QQN

            ! direction away from boundary
            invDir = scheme%layout%fStencil%cxDirInv( iDir )

             ! compute force only if current direction has boundary
            bcID = int( fPtr%solverData%geometry%boundary              &
              &             %boundary_ID( iMeshDir(iDir), posInBndID ) )

            if (bcID > 0) then

              ! element has qVal
              if ( btest( tree%elemPropertyBits( elemPos(iElem) ), &
                &         prp_hasQval ) ) then
                posInQVal  = fPtr%solverData%geometry%posInQVal(elemPos(iElem))
                if (posInQVal > 0) then
                  qVal = fPtr%solverData%geometry%boundary    &
                    &        %qVal( iMeshDir(iDir), posInQVal )
                else
                  qVal = 0.5_rk
                end if
              else
                qVal = 0.5_rk
              end if ! element has qVal

              ! if no intersection in certain direction use default qVal=0.5
              ! this happens only if element has more than one boundary and
              ! one has qVal and other not
              if (qVal <= 0) qVal = 0.5_rk

              ! @todo: can be done in the same way as bc routine
              if ( qVal >= 0.5_rk ) then
                fOut  = f(   iDir )
                fIn   = f( invDir )
                fIntp = (1._rk - 0.5_rk/qVal) * fIn + (0.5_rk/qVal) * fOut
              else
                fOutNeigh = scheme%state( iLevel )%val(                      &
&?FETCH?(iDir,1,statePos,QQ,varSys%nScalars,nSize,scheme%pdf(iLevel)%neigh), &
          & scheme%pdf( iLevel )%nNext )
                fOut  = f( iDir )
                fIntp = (1._rk - 2._rk*qVal) * fOutNeigh + (2._rk*qVal) * fOut
              end if ! qVal >= or < 0.5_rk
              force = force + scheme%layout%fStencil%cxDirRK(:,iDir) &
                &           * ( fOut + fIntp )

            end if ! bcID>0
          end do !iDir
        end if ! posInBndID > 0

      end if ! element has_boundary

      res( (iElem-1)*fun%nComponents+1: iElem*fun%nComponents ) = force

    end do ! iElem

    deallocate( iMeshDir )
    deallocate( f )
    deallocate( varPos )

  end subroutine deriveBndForce
! ****************************************************************************** !


! ****************************************************************************** !
  !> Calculate momentum from density and velocity stored in auxField
  !!
  !! The interface has to comply to the abstract interface
  !! [[tem_varSys_module:tem_varSys_proc_element]].
  !!
?? copy :: get_element_headtxt(deriveMomentum)
    ! ---------------------------------------------------------------------------
    !> 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 density in auxField array
        dens_pos = varSys%method%val( fun%input_varPos(1) ) &
          &                     %auxField_varPos(1)

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

        ! density
        rho = auxField(iLevel)%val( elemOff + dens_pos )

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

      end do ! iElem
    end associate

  end subroutine deriveMomentum
! ****************************************************************************** !


! ****************************************************************************** !
  !> 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(deriveEquil)
    ! ---------------------------------------------------------------------------
    !> Function pointer to perform specific operation.
    procedure(mus_derive_fromPDF), pointer :: fnCalcPtr
    ! ---------------------------------------------------------------------------
    fnCalcPtr => mus_deriveEquil

    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 deriveEquil
! ****************************************************************************** !


! ****************************************************************************** !
  !> 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(deriveNonEquil)
    ! ---------------------------------------------------------------------------
    !> Function pointer to perform specific operation.
    procedure(mus_derive_fromPDF), pointer :: fnCalcPtr
    ! ---------------------------------------------------------------------------
    fnCalcPtr => mus_deriveNonEquil

    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 deriveNonEquil
! ****************************************************************************** !


! ****************************************************************************** !
  !> 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(deriveStrainRate)
    ! ---------------------------------------------------------------------------
    !> Function pointer to perform specific operation.
    procedure(mus_derive_fromPDF), pointer :: fnCalcPtr
    ! ---------------------------------------------------------------------------
    fnCalcPtr => mus_deriveStrainRate

    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 deriveStrainRate
! ****************************************************************************** !


! ****************************************************************************** !
  !> author: Jiaxing Qi
  !! Calculate the shear rate
  !!
  !! The Shear Rate is defined as
  !! \[
  !!  \dot{\gamma} = \sqrt{ 2D_{II} }
  !! \]
  !! where \( D_{II} \) is the second invariant of the strain rate tensor and
  !! defined as
  !! \[
  !!    D_{II} = \sum^{l}_{\alpha,\beta=l} S_{\alpha\beta} S_{\alpha\beta}
  !! \]
  !! where \( S_{\alpha\beta} \) is the strain rate tensor.
  !!
?? copy :: get_element_headtxt(deriveShearRate)
    ! ---------------------------------------------------------------------------
    integer :: posStrain, iElem
    real(kind=rk) :: strain(6*nElems)
    ! ---------------------------------------------------------------------------

    posStrain = fun%input_varPos(1)

    call varSys%method%val(posStrain)%get_element( varSys  = varSys,  &
      &                                            elemPos = elemPos, &
      &                                            time    = time,    &
      &                                            tree    = tree,    &
      &                                            nElems  = nElems,  &
      &                                            nDofs   = nDofs,   &
      &                                            res     = strain   )

    do iElem = 1, nElems

      res( iElem ) = getShearRate( strain( (iElem-1)*6+1:iElem*6 ) )

    end do ! iElem

  end subroutine deriveShearRate
! ****************************************************************************** !

! ****************************************************************************** !
  !> author: Jiaxing Qi
  !! Calculate the deviatoric shear stress for Newtonian fluid
  !! (exclude pressure) (no mixtures).\n
  !! Shear Stress depends on variable: nonEquilibirium
  !!
  !! The formula is
  !! \[
  !!  \tau_{\alpha \beta}=
  !!    -(1-\frac{\omega}{2}) \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-equilibirium density.\n
  !! For more information, please refer to:\n
  !! Krueger T, Varnik F, Raabe D. Shear stress in lattice Boltzmann
  !! simulations. Physical Review E. 2009;79(4):1-14.\n
  !! Thus this variable is dependent on another variable: nonEquilibirium.
  !!
  !! 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
  !!
?? copy :: get_element_headtxt(deriveShearStress)
    ! ---------------------------------------------------------------------------
    type(mus_varSys_data_type), pointer :: fPtr
    type(mus_scheme_type), pointer :: scheme
    integer :: posNonEq, nCompsNonEq, iLevel, iElem
    real(kind=rk) :: omega
    real(kind=rk), allocatable :: nonEq(:)
    real(kind=rk), allocatable :: tau(:)
    ! ---------------------------------------------------------------------------
    call C_F_POINTER( fun%method_Data, fPtr )
    scheme => fPtr%solverData%scheme

    posNonEq = fun%input_varPos(1)
    nCompsNonEq = varSys%method%val( posNonEq )%nComponents
    allocate( nonEq( nCompsNonEq * nElems ) )
    allocate(   tau( fun%nComponents ) )

    ! calculate nonEq for all elements
    call varSys%method%val(posNonEq)%get_element(  varSys  = varSys,  &
      &                                            elemPos = elemPos, &
      &                                            time    = time,    &
      &                                            tree    = tree,    &
      &                                            nElems  = nElems,  &
      &                                            nDofs   = nDofs,   &
      &                                            res     = nonEq    )

    do iElem = 1, nElems

      iLevel = tem_levelOf( tree%treeID( elemPos(iElem) ) )
      ! get the correct omega value
      omega = scheme%field(1)%fieldProp%fluid%viscKine%omLvl(iLevel)%val(iElem)

      ! compute second moment of nonEq
      tau(:) = secondMom( cxcx = scheme%layout%fStencil%cxcx(:,:), &
        &                 f    = nonEq( (iElem-1)*nCompsNonEq+1    &
        &                              :iElem*nCompsNonEq ),       &
        &                 QQ   = scheme%layout%fStencil%QQ         )

      res( (iElem-1)*fun%nComponents+1: iElem*fun%nComponents ) = &
        &                   tau(:) * ( 0.5_rk * omega - 1._rk )

    end do ! iElem

    deallocate( nonEq )
    deallocate( tau )

  end subroutine deriveShearStress
! ****************************************************************************** !


! ****************************************************************************** !
  !> author: Jiaxing Qi
  !! Calculate the wall shear stress (WSS) of a given element with the given
  !! input
  !!
  !! The wall shear stress is shear stress exerted on the wall. Since it is well
  !! known that shear stress gets its maximum value on the wall, here we can
  !! directly calculate the principle stress of the stress tensor instead of
  !! multiplying the stress tensor by normal vector of the plane.
  !! The principle stresses are just the eigenvalues of the stress tensor. To
  !! get those eigenvalues, we need solve the characteristic equation:\n
  !! \( {\tau}^3 + a_{2}{\tau}^2 + a_1{\tau} + a_0 = 0 \)
  !! where \( a_2 = \tau_{ii} \)
  !! \( a_1 = (\tau_{ii}\tau_{jj} - \tau_{ij}\tau_{ji}) / 2 \)
  !! \( a_0 = det(\tau_{ij}) \). Here Einstein notation is used.
  !!
?? copy :: get_element_headtxt(deriveWSS3D)
    ! ---------------------------------------------------------------------------
    ! local variables
    integer :: iElem, shear_pos
    real(kind=rk) :: shear(6*nElems) ! shear stress of all
    ! ---------------------------------------------------------------------------

    ! wall shear stress is depend on shear stress. get its position in globSys
    shear_pos = fun%input_varPos(1)

    ! first compute shear stress, save it to array tau
    ! the shear tensor has the following order: x, y, z, xy, xz, yz
    call varSys%method%val(shear_pos)%get_element( varSys  = varSys,  &
      &                                            elemPos = elemPos, &
      &                                            time    = time,    &
      &                                            tree    = tree,    &
      &                                            nElems  = nElems,  &
      &                                            nDofs   = nDofs,   &
      &                                            res     = shear    )

    do iElem = 1, nElems

      res( iElem ) = getWSS( shear( (iElem-1)*6+1 : iElem*6 ) )

    end do ! iElem

  end subroutine deriveWSS3D
! ****************************************************************************** !

! ****************************************************************************** !
?? copy :: get_element_headtxt(deriveWSS2D)
    ! ---------------------------------------------------------------------------
    ! local variables
    integer :: iElem, shear_pos
    real(kind=rk) :: a0, a1
    real(kind=rk) :: tau(6*nElems) ! shear stress of all
    ! ---------------------------------------------------------------------------

    ! wall shear stress is depend on shear stress. get its position in globSys
    shear_pos = fun%input_varPos(1)

    ! first compute shear stress, save it to array tau
    ! the shear tensor has the following order: x, y, z, xy, xz, yz
    call varSys%method%val(shear_pos)%get_element( varSys  = varSys,  &
      &                                            elemPos = elemPos, &
      &                                            time    = time,    &
      &                                            tree    = tree,    &
      &                                            nElems  = nElems,  &
      &                                            nDofs   = nDofs,   &
      &                                            res     = tau      )


    do iElem = 1, nElems

      a1 = - tau((iElem-1)*6+1) - tau((iElem-1)*6+2)
      a0 =   tau((iElem-1)*6+1) * tau((iElem-1)*6+2) &
        &  - tau((iElem-1)*6+3) * tau((iElem-1)*6+3)

      res(iElem) = sqrt( a1 * a1 - 4._rk * a0 ) * div1_2

    end do ! iElem

  end subroutine deriveWSS2D
! ****************************************************************************** !


! ****************************************************************************** !
  !> Calculate the shear stress magnitude of a given element number with the
  !! given
  !!
  !! The interface has to comply to the abstract interface
  !! [[tem_varSys_module:tem_varSys_proc_element]].
  !!
?? copy :: get_element_headtxt(deriveShearMag)
    ! ---------------------------------------------------------------------------
    integer :: shear_pos, iElem
    real(kind=rk) :: shear(6*nElems) ! shear stress of all
    real(kind=rk) :: tau(6)
    ! ---------------------------------------------------------------------------

    shear_pos    = fun%input_varPos(1)

    ! first compute shear stress, save it to array tau
    ! the shear tensor has the following order: x, y, z, xy, xz, yz
    call varSys%method%val(shear_pos)%get_element( varSys  = varSys,  &
      &                                            elemPos = elemPos, &
      &                                            time    = time,    &
      &                                            tree    = tree,    &
      &                                            nElems  = nElems,  &
      &                                            nDofs   = nDofs,   &
      &                                            res     = shear    )


    do iElem = 1, nElems

      tau = shear( (iElem-1)*6+1 : iElem*6 )

      ! Von mises-criterion (see wikipedia article on von Mises yield criterion)
      ! wss = (tau_xx-tau_yy)**2 + (tau_yy-tau_zz)**2 + (tau_xx-tau_zz)**2  &
      ! &   + 6._rk*( tau_xy*tau_xy + tau_yz*tau_yz + tau_xz*tau_xz)
      ! wss = sqrt(0.5_rk*wss )
      res( iElem ) = sqrt( 0.5_rk*(   (tau(1)-tau(2))**2 &
        &                           + (tau(2)-tau(3))**2 &
        &                           + (tau(1)-tau(3))**2 &
        &                           + 6._rk * (   tau(4)*tau(4) &
        &                                       + tau(5)*tau(5) &
        &                                       + tau(6)*tau(6) )))

    end do ! iElem

  end subroutine deriveShearMag
! ****************************************************************************** !


! ****************************************************************************** !
  !> Calculate kinetic energy from density and velocity in auxField
  !! This routine sets the function Pointer for kinetic energy 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(deriveKE)
    ! ---------------------------------------------------------------------------
    !> Function pointer to perform specific operation.
    integer :: statePos, iElem, iLevel, elemOff
    type(mus_varSys_data_type), pointer :: fPtr
    integer :: dens_pos, vel_pos(3)
    real(kind=rk) :: rho, 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 density in auxField array
        dens_pos = varSys%method%val( fun%input_varPos(1) ) &
          &                     %auxField_varPos(1)

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

        ! density
        rho = auxField(iLevel)%val( elemOff + dens_pos )

        ! 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 * rho
      end do ! iElem
    end associate

  end subroutine deriveKE
! ****************************************************************************** !

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


! ****************************************************************************** !
   !> Derive external force variable defined as a source term.
   !! It evaluates spacetime function defined in lua file for force variable
   !! and convert it to state value which is to be added to the state
?? copy :: get_element_headtxt(derive_force)
    ! ---------------------------------------------------------------------------
    call tem_abort('Not implemented')
    ! ---------------------------------------------------------------------------
  end subroutine derive_force
! ****************************************************************************** !


! ************************************************************************** !
   !> Derive external force variable defined as a source term.
   !! It evaluates spacetime function defined in lua file for force variable
   !! and convert it to state value which is to be added to the state
   !! @todo KM: Not use we need seperate force for incompressible fluid model
?? copy :: get_element_headtxt(derive_force1stOrd)
    ! ---------------------------------------------------------------------- !
    call tem_abort('Not implemented')
    ! ---------------------------------------------------------------------- !
  end subroutine derive_force1stOrd
! ************************************************************************** !


! **************************************************************************** !
  !> 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.
  !! 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_absorbLayer( 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_absorbLayer
! **************************************************************************** !


! ****************************************************************************** !
  !> Update state with source variable "force".
  !! Force term used here is from:
  !! "Discrete lattice effects on the forcing term in the lattice Boltzmann
  !!  method", Zhaoli Guo, Chugung Zheng and Baochang Shi.
  !! In the paper, use force term is referred as Method 2 as:
  !! \[ F_i = w_i( (\vec{e}_i-\vec{u}*)/cs2 +
  !!       (\vec{e}_i \cdot \vec{u}*)\vec{e}_i/cs4) \cdot \vec{F} \]
  !! Force must be defined as body force per unit volume
  !! KM: If this force formula is used then velocity needs to be
  !! computed as u = \sum c_i f_i + \vec{F}/2
  !!
  !! Simuilar to derive routine but it updates the state whereas derive
  !! is used for tracking.
  !!
  !! 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_force( 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(:)
    ! -------------------------------------------------------------------- !
    type(mus_varSys_data_type), pointer :: fPtr
    type(mus_scheme_type), pointer :: scheme
    real(kind=rk) :: forceField(fun%elemLvl(iLevel)%nElems*3)
    real(kind=rk) :: G(3), velocity(3), ucx, uMinusCX(3), forceTerm
    integer :: nElems, iElem, iDir, QQ, nScalars, posInTotal, statePos, elemOff
    integer :: vel_pos(3)
    real(kind=rk) :: omega, omega_fac
    ! ---------------------------------------------------------------------------
    ! convert c pointer to solver type fortran pointer
    call c_f_pointer( varSys%method%val( fun%srcTerm_varPos )%method_data, &
      &               fPtr )
    scheme => fPtr%solverData%scheme

    ! Number of elements to apply source terms
    nElems = fun%elemLvl(iLevel)%nElems

    ! Get force which is refered in config file either its
    ! spacetime variable or operation variable
    call varSys%method%val(fun%data_varPos)%get_valOfIndex( &
      & varSys  = varSys,                                   &
      & time    = time,                                     &
      & iLevel  = iLevel,                                   &
      & idx     = fun%elemLvl(iLevel)%idx(1:nElems),        &
      & nVals   = nElems,                                   &
      & res     = forceField                                )

    ! convert physical to lattice
    forceField = forceField / fPtr%solverData%physics%fac(iLevel)%body_force

    ! constant parameter
    QQ = scheme%layout%fStencil%QQ
    nScalars = varSys%nScalars
    ! Position of velocity variable in auxField
    vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:3)

!$omp parallel do schedule(static), private( posInTotal, velocity, G, ucx, uMinusCx, forceTerm, omega )
    do iElem = 1, nElems
      posInTotal = fun%elemLvl(iLevel)%posInTotal(iElem)

      ! element offset
      elemoff = (posInTotal-1)*varSys%nAuxScalars
      ! obtain velocity from auxField
      velocity(1) = auxField(elemOff + vel_pos(1))
      velocity(2) = auxField(elemOff + vel_pos(2))
      velocity(3) = auxField(elemOff + vel_pos(3))

      ! force field on current element
      ! For incompressible model: this forceField should be divided by rho0.
      ! Since rho0 =1, this term is also valid for incompressible model
      G = forceField((iElem-1)*3+1 : iElem*3)

      ! get the correct omega value
      omega = scheme%field(1)%fieldProp%fluid%viscKine              &
        &                              %omLvl(iLevel)%val(posInTotal)
      omega_fac = 1.0_rk - omega * 0.5_rk

      ! force term:
      ! F_i = w_i( (\vec{e}_i-\vec{u}*)/cs2 +
      !       (\vec{e}_i \cdot \vec{u}*)\vec{e}_i/cs4) \cdot \vec{F}
      do iDir = 1, QQ
        ucx = dot_product( scheme%layout%fStencil%cxDirRK(:, iDir), &
          &                velocity )
        uMinusCx = scheme%layout%fStencil%cxDirRK(:, iDir) - velocity

        forceTerm = dot_product( uMinusCx * cs2inv               &
          &       + ucx * scheme%layout%fStencil%cxDirRK(:,iDir) &
          &       * cs4inv, G )

        ! position in state array
        statePos = ?SAVE?(iDir, 1, posInTotal, QQ, nScalars, nPdfSize, neigh)
        ! update outstate
        outState(statePos) = outState(statePos)                         &
          & + omega_fac * scheme%layout%weight( iDir ) * forceTerm

      end do

    end do !iElem
!$omp end parallel do

  end subroutine applySrc_force
! ****************************************************************************** !

! ****************************************************************************** !
  !> Update state with source variable "force" for MRT collision model.
  !! Force term used here is from:
  !! Chai, Z., & Zhao, T. (2012). Effect of the forcing term in the
  !! multiple-relaxation-time lattice Boltzmann equation on the shear stress
  !! or the strain rate tensor. Physical Review E, 86(1), 1–11.
  !! Force term for MRT is
  !! \[ \bar{F} = M^-1 (I-0.5 S) M F' \] and
  !! \[ F'_i = w_i( (\vec{e}_i-\vec{u}*)/cs2 +
  !!       (\vec{e}_i \cdot \vec{u}*)\vec{e}_i/cs4) \cdot \vec{F} \]
  !!
  !! \vec{F} is the force that must be defined as body force per unit volume
  !!
  !! Simuilar to derive routine but it updates the state whereas derive
  !! is used for tracking.
  !!
  !! 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_force_MRT( 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(:)
    ! -------------------------------------------------------------------- !
    type(mus_varSys_data_type), pointer :: fPtr
    type(mus_scheme_type), pointer :: scheme
    real(kind=rk) :: forceField(fun%elemLvl(iLevel)%nElems*3)
    real(kind=rk) :: G(3), velocity(3), ucx, uMinusCX(3)
    integer :: nElems, iElem, iDir, QQ, nScalars, posInTotal, statePos, elemOff
    integer :: vel_pos(3)
    real(kind=rk) :: omegaKine, omegaBulk, discForce
    !real(kind=rk), allocatable :: forceTerm(:), momForce(:), s_mrt(:)
    !real(kind=rk), allocatable :: mInvXOmega(:,:)
    ! Allocate with maximum QQ to avoid allocation and deallocation
    real(kind=rk) :: forceTerm(27), momForce(27), s_mrt(27)
    real(kind=rk) :: mInvXOmega(27,27)
    ! ---------------------------------------------------------------------------
    ! convert c pointer to solver type fortran pointer
    call c_f_pointer( varSys%method%val( fun%srcTerm_varPos )%method_data, &
      &               fPtr )
    scheme => fPtr%solverData%scheme

    ! Number of elements to apply source terms
    nElems = fun%elemLvl(iLevel)%nElems

    ! Get force which is refered in config file either its
    ! spacetime variable or operation variable
    call varSys%method%val(fun%data_varPos)%get_valOfIndex( &
      & varSys  = varSys,                                   &
      & time    = time,                                     &
      & iLevel  = iLevel,                                   &
      & idx     = fun%elemLvl(iLevel)%idx(1:nElems),        &
      & nVals   = nElems,                                   &
      & res     = forceField                                )

    ! convert physical to lattice
    forceField = forceField / fPtr%solverData%physics%fac(iLevel)%body_force

    ! constant parameter
    QQ = scheme%layout%fStencil%QQ
!    allocate(s_mrt(QQ))
!    allocate(forceTerm(QQ))
!    allocate(momForce(QQ))
!    allocate(mInvXOmega(QQ,QQ))

    nScalars = varSys%nScalars
    ! Position of velocity variable in auxField
    vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:3)

    omegaBulk = scheme%field(1)%fieldProp%fluid%omegaBulkLvl(iLevel)

!$omp parallel do schedule(static), private( posInTotal, velocity, G, ucx, uMinusCx, forceTerm )
    do iElem = 1, nElems
      posInTotal = fun%elemLvl(iLevel)%posInTotal(iElem)

      ! element offset
      elemoff = (posInTotal-1)*varSys%nAuxScalars
      ! obtain velocity from auxField
      velocity(1) = auxField(elemOff + vel_pos(1))
      velocity(2) = auxField(elemOff + vel_pos(2))
      velocity(3) = auxField(elemOff + vel_pos(3))

      ! force field on current element
      ! For incompressible model: this forceField should be divided by rho0.
      ! Since rho0 =1, this term is also valid for incompressible model
      G = forceField((iElem-1)*3+1 : iElem*3)

      ! get the correct omega value
      omegaKine = scheme%field(1)%fieldProp%fluid%viscKine          &
        &                              %omLvl(iLevel)%val(posInTotal)
      ! MRT omegas
      ! overwrite omegaKine term in the element loop
      s_mrt = scheme%field(1)%fieldProp%fluid                              &
        &           %mrtPtr(omegaKine=omegaKine, omegaBulk=omegaBulk, QQ=QQ)

      ! M^-1 * (I-0.5 S)
      s_mrt = 1.0_rk - 0.5_rk * s_mrt
      do iDir = 1, QQ
        mInvXOmega(1:QQ,iDir) = scheme%layout%moment%toPDF%A(1:QQ,iDir) &
          &                   * s_mrt(iDir)
      end do

      ! force term:
      ! F_i = w_i( (\vec{e}_i-\vec{u}*)/cs2 +
      !       (\vec{e}_i \cdot \vec{u}*)\vec{e}_i/cs4) \cdot \vec{F}
      do iDir = 1, QQ
        ucx = dot_product( scheme%layout%fStencil%cxDirRK(:, iDir), &
          &                velocity )
        uMinusCx = scheme%layout%fStencil%cxDirRK(:, iDir) - velocity

        forceTerm(iDir) = scheme%layout%weight(iDir)             &
          &       * dot_product( uMinusCx * cs2inv               &
          &       + ucx * scheme%layout%fStencil%cxDirRK(:,iDir) &
          &       * cs4inv, G )
      end do

      ! Force moments: M * F
      !do iDir = 1, QQ
      !  momForce(iDir) = sum(scheme%layout%moment%toMoments%A(iDir,:) * forceTerm)
      !end do
      momForce = matmul( scheme%layout%moment%toMoments%A(1:QQ, 1:QQ), &
        &                forceTerm(1:QQ) )

      !discForce = matmul( omegaTerm, forceTerm )
      do iDir = 1, QQ
        ! discrete force
        ! \bar{F} =  M^-1 (I-S/2) M F
        discForce = dot_product(mInvXOmega(iDir,1:QQ),  momForce(1:QQ))
        ! position in state array
        statePos = ?SAVE?(iDir, 1, posInTotal, QQ, nScalars, nPdfSize, neigh)
        ! update outstate
        outState(statePos) = outState(statePos) + discForce
      end do

    end do !iElem
!$omp end parallel do

    !deallocate(s_mrt)
    !deallocate(forceTerm)
    !deallocate(momForce)
    !deallocate(mInvXOmega)
  end subroutine applySrc_force_MRT
! ****************************************************************************** !


! ****************************************************************************** !
  !> Update state with source variable "force" for d3q19 MRT collision model.
  !! Force term used here is from:
  !! Chai, Z., & Zhao, T. (2012). Effect of the forcing term in the
  !! multiple-relaxation-time lattice Boltzmann equation on the shear stress
  !! or the strain rate tensor. Physical Review E, 86(1), 1–11.
  !! Force term for MRT is
  !! \[ \bar{F} = M^-1 (I-0.5 S) M F' \] and
  !! \[ F'_i = w_i( (\vec{e}_i-\vec{u}*)/cs2 +
  !!       (\vec{e}_i \cdot \vec{u}*)\vec{e}_i/cs4) \cdot \vec{F} \]
  !!
  !! Force must be defined as body force per unit volume
  !!
  !! Simuilar to derive routine but it updates the state whereas derive
  !! is used for tracking.
  !!
  !! 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_force_MRT_d3q19( 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(:)
    ! -------------------------------------------------------------------- !
    type(mus_varSys_data_type), pointer :: fPtr
    type(mus_scheme_type), pointer :: scheme
    real(kind=rk) :: forceField(fun%elemLvl(iLevel)%nElems*3)
    real(kind=rk) :: G(3), velocity(3), ucx, uMinusCX(3)
    integer :: nElems, iElem, iDir, QQ, nScalars, posInTotal, statePos, elemOff
    integer :: vel_pos(3)
    real(kind=rk) :: omegaKine, omegaBulk, discForce
    real(kind=rk) :: forceTerm(19), momForce(19), s_mrt(19)
    real(kind=rk) :: mInvXOmega(19,19)
    ! ---------------------------------------------------------------------------
    ! convert c pointer to solver type fortran pointer
    call c_f_pointer( varSys%method%val( fun%srcTerm_varPos )%method_data, &
      &               fPtr )
    scheme => fPtr%solverData%scheme

    ! Number of elements to apply source terms
    nElems = fun%elemLvl(iLevel)%nElems

    ! Get force which is refered in config file either its
    ! spacetime variable or operation variable
    call varSys%method%val(fun%data_varPos)%get_valOfIndex( &
      & varSys  = varSys,                                   &
      & time    = time,                                     &
      & iLevel  = iLevel,                                   &
      & idx     = fun%elemLvl(iLevel)%idx(1:nElems),        &
      & nVals   = nElems,                                   &
      & res     = forceField                                )

    ! convert physical to lattice
    forceField = forceField / fPtr%solverData%physics%fac(iLevel)%body_force

    ! constant parameter
    QQ = 19

    nScalars = varSys%nScalars
    ! Position of velocity variable in auxField
    vel_pos = varSys%method%val(derVarPos(1)%velocity)%auxField_varPos(1:3)

    omegaBulk = scheme%field(1)%fieldProp%fluid%omegaBulkLvl(iLevel)
    ! MRT omegas
    ! overwrite omegaKine term in the element loop
    ! KM: For incompressible model: omegaBulk is unused in mrtPtr
    s_mrt = scheme%field(1)%fieldProp%fluid                           &
      &           %mrtPtr(omegaKine=1.0_rk, omegaBulk=omegaBulk, QQ=QQ)

    ! F = M^-1 (I-0.5 S) M F
    ! (I-0.5 S) - omega for force term
    s_mrt = 1.0_rk - 0.5_rk * s_mrt

!$omp parallel do schedule(static), private( posInTotal, velocity, G, ucx, uMinusCx, forceTerm )
    do iElem = 1, nElems
      posInTotal = fun%elemLvl(iLevel)%posInTotal(iElem)

      ! element offset
      elemoff = (posInTotal-1)*varSys%nAuxScalars
      ! obtain velocity from auxField
      velocity(1) = auxField(elemOff + vel_pos(1))
      velocity(2) = auxField(elemOff + vel_pos(2))
      velocity(3) = auxField(elemOff + vel_pos(3))

      ! force field on current element
      ! For incompressible model: this forceField should be divided by rho0.
      ! Since rho0 =1, this term is also valid for incompressible model
      G = forceField((iElem-1)*3+1 : iElem*3)

      ! get the correct omega value
      omegaKine = scheme%field(1)%fieldProp%fluid%viscKine          &
        &                              %omLvl(iLevel)%val(posInTotal)
      ! MRT omegas
      ! overwrite omegaKine term in the element loop
      s_mrt(10) = 1.0_rk - 0.5_rk * omegaKine
      s_mrt(12) = s_mrt(10)
      s_mrt(14) = s_mrt(10)
      s_mrt(15) = s_mrt(10)
      s_mrt(16) = s_mrt(10)

      ! M^-1 (1-0.5 S)
      do iDir = 1, QQ
        mInvXOmega(:,iDir) = scheme%layout%moment%toPDF%A(:,iDir) * s_mrt(iDir)
      end do

      ! force term:
      ! F_i = w_i( (\vec{e}_i-\vec{u}*)/cs2 +
      !       (\vec{e}_i \cdot \vec{u}*)\vec{e}_i/cs4) \cdot \vec{F}
      do iDir = 1, QQ
        ucx = dot_product( scheme%layout%fStencil%cxDirRK(:, iDir), &
          &                velocity )
        uMinusCx = scheme%layout%fStencil%cxDirRK(:, iDir) - velocity

        forceTerm(iDir) = scheme%layout%weight(iDir)             &
          &       * dot_product( uMinusCx * cs2inv               &
          &       + ucx * scheme%layout%fStencil%cxDirRK(:,iDir) &
          &       * cs4inv, G )
      end do

      ! Force moments: M * F
      !do iDir = 1, QQ
      !  momForce(iDir) = sum(scheme%layout%moment%toMoments%A(iDir,:) * forceTerm)
      !end do
      momForce = matmul(scheme%layout%moment%toMoments%A(1:QQ,1:QQ), &
        &               forceTerm(1:QQ))

      do iDir = 1, QQ
        ! discrete force
        ! \bar{F} =  M^-1 (I-S/2) M F
        discForce = sum(mInvXOmega(iDir,1:QQ) * momForce(1:QQ))
        ! position in state array
        statePos = ?SAVE?(iDir, 1, posInTotal, QQ, nScalars, nPdfSize, neigh)
        ! update outstate
        outState(statePos) = outState(statePos) + discForce
      end do

    end do !iElem
!$omp end parallel do

  end subroutine applySrc_force_MRT_d3q19
! ****************************************************************************** !

! ************************************************************************** !
  !> Update state with source variable "force_1stOrd"
  !! Force term used here is from:
  !! "A D3Q27 multiple-relaxation-time lattice Boltzmann method for
  !! turbulent flows", K. Suga, Y. Kuwata, K. Takashima, R. Chikasue
  !!
  !! \[ F_i = w_i/c_s^2 ( \vec{e}_i \cdot \vec{F} ) \]
  !! Force must be defined as body force per unit volume
  !! This force term can be applied for both compressible and incompressible
  !! LBM models
  !!
  !! Similar to derive routine but it updates the state whereas derive
  !! is used for tracking
  !!
  !! 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_force1stOrd( 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(:)
    ! -------------------------------------------------------------------- !
    type(mus_varSys_data_type), pointer :: fPtr
    type(mus_scheme_type), pointer :: scheme
    real(kind=rk) :: forceField(fun%elemLvl(iLevel)%nElems*3)
    real(kind=rk) :: FF_elem(3)
    real(kind=rk) :: forceTerm
    integer :: nElems, iElem, iDir, QQ, nScalars
    integer :: posInTotal
    ! ---------------------------------------------------------------------- !
    ! convert c pointer to solver type fortran pointer
    call c_f_pointer( varSys%method%val( fun%srcTerm_varPos )%method_data, &
      &               fPtr )
    scheme => fPtr%solverData%scheme

    ! Number of elements to apply source terms
    nElems = fun%elemLvl(iLevel)%nElems

    ! Get force which is refered in config file either its
    ! spacetime variable or operation variable
    call varSys%method%val(fun%data_varPos)%get_valOfIndex( &
      & varSys  = varSys,                                   &
      & time    = time,                                     &
      & iLevel  = iLevel,                                   &
      & idx     = fun%elemLvl(iLevel)%idx(1:nElems),        &
      & nVals   = nElems,                                   &
      & res     = forceField                                )

!write(dbgUnit(1),*) 'ApplySrc_force1stOrdIncomp'
!    do iElem = 1, nElems
!      posInTotal = fun%elemLvl(iLevel)%posInTotal(iElem)
!      FF_elem = forceField((iElem-1)*3+1 : iElem*3)
!write(dbgUnit(1),*) 'treeID ', scheme%levelDesc(iLevel)%total(posInTotal) &
!  & , 'force ', FF_elem
!    end do

    ! convert physical to lattice
    forceField = forceField / fPtr%solverData%physics%fac(iLevel)%body_force

    ! constant parameter
    QQ = scheme%layout%fStencil%QQ
    nScalars = varSys%nScalars

    do iElem = 1, nElems
      ! to access level wise state array
      posInTotal = fun%elemLvl(iLevel)%posInTotal(iElem)

      ! force field on current element
      FF_elem = forceField((iElem-1)*3+1 : iElem*3)

      ! force term:
      ! F_i = w_i/cs2 (\vec{e}_i \cdot \vec{F}
      do iDir = 1, QQ

        forceTerm = dot_product( scheme%layout%fStencil    &
          &                      %cxDirRK(:,iDir), FF_elem )

        outState( ?SAVE?(iDir,1,posInTotal,QQ,nScalars,nPdfSize,neigh) ) &
          & = outState(                                                  &
          & ?SAVE?(iDir,1,posInTotal,QQ,nScalars,nPdfSize,neigh) )       &
          & + scheme%layout%weight( iDir ) * cs2inv * forceTerm

      end do

    end do !iElem
  end subroutine applySrc_force1stOrd
! ************************************************************************** !


! ****************************************************************************** !
  !> Calculate wss from shear stress (tau)
  !! tau: x, y, z, xy, yz, xz
  pure function getWSS( tau ) result( wss )
    ! ---------------------------------------------------------------------------
    real(kind=rk), intent(in) :: tau(6)
    real(kind=rk) :: a0, a1, a2
    real(kind=rk) :: q_stress, q_stress_cube, q_term, r_stress, d_stress
    real(kind=rk) :: theta_stress, cos_term, shear_1, shear_2, shear_3
    real(kind=rk) :: wss
    ! ---------------------------------------------------------------------------

      ! tau_x + tau_y + tau_z
      a2 = -1._rk * ( tau(1) + tau(2) + tau(3) )

      ! tau_x * tau_y + tau_y * tau_z + tau_x * tau_z - tau_xy^2 - tau_xz^2 - tau_yz^2
      a1 =    tau(1) * tau(2) + tau(2) * tau(3) + tau(1) * tau(3)              &
        &   - tau(4) * tau(4) - tau(5) * tau(5) - tau(6) * tau(6)

      a0 = -1._rk * (   tau(1) * tau(2) * tau(3)                               &
        &             + 2._rk * tau(4) * tau(5) * tau(6)                       &
        &             - tau(1) * tau(6) * tau(6)                               &
        &             - tau(2) * tau(5) * tau(5)                               &
        &             - tau(3) * tau(4) * tau(4) )

      q_stress = div1_9 * ( 3._rk * a1 - a2 * a2 )
      q_stress_cube = q_stress * q_stress * q_stress
      if ( q_stress < 0._rk ) then
        q_term = 2._rk * sqrt( -q_stress )
      else
        q_term = 0._rk
      end if
      r_stress = div1_54 * ( 9._rk*a2*a1 - 27._rk*a0 - 2._rk*a2*a2*a2 )
      d_stress = q_stress_cube + r_stress * r_stress

      if ( d_stress < 0._rk ) then
        ! the solutions are real and unequal
        theta_stress = div1_3 * acos( r_stress / sqrt( -q_stress_cube ))

        cos_term = div3_4 * q_term * cos( theta_stress )

        ! three values of maximum shear stress
        shear_1 = abs( sqrt3 * div1_2 * q_term * sin( theta_stress ))

        shear_2 = abs( cos_term - shear_1*div1_2 )

        shear_3 = abs( cos_term + shear_1*div1_2 )

        ! choose the max shear as Wall Shear Stress(wss)
        wss = max(shear_1, shear_2, shear_3)
      else
        wss = 0._rk
      end if
  end function getWSS
! ****************************************************************************** !


! ****************************************************************************** !
  !> This routine computes equilbrium from density and velocity
  !! This must comply with interface in mus_variable_module derive_FromMacro
  !!
  !! 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 deriveEquil_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)      &
          & * ( 1._rk + cs2inv * ucx + ucx * ucx * t2cs4inv - usq )
      enddo

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

! ************************************************************************** !
  !> 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 deriveEquil_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 * ( 1.0_rk &
          & + cs2inv * ucx + ucx * ucx * t2cs4inv - usq )
      enddo

    end do
  end subroutine deriveEquil_fromAux
! ************************************************************************** !

! ****************************************************************************** !
  !> This routine computes density 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 deriveRho_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
    ! ---------------------------------------------------------------------------

    res = 0.0_rk
    do iElem = 1, nElems
      do iDir = 1, layout%fStencil%QQ
        res( iElem ) = res( iElem ) + state( iDir+(iElem-1)*varSys%nScalars )
      end do
    end do

  end subroutine deriveRho_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 deriveVel_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) :: rho, pdf( layout%fStencil%QQ )
    ! ---------------------------------------------------------------------------

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

  end subroutine deriveVel_FromState
! ****************************************************************************** !

! ****************************************************************************** !
  !> This routine computes velocity from pre collision state array using Fetch
  !!
  !! 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 deriveVel_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) :: rho, 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
      ! density
      rho = sum( pdf )
      ! 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) / rho
    end do

  end subroutine deriveVel_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 deriveAux_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 ), rho, inv_rho
    ! ------------------------------------------------------------------------ !
    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
      rho = sum( pdf )
      auxField(elemOff+dens_pos) = rho

      ! velocity
      inv_rho = 1.0_rk/rho
      auxField(elemOff+vel_pos(1)) = sum( pdf * stencil%cxDirRK(1,:) ) * inv_rho
      auxField(elemOff+vel_pos(2)) = sum( pdf * stencil%cxDirRK(2,:) ) * inv_rho
      auxField(elemOff+vel_pos(3)) = sum( pdf * stencil%cxDirRK(3,:) ) * inv_rho
    end do

  end subroutine deriveAux_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 deriveEq_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) :: pdf( layout%fStencil%QQ )
    real(kind=rk) :: rho, vel(3), usq, ucx
    ! ---------------------------------------------------------------------------

    QQ = layout%fStencil%QQ

    do iElem = 1, nElems

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

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

      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  &
          &             * ( 1._rk + cs2inv * ucx + ucx * ucx * t2cs4inv - usq )
      enddo

    end do

  end subroutine deriveEq_FromState
! ****************************************************************************** !

! ****************************************************************************** !
?? copy :: get_element_headtxt(deriveMomentumChange)
    ! ---------------------------------------------------------------------------
    integer :: statePos, iElem, iComp, iLevel, varPos
    type(mus_varSys_data_type), pointer :: fPtr
    integer :: pdfPos, nCompsPDF
    real(kind=rk), allocatable :: f(:)
    integer :: nSize
    ! ---------------------------------------------------------------------------
    call C_F_POINTER( fun%method_Data, fPtr )

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

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

      ! transfer PDF to array f
      do iComp = 1, nCompsPDF
        varPos = varSys%method%val(pdfPos)%state_varPos(iComp)
        f( iComp ) =  &
          & fPtr%solverData%scheme%state( iLevel )%val(       &
          ! position of this state variable in the state array
          & ?IDX?( varPos, statePos, varSys%nScalars, nSize), &
          & fPtr%solverData%scheme%pdf( iLevel )%nNext )
      end do !iComp

      ! only count on directions that insect with BC

    end do ! iElem

    deallocate( f )

  end subroutine deriveMomentumChange
! ****************************************************************************** !


! ****************************************************************************** !
  !> Initiates the calculation of moment for 2D
  !! This routine sets the function Pointer for moment for 2D 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]].
  !!
  !> For 2D only!
?? copy :: get_element_headtxt(deriveMoment)
    ! ---------------------------------------------------------------------------
    !> Function pointer to perform specific operation.
    procedure(mus_derive_fromPDF), pointer :: fnCalcPtr
    ! ---------------------------------------------------------------------------
    fnCalcPtr => mus_deriveMoment

    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 deriveMoment
! ****************************************************************************** !

  ! ************************************************************************* !
  !         Subroutines with common interface for values from index           !
  ! ************************************************************************* !

! ****************************************************************************** !
  !> Initiates the calculation of density.
  !! This routine sets the function Pointer for density 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 deriveDensity_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_derivedensity

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

  end subroutine deriveDensity_fromIndex
! ****************************************************************************** !


! ****************************************************************************** !
  !> 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 deriveVelocity_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_derivevelocity

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

  end subroutine deriveVelocity_fromIndex
! ****************************************************************************** !


! ****************************************************************************** !
  !> Initiates the calculation of pressure.
  !! This routine sets the function Pointer for pressure 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 derivePressure_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(:)
    ! ---------------------------------------------------------------------------
    integer :: dens_pos
    type(mus_varSys_data_type), pointer     :: fPtr
    ! ---------------------------------------------------------------------------
    !convert pointer from C to Fotran
    call C_F_POINTER( fun%method_Data, fPtr )
    ! position of density in glob system
    dens_pos = fun%input_varPos(1)
    ! get density variable from auxField
    call varSys%method%val( dens_pos )%get_ValOfIndex( &
      &     varSys = varSys,                           &
      &     time   = time,                             &
      &     iLevel = iLevel,                           &
      &     idx    = fPtr%opData%input_pntIndex(1)     &
      &              %indexLvl(iLevel)%val( idx(:) ),  &
      &     nVals  = nVals,                            &
      &     res    = res                               )

    ! convert density to pressure
    res(1:nVals) = ( res(1:nVals) - rho0 ) * cs2

  end subroutine derivePressure_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 deriveEquil_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_deriveEquil

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

  end subroutine deriveEquil_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 deriveNonEquil_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_deriveNonEquil

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

  end subroutine deriveNonEquil_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 deriveKe_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 :: dens_pos, vel_pos, iVal
    real(kind=rk) :: mass_dens(nVals), 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 mass density values for IDX
    dens_pos = fun%input_varPos(1)
    call varSys%method%val( dens_pos )%get_ValOfIndex( &
      &     varSys = varSys,                           &
      &     time   = time,                             &
      &     iLevel = iLevel,                           &
      &     idx    = fPtr%opData%input_pntIndex(1)     &
      &              %indexLvl(iLevel)%val( idx(:) ),  &
      &     nVals  = nVals,                            &
      &     res    = mass_dens                         )

    ! get velocity values for IDX
    vel_pos = fun%input_varPos(2)
    call varSys%method%val( vel_pos )%get_ValOfIndex( &
      &     varSys = varSys,                          &
      &     time   = time,                            &
      &     iLevel = iLevel,                          &
      &     idx    = fPtr%opData%input_pntIndex(2)    &
      &              %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 * mass_dens(iVal)
    end do

  end subroutine deriveKe_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 deriveStrainRate_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_deriveStrainRate

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

  end subroutine deriveStrainRate_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 deriveMomentum_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 :: dens_pos, vel_pos, iVal
    real(kind=rk) :: mass_dens(nVals)
    ! ---------------------------------------------------------------------------

?? 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 mass density values for IDX
    dens_pos = fun%input_varPos(1)
    call varSys%method%val( dens_pos )%get_ValOfIndex( &
      &     varSys = varSys,                           &
      &     time   = time,                             &
      &     iLevel = iLevel,                           &
      &     idx    = fPtr%opData%input_pntIndex(1)     &
      &              %indexLvl(iLevel)%val( idx(:) ),  &
      &     nVals  = nVals,                            &
      &     res    = mass_dens                         )

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

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

  end subroutine deriveMomentum_fromIndex
! ****************************************************************************** !


  ! ************************************************************************* !
  !          Subroutines with common interface called by fnCalcPtr            !
  ! ************************************************************************* !

! ****************************************************************************** !
  !> Calculate the density of a given set of elements (sum up all links).
  !! This routine is used to compute density for all scheme kinds
  !! For multispecies, it can compute both species density and mixture density
  !!
  !! The interface has to comply to the abstract interface
  !! [[mus_varSys_module:mus_derive_fromPDF]].
  !!
  recursive subroutine mus_derivedensity(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
    ! ---------------------------------------------------------------------------
    real(kind=rk), allocatable                :: tmpPDF(:)
    real(kind=rk)                             :: dens
    integer                                   :: iDep, iVal, iComp
    integer                                   :: pdfPos, nCompPDF
    ! ---------------------------------------------------------------------------
    pdfPos = fun%input_varPos(1)
    nCompPDF = varSys%method%val(pdfPos)%nComponents
    allocate( tmpPDF( nCompPDF ) )
    res = 0.0_rk

    do iVal = 1, nVals
      dens = 0.0_rk
      do iDep = 1, fun%nInputs
        tmpPDF = pdf(                                &
          &  (iDep-1)*nVals + (iVal-1)*nCompPDF + 1  &
          &  :                                       &
          &  (iDep-1)*nVals + iVal*nCompPDF          )
        do iComp = 1, nCompPDF
          dens = dens + tmpPDF( iComp )
        end do !iComp
      end do !iDep
      res( iVal ) = dens
    end do !iVal
    deallocate( tmpPDF )
  end subroutine mus_deriveDensity
! ****************************************************************************** !


! ****************************************************************************** !
  !> Calculate the velocity of a given element number with the given input
  !! vector (sum up all values)
  !!
  !! The interface has to comply to the abstract interface
  !! [[mus_varSys_module:mus_derive_fromPDF]].
  !!
  recursive subroutine mus_deriveVelocity(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)                             :: densInv
    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
        densInv = 1._rk / sum( tmpPDF )
        res( (iVal-1)*fun%nComponents + iComp : iVal*fun%nComponents ) = &
          &  sum(tmpPDF * scheme%layout%fStencil%cxDirRK(iComp,:))       &
          &  * densInv
      end do ! iComp
    end do ! iVal
    deallocate( tmpPDF )
  end subroutine mus_deriveVelocity
! ****************************************************************************** !


! ****************************************************************************** !
  !> Calculate the equlibrium of given elements with the given input state
  !! array.
  !!
  !! The interface has to comply to the abstract interface
  !! [[mus_varSys_module:mus_derive_fromPDF]].
  !!
  recursive subroutine mus_deriveEquil(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), densInv
    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 )
      ! computes density and velocity
      dens   = sum(tmpPDF)
      densInv = 1.0_rk / dens
      vel(1) = sum(tmpPDF * scheme%layout%fStencil%cxDirRK(1,:)) &
        &    * densInv
      vel(2) = sum(tmpPDF * scheme%layout%fStencil%cxDirRK(2,:)) &
        &    * densInv
      vel(3) = sum(tmpPDF * scheme%layout%fStencil%cxDirRK(3,:)) &
        &    * densInv

      do iDir = 1, scheme%layout%fStencil%QQ
?? copy :: EqDir_comp(fEq(iDir), dens, vel, scheme%layout, iDir, cs2inv, rk)
      end do

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

  end subroutine mus_deriveEquil
! ****************************************************************************** !


! ****************************************************************************** !
  !> Calculate the kinetic energy as
  !!  E = 1/2 rho ( ux^2 + uy^2 +uz^2 )
  !!
  !! The interface has to comply to the abstract interface
  !! [[mus_varSys_module:mus_derive_fromPDF]].
  !!
  recursive subroutine mus_deriveKE(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)                             :: dens, 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 )
      dens   = sum(tmpPDF)
      vel(1) = sum(tmpPDF * scheme%layout%fStencil%cxDirRK(1,:)) / dens
      vel(2) = sum(tmpPDF * scheme%layout%fStencil%cxDirRK(2,:)) / dens
      vel(3) = sum(tmpPDF * scheme%layout%fStencil%cxDirRK(3,:)) / dens
      res( iVal ) = sum( vel(:)*vel(:) ) * 0.5_rk * dens
    end do !iVal
    deallocate( tmpPDF )

  end subroutine mus_deriveKE
! ****************************************************************************** !


! ****************************************************************************** !
  !> Calculate the momentum of a given element number with the given input
  !!        vector (sum up all values)
  !!
  !! The interface has to comply to the abstract interface
  !! [[mus_varSys_module:mus_derive_fromPDF]].
  !!
  recursive subroutine mus_deriveMomentum(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                                   :: 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 )
      res((iVal-1)*3+1) = sum( tmpPDF * scheme%layout%fStencil%cxDirRK(1,:) )
      res((iVal-1)*3+2) = sum( tmpPDF * scheme%layout%fStencil%cxDirRK(2,:) )
      res((iVal-1)*3+3) = sum( tmpPDF * scheme%layout%fStencil%cxDirRK(3,:) )
    end do !iVal
    deallocate( tmpPDF )

  end subroutine mus_deriveMomentum
! ****************************************************************************** !


! ****************************************************************************** !
  !> For 2D only!
  !!
  !! The interface has to comply to the abstract interface
  !! [[mus_varSys_module:mus_derive_fromPDF]].
  !!
  recursive subroutine mus_deriveMoment(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
    ! ---------------------------------------------------------------------------
    real(kind=rk), allocatable                :: tmpPDF(:)
    integer                                   :: pdfPos, nCompsPDF
    integer                                   :: iVal
    ! ---------------------------------------------------------------------------
    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 )

      res( (iVal-1)*9+1 ) = sum( tmpPDF ) - 1.0_rk
      res( (iVal-1)*9+2 ) = - sum(tmpPDF(1:4)) + 2.0_rk &
        &                   * sum(tmpPDF(5:8)) - 4.0_rk &
        &                   * tmpPDF(9) + 2.0_rk
      res( (iVal-1)*9+3 ) = - 2.0_rk * sum(tmpPDF(1:4)) &
        &                   + sum(tmpPDF(5:8)) + 4.0_rk &
        &                   * tmpPDF(9) - 1.0_rk
      res( (iVal-1)*9+4 ) = tmpPDF(3) + tmpPDF(7) + tmpPDF(8) &
        &                   - tmpPDF(1) - tmpPDF(5) - tmpPDF(6)
      res( (iVal-1)*9+6 ) = tmpPDF(4) + tmpPDF(6) + tmpPDF(8) &
        &                   - tmpPDF(2) - tmpPDF(5) - tmpPDF(7)
      res( (iVal-1)*9+5 ) = 2.0_rk * ( tmpPDF(1) - tmpPDf(3) ) - tmpPDF(5) &
        &                   - tmpPDF(6) + tmpPDF(7) + tmpPDF(8)
      res( (iVal-1)*9+7 ) = 2.0_rk * ( tmpPDF(2) - tmpPDF(4) ) - tmpPDf(5) &
        &                   + tmpPDF(6) - tmpPDF(7) + tmpPDF(8)
      res( (iVal-1)*9+8 ) = tmpPDF(1) - tmpPDF(2) + tmpPDF(3) - tmpPDF(4)
      res( (iVal-1)*9+9 ) = tmpPDF(5) - tmpPDF(6) - tmpPDF(7) + tmpPDF(8)
    end do ! iVal
    deallocate( tmpPDF )

  end subroutine mus_deriveMoment
! ****************************************************************************** !


! ****************************************************************************** !
  recursive subroutine mus_derivePressure(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
    ! ---------------------------------------------------------------------------
    real(kind=rk), allocatable                :: tmpPDF(:)
    integer                                   :: pdfPos, nCompsPDF
    integer                                   :: iVal
    ! ---------------------------------------------------------------------------
    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 )
      ! convert pdf to density(sum of tmpPDF) and then density into pressure
      res(iVal) = ( sum(tmpPDF) - 1.0_rk )  * cs2
    end do !iVal

  end subroutine mus_derivePressure
! ****************************************************************************** !


! ****************************************************************************** !
  !> Calculate the Non-Equlibrium
  !!
  recursive subroutine mus_deriveNonEquil(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,:)) / dens
      vel(2) = sum(tmpPDF * scheme%layout%fStencil%cxDirRK(2,:)) / dens
      vel(3) = sum(tmpPDF * scheme%layout%fStencil%cxDirRK(3,:)) / dens

      ! computes equilibrium
      do iDir = 1, scheme%layout%fStencil%QQ
?? copy :: EqDir_comp(fEq(iDir), dens, vel, scheme%layout, iDir, 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_deriveNonEquil
! ****************************************************************************** !


! ****************************************************************************** !
  !> author: Jiaxing Qi
  !! Calculate the strain rate ( or rate of strain, or rate of deformation)
  !!
  !! The equation is:
  !! \[
  !!  \tau_{\alpha \beta}=
  !!    -\frac{3\omega}{2\rho} \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_deriveStrainRate(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                :: nonEq(:)
    real(kind=rk), allocatable                :: tau(:)
    real(kind=rk), allocatable                :: fEq(:)
    real(kind=rk)                             :: dens, vel(3), omega
    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( nonEq( nCompsPDF ) )
    allocate( fEq( 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,:)) / dens
      vel(2) = sum(tmpPDF * scheme%layout%fStencil%cxDirRK(2,:)) / dens
      vel(3) = sum(tmpPDF * scheme%layout%fStencil%cxDirRK(3,:)) / dens

      ! computes equilibrium
      do iDir = 1, scheme%layout%fStencil%QQ
?? copy :: EqDir_comp(fEq(iDir), dens, vel, scheme%layout, iDir, 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( scheme%layout%fStencil%cxcx(:,:), &
        &                 nonEq(:), &
        &                 scheme%layout%fStencil%QQ )

      ! S_{ij} = - (omega/(2*rho*cs2))* sum_k(c_{ki} c{kj} f^neq_k)
      res( (iVal-1)*fun%nComponents + 1 : iVal*fun%nComponents ) &
        &  = tau(:) * (-1.5_rk) * omega / dens
    end do !iVal

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

  end subroutine mus_deriveStrainRate
! ****************************************************************************** !

end module mus_derQuan_module
