! Copyright (c) 2016-2017, 2019, 2021 Kannan Masilamani <kannan.masilamani@uni-siegen.de>
! Copyright (c) 2016-2017 Verena Krupp <verena.krupp@uni-siegen.de>
! Copyright (c) 2016-2020 Peter Vitt <peter.vitt2@uni-siegen.de>
! Copyright (c) 2016-2017 Tobias Schneider <tobias1.schneider@student.uni-siegen.de>
! Copyright (c) 2016 Jiaxing Qi <jiaxing.qi@uni-siegen.de>
! Copyright (c) 2016 Kartik Jain <kartik.jain@uni-siegen.de>
! Copyright (c) 2016-2017, 2019, 2021 Harald Klimach <harald.klimach@uni-siegen.de>
! Copyright (c) 2016 Jana Gericke <jana.gericke@student.uni-siegen.de>
! Copyright (c) 2017 Michael Gaida  <michael.gaida@student.uni-siegen.de>
! Copyright (c) 2017 Neda Ebrahimi Pour <neda.epour@uni-siegen.de>
! Copyright (c) 2018 Robin Weihe <robin.weihe@student.uni-siegen.de>
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions are met:
!
! 1. Redistributions of source code must retain the above copyright notice, this
! list of conditions and the following disclaimer.
!
! 2. Redistributions in binary form must reproduce the above copyright notice,
! this list of conditions and the following disclaimer in the documentation
! and/or other materials provided with the distribution.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "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 THE COPYRIGHT HOLDER 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.
! Copyright (c) 2013 Harald Klimach <harald.klimach@uni-siegen.de>
! Copyright (c) 2013-2014 Nikhil Anand <nikhil.anand@uni-siegen.de>
! Copyright (c) 2014, 2016 Kannan Masilamani <kannan.masilamani@uni-siegen.de>
! Copyright (c) 2015, 2018, 2020 Peter Vitt <peter.vitt2@uni-siegen.de>
! Copyright (c) 2016 Verena Krupp <verena.krupp@uni-siegen.de>
! Copyright (c) 2016 Tobias Schneider <tobias1.schneider@student.uni-siegen.de>
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions are met:
!
! 1. Redistributions of source code must retain the above copyright notice, this
! list of conditions and the following disclaimer.
!
! 2. Redistributions in binary form must reproduce the above copyright notice,
! this list of conditions and the following disclaimer in the documentation
! and/or other materials provided with the distribution.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "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 THE COPYRIGHT HOLDER 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.
!--------------------------------------------
!    A O S - Array of structures layout new
!-------------------------------------------
! Access to get_point value output
! Access to get_element value output

! *************************************************************************** !
!> This module provides a mechanism to define new variables in the simulation
!! by applying some operator on existing variables.
!!
module tem_operation_var_module
  use, intrinsic :: iso_c_binding,  only: c_ptr, c_loc, c_f_pointer, &
    &                                     c_null_ptr

  ! include treelm modules
  use env_module,                   only: rk, long_k, labelLen
  use tem_varSys_module,            only: tem_varSys_type,               &
    &                                     tem_varSys_check_inArgs,       &
    &                                     tem_get_element_chunk,         &
    &                                     tem_get_point_chunk,           &
    &                                     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_solverData_evalElem_type
  use tem_variable_module,          only: tem_variable_type
  use tem_varMap_module,            only: tem_varMap_type, tem_create_varMap
  use tem_time_module,              only: tem_time_type
  use treelmesh_module,             only: treelmesh_type
  use tem_logging_module,           only: logUnit, llerror
  use tem_aux_module,               only: tem_abort
  use tem_dyn_array_module,         only: PositionOfVal
  use tem_grow_array_module,        only: append, init, &
    &                                     truncate, grw_intArray_type, &
    &                                     grw_labelArray_type
  use tem_operation_module,         only: tem_varSys_op_data_type, &
    &                                     tem_indexLvl_type
  use tem_reduction_transient_module, only: tem_reduction_transient_type,    &
    &                                       tem_reduction_transient_init,    &
    &                                       tem_reduction_transient_reset,   &
    &                                       tem_reduction_transient_update,  &
    &                                       tem_reduction_transient_getElement
  use tem_logical_operation_var_module,                                        &
    &                              only: evalLogicalAnd_forPoint,              &
    &                                    evalLogicalOr_forPoint,               &
    &                                    evalLogicalGreater_forPoint,          &
    &                                    evalLogicalGreaterOrEqual_forPoint,   &
    &                                    evalLogicalLess_forPoint,             &
    &                                    evalLogicalLessOrEqual_forPoint,      &
    &                                    evalLogicalEqual_forPoint,            &
    &                                    evalLogicalNotEqual_forPoint,         &
    &                                    evalLogicalAnd_forElement,            &
    &                                    evalLogicalOr_forElement,             &
    &                                    evalLogicalGreater_forElement,        &
    &                                    evalLogicalGreaterOrEqual_forElement, &
    &                                    evalLogicalLess_forElement,           &
    &                                    evalLogicalLessOrEqual_forElement,    &
    &                                    evalLogicalEqual_forElement,          &
    &                                    evalLogicalNotEqual_forElement,       &
    &                                    evalLogicalAnd_fromIndex,             &
    &                                    evalLogicalOr_fromIndex,              &
    &                                    evalLogicalGreater_fromIndex,         &
    &                                    evalLogicalGreaterOrEqual_fromIndex,  &
    &                                    evalLogicalLess_fromIndex,            &
    &                                    evalLogicalLessOrEqual_fromIndex,     &
    &                                    evalLogicalEqual_fromIndex,           &
    &                                    evalLogicalNotEqual_fromIndex

  implicit none

  private

  public :: tem_divideVecByScal_forPoint
  public :: tem_divideVecByScal_fromIndex
  public :: tem_varSys_append_operVar
  public :: tem_opVar_setupIndices
  public :: tem_opVar_fill_inputIndex
  public :: tem_varSys_op_data_type
  public :: tem_get_new_varSys_data_ptr
  public :: tem_free_varSys_data_ptr
  public :: tem_evalMag_forPoint, tem_evalMag_forElement, tem_evalMag_fromIndex
  public :: tem_evalAdd_forPoint, tem_evalAdd_forElement, tem_evalAdd_fromIndex
  public :: tem_opVar_getParams, tem_opVar_setParams

  public :: tem_opVar_reduction_transient_init
  public :: tem_opVar_reduction_transient_update


contains


  ! ************************************************************************** !
  !> Routine to get a pointer to a new instance of method_data for an operation
  !! variable
  function tem_get_new_varSys_data_ptr(solver_bundle) result(resPtr)
    ! ---------------------------------------------------------------------- !
    !> Pointer to the newly created instance.
    type(c_ptr) ::  resPtr
    !> Optional solver data to store in tem_varSys_op_data_type
    type(c_ptr), optional, intent(in) :: solver_bundle
    ! ---------------------------------------------------------------------- !
    !> Local variable to allocate a new instance.
    type(tem_varSys_op_data_type), pointer :: res
    ! ---------------------------------------------------------------------- !

    allocate(res)
    if (present(solver_bundle)) res%solver_bundle = solver_bundle
    resPtr = c_loc(res)

  end function tem_get_new_varSys_data_ptr
  ! ************************************************************************** !

  ! ************************************************************************** !
  !> Free a method data structure again.
  subroutine tem_free_varSys_data_ptr(vardat_ptr)
    ! ---------------------------------------------------------------------- !
    !> Data pointer to free
    type(c_ptr), intent(inout) :: vardat_ptr
    ! ---------------------------------------------------------------------- !
    type(tem_varSys_op_data_type), pointer :: vardat
    ! ---------------------------------------------------------------------- !

    call c_f_pointer(vardat_ptr, vardat)
    deallocate(vardat)
    vardat_ptr = c_null_ptr

  end subroutine tem_free_varSys_data_ptr
  ! ************************************************************************** !


  ! ************************************************************************** !
  !> subroutine to add the variables from the input lua script to the varsys
  subroutine tem_varSys_append_operVar( operVar, varSys, pos,     &
    &                                   solverData_evalElem       )
    !--------------------------------------------------------------------------
    !> variables defined in the lua file
    type(tem_variable_type), intent(in)           :: operVar

    !> global variable system to which operVar to be appended
    type(tem_varSys_type), intent(inout)            :: varSys

    !> Position of the variable in the system.
    integer, optional, intent(out) :: pos

    !> A setter routine that allows the caller to define routine for the
    !! construction of an element representation.
    type(tem_varSys_solverData_evalElem_type), &
      &  optional, intent(in) :: solverData_evalElem
    !--------------------------------------------------------------------------
    integer :: addedPos
    integer :: nComps, nInputs
    integer, allocatable :: inPos(:)
    logical :: wasAdded
    character(len=labelLen), allocatable :: input_varname(:)
    integer, allocatable :: input_varIndex(:)
    logical :: isSatisfied
    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
    type(tem_varSys_op_data_type), pointer :: opData
    !--------------------------------------------------------------------------
    nullify(get_point, get_element, set_params, get_params, setup_indices, &
      &     get_valOfIndex)

    nComps = operVar%nComponents

    nInputs = size(operVar%input_varname)
    allocate(input_varname(nInputs))
    input_varname = operVar%input_varname

    ! check if input_varnames satisfy requirements for operType
    ! and correct user defined nComps if it does not match operType
    call check_opVar_prerequisites( operType      = operVar%operType, &
      &                             nInputs       = nInputs,          &
      &                             input_varname = input_varname,    &
      &                             varSys        = varSys,           &
      &                             nComps        = nComps,           &
      &                             isSatisfied   = isSatisfied       )

    !> If not satisfied then it is not possible to append current variable
    !! to varSys
    if (.not. isSatisfied) then
      write(logUnit(1),*) 'WARNING: input varnames does not satisfy'
      write(logUnit(1),*) '         requirements for operType '&
        &               //trim(operVar%operType)
      write(logUnit(1),*) 'Variable: "'//trim(operVar%label) &
        &              //'" is not appended.'
      return
    end if

    ! for operation variables, set_params, get_params and setup_indices are same
    ! since they send information to depend variables
    set_params => tem_opVar_setParams
    get_params => tem_opVar_getParams
    setup_indices => tem_opVar_setupIndices
    ! Get method data container to store indices for getValOfIndex
    ! Overwrite this method data with solver method data if operation
    ! is solver-specific
    method_data = tem_get_new_varSys_data_ptr()
    call C_F_POINTER(method_data, opData)

    select case(trim(operVar%operType))
    ! magnitude, division, multiplication are solver specific
    ! so set using getEvalFuncionsCallback
    case( 'difference' )
      get_element => evalDiff_forElement
      get_point => evalDiff_forPoint
      get_valOfIndex => evalDiff_fromIndex

    case( 'rel_difference' )
      get_element => evalRelDiff_forElement
      get_point => evalRelDiff_forPoint
      get_valOfIndex => evalRelDiff_fromIndex

    case( 'addition' )
      get_element => tem_evalAdd_forElement
      get_point => tem_evalAdd_forPoint
      get_valOfIndex => tem_evalAdd_fromIndex

    case ('multiplication')
      get_point => evalMultiply_forPoint
      get_element => evalMultiply_forElement
      get_valOfindex => evalMultiply_fromIndex

    case( 'multiply_scalar_times_vector' )
      get_point => tem_multiplyScalTimesVec_forPoint
      get_element => multiplyScalTimesVec_forElement
      get_valOfindex => tem_multiplyScalTimesVec_fromIndex

    case( 'division', 'div' )
      get_point => division_forPoint
      get_element => division_forElement
      get_valOfindex => division_fromIndex

    case( 'divide_vector_by_scalar' )
      get_point => tem_divideVecByScal_forPoint
      get_element => divideVecByScal_forElement
      get_valOfindex => tem_divideVecByScal_fromIndex

    case( 'gradient', 'grad', 'gradientX', 'gradX','gradientY','gradY', &
      &   'gradientZ', 'gradZ'                                          )
    ! Pointers set by the solvers using opVar_setter callback, see below

    case( 'magnitude' )
      get_point => tem_evalMag_forPoint
      get_element => tem_evalMag_forElement
      get_valOfindex => tem_evalMag_fromIndex

    case( 'extract' )
      allocate( input_varIndex(size(operVar%input_varIndex)) )
      input_varIndex = operVar%input_varIndex

      get_element => extract_forElement
      get_point => extract_forPoint
      get_valOfIndex => extract_fromIndex

    case( 'combine' )
      get_element => combine_forElement
      get_point => combine_forPoint
      get_valOfIndex => combine_fromIndex

    case( 'greater_than', 'gt', '>' )
      get_point => evalLogicalGreater_forPoint
      get_element => evalLogicalGreater_forElement
      get_valOfIndex => evalLogicalGreater_fromIndex

    case( 'greater_than_or_equal', 'ge', '>=' )
      get_point => evalLogicalGreaterOrEqual_forPoint
      get_element => evalLogicalGreaterOrEqual_forElement
      get_valOfIndex => evalLogicalGreaterOrEqual_fromIndex

    case( 'less_than', 'lt', '<' )
      get_point => evalLogicalLess_forPoint
      get_element => evalLogicalLess_forElement
      get_valOfIndex => evalLogicalLess_fromIndex

    case( 'less_than_or_equal', 'le', '<=' )
      get_point => evalLogicalLessOrEqual_forPoint
      get_Element => evalLogicalLessOrEqual_forElement
      get_valOfIndex => evalLogicalLessOrEqual_fromIndex

    case( 'equal', 'eq', '=' )
      get_point => evalLogicalEqual_forPoint
      get_Element => evalLogicalEqual_forElement
      get_valOfIndex => evalLogicalEqual_fromIndex

    case( 'not_equal', 'ne', '/=' )
      get_point => evalLogicalNotEqual_forPoint
      get_Element => evalLogicalNotEqual_forElement
      get_valOfIndex => evalLogicalNotEqual_fromIndex

    case( 'and' )
      get_point => evalLogicalAnd_forPoint
      get_Element => evalLogicalAnd_forElement
      get_valOfIndex => evalLogicalAnd_fromIndex

    case( 'or' )
      get_point => evalLogicalOr_forPoint
      get_Element => evalLogicalOr_forElement
      get_valOfIndex => evalLogicalOr_fromIndex

    case('reduction_transient')
      opData%redTrans%config = operVar%redTransConfig
      get_point => reductionTransient_forPoint
      get_Element => reductionTransient_forElement
      get_valOfIndex => reductionTransient_fromIndex

    case default
      if (.not. (associated(get_point)      &
        &   .or. associated(get_element)    &
        &   .or. associated(set_params)     &
        &   .or. associated(get_params)     &
        &   .or. associated(setup_indices)  &
        &   .or. associated(get_valOfIndex))) then
        write(logUnit(3),*) 'operType: '                 &
          & // trim(operVar%operType)                    &
          & // ' not supported. Variable is not appended.'
        return ! go to next variable
      end if
    end select

    ! Workaround for Intel 15 compiler
    if ( .not. allocated(input_varname) ) then
      allocate( input_varname(0) )
    end if

    if ( .not. allocated(input_varIndex) ) then
      allocate( input_varIndex(0) )
    end if

    ! append variable to varSys
    call tem_varSys_append_derVar(          &
      &  me             = varSys,           &
      &  varName        = operVar%label,    &
      &  operType       = operVar%operType, &
      &  nComponents    = nComps,           &
      &  input_varname  = input_varname,    &
      &  input_varIndex = input_varIndex,   &
      &  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
      if (present(solverData_evalElem)) then
        ! If an solverData_evalElem operation is provided, override the
        ! get_element pointer and use the provided setter solverData_evalElem
        ! instead to define the get_element routine.
        call solverData_evalElem%opVar_setter(varSys%method%val(addedPos))
      end if
      write(logUnit(9),*) 'Successfully appended variable "' &
        & // trim(operVar%label)// '" to the variable system'
    else if (addedpos < 1) then
      write(logUnit(7),*) 'WARNING: variable '//trim(operVar%label)// &
        &                 ' is not added to variable system'
    end if

    if (present(pos)) pos = addedPos

    ! deallocate here to be allocated for next variable
    deallocate(input_varname)
    if (allocated(input_varIndex)) deallocate(input_varIndex)
    if (allocated(inPos)) deallocate(inPos)

  end subroutine tem_varSys_append_operVar
  ! ************************************************************************** !


  ! ************************************************************************** !
  !> This subroutine checks whether input variables satisfy requirements for
  !! opertype.
  !!
  !! For example: nComponents, number of inputs etc.
  !! If input_varname not found in varSys then this function returns false.
  !! If user defined nComps = -1 then nComps is set according to operType.
  !!
  !! This subroutine checks operations used in both treelm and solvers
  subroutine check_opVar_prerequisites( operType, nInputs, input_varname, &
    &                                   varSys, nComps, isSatisfied )
    !--------------------------------------------------------------------------!
    !> Operation type
    character(len=*), intent(in) :: operType
    !> Number of inputs
    integer, intent(in) :: nInputs
    !> Input varnames for current operation
    character(len=*), intent(in) :: input_varname(nInputs)
    !> Variable system to look for input_varname
    type(tem_varSys_type), intent(in) :: varSys
    !> Number of components defined for operation variable.
    !! If nComps == -1 then current nComps is set here
    integer, intent(inout) :: nComps
    !> true if all requirements for opertype are satisfied
    logical, intent(out) :: isSatisfied
    !--------------------------------------------------------------------------!
    integer :: iIn, total_input_nComps
    integer :: inpos(nInputs), input_nComps(nInputs)
    !--------------------------------------------------------------------------!
    write(logUnit(7),*) 'Checking prerequisites for opertype: '//trim(operType)
    isSatisfied = .true.

    ! Position of input variables in varSys
    do iIn = 1, nInputs
      inPos(iIn) = PositionofVal(varSys%varname, input_varname(iIn))
    end do

    if ( all(inPos > 0) ) then
      input_nComps(:) = varSys%method%val(inPos(:))%nComponents
    else
      isSatisfied = .false.
      write(logUnit(1),*) 'Error: input varnames not found in varsys'
      return
    end if

    select case (trim(operType))
    case ('addition', 'difference', 'rel_difference', 'multiplication', &
      &   'division')
      ! operations which require two inputs and both to have same number of
      ! components
      if (size(input_varname) /= 2) then
        write(logUnit(1),*) 'Error: In operation type: '//trim(operType)
        write(logUnit(1),*) 'Number of input_varname /= 2'
        call tem_abort()
      end if

      if (nComps == -1) then
        nComps = input_nComps(1)
        write(logUnit(7),*) 'INFO: nComponents is not defined by user.'
        write(logUnit(7),*)  '     nComps set to nComps of 1st input variable'
      end if

      if (.not. all(input_nComps == nComps)) then
        write(logUnit(1),*) 'Error: nComps of operation variable does not ' &
          &               //'match with input variables nComps'
        write(logUnit(1),*) 'Input nComps: ', input_nComps
        call tem_abort()
      end if

    case ('multiply_scalar_times_vector')
      ! operations which require two inputs and 1st input must be a scalar
      if (size(input_varname) /= 2) then
        write(logUnit(1),*) 'Error: In operation type: '//trim(operType)
        write(logUnit(1),*) 'Number of input_varname /= 2'
        call tem_abort()
      end if

      ! 1st input variable must be scalar
      if (input_nComps(1) /= 1) then
        write(logUnit(1),*) 'Error: In operation type: '//trim(operType)
        write(logUnit(1),*) '1st input variable is not a scalar'
        call tem_abort()
      end if

      if (nComps == -1 .or. nComps /= input_nComps(2)) then
        write(logUnit(1),*) 'Warning: nComps for operVar is wrong.'
        write(logUnit(1),*) '         Setting nComps =', input_nComps(2)
        nComps = input_nComps(2)
      end if

    case ('divide_vector_by_scalar')
      ! operations which require two inputs and 1st input must be a scalar
      if (size(input_varname) /= 2) then
        write(logUnit(1),*) 'Error: In operation type: '//trim(operType)
        write(logUnit(1),*) 'Number of input_varname /= 2'
        call tem_abort()
      end if

      ! 2nd input variable must be scalar
      if (input_nComps(2) /= 1) then
        write(logUnit(1),*) 'Error: In operation type: '//trim(operType)
        write(logUnit(1),*) '1st input variable is not a scalar'
        call tem_abort()
      end if

      if (nComps == -1 .or. nComps /= input_nComps(1)) then
        write(logUnit(1),*) 'Warning: nComps for operVar is wrong.'
        write(logUnit(1),*) '         Setting nComps =', input_nComps(1)
        nComps = input_nComps(1)
      end if

    case ('magnitude')
      ! operations which require only one input and ncomponents must be 1
      if (nInputs /= 1 ) then
        write(logUnit(1),*) 'Error: In operation type: '//trim(operType)
        write(logUnit(1),*) 'Number of input_varname /= 1'
        call tem_abort()
      end if

      if ( nComps /= 1 ) then
        write(logUnit(1),*) 'Warning: nComps /= 1. Setting nComps = 1'
        nComps = 1
      end if

    case ('meansquare')
      ! operations which require only one input and ncomps == input_ncomps
      if (nInputs /= 1 ) then
        write(logUnit(1),*) 'Error: In operation type: '//trim(operType)
        write(logUnit(1),*) 'Number of input_varname /= 1'
        call tem_abort()
      end if

      if ( nComps /= input_nComps(1) ) then
        write(logUnit(1),*) 'Warning: nComps /= input, should be:', &
          &                 input_nComps(1)
        write(logUnit(1),*) 'but is:', nComps
        write(logUnit(1),*) 'Setting it to ', input_nComps(1)
        nComps = input_nComps(1)
      end if

    case ('locall2mean','deviation')
      ! operations which require only one input and ncomps == input_ncomps
      if (nInputs /= 1 ) then
        write(logUnit(1),*) 'Error: In operation type: '//trim(operType)
        write(logUnit(1),*) 'Number of input_varname /= 1'
        call tem_abort()
      end if

      if ( nComps /= input_nComps(1) ) then
        write(logUnit(1),*) 'Warning: nComps /= input, should be:', &
          &                 input_nComps(1)
        write(logUnit(1),*) 'but is:', nComps
        write(logUnit(1),*) 'Setting it to ', input_nComps(1)
        nComps = input_nComps(1)
      end if


    case ('extract', 'gradient', 'gradientX', 'gradientY', 'gradientZ', &
      &   'reduction_transient'                                         )
      ! operations which require only one input and ncomponents can be >= 1
      if (nInputs /= 1 ) then
        write(logUnit(1),*) 'Error: In operation type: '//trim(operType)
        write(logUnit(1),*) 'Number of input_varname /= 1'
        call tem_abort()
      end if

      if (nComps == -1) then
        write(logUnit(1),*) 'Error: nComps for operVar is not defined.'
        call tem_abort()
      end if

    case ('combine')
      ! operations which require more than two inputs and
      ! nComps is sum of nComps of input vars
      total_input_nComps = sum(input_nComps)
      if (nComps /= total_input_nComps) then
        write(logUnit(llerror),*) 'Warning: In variable operation combine'
        write(logUnit(llerror),*) '  User defined nComps ' &
          &                    // '/= sum(input_nComponents).'
        write(logUnit(llerror),*) &
          &  '  So, setting nComps = sum(input_nComponents)'
        write(logUnit(llerror),*) '  nComps: ', total_input_nComps
        nComps = total_input_nComps
      end if

    case( 'greater_than', 'gt', '>',         &
      & 'greater_than_or_equal', 'ge', '>=', &
      & 'less_than', 'lt', '<',              &
      & 'less_than_or_equal', 'le', '<=',    &
      & 'equal', 'eq', '=',                  &
      & 'not_equal', 'ne', '/=',             &
      & 'and',                               &
      & 'or'                                 )
      if (size(input_varname) /= 2) then
        write(logUnit(1),*) 'Error: In operation type: '//trim(operType)
        write(logUnit(1),*) 'Number of input_varname /= 2'
        call tem_abort()
      end if

    case default
      write(logUnit(1),*) 'ERROR: operType: ' // trim(operType) &
        & // ' not supported. Variable is not appended.'
      call tem_abort()
    end select

  end subroutine check_opVar_prerequisites
  ! ************************************************************************** !

  ! ************************************************************************** !
  !> Initialize time reduction operation variable
  !! Loop over all variable in varSys and allocate redTrans%val for
  !! reduction_transient operation variable with nElems
  subroutine tem_opVar_reduction_transient_init(varSys, tree, redTransVarMap,&
    &                                           nDofs, time)
    ! -------------------------------------------------------------------------
    !> Global variable system
    type(tem_varSys_type), intent(in) :: varSys
    !> treelmesh_type
    type(treelmesh_type), intent(in) :: tree
    !> position of time reduction variable in varSys
    type(tem_varMap_type), intent(out) :: redTransVarMap
    !> Solver nDegrees of freedom
    integer, intent(in), optional :: nDofs
    !> Current time
    type(tem_time_type), intent(in) :: time
    ! -------------------------------------------------------------------------
    type(tem_varSys_op_data_type), pointer :: opData
    integer :: iVar, iElem, varPos, nDofs_loc, posDepVar, nCompMax, idxMax
    integer :: nRedVars
    type(grw_labelArray_type) :: redTransVarName
    integer :: elemPos(tree%nElems)
    real(kind=rk), allocatable :: input_varRes(:)
    ! -------------------------------------------------------------------------
    if (present(nDofs)) then
      nDofs_loc = nDofs
    else
      nDofs_loc = 1
    end if

    ! Gather list of variable names which has reduction_transient operation
    call init(redTransVarName)
    do iVar = 1, varSys%varName%nVals
      if (trim(varSys%method%val(iVar)%operType) == 'reduction_transient') then
        call append(me=redTransVarName, val=varSys%varName%val(iVar))
      end if
    end do

    ! create varMap to store position of reduction_transient variable in varSys
    call tem_create_varMap(                                                 &
      &             varName = redTransVarName%val(1:redTransVarName%nVals), &
      &             varSys  = varSys,                                       &
      &             varMap  = redTransVarMap                                )

    elemPos(1:tree%nElems) = (/ (iElem, iElem=1, tree%nElems) /)
    nRedVars =  redTransVarMap%varPos%nVals
    nCompMax = maxval(varSys%method%val(redTransVarMap%varPos%val(1:nRedVars)) &
      &                            %nComponents)
    allocate(input_varRes(tree%nElems*nCompMax*nDofs_loc))

    ! Initialize time reduction
    do iVar = 1, redTransVarMap%varPos%nVals
      varPos = redTransVarmap%varPos%val(iVar)
      call C_F_POINTER(varSys%method%val(varPos)%method_data, opData)

      call tem_reduction_transient_init(                            &
        &      me          = opData%redTrans,                       &
        &      nElems      = tree%nElems,                           &
        &      nComponents = varSys%method%val(varPos)%nComponents, &
        &      nDofs       = nDofs_loc                              )
      ! Fill last
      posDepVar = varSys%method%val(varPos)%input_varPos(1)
      call varSys%method%val(posDepVar)%get_element(                        &
        &                                   varSys  = varSys,               &
        &                                   elemPos = elemPos,              &
        &                                   time    = time,                 &
        &                                   tree    = tree,                 &
        &                                   nElems  = tree%nElems,          &
        &                                   nDofs   = nDofs_loc,            &
        &                                   res     = input_varRes(:)       )

      idxMax = opData%redTrans%nEntries
      opData%redTrans%val(:, opData%redTrans%last) = input_varRes(1:idxMax)
    end do

  end subroutine tem_opVar_reduction_transient_init
  ! ************************************************************************** !

  ! ************************************************************************** !
  !> Update all time reduction operation variables for entire domain
  subroutine tem_opVar_reduction_transient_update(redTransVarPos, varSys, tree,&
    &                                             time)
    ! ---------------------------------------------------------------------- !
    !> Position of time reduction variables in varSys
    integer, intent(in) :: redTransVarPos(:)
    !> Global Variable system
    type(tem_varSys_type), intent(in) :: varSys
    !> treelmesh_type
    type(treelmesh_type), intent(in) :: tree
    !> Current time
    type(tem_time_type), intent(in) :: time
    ! ---------------------------------------------------------------------- !
    integer :: elemPos(tree%nElems)
    real(kind=rk), allocatable :: input_varRes(:)
    integer :: iVar, iElem, posDepVar, nCompMax, idxMax, nDofs
    type(tem_varSys_op_data_type), pointer :: opData
    ! ---------------------------------------------------------------------- !

    ! Only need to do anything here if at least one variable with reduction
    ! is to be computed.
    if (size(redTransVarPos) < 1) RETURN

    ! nDofs of solver are stored in opData%redTrans which is same for all
    ! variables so get from any reduction_transient variable
    call C_F_Pointer(varSys%method%val(redTransVarPos(1))%method_data, &
      &              opData                                            )
    nDofs = opData%redTrans%nDofs

    elemPos(1:tree%nElems) = (/ (iElem, iElem=1, tree%nElems) /)

    nCompMax = maxval(varSys%method%val(redTransVarPos(:))%nComponents)
    allocate(input_varRes(tree%nElems*nCompMax*nDofs))

    do iVar = 1, size(redTransVarPos)
      call C_F_Pointer(varSys%method%val(redTransVarPos(iVar))%method_data, &
        &              opData                                               )

      posDepVar = varSys%method%val(redTransVarPos(iVar))%input_varPos(1)
      call varSys%method%val(posDepVar)%get_element(                        &
        &                                   varSys  = varSys,               &
        &                                   elemPos = elemPos,              &
        &                                   time    = time,                 &
        &                                   tree    = tree,                 &
        &                                   nElems  = tree%nElems,          &
        &                                   nDofs   = nDofs,                &
        &                                   res     = input_varRes(:)       )

      idxMax = opData%redTrans%nEntries

      call tem_reduction_transient_update(me     = opData%redTrans,       &
        &                                 res    = input_varRes(1:idxMax) )
    end do
    deallocate(input_varRes)
  end subroutine tem_opVar_reduction_transient_update
  ! ************************************************************************** !


  ! ************************************************************************** !
  !> Evaluate the function pointers of the dependent variables,
  !! and then calculate the difference between these two. ( scalar or vector )
  !! In lua file, first define new variable with varType operation kind as
  !! "difference" and provide two dependent variable via input_varname.
  !! If input_varname variable is not part of predefined solver variables then
  !! add also that variable via variable table for example spacetime function
  !! variable.
  !! For example: Define a variable called difference, which depend on density
  !! and spacetime. one can get an error between simulation
  !! results and analytical solution.
  !!
  !! \verbatim
  !! -- in lua file, one can define as following:
  !! variable = {{
  !!   name = 'dens_reference',
  !!   ncomponents = 1,
  !!   vartype = "st_fun",
  !!   st_fun =  luaFun
  !!   },
  !!   {
  !!   name = 'dens_difference',
  !!   ncomponents = 1,
  !!   vartype = "operation",
  !!   operation = {kind='difference', input_varname={density, dens_reference}}
  !!   }
  !!   ...
  !!  }
  !! tracking = {
  !!   variable = {'dens_difference'},
  !!   folder = 'tracking/',
  !!   shape = {kind = 'canoND', object = {origin = {3.0,3.1,3.0} } },
  !!   format = 'ascii',
  !!   time = {min = 0, max = tmax, interval = 1},
  !!  }
  !! \endverbatim
  !!
  !! The interface has to comply to the abstract interface
  !! tem_varSys_module#tem_varSys_proc_element.
  !!
  recursive subroutine evalDiff_forElement( fun, varsys, elempos, time, tree, &
    &                                       nElems, nDofs, 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

    !> Position of the TreeID of the element to get the variable for in the
    !! global treeID list.
    integer, intent(in) :: elempos(:)

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

    !> global treelm mesh info
    type(treelmesh_type), intent(in) :: tree

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

    !> Number of degrees of freedom within an element.
    integer, intent(in) :: nDofs

    !> Resulting values for the requested variable.
    !!
    !! Linearized array dimension:
    !! (n requested entries) x (nComponents of this variable)
    !! x (nDegrees of freedom)
    !! Access: (iElem-1)*fun%nComponents*nDofs +
    !!         (iDof-1)*fun%nComponents + iComp
    real(kind=rk), intent(out) :: res(:)
    ! ---------------------------------------------------------------------- !
    integer :: iDep, posDepVar
    ! first dimension:nElems,
    ! second dimension:size of difference variable in input_varname
    real(kind=rk), allocatable :: input_varRes(:,:)
    ! ---------------------------------------------------------------------- !
    integer :: nTotal
    ! ---------------------------------------------------------------------- !

    nTotal = nElems*nDofs*fun%nComponents

    ! nInputs must be two
    allocate(input_varRes(nTotal, fun%nInputs))

    do iDep = 1, fun%nInputs
      ! get position of dependent var in varSys
      posDepVar = fun%input_varPos(iDep)

      ! derive dependent variable
      call varSys%method%val(posDepVar)%get_element(                        &
        &                                   varSys  = varSys,               &
        &                                   elemPos = elemPos,              &
        &                                   time    = time,                 &
        &                                   tree    = tree,                 &
        &                                   nElems  = nElems,               &
        &                                   nDofs   = nDofs,                &
        &                                   res     = input_varRes(:, iDep) )

    end do

    res(:nTotal) = input_varRes(:, 1) - input_varRes(:, 2)

    deallocate(input_varRes)

  end subroutine evalDiff_forElement
  ! ************************************************************************** !


  ! ************************************************************************** !
  !> Same as evalDiff_forElement except it evaluate diff from points
  !!
  !! The interface has to comply to the abstract interface
  !! tem_varSys_module#tem_varSys_proc_point
  recursive subroutine evalDiff_forPoint( fun, varsys, point, time, tree, &
    &                                     nPnts, 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

    !> Three-dimensional coordinates at which the variable should be
    !! evaluated. Only useful for variables provided as space-time functions.
    real(kind=rk), intent(in) :: point(:,:)

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

    !> global treelm mesh info
    type(treelmesh_type), intent(in) :: tree

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

    !> 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 :: iDep, posDepVar
    integer :: nTotal
    ! first dimension:nElems,
    ! second dimension:size of difference variable in input_varname
    real(kind=rk), allocatable :: input_varRes(:,:)
    ! ---------------------------------------------------------------------- !

    nTotal = nPnts*fun%nComponents

    ! nInputs must be two
    allocate(input_varRes(nPnts*fun%nComponents, fun%nInputs))

    do iDep = 1, fun%nInputs
      ! get position of dependent var in varSys
      posDepVar = fun%input_varPos(iDep)

      ! derive dependent variable
      call varSys%method%val(posDepVar)%get_point(                         &
        &                                   varSys = varSys,               &
        &                                   point  = point,                &
        &                                   time   = time,                 &
        &                                   tree   = tree,                 &
        &                                   nPnts  = nPnts,                &
        &                                   res    = input_varRes(:, iDep) )

    end do

    res(:nTotal) = input_varRes(:, 1) - input_varRes(:, 2)

    deallocate(input_varRes)

  end subroutine evalDiff_forPoint
  ! ************************************************************************** !


  ! ************************************************************************** !
  !> Same as evalDiff_forPoint except it evaluate diff from points via index
  !! which are setup before hand
  !!
  !! The interface has to comply to the abstract interface
  !! tem_varSys_module#tem_varSys_proc_point
  recursive subroutine evalDiff_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: nVals
    integer, intent(in) :: idx(:)

    !> With idx as start index in contiguous memory,
    !! idxLength defines length of each contiguous memory
    !! Size: 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 :: iDep, posDepVar
    ! first dimension:nElems,
    ! second dimension:size of difference variable in input_varname
    real(kind=rk), allocatable :: input_varRes(:,:)
    ! ---------------------------------------------------------------------- !
    type(tem_varSys_op_data_type), pointer :: opData
    integer :: nTotal
    ! ---------------------------------------------------------------------- !

    call C_F_POINTER( fun%method_Data, opData )

    ! check if number of index are the same as number of values asked for
    call tem_varSys_check_inArgs( fun, varSys, time, iLevel, idx, idxLen, &
      &    nVals, label = 'evalDiff_fromIndex'                            )

    nTotal = nVals*fun%nComponents
    ! nInputs must be two
    allocate(input_varRes(nTotal, fun%nInputs))

    do iDep = 1, fun%nInputs
      ! get position of dependent var in varSys
      posDepVar = fun%input_varPos(iDep)

      ! derive dependent variable for indices
      !rent derive dependent variable
      call varSys%method%val(posDepVar)%get_valOfIndex( &
          & varSys  = varSys,                           &
          & time    = time,                             &
          & iLevel  = iLevel,                           &
          & idx     = opData%input_pntIndex(iDep)       &
          &           %indexLvl(iLevel)%val( idx(:) ),  &
          & nVals   = nVals,                            &
          & res     = input_varRes(:,iDep)              )
    end do

    res(:nTotal) = input_varRes(:, 1) - input_varRes(:, 2)

    deallocate(input_varRes)

  end subroutine evalDiff_fromIndex
  ! ************************************************************************** !


  ! ************************************************************************** !
  !> Same as evalDiff but output of 2nd dependent variable is used to compute
  !! relative difference.
  !!
  !! The interface has to comply to the abstract interface
  !! tem_varSys_module#tem_varSys_proc_element.
  !!
  recursive subroutine evalRelDiff_forElement( fun, varsys, elempos, time, &
    &                                          tree, nElems, nDofs, 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

    !> Position of the TreeID of the element to get the variable for in the
    !! global treeID list.
    integer, intent(in) :: elempos(:)

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

    !> global treelm mesh info
    type(treelmesh_type), intent(in) :: tree

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

    !> Number of degrees of freedom within an element.
    integer, intent(in) :: nDofs

    !> Resulting values for the requested variable.
    !!
    !! Linearized array dimension:
    !! (n requested entries) x (nComponents of this variable)
    !! x (nDegrees of freedom)
    !! Access: (iElem-1)*fun%nComponents*nDofs +
    !!         (iDof-1)*fun%nComponents + iComp
    real(kind=rk), intent(out) :: res(:)
    ! ---------------------------------------------------------------------- !
    integer :: iDep, posDepVar
    integer :: nTotal
    ! first dimension:nElems,
    ! second dimension:size of difference variable in input_varname
    real(kind=rk), allocatable :: input_varRes(:,:)
    ! ---------------------------------------------------------------------- !

    nTotal = nElems*nDofs*fun%nComponents
    ! nInputs must be two
    allocate(input_varRes(nTotal, fun%nInputs))

    if (nDofs > 1) then
      write(*,*) 'TODO: evalreldiff does not work for polynomial data yet'
      write(*,*) '      It makes use of a division, which can not directly'
      write(*,*) '      be done in modal space!'
      write(*,*) ''
      write(*,*) 'Need to replace this routine in Ateles!'
      write(*,*) 'Stopping now.'
      call tem_abort()
    end if

    do iDep = 1, fun%nInputs
      ! get position of dependent var in varSys
      posDepVar = fun%input_varPos(iDep)

      ! derive dependent variable
      call varSys%method%val(posDepVar)%get_element(                        &
        &                                   varSys  = varSys,               &
        &                                   elemPos = elemPos,              &
        &                                   time    = time,                 &
        &                                   tree    = tree,                 &
        &                                   nElems  = nElems,               &
        &                                   nDofs   = nDofs,                &
        &                                   res     = input_varRes(:, iDep) )

    end do

    res(:nTotal) = ( input_varRes(:, 1) - input_varRes(:, 2) ) &
      &          / input_varRes(:, 2)

    deallocate(input_varRes)

  end subroutine evalRelDiff_forElement
  ! ************************************************************************** !


  ! ************************************************************************** !
  !> Same as evalRelDiff_forElement except it evaluate relatice diff
  !! from points
  !!
  !! The interface has to comply to the abstract interface
  !! tem_varSys_module#tem_varSys_proc_point.
  recursive subroutine evalRelDiff_forPoint( fun, varsys, point, time, tree, &
    &                                        nPnts, 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

    !> Three-dimensional coordinates at which the variable should be
    !! evaluated. Only useful for variables provided as space-time functions.
    real(kind=rk), intent(in) :: point(:,:)

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

    !> global treelm mesh info
    type(treelmesh_type), intent(in) :: tree

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

    !> 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 :: iDep, posDepVar
    integer :: nTotal
    ! first dimension:nElems,
    ! second dimension:size of difference variable in input_varname
    real(kind=rk), allocatable :: input_varRes(:,:)
    ! ---------------------------------------------------------------------- !

    nTotal = nPnts*fun%nComponents
    ! nInputs must be two
    allocate(input_varRes(nTotal, fun%nInputs))

    do iDep = 1, fun%nInputs
      ! get position of dependent var in varSys
      posDepVar = fun%input_varPos(iDep)

      ! derive dependent variable
      call varSys%method%val(posDepVar)%get_point(                         &
        &                                   varSys = varSys,               &
        &                                   point  = point,                &
        &                                   time   = time,                 &
        &                                   tree   = tree,                 &
        &                                   nPnts  = nPnts,                &
        &                                   res    = input_varRes(:, iDep) )

    end do

    res(:nTotal) = ( input_varRes(:, 1) - input_varRes(:, 2) ) &
      &            / input_varRes(:, 2)

    deallocate(input_varRes)

  end subroutine evalRelDiff_forPoint
  ! ************************************************************************** !


  ! ************************************************************************** !
  !> Same as evalRelDiff_forElement except it evaluate relative diff
  !! from points via index which are setup before
  !!
  !!
  recursive subroutine evalRelDiff_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: nVals
    integer, intent(in) :: idx(:)

    !> With idx as start index in contiguous memory,
    !! idxLength defines length of each contiguous memory
    !! Size: 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 :: iDep, posDepVar
    ! first dimension:nElems,
    ! second dimension:size of difference variable in input_varname
    real(kind=rk), allocatable :: input_varRes(:,:)
    ! ---------------------------------------------------------------------- !
    type(tem_varSys_op_data_type), pointer :: opData
    integer :: nTotal
    ! ---------------------------------------------------------------------- !

    call C_F_POINTER( fun%method_Data, opData )

    nTotal = nVals*fun%nComponents

    call tem_varSys_check_inArgs( fun, varSys, time, iLevel, idx, idxLen, &
      &    nVals, label = 'evalRelDiff_fromIndex'                         )

    ! nInputs must be two
    allocate(input_varRes(nTotal, fun%nInputs))

    do iDep = 1, fun%nInputs
      ! get position of dependent var in varSys
      posDepVar = fun%input_varPos(iDep)

      ! derive dependent variable for indices
      call varSys%method%val(posDepVar)%get_valOfIndex( &
          & varSys  = varSys,                           &
          & time    = time,                             &
          & iLevel  = iLevel,                           &
          & idx     = opData%input_pntIndex(iDep)       &
          &           %indexLvl(iLevel)%val( idx(:) ),  &
          & nVals   = nVals,                            &
          & res     = input_varRes(:,iDep)              )
    end do

    res(:nTotal) = ( input_varRes(:, 1) - input_varRes(:, 2) ) &
      &           / input_varRes(:, 2)

    deallocate(input_varRes)

  end subroutine evalRelDiff_fromIndex
  ! ************************************************************************** !


  ! ************************************************************************** !
  !> Evaluate the function pointers of the dependent variables,
  !! and then calculate the addition or arbinary number of input variables
  !! with all input variables have same nComponents. ( scalar or vector )
  !! In lua file, first define new variable with varType operation kind as
  !! "addition" and provide two dependent variable via input_varname.
  !! If input_varname variable is not part of predefined solver variables then
  !! add also that variable via variable table for example spacetime function
  !! variable.
  !!
  !! \verbatim
  !! -- in lua file, one can define as following:
  !! variable = {{
  !!   name = 'dens_reference',
  !!   ncomponents = 1,
  !!   vartype = "st_fun",
  !!   st_fun =  luaFun
  !!   },
  !!   {
  !!   name = 'dens_difference',
  !!   ncomponents = 1,
  !!   vartype = "operation",
  !!   operation = {kind='addition', input_varname={density, dens_reference}}
  !!   }
  !!   ...
  !!  }
  !! \endverbatim
  !!
  !! The interface has to comply to the abstract interface
  !! tem_varSys_module#tem_varSys_proc_element.
  !!
  recursive subroutine tem_evalAdd_forElement( fun, varsys, elempos, time, &
    &                                          tree, nElems, nDofs, 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

    !> Position of the TreeID of the element to get the variable for in the
    !! global treeID list.
    integer, intent(in) :: elempos(:)

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

    !> global treelm mesh info
    type(treelmesh_type), intent(in) :: tree

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

    !> Number of degrees of freedom within an element.
    integer, intent(in) :: nDofs

    !> Resulting values for the requested variable.
    !!
    !! Linearized array dimension:
    !! (n requested entries) x (nComponents of this variable)
    !! x (nDegrees of freedom)
    !! Access: (iElem-1)*fun%nComponents*nDofs +
    !!         (iDof-1)*fun%nComponents + iComp
    real(kind=rk), intent(out) :: res(:)
    ! ---------------------------------------------------------------------- !
    integer :: iDep, posDepVar
    real(kind=rk), allocatable :: input_varRes(:)
    integer :: nTotal
    ! ---------------------------------------------------------------------- !

    nTotal = nElems*nDofs*fun%nComponents
    allocate(input_varRes(nTotal))

    res = 0.0_rk
    do iDep = 1, fun%nInputs
      ! get position of dependent var in varSys
      posDepVar = fun%input_varPos(iDep)

      ! derive dependent variable
      call varSys%method%val(posDepVar)%get_element(               &
        &                                   varSys  = varSys,      &
        &                                   elemPos = elemPos,     &
        &                                   time    = time,        &
        &                                   tree    = tree,        &
        &                                   nElems  = nElems,      &
        &                                   nDofs   = nDofs,       &
        &                                   res     = input_varRes )

      res(:nTotal) = res(:nTotal) + input_varRes
    end do

    deallocate(input_varRes)

  end subroutine tem_evalAdd_forElement
  ! ************************************************************************** !


  ! ************************************************************************** !
  !> Same as tem_evalAdd_forElement except it evaluate addition from points
  !!
  !! The interface has to comply to the abstract interface
  !! tem_varSys_module#tem_varSys_proc_point
  recursive subroutine tem_evalAdd_forPoint( fun, varsys, point, time, tree, &
    &                                        nPnts, 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

    !> Three-dimensional coordinates at which the variable should be
    !! evaluated. Only useful for variables provided as space-time functions.
    real(kind=rk), intent(in) :: point(:,:)

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

    !> global treelm mesh info
    type(treelmesh_type), intent(in) :: tree

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

    !> 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 :: iDep, posDepVar
    real(kind=rk) :: input_varRes(nPnts*fun%nComponents)
    ! ---------------------------------------------------------------------- !
    res = 0.0_rk
    do iDep = 1, fun%nInputs
      ! get position of dependent var in varSys
      posDepVar = fun%input_varPos(iDep)

      ! derive dependent variable
      call varSys%method%val(posDepVar)%get_point( varSys = varSys,      &
        &                                          point  = point,       &
        &                                          time   = time,        &
        &                                          tree   = tree,        &
        &                                          nPnts  = nPnts,       &
        &                                          res    = input_varRes )

      res(:nPnts*fun%nComponents) = res(:nPnts*fun%nComponents) &
        &                         + input_varRes
    end do

  end subroutine tem_evalAdd_forPoint
  ! ************************************************************************** !


  ! ************************************************************************** !
  !> Same as tem_evalAdd_forPoint except it evaluate addition from points via
  !! index which are setup before hand
  !!
  !! The interface has to comply to the abstract interface
  !! tem_varSys_module#tem_varSys_proc_point
  recursive subroutine tem_evalAdd_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: nVals
    integer, intent(in) :: idx(:)

    !> With idx as start index in contiguous memory,
    !! idxLength defines length of each contiguous memory
    !! Size: 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 :: iDep, posDepVar
    ! first dimension:nElems,
    ! second dimension:size of difference variable in input_varname
    real(kind=rk) :: input_varRes(nVals * fun%nComponents)
    ! ---------------------------------------------------------------------- !
    type(tem_varSys_op_data_type), pointer :: opData
    ! ---------------------------------------------------------------------- !

    call C_F_POINTER( fun%method_Data, opData )

    call tem_varSys_check_inArgs( fun, varSys, time, iLevel, idx, idxLen, &
      &    nVals, label = 'evalAdd_fromIndex'                             )

    res = 0.0_rk
    do iDep = 1, fun%nInputs
      ! get position of dependent var in varSys
      posDepVar = fun%input_varPos(iDep)

      ! derive dependent variable for indices
      !rent derive dependent variable
      call varSys%method%val(posDepVar)%get_valOfIndex( &
          & varSys  = varSys,                           &
          & time    = time,                             &
          & iLevel  = iLevel,                           &
          & idx     = opData%input_pntIndex(iDep)       &
          &           %indexLvl(iLevel)%val( idx(:) ),  &
          & nVals   = nVals,                            &
          & res     = input_varRes                      )
      res(:nVals*fun%nComponents) = res(:nVals*fun%nComponents) + input_varRes
    end do

  end subroutine tem_evalAdd_fromIndex
  ! ************************************************************************** !


  ! ************************************************************************** !
  !> Evaluate magnitude of any vectorial variable.
  !! In lua file, first define new variable with varType operation kind as
  !! "magnitude" and provide name of the variable from which magnitude
  !! to be derived in input_varname.
  !! If input_varname variable is not part of predefined solver variables then
  !! add also that variable via variable table.
  !!
  !! \verbatim
  !! -- in lua file, one can define as following:
  !! variable = {{
  !!   name = 'velMag',
  !!   ncomponents = 1,
  !!   vartype = "operation",
  !!   operation = {kind='magnitude',input_varname={'velocity'}}
  !!   },
  !!  }
  !! tracking = {
  !!   variable = {'velMag'},
  !!   folder = 'tracking/',
  !!   shape = {kind = 'canoND', object = {origin = {3.0,3.1,3.0} } },
  !!   format = 'ascii',
  !!   time = {min = 0, max = tmax, interval = 1},
  !!  }
  !! \endverbatim
  !!
  !! The interface has to comply to the abstract interface
  !! tem_varSys_module#tem_varSys_proc_element.
  !!
  recursive subroutine tem_evalMag_forElement( fun, varsys, elempos, time, &
    &                                          tree, nElems, nDofs, 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

    !> Position of the TreeID of the element to get the variable for in the
    !! global treeID list.
    integer, intent(in) :: elempos(:)

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

    !> global treelm mesh info
    type(treelmesh_type), intent(in) :: tree

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

    !> Number of degrees of freedom within an element.
    integer, intent(in) :: nDofs

    !> Resulting values for the requested variable.
    !!
    !! Linearized array dimension:
    !! (n requested entries) x (nComponents of this variable)
    !! x (nDegrees of freedom)
    !! Access: (iElem-1)*fun%nComponents*nDofs +
    !!         (iDof-1)*fun%nComponents + iComp
    real(kind=rk), intent(out) :: res(:)
    ! ---------------------------------------------------------------------- !
    integer :: iElem, depVar_pos, depVar_nComps, offset, iDof
    real(kind=rk), allocatable :: input_varRes(:)
    ! ---------------------------------------------------------------------- !

    ! get the position of dependent variable, assuming only one
    depVar_pos = fun%input_varPos(1)

    ! get the nComp of dependent variable
    depVar_nComps = varSys%method%val( depVar_pos )%nComponents

    ! allocate array to save results
    allocate(input_varRes( nElems * nDofs * depVar_nComps ))

    ! derive dependent variable
    call varSys%method%val(depVar_pos)%get_element(                 &
      &                                   varSys  = varSys,         &
      &                                   elemPos = elemPos,        &
      &                                   time    = time,           &
      &                                   tree    = tree,           &
      &                                   nElems  = nElems,         &
      &                                   nDofs   = nDofs,          &
      &                                   res     = input_varRes(:) )

    ! compute magnitude
    ! assuming fun%nComponents = 1
    input_varRes = input_varRes * input_varRes
    do iElem = 1, nElems
      do iDof = 1, nDofs
        offset = (( ielem-1)* depvar_ncomps* ndofs+( idof-1)* depvar_ncomps+0)
        res((( ielem-1)* 1* ndofs+( idof-1)* 1+1) )                    &
          & = sqrt(sum(input_varRes(offset+1 : offset+depVar_nComps)))
      end do
    end do

    deallocate(input_varRes)

  end subroutine tem_evalMag_forElement
  ! ************************************************************************** !


  ! ************************************************************************** !
  !> Same as evalMag_forElement except it evaluate magnitude from points
  !!
  !! The interface has to comply to the abstract interface
  !! tem_varSys_module#tem_varSys_proc_element.
  recursive subroutine tem_evalMag_forPoint( fun, varsys, point, time, tree, &
    &                                        nPnts, 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

    !> Three-dimensional coordinates at which the variable should be
    !! evaluated. Only useful for variables provided as space-time functions.
    real(kind=rk), intent(in) :: point(:,:)

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

    !> global treelm mesh info
    type(treelmesh_type), intent(in) :: tree

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

    !> 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 :: iPnt, depVar_pos, depVar_nComps, offset
    real(kind=rk), allocatable :: input_varRes(:)
    ! ---------------------------------------------------------------------- !

    ! get the position of dependent variable, assuming only one
    depVar_pos = fun%input_varPos(1)

    ! get the nComp of dependent variable
    depVar_nComps = varSys%method%val( depVar_pos )%nComponents

    ! allocate array to save results
    allocate(input_varRes( nPnts * depVar_nComps ))

    ! derive dependent variable
    call varSys%method%val(depVar_pos)%get_point(                  &
      &                                   varSys = varSys,         &
      &                                   point  = point,          &
      &                                   time   = time,           &
      &                                   tree   = tree,           &
      &                                   nPnts  = nPnts,          &
      &                                   res    = input_varRes(:) )

    ! compute magnitude
    ! assuming fun%nComponents = 1
    input_varRes = input_varRes * input_varRes
    do iPnt = 1, nPnts
      offset = (( ipnt-1)* depvar_ncomps+0)
      res(iPnt) = sqrt(sum(input_varRes(offset+1 : offset+depVar_nComps)))
    end do

    deallocate(input_varRes)

  end subroutine tem_evalMag_forPoint
  ! ************************************************************************** !


  ! ************************************************************************** !
  !> Same as evalMag_forPoint except it evaluate magnitude from points via
  !! indices which need to be setup before
  !!
  !! The interface has to comply to the abstract interface
  !! tem_varSys_module#tem_varSys_proc_element.
  recursive subroutine tem_evalMag_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: nVals
    integer, intent(in) :: idx(:)

    !> With idx as start index in contiguous memory,
    !! idxLength defines length of each contiguous memory
    !! Size: 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 :: iVal, depVar_pos, depVar_nComps, offset
    real(kind=rk), allocatable :: input_varRes(:)
    ! ---------------------------------------------------------------------- !
    type(tem_varSys_op_data_type), pointer :: opData
    ! ---------------------------------------------------------------------- !

    call C_F_POINTER( fun%method_Data, opData )

    call tem_varSys_check_inArgs( fun, varSys, time, iLevel, idx, idxLen, &
      &    nVals, label = 'evalMag_fromIndex'                             )

    ! get the position of dependent variable, assuming only one
    depVar_pos = fun%input_varPos(1)

    ! get the nComp of dependent variable
    depVar_nComps = varSys%method%val( depVar_pos )%nComponents

    ! allocate array to save results
    allocate(input_varRes( nVals * depVar_nComps ))

    ! derive dependent variable
    call varSys%method%val(depVar_pos)%get_valOfIndex(            &
      &                varSys  = varSys,                          &
      &                time    = time,                            &
      &                iLevel  = iLevel,                          &
      &                idx     = opData%input_pntIndex(1)         &
      &                          %indexLvl(iLevel)%val( idx(:) ), &
      &                nVals   = nVals,                           &
      &                res     = input_varRes(:)                  )

    ! compute magnitude
    ! assuming fun%nComponents = 1
    input_varRes = input_varRes * input_varRes
    do iVal = 1, nVals
      offset = (( ival-1)* depvar_ncomps+0)
      res(iVal) = sqrt(sum(input_varRes(offset+1 : offset+depVar_nComps)))
    end do

    deallocate(input_varRes)

  end subroutine tem_evalMag_fromIndex
  ! ************************************************************************** !


  ! ************************************************************************** !
  !> Routine to multiply variables if all variables have same number of
  !! components.
  !!
  !! \verbatim
  !! -- in lua file, one can define as following:
  !! variable = {{
  !!   name = 'coeff',
  !!   ncomponents = 1,
  !!   vartype = "st_fun",
  !!   st_fun =  0.25
  !!   },
  !!   {
  !!   name = 'newVel',
  !!   ncomponents = 1,
  !!   vartype = "operation",
  !!   operation = {kind='multiplication',
  !!                input_varname={coeff, vel_mag}}
  !!   }
  !!   ...
  !!  }
  !! \endverbatim
  !!
  !! The interface has to comply to the abstract interface
  !! tem_varSys_module#tem_varSys_proc_element.
  !!
  recursive subroutine evalMultiply_forElement( fun, varsys, elempos, time, &
    &                                           tree, nElems, nDofs, 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

    !> Position of the TreeID of the element to get the variable for in the
    !! global treeID list.
    integer, intent(in) :: elempos(:)

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

    !> global treelm mesh info
    type(treelmesh_type), intent(in) :: tree

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

    !> Number of degrees of freedom within an element.
    integer, intent(in) :: nDofs

    !> Resulting values for the requested variable.
    !!
    !! Linearized array dimension:
    !! (n requested entries) x (nComponents of this variable)
    !! x (nDegrees of freedom)
    !! Access: (iElem-1)*fun%nComponents*nDofs +
    !!         (iDof-1)*fun%nComponents + iComp
    real(kind=rk), intent(out) :: res(:)
    ! ---------------------------------------------------------------------- !
    integer :: iDep, posDepVar
    integer :: nTotal
    real(kind=rk), allocatable :: input_varRes(:)
    ! ---------------------------------------------------------------------- !

    if (nDofs > 1) then
      write(*,*) 'TODO: evalmultiply does not work for polynomial data yet'
      write(*,*) '      It makes use of a multiplication, which can not'
      write(*,*) '      directly be done that simple in modal space!'
      write(*,*) ''
      write(*,*) 'Need to replace this routine in Ateles!'
      write(*,*) 'Stopping now.'
      call tem_abort()
    end if

    nTotal = nElems*nDofs*fun%nComponents
    allocate( input_varRes(nTotal) )

    res = 1.0_rk

    do iDep = 1, fun%nInputs
      ! get position of dependent var in varSys
      posDepVar = fun%input_varPos(iDep)

      ! derive dependent variable
      call varSys%method%val(posDepVar)%get_element(               &
        &                                   varSys  = varSys,      &
        &                                   elemPos = elemPos,     &
        &                                   time    = time,        &
        &                                   tree    = tree,        &
        &                                   nElems  = nElems,      &
        &                                   nDofs   = nDofs,       &
        &                                   res     = input_varRes )

      res(:nTotal) = res(:nTotal) * input_varRes

    end do

    deallocate(input_varRes)

  end subroutine evalMultiply_forElement
  ! ************************************************************************** !


  ! ************************************************************************** !
  !> Same as evalMultiply_forElement except it evaluate it multiply values
  !! from points
  !!
  !! The interface has to comply to the abstract interface
  !! tem_varSys_module#tem_varSys_proc_point.
  recursive subroutine evalMultiply_forPoint( fun, varsys, point, time, tree, &
    &                                         nPnts, 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

    !> Three-dimensional coordinates at which the variable should be
    !! evaluated. Only useful for variables provided as space-time functions.
    real(kind=rk), intent(in) :: point(:,:)

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

    !> global treelm mesh info
    type(treelmesh_type), intent(in) :: tree

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

    !> 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 :: iDep, posDepVar
    integer :: nTotal
    real(kind=rk), allocatable :: input_varRes(:)
    ! ---------------------------------------------------------------------- !

    nTotal = nPnts*fun%nComponents
    allocate( input_varRes(nTotal) )

    res = 1.0_rk

    do iDep = 1, fun%nInputs
      ! get position of dependent var in varSys
      posDepVar = fun%input_varPos(iDep)

      ! derive dependent variable
      call varSys%method%val(posDepVar)%get_point(                &
        &                                   varSys = varSys,      &
        &                                   point  = point,       &
        &                                   time   = time,        &
        &                                   tree   = tree,        &
        &                                   nPnts  = nPnts,       &
        &                                   res    = input_varRes )

      res(:nTotal) = res(:nTotal) * input_varRes

    end do

    deallocate(input_varRes)

  end subroutine evalMultiply_forPoint
  ! ************************************************************************** !


  ! ************************************************************************** !
  !> Same as evalMultiply_forPoint except it multiply values from points via
  !! indices which are setup before
  !!
  !! The interface has to comply to the abstract interface
  !! tem_varSys_module#tem_varSys_proc_getvalofindex.
  recursive subroutine evalMultiply_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: nVals
    integer, intent(in) :: idx(:)

    !> With idx as start index in contiguous memory,
    !! idxLength defines length of each contiguous memory
    !! Size: 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 :: iDep, posDepVar
    real(kind=rk), allocatable :: input_varRes(:)
    ! ---------------------------------------------------------------------- !
    type(tem_varSys_op_data_type), pointer :: opData
    integer :: nTotal
    ! ---------------------------------------------------------------------- !

    call C_F_POINTER( fun%method_Data, opData )

    call tem_varSys_check_inArgs( fun, varSys, time, iLevel, idx, idxLen, &
      &    nVals, label = 'evalMultiply_fromIndex'                        )

    nTotal = nVals*fun%nComponents
    allocate( input_varRes(nTotal) )

    res = 1.0_rk

    do iDep = 1, fun%nInputs
      ! get position of dependent var in varSys
      posDepVar = fun%input_varPos(iDep)

      ! derive dependent variable
      call varSys%method%val(posDepVar)%get_valOfIndex( &
          & varSys  = varSys,                           &
          & time    = time,                             &
          & iLevel  = iLevel,                           &
          & idx     = opData%input_pntIndex(iDep)       &
          &           %indexLvl(iLevel)%val( idx(:) ),  &
          & nVals   = nVals,                            &
          & res     = input_varRes                      )

      res(:nTotal) = res(:nTotal) * input_varRes

    end do

    deallocate(input_varRes)

  end subroutine evalMultiply_fromIndex
  ! ************************************************************************** !


  ! ***************************************************************************!
  !> Routine to divide variables if all variables have same number of
  !! components.
  !!
  !! \verbatim
  !! -- in lua file, one can define as following:
  !! variable = {{
  !!   name = 'coeff',
  !!   ncomponents = 3,
  !!   vartype = "st_fun",
  !!   st_fun =  0.25
  !!   },
  !!   {
  !!   name = 'newVel',
  !!   ncomponents = 1,
  !!   vartype = "operation",
  !!   operation = {kind='division',
  !!                input_varname={velocity, coeff}} -- numerator, denominator
  !!   }
  !!   ...
  !!  }
  !! \endverbatim
  !!
  !! The interface has to comply to the abstract interface
  !! tem_varSys_module#tem_varSys_proc_element.
  !!
  recursive subroutine division_forElement( fun, varsys, elempos, time, tree, &
    &                                       nElems, nDofs, 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

    !> Position of the TreeID of the element to get the variable for in the
    !! global treeID list.
    integer, intent(in) :: elempos(:)

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

    !> global treelm mesh info
    type(treelmesh_type), intent(in) :: tree

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

    !> Number of degrees of freedom within an element.
    integer, intent(in) :: nDofs

    !> Resulting values for the requested variable.
    !!
    !! Linearized array dimension:
    !! (n requested entries) x (nComponents of this variable)
    !! x (nDegrees of freedom)
    !! Access: (iElem-1)*fun%nComponents*nDofs +
    !!         (iDof-1)*fun%nComponents + iComp
    real(kind=rk), intent(out) :: res(:)
    ! -------------------------------------------------------------------- !
    real(kind=rk) :: divisor(nElems*nDofs*fun%nComponents)
    real(kind=rk) :: dividend(nElems*nDofs*fun%nComponents)
    integer :: iComp, iDof, idx, iElem
    ! -------------------------------------------------------------------- !

    call varSys%method%val(fun%input_varPos(1))%get_element( &
      & varSys  = varSys,                                    &
      & elemPos = elemPos,                                   &
      & time    = time,                                      &
      & tree    = tree,                                      &
      & nElems  = nElems,                                    &
      & nDofs   = nDofs,                                     &
      & res     = dividend                                   )

    call varSys%method%val(fun%input_varPos(2))%get_element( &
      & varSys  = varSys,                                    &
      & elemPos = elemPos,                                   &
      & time    = time,                                      &
      & tree    = tree,                                      &
      & nElems  = nElems,                                    &
      & nDofs   = nDofs,                                     &
      & res     = divisor                                    )

    do iElem = 1, nElems
      do iDof = 1, nDofs
        do iComp = 1, fun%nComponents
          idx = (( ielem-1)* fun%ncomponents* ndofs+( idof-1)* fun%ncomponents+icomp)
          res(idx) = dividend(idx) / divisor(idx)
        end do
      end do
    end do

  end subroutine division_forElement
  ! ***************************************************************************!


  ! ************************************************************************** !
  !> Same as division_forElement except it evaluate it division values
  !! from points
  !!
  !! The interface has to comply to the abstract interface
  !! tem_varSys_module#tem_varSys_proc_point.
  recursive subroutine division_forPoint( fun, varsys, point, time, tree, &
    &                                     nPnts, 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

    !> Three-dimensional coordinates at which the variable should be
    !! evaluated. Only useful for variables provided as space-time functions.
    real(kind=rk), intent(in) :: point(:,:)

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

    !> global treelm mesh info
    type(treelmesh_type), intent(in) :: tree

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

    !> 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(:)
    ! -------------------------------------------------------------------- !
    real(kind=rk) :: divisor(nPnts*fun%nComponents)
    real(kind=rk) :: dividend(nPnts*fun%nComponents)
    integer :: iComp, iPnt, idx
    ! -------------------------------------------------------------------- !

    call varSys%method%val(fun%input_varPos(1))%get_point( &
      & varSys  = varSys,                                  &
      & point   = point,                                   &
      & time    = time,                                    &
      & tree    = tree,                                    &
      & nPnts   = nPnts,                                   &
      & res     = dividend                                 )

    call varSys%method%val(fun%input_varPos(2))%get_point( &
      & varSys  = varSys,                                  &
      & point   = point,                                   &
      & time    = time,                                    &
      & tree    = tree,                                    &
      & nPnts   = nPnts,                                   &
      & res     = divisor                                  )

    do iPnt = 1, nPnts
      do iComp = 1, fun%nComponents
        idx = (( ipnt-1)* fun%ncomponents+icomp)
        res(idx) = dividend(idx) / divisor(idx)
      end do
    end do

  end subroutine division_forPoint
  ! ***************************************************************************!


  ! ************************************************************************** !
  !> Same as division_forElement except it division values from points via
  !! indices which are setup before
  !!
  !! The interface has to comply to the abstract interface
  !! tem_varSys_module#tem_varSys_proc_getvalofindex.
  recursive subroutine division_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: nVals
    integer, intent(in) :: idx(:)

    !> With idx as start index in contiguous memory,
    !! idxLength defines length of each contiguous memory
    !! Size: 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(tem_varSys_op_data_type), pointer :: opData
    real(kind=rk) :: divisor(nVals*fun%nComponents)
    real(kind=rk) :: dividend(nVals*fun%nComponents)
    integer :: iComp, iVal
    ! -------------------------------------------------------------------- !
    call C_F_POINTER( fun%method_Data, opData )

    ! get the nodal values for the 2 inputs
    !>TODO make it working for idxLen and contiguous access of index array
    call tem_varSys_check_inArgs( fun, varSys, time, iLevel, idx, idxLen, &
      &    nVals, label = 'division_fromIndex'                            )

    call varSys%method%val(fun%input_varPos(1))%get_valOfIndex( &
      & varSys  = varSys,                                       &
      & time    = time,                                         &
      & iLevel  = iLevel,                                       &
      & idx     = opData%input_pntIndex(1)                      &
      &           %indexLvl(iLevel)%val( idx(:) ),              &
      & nVals   = nVals,                                        &
      & res     = dividend                                      )

    call varSys%method%val(fun%input_varPos(2))%get_valOfIndex( &
      & varSys  = varSys,                                       &
      & time    = time,                                         &
      & iLevel  = iLevel,                                       &
      & idx     = opData%input_pntIndex(2)                      &
      &           %indexLvl(iLevel)%val( idx(:) ),              &
      & nVals   = nVals,                                        &
      & res     = divisor                                       )

    do iVal = 1, nVals
      do iComp = 1, fun%nComponents
        res((iVal-1)*fun%nComponents + iComp) =         &
          & dividend((iVal-1)*fun%nComponents + iComp)  &
          & / divisor((iVal-1)*fun%nComponents + iComp)
      end do
    end do

  end subroutine division_fromIndex
  ! ***************************************************************************!


  ! ***************************************************************************!
  !> Routine to multiply scalat times vector
  !!
  !! \verbatim
  !! -- in lua file, one can define as following:
  !! variable = {{
  !!   name = 'coeff',
  !!   ncomponents = 1,
  !!   vartype = "st_fun",
  !!   st_fun =  0.25
  !!   },
  !!   {
  !!   name = 'newVel',
  !!   ncomponents = 3,
  !!   vartype = "operation",
  !!   operation = {kind='multiply_scalar_times_vector',
  !!                input_varname={coeff, velocity }} -- scalar, vector
  !!   }
  !!   ...
  !!  }
  !! \endverbatim
  !!
  !! The interface has to comply to the abstract interface
  !! tem_varSys_module#tem_varSys_proc_element.
  !!
  recursive subroutine multiplyScalTimesVec_forElement( fun, varsys, elempos, &
    &                                                   time, tree, nElems,   &
    &                                                   nDofs, 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

    !> Position of the TreeID of the element to get the variable for in the
    !! global treeID list.
    integer, intent(in) :: elempos(:)

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

    !> global treelm mesh info
    type(treelmesh_type), intent(in) :: tree

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

    !> Number of degrees of freedom within an element.
    integer, intent(in) :: nDofs

    !> Resulting values for the requested variable.
    !!
    !! Linearized array dimension:
    !! (n requested entries) x (nComponents of this variable)
    !! x (nDegrees of freedom)
    !! Access: (iElem-1)*fun%nComponents*nDofs +
    !!         (iDof-1)*fun%nComponents + iComp
    real(kind=rk), intent(out) :: res(:)
    ! -------------------------------------------------------------------- !
    real(kind=rk) :: scalar(nElems*nDofs)
    real(kind=rk) :: vector(nElems*nDofs*fun%nComponents)
    integer :: iComp, iElem, iDof
    ! -------------------------------------------------------------------- !

    call varSys%method%val(fun%input_varPos(1))%get_element( &
      & varSys  = varSys,                                    &
      & elemPos = elemPos,                                   &
      & time    = time,                                      &
      & tree    = tree,                                      &
      & nElems  = nElems,                                    &
      & nDofs   = nDofs,                                     &
      & res     = scalar                                     )

    call varSys%method%val(fun%input_varPos(2))%get_element( &
      & varSys  = varSys,                                    &
      & elemPos = elemPos,                                   &
      & time    = time,                                      &
      & tree    = tree,                                      &
      & nElems  = nElems,                                    &
      & nDofs   = nDofs,                                     &
      & res     = vector                                     )

    do iElem = 1, nElems
      do iDof = 1, nDofs
        do iComp = 1, fun%nComponents
          res((( ielem-1)* fun%ncomponents* ndofs+( idof-1)* fun%ncomponents+icomp)) =        &
            & scalar((iElem-1)*nDofs + iDof)                                  &
            & * vector((( ielem-1)* fun%ncomponents* ndofs+( idof-1)* fun%ncomponents+icomp))
        end do
      end do
    end do

  end subroutine multiplyScalTimesVec_forElement
  ! ***************************************************************************!


  ! ***************************************************************************!
  !> Same as multiplyScalTimesVec_forElement except it multiply values for
  !! given corrdinate points
  !!
  !! The interface has to comply to the abstract interface
  !! tem_varSys_module#tem_varSys_proc_getvalofindex.
  recursive subroutine tem_multiplyScalTimesVec_forPoint( fun, varsys, point, &
    &                                                     time, tree, nPnts,  &
    &                                                     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

    !> Three-dimensional coordinates at which the variable should be
    !! evaluated. Only useful for variables provided as space-time functions.
    real(kind=rk), intent(in) :: point(:,:)

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

    !> global treelm mesh info
    type(treelmesh_type), intent(in) :: tree

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

    !> 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(:)
    !--------------------------------------------------------------------------!
    real(kind=rk) :: scalar(nPnts)
    real(kind=rk) :: vector(nPnts*fun%nComponents)
    integer :: iComp, iPnt
    !--------------------------------------------------------------------------!

    call varSys%method%val(fun%input_varPos(1))%get_point( &
      & varSys  = varSys,                                  &
      & point   = point,                                   &
      & time    = time,                                    &
      & tree    = tree,                                    &
      & nPnts   = nPnts,                                   &
      & res     = scalar                                   )

    call varSys%method%val(fun%input_varPos(2))%get_point( &
      & varSys  = varSys,                                  &
      & point   = point,                                   &
      & time    = time,                                    &
      & tree    = tree,                                    &
      & nPnts   = nPnts,                                   &
      & res     = vector                                   )

    do iPnt = 1, nPnts
      do iComp = 1, fun%nComponents
        res((( ipnt-1)* fun%ncomponents+icomp)) =                   &
          & scalar(iPnt) * vector((( ipnt-1)* fun%ncomponents+icomp))
      end do
    end do

  end subroutine tem_multiplyScalTimesVec_forPoint
  ! ***************************************************************************!


  ! ***************************************************************************!
  !> Same as multiplyScalTimesVec_fromElement except it multiply values
  !! from points via indices which are setup before
  !!
  !! The interface has to comply to the abstract interface
  !! tem_varSys_module#tem_varSys_proc_getvalofindex.
  recursive subroutine tem_multiplyScalTimesVec_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: nVals
    integer, intent(in) :: idx(:)

    !> With idx as start index in contiguous memory,
    !! idxLength defines length of each contiguous memory
    !! Size: 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(tem_varSys_op_data_type), pointer :: opData
    real(kind=rk) :: scalar(nVals)
    real(kind=rk) :: vector(nVals*fun%nComponents)
    integer :: iComp, iVal
    !--------------------------------------------------------------------------!
    write(logUnit(10),*) 'Get the values of indices for derived variable',&
      &                  'by multiplyScalTimesVec ',                      &
      &                  trim(varSys%varname%val(fun%myPos))
    call C_F_POINTER( fun%method_Data, opData )

    ! get the nodal values for the 2 inputs
    !>TODO make it working for idxLen and contiguous access of index array
    call tem_varSys_check_inArgs( fun, varSys, time, iLevel, idx, idxLen, &
      &    nVals, label = 'tem_multiplyScalTimesVec_fromIndex'            )

    call varSys%method%val(fun%input_varPos(1))%get_valOfIndex( &
      & varSys  = varSys,                                       &
      & time    = time,                                         &
      & iLevel  = iLevel,                                       &
      & idx     = opData%input_pntIndex(1)                      &
      &           %indexLvl(iLevel)%val( idx(:) ),              &
      & nVals   = nVals,                                        &
      & res     = scalar                                        )

    call varSys%method%val(fun%input_varPos(2))%get_valOfIndex( &
      & varSys  = varSys,                                       &
      & time    = time,                                         &
      & iLevel  = iLevel,                                       &
      & idx     = opData%input_pntIndex(2)                      &
      &           %indexLvl(iLevel)%val( idx(:) ),              &
      & nVals   = nVals,                                        &
      & res     = vector                                        )

    do iVal = 1, nVals
      do iComp = 1, fun%nComponents
        res((iVal-1)*fun%nComponents + iComp) =                   &
          & scalar(iVal) * vector((iVal-1)*fun%nComponents + iComp)
      end do
    end do
  end subroutine tem_multiplyScalTimesVec_fromIndex
  ! ***************************************************************************!


  ! ***************************************************************************!
  !> Routine to divide vector by scalar
  !!
  !! \verbatim
  !! -- in lua file, one can define as following:
  !! variable = {{
  !!   name = 'coeff',
  !!   ncomponents = 1,
  !!   vartype = "st_fun",
  !!   st_fun =  0.25
  !!   },
  !!   {
  !!   name = 'newVel',
  !!   ncomponents = 3,
  !!   vartype = "operation",
  !!   operation = {kind='division',
  !!                input_varname={velocity, coeff}} -- numerator, denominator
  !!   }
  !!   ...
  !!  }
  !! \endverbatim
  !!
  !! The interface has to comply to the abstract interface
  !! tem_varSys_module#tem_varSys_proc_element.
  !!
  recursive subroutine divideVecByScal_forElement( fun, varsys, elempos, time, &
    &                                              tree, nElems, nDofs, 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

    !> Position of the TreeID of the element to get the variable for in the
    !! global treeID list.
    integer, intent(in) :: elempos(:)

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

    !> global treelm mesh info
    type(treelmesh_type), intent(in) :: tree

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

    !> Number of degrees of freedom within an element.
    integer, intent(in) :: nDofs

    !> Resulting values for the requested variable.
    !!
    !! Linearized array dimension:
    !! (n requested entries) x (nComponents of this variable)
    !! x (nDegrees of freedom)
    !! Access: (iElem-1)*fun%nComponents*nDofs +
    !!         (iDof-1)*fun%nComponents + iComp
    real(kind=rk), intent(out) :: res(:)
    !--------------------------------------------------------------------------!
    real(kind=rk) :: divisor(nElems*nDofs)
    real(kind=rk) :: dividend(nElems*nDofs*fun%nComponents)
    integer :: iComp, iElem, iDof
    !--------------------------------------------------------------------------!

    call varSys%method%val(fun%input_varPos(1))%get_element( &
      & varSys  = varSys,                                    &
      & elemPos = elemPos,                                   &
      & time    = time,                                      &
      & tree    = tree,                                      &
      & nElems  = nElems,                                    &
      & nDofs   = nDofs,                                     &
      & res     = dividend                                   )

    call varSys%method%val(fun%input_varPos(2))%get_element( &
      & varSys  = varSys,                                    &
      & elemPos = elemPos,                                   &
      & time    = time,                                      &
      & tree    = tree,                                      &
      & nElems  = nElems,                                    &
      & nDofs   = nDofs,                                     &
      & res     = divisor                                    )

    do iElem = 1, nElems
      do iDof = 1, nDofs
        do iComp = 1, fun%nComponents
          res((( ielem-1)* fun%ncomponents* ndofs+( idof-1)* fun%ncomponents+icomp)) =        &
            & dividend((( ielem-1)* fun%ncomponents* ndofs+( idof-1)* fun%ncomponents+icomp)) &
            & / divisor((iElem-1)*nDofs + iDof)
        end do
      end do
    end do

  end subroutine divideVecByScal_forElement
  ! ***************************************************************************!


  ! ***************************************************************************!
  !> Same as divideVecByScal_forElement except it divides values from points
  !!
  !! The interface has to comply to the abstract interface
  !! tem_varSys_module#tem_varSys_proc_getvalofindex.
  recursive subroutine tem_divideVecByScal_forPoint( fun, varsys, point, time, &
    &                                                tree, nPnts, 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

    !> Three-dimensional coordinates at which the variable should be
    !! evaluated. Only useful for variables provided as space-time functions.
    real(kind=rk), intent(in) :: point(:,:)

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

    !> global treelm mesh info
    type(treelmesh_type), intent(in) :: tree

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

    !> 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(:)
    !--------------------------------------------------------------------------!
    real(kind=rk) :: divisor(nPnts)
    real(kind=rk) :: dividend(nPnts*fun%nComponents)
    integer :: iComp, iPnt
    !--------------------------------------------------------------------------!

    call varSys%method%val(fun%input_varPos(1))%get_point( &
      & varSys  = varSys,                                  &
      & point   = point,                                   &
      & time    = time,                                    &
      & tree    = tree,                                    &
      & nPnts   = nPnts,                                   &
      & res     = dividend                                 )

    call varSys%method%val(fun%input_varPos(2))%get_point( &
      & varSys  = varSys,                                  &
      & point   = point,                                   &
      & time    = time,                                    &
      & tree    = tree,                                    &
      & nPnts   = nPnts,                                   &
      & res     = divisor                                  )

    do iPnt = 1, nPnts
      do iComp = 1, fun%nComponents
        res((( ipnt-1)* fun%ncomponents+icomp)) =                      &
          & dividend((( ipnt-1)* fun%ncomponents+icomp)) / divisor(iPnt)
      end do
    end do

  end subroutine tem_divideVecByScal_forPoint
  ! ***************************************************************************!


  ! ***************************************************************************!
  !> Same as divideVecByScal_fromElement except it divide values from points
  !! via indices which are setup before
  !!
  !! The interface has to comply to the abstract interface
  !! tem_varSys_module#tem_varSys_proc_getvalofindex.
  recursive subroutine tem_divideVecByScal_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: nVals
    integer, intent(in) :: idx(:)

    !> With idx as start index in contiguous memory,
    !! idxLength defines length of each contiguous memory
    !! Size: 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(tem_varSys_op_data_type), pointer :: opData
    real(kind=rk) :: divisor(nVals)
    real(kind=rk) :: dividend(nVals*fun%nComponents)
    integer :: iComp, iVal
    !--------------------------------------------------------------------------!
    write(logUnit(4),*) 'Get the values of indices for derived variable',&
      &                  'by divideVecByScal ',                          &
      &                  trim(varSys%varname%val(fun%myPos))
    call C_F_POINTER( fun%method_Data, opData )

    ! get the nodal values for the 2 inputs
    !>TODO make it working for idxLen and contiguous access of index array
    call tem_varSys_check_inArgs( fun, varSys, time, iLevel, idx, idxLen, &
      &    nVals, label = 'tem_divideVecByScal_fromIndex'            )

    call varSys%method%val(fun%input_varPos(1))%get_valOfIndex( &
      & varSys  = varSys,                                       &
      & time    = time,                                         &
      & iLevel  = iLevel,                                       &
      & idx     = opData%input_pntIndex(1)                      &
      &           %indexLvl(iLevel)%val( idx(:) ),              &
      & nVals   = nVals,                                        &
      & res     = dividend                                      )

    call varSys%method%val(fun%input_varPos(2))%get_valOfIndex( &
      & varSys  = varSys,                                       &
      & time    = time,                                         &
      & iLevel  = iLevel,                                       &
      & idx     = opData%input_pntIndex(2)                      &
      &           %indexLvl(iLevel)%val( idx(:) ),              &
      & nVals   = nVals,                                        &
      & res     = divisor                                       )

    do iVal = 1, nVals
      do iComp = 1, fun%nComponents
        res((iVal-1)*fun%nComponents + iComp) =                      &
          & dividend((iVal-1)*fun%nComponents + iComp) / divisor(iVal)
      end do
    end do
  end subroutine tem_divideVecByScal_fromIndex
  ! ***************************************************************************!


  ! ************************************************************************** !
  !> Extract component index of any vectorial variable.
  !! In lua file, first define new variable with varType operation kind as
  !! "extract" and provide name of the variable from which to extract
  !! component index via input_varname (it must be single variable) and
  !! index to extract via input_varIndex.
  !! If input_varname variable is not part of predefined solver variables then
  !! add also that variable via variable table.
  !!
  !! \verbatim
  !! -- in lua file, one can define as following:
  !! variable = {{
  !!   name = 'vel_y',
  !!   ncomponents = 1,
  !!   vartype = "operation",
  !!   operation = {kind='extract',
  !!                input_varname={'velocity'},
  !!                input_varindex = {2}
  !!               }
  !!   },
  !!  }
  !! \endverbatim
  !!
  !! The interface has to comply to the abstract interface
  !! tem_varSys_module#tem_varSys_proc_element.
  !!
  recursive subroutine extract_forElement( fun, varsys, elempos, time, tree, &
    &                                      nElems, nDofs, 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

    !> Position of the TreeID of the element to get the variable for in the
    !! global treeID list.
    integer, intent(in) :: elempos(:)

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

    !> global treelm mesh info
    type(treelmesh_type), intent(in) :: tree

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

    !> Number of degrees of freedom within an element.
    integer, intent(in) :: nDofs

    !> Resulting values for the requested variable.
    !!
    !! Linearized array dimension:
    !! (n requested entries) x (nComponents of this variable)
    !! x (nDegrees of freedom)
    !! Access: (iElem-1)*fun%nComponents*nDofs +
    !!         (iDof-1)*fun%nComponents + iComp
    real(kind=rk), intent(out) :: res(:)
    ! ---------------------------------------------------------------------- !
    integer :: iElem, depVar_pos, depVar_nComps, iDof, comp_index, iComp
    real(kind=rk), allocatable :: input_varRes(:)
    ! ---------------------------------------------------------------------- !

    ! get the position of dependent variable, assuming only one
    depVar_pos = fun%input_varPos(1)

    ! get the nComp of dependent variable
    depVar_nComps = varSys%method%val( depVar_pos )%nComponents

    ! allocate array to save results
    allocate(input_varRes( nElems * nDofs * depVar_nComps ))

    ! derive dependent variable
    call varSys%method%val(depVar_pos)%get_element(              &
      &                                   varSys  = varSys,      &
      &                                   elemPos = elemPos,     &
      &                                   time    = time,        &
      &                                   tree    = tree,        &
      &                                   nElems  = nElems,      &
      &                                   nDofs   = nDofs,       &
      &                                   res     = input_varRes )

    ! Extract component index from input_varRes
    do iElem = 1, nElems
      do iDof = 1, nDofs
        do iComp = 1, fun%nComponents
          comp_index = fun%input_varIndex(iComp)
          res((( ielem-1)* fun%ncomponents* ndofs+( idof-1)* fun%ncomponents+icomp) )  &
            & = input_varRes(                                          &
            & (( ielem-1)* depvar_ncomps* ndofs+( idof-1)* depvar_ncomps+comp_index))
        end do
      end do
    end do

    deallocate(input_varRes)

  end subroutine extract_forElement
  ! ************************************************************************** !


  ! ************************************************************************** !
  !> Same as extract_forElement except it extract from points
  !!
  !! The interface has to comply to the abstract interface
  !! tem_varSys_module#tem_varSys_proc_element.
  !!
  recursive subroutine extract_forPoint( fun, varsys, point, time, tree, &
    &                                    nPnts, 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

    !> Three-dimensional coordinates at which the variable should be
    !! evaluated. Only useful for variables provided as space-time functions.
    real(kind=rk), intent(in) :: point(:,:)

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

    !> global treelm mesh info
    type(treelmesh_type), intent(in) :: tree

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

    !> 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 :: iPnt, depVar_pos, depVar_nComps, iComp, comp_index
    real(kind=rk), allocatable :: input_varRes(:)
    ! ---------------------------------------------------------------------- !

    ! get the position of dependent variable, assuming only one
    depVar_pos = fun%input_varPos(1)

    ! get the nComp of dependent variable
    depVar_nComps = varSys%method%val( depVar_pos )%nComponents

    ! allocate array to save results
    allocate(input_varRes( nPnts * depVar_nComps ))

    ! derive dependent variable
    call varSys%method%val(depVar_pos)%get_point(                &
      &                                   varSys  = varSys,      &
      &                                   point   = point,       &
      &                                   time    = time,        &
      &                                   tree    = tree,        &
      &                                   nPnts   = nPnts,       &
      &                                   res     = input_varRes )

    ! Extract component index from input_varRes
    do iPnt = 1, nPnts
      do iComp = 1, fun%nComponents
        comp_index = fun%input_varIndex(iComp)
        res((( ipnt-1)* fun%ncomponents+icomp) )                  &
          & = input_varRes( (( ipnt-1)* depvar_ncomps+comp_index) )
      end do
    end do

    deallocate(input_varRes)

  end subroutine extract_forPoint
  ! ************************************************************************** !


  ! ************************************************************************** !
  !> Same as extract_from Point except it extract from points via indices
  !! which are setup before
  !!
  !! The interface has to comply to the abstract interface
  !! tem_varSys_module#tem_varSys_proc_getvalofIndex.
  !!
  recursive subroutine extract_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: nVals
    integer, intent(in) :: idx(:)

    !> With idx as start index in contiguous memory,
    !! idxLength defines length of each contiguous memory
    !! Size: 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 :: iVal, depVar_pos, depVar_nComps, iComp, comp_index
    real(kind=rk), allocatable :: input_varRes(:)
    ! ---------------------------------------------------------------------- !
    type(tem_varSys_op_data_type), pointer :: opData
    ! ---------------------------------------------------------------------- !

    call C_F_POINTER( fun%method_Data, opData )

    call tem_varSys_check_inArgs( fun, varSys, time, iLevel, idx, idxLen, &
      &    nVals, label = 'extract_fromIndex'                             )

    ! get the position of dependent variable, assuming only one
    depVar_pos = fun%input_varPos(1)

    ! get the nComp of dependent variable
    depVar_nComps = varSys%method%val( depVar_pos )%nComponents

    ! allocate array to save results
    allocate(input_varRes( nVals * depVar_nComps ))
    ! derive dependent variable
    call varSys%method%val(depVar_pos)%get_valOfIndex( &
        & varSys  = varSys,                            &
        & time    = time,                              &
        & iLevel  = iLevel,                            &
        & idx     = opData%input_pntIndex(1)           &
        &           %indexLvl(iLevel)%val( idx(:) ),   &
        & nVals   = nVals,                             &
        & res     = input_varRes                       )

    ! Extract component index from input_varRes
    do iVal = 1, nVals
      do iComp = 1, fun%nComponents
        comp_index = fun%input_varIndex(iComp)
        res((( ival-1)* fun%ncomponents+icomp) )                  &
          & = input_varRes( (( ival-1)* depvar_ncomps+comp_index) )
      end do
    end do

    deallocate(input_varRes)

  end subroutine extract_fromIndex
  ! ************************************************************************** !


  ! ************************************************************************** !
  !> Combine multiple variables into single variable with nComponent of
  !! output variable as sum of all input variables nComponents.
  !! In lua file, first define new variable with varType operation kind as
  !! "combine" and provide name of the variable from which to extract
  !! component index via input_varname (it must be single variable) and
  !! index to combine via input_varIndex.
  !! If input_varname variable is not part of predefined solver variables then
  !! add also that variable via variable table.
  !!
  !! \verbatim
  !! -- in lua file, one can define as following:
  !! variable = {{
  !!   name = 'dens_and_vel',
  !!   ncomponents = 4,
  !!   vartype = "operation",
  !!   operation = {kind='combine',
  !!                input_varname={'density','velocity'}
  !!               }
  !!   },
  !!  }
  !! \endverbatim
  !!
  !! The interface has to comply to the abstract interface
  !! tem_varSys_module#tem_varSys_proc_element.
  !!
  recursive subroutine combine_forElement( fun, varsys, elempos, time, tree, &
    &                                      nElems, nDofs, 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

    !> Position of the TreeID of the element to get the variable for in the
    !! global treeID list.
    integer, intent(in) :: elempos(:)

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

    !> global treelm mesh info
    type(treelmesh_type), intent(in) :: tree

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

    !> Number of degrees of freedom within an element.
    integer, intent(in) :: nDofs

    !> Resulting values for the requested variable.
    !!
    !! Linearized array dimension:
    !! (n requested entries) x (nComponents of this variable)
    !! x (nDegrees of freedom)
    !! Access: (iElem-1)*fun%nComponents*nDofs +
    !!         (iDof-1)*fun%nComponents + iComp
    real(kind=rk), intent(out) :: res(:)
    ! ---------------------------------------------------------------------- !

    ! This routine combines multiple variable together
    call tem_get_element_chunk( varSys  = varSys,           &
      &                         varPos  = fun%input_varPos, &
      &                         elemPos = elemPos,          &
      &                         time    = time,             &
      &                         tree    = tree,             &
      &                         nElems  = nElems,           &
      &                         nDofs   = nDofs,            &
      &                         res     = res               )

  end subroutine combine_forElement
  ! ************************************************************************** !


  ! ************************************************************************** !
  !> same as combine_fromelement except it extract from points
  !!
  !! the interface has to comply to the abstract interface
  !! tem_varsys_module#tem_varsys_proc_element.
  !!
  recursive subroutine combine_forPoint( fun, varsys, point, time, tree, &
    &                                    nPnts, 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

    !> Three-dimensional coordinates at which the variable should be
    !! evaluated. Only useful for variables provided as space-time functions.
    real(kind=rk), intent(in) :: point(:,:)

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

    !> global treelm mesh info
    type(treelmesh_type), intent(in) :: tree

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

    !> 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(:)
    ! ---------------------------------------------------------------------- !
    ! This routine combines multiple variable together
    call tem_get_point_chunk( varSys = varSys,           &
      &                       varPos = fun%input_varPos, &
      &                       point  = point,            &
      &                       time   = time,             &
      &                       tree   = tree,             &
      &                       nPnts  = nPnts,            &
      &                       res    = res               )

  end subroutine combine_forPoint
  ! ************************************************************************** !


  ! ************************************************************************** !
  !> Same as combine_from Point except it combine from points via indices
  !! which are setup before
  !!
  !! The interface has to comply to the abstract interface
  !! tem_varSys_module#tem_varSys_proc_getvalofIndex.
  !!
  recursive subroutine combine_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: nVals
    integer, intent(in) :: idx(:)

    !> With idx as start index in contiguous memory,
    !! idxLength defines length of each contiguous memory
    !! Size: 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 :: iVal, iDep
    integer :: maxComponents, compOff, dep_nComps, depVar_pos
    real(kind=rk), allocatable :: input_varRes(:)
    integer :: e_start, t_start, res_size
    type(tem_varSys_op_data_type), pointer :: opData
    ! ---------------------------------------------------------------------- !
    call C_F_POINTER( fun%method_Data, opData )

    call tem_varSys_check_inArgs( fun, varSys, time, iLevel, idx, idxLen, &
      &    nVals, label = 'combine_fromIndex'                             )

    ! Need to obtain the data variable for variable, and store it in an
    ! intermediate array, because all components should be put together in the
    ! res array.
    ! The temporary array therefore needs to be sufficiently large to store the
    ! maximal number of components.
    maxComponents = maxval(varSys%method%val(fun%input_varPos(:))%nComponents)

    ! Using a temporary array to store the variables and transfer them to res
    ! in the correct ordering afterwards.
    allocate(input_varRes(nVals*maxComponents))

    compOff = 0
    do iDep = 1, fun%nInputs

      ! get the position of dependent variable
      depVar_pos = fun%input_varPos(iDep)

      ! get the number of components for variable iVar
      dep_nComps = varSys%method%val(depVar_pos)%nComponents

      ! get the size of the needed part of the res array
      res_size = nVals * dep_nComps

      ! derive dependent variable
      call varSys%method%val(depVar_pos)%get_valOfIndex( &
          & varSys  = varSys,                            &
          & time    = time,                              &
          & iLevel  = iLevel,                            &
          & idx     = opData%input_pntIndex(iDep)        &
          &           %indexLvl(iLevel)%val( idx(:) ),   &
          & nVals   = nVals,                             &
          & res     = input_varRes(:res_size)            )

      ! copy the information to the right positions in the result array
      ! res contains results for all variables,
      ! input_varRes is only for one variable
      do iVal = 1, nVals
        e_start = (iVal-1)*fun%nComponents + compOff
        t_start = (iVal-1)*dep_nComps
        res( (e_start+1) : (e_start+dep_nComps) )               &
          &  = input_varRes( t_start + 1 : t_start + dep_nComps )
      end do
      ! Increase the component offset for the next variables.
      compOff = compOff + dep_nComps
    end do !iDep

    deallocate(input_varRes)

  end subroutine combine_fromIndex
  ! ************************************************************************** !

  ! ************************************************************************** !
  !> Routine to return time reduction variables
  !!
  !! \verbatim
  !! -- in lua file, one can define as following:
  !! variable = {
  !!   {
  !!     name = 'press_timeavg',
  !!     ncomponents = 1,
  !!     vartype = "operation",
  !!     operation = {
  !!       kind='reduction_transient',
  !!       input_varname={'pressure'},
  !!       reduction_transient = {kind = 'average', nrecord = 1000}
  !!     }
  !!   }
  !!   ...
  !!  }
  !! \endverbatim
  !!
  !! The interface has to comply to the abstract interface
  !! tem_varSys_module#tem_varSys_proc_element.
  !!
  recursive subroutine reductionTransient_forElement( fun, varsys, elempos, &
    &                                                 time, tree, nElems,   &
    &                                                 nDofs, 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

    !> Position of the TreeID of the element to get the variable for in the
    !! global treeID list.
    integer, intent(in) :: elempos(:)

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

    !> global treelm mesh info
    type(treelmesh_type), intent(in) :: tree

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

    !> Number of degrees of freedom within an element.
    integer, intent(in) :: nDofs

    !> Resulting values for the requested variable.
    !!
    !! Linearized array dimension:
    !! (n requested entries) x (nComponents of this variable)
    !! x (nDegrees of freedom)
    !! Access: (iElem-1)*fun%nComponents*nDofs +
    !!         (iDof-1)*fun%nComponents + iComp
    real(kind=rk), intent(out) :: res(:)
    ! ------------------------------------------------------------------------ !
    type(tem_varSys_op_data_type), pointer :: opData
    ! ------------------------------------------------------------------------ !
    !call C_F_POINTER( fun%method_Data, opData )
    ! Avoid warning about unused varsys dummy argument
    call C_F_POINTER( varsys%method%val(fun%mypos)%method_Data, opData )

    if (time%sim < 0.0_rk) then
      write(logunit(10),*) 'Avoid unused argument warning for time'
      write(logunit(10),*) 'tree%nElems:', tree%nElems
    end if

    res = 0.0_rk
    call tem_reduction_transient_getElement( me      = opData%redTrans, &
      &                                      elemPos = elemPos,         &
      &                                      nElems  = nElems,          &
      &                                      nDofs   = nDofs,           &
      &                                      res     = res              )

  end subroutine reductionTransient_forElement
  ! ************************************************************************** !


  ! ************************************************************************** !
  !> Same as reductionTransient_forElement except it evaluate it multiply values
  !! from points
  !!
  !! The interface has to comply to the abstract interface
  !! tem_varSys_module#tem_varSys_proc_point.
  recursive subroutine reductionTransient_forPoint( fun, varsys, point, time, &
    &                                               tree, nPnts, 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

    !> Three-dimensional coordinates at which the variable should be
    !! evaluated. Only useful for variables provided as space-time functions.
    real(kind=rk), intent(in) :: point(:,:)

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

    !> global treelm mesh info
    type(treelmesh_type), intent(in) :: tree

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

    !> 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(:)
    ! ---------------------------------------------------------------------- !
    res = 0.0_rk
    write(logunit(1),*) 'Variable: ', trim(varsys%varname%val(fun%mypos))
    write(logunit(1),*) 'nPnts=', nPnts
    write(logunit(1),*) 'tree%nElems=', tree%nElems
    write(logunit(1),*) 'time%sim=', time%sim
    write(logunit(1),*) 'size(point)=', size(point)
    flush(logunit(1))
    call tem_abort('Reduction_transient for Point is not implemented yet')

  end subroutine reductionTransient_forPoint
  ! ************************************************************************** !


  ! ************************************************************************** !
  !> Same as reductionTransient_forPoint except it multiply values from points
  !! via indices which are setup before
  !!
  !! The interface has to comply to the abstract interface
  !! tem_varSys_module#tem_varSys_proc_getvalofindex.
  recursive subroutine reductionTransient_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: nVals
    integer, intent(in) :: idx(:)

    !> With idx as start index in contiguous memory,
    !! idxLength defines length of each contiguous memory
    !! Size: 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(:)
    ! ---------------------------------------------------------------------- !
    res = 0.0_rk
    call tem_varSys_check_inArgs( fun, varSys, time, iLevel, idx, idxLen, &
      &    nVals, label = 'reductionTransient_fromIndex'                  )
    call tem_abort('Reduction_transient from Index is not implemented yet')
  end subroutine reductionTransient_fromIndex
  ! ************************************************************************** !


  ! ************************************************************************** !
  !> This subroutine call set_params of input_variable
  !!
  !! the interface has to comply to the abstract interface
  !! tem_varsys_module#tem_varsys_proc_setParams.
  !!
  recursive subroutine tem_opVar_setParams(fun, varSys, instring)
    ! ---------------------------------------------------------------------- !
    !> 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

    !> Input string with parameter to set in method_data
    character(len=*), intent(in) :: instring
    ! ---------------------------------------------------------------------- !
    integer :: iDep, posDepVar
    ! ---------------------------------------------------------------------- !
    do iDep = 1, fun%nInputs
      ! get position of dependent var in varSys
      posDepVar = fun%input_varPos(iDep)

      ! set params in dependent variable
      call varSys%method%val(posDepVar)%set_params( varSys   = varSys,  &
        &                                           inString = inString )
    end do

  end subroutine tem_opVar_setParams
  ! ************************************************************************** !


  ! ************************************************************************** !
  !> This subroutine call get_params of input_variable
  !!
  !! the interface has to comply to the abstract interface
  !! tem_varsys_module#tem_varsys_proc_setParams.
  !!
  recursive subroutine tem_opVar_getParams( fun, varSys, instring, outstring )
    ! ---------------------------------------------------------------------- !
    !> 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

    !> Input string with parameter to set in method_data
    character(len=*), intent(in) :: instring

    !> Output string with requested parameter value from method_data
    character(len=*), intent(out) :: outstring
    ! ---------------------------------------------------------------------- !
    integer :: iDep, posDepVar
    ! ---------------------------------------------------------------------- !
    do iDep = 1, fun%nInputs
      ! get position of dependent var in varSys
      posDepVar = fun%input_varPos(iDep)

      ! get params from dependent variable
      call varSys%method%val(posDepVar)%get_params( varSys    = varSys,   &
        &                                           inString  = inString, &
        &                                           outString = outString )

      ! if outString is filled by any dependent variable then exit loop
      if (len(trim(outString)) > 0) then
        outString = trim(outString)//'_oper'
        exit
      end if
    end do

  end subroutine tem_opVar_getParams
  ! ************************************************************************** !


  ! ************************************************************************** !
  !> This subroutine call setup indices of input_variable
  !!
  !! the interface has to comply to the abstract interface
  !! tem_varsys_module#tem_varsys_proc_setupIndices.
  !!
  recursive subroutine tem_opVar_setupIndices( fun, varSys, point, offset_bit, &
    &                                          iLevel, tree, nPnts, idx        )
    ! ---------------------------------------------------------------------- !
    !> 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

    !> List of space coordinate points to store as growing array in
    !! method_data
    real(kind=rk), intent(in) :: point(:,:)

    !> Offset bit encoded as character for every point.
    !!
    !! Offset integer coord(3) is converted into a character with
    !! offset_bit = achar( (coord(1)+1) + (coord(2)+1)*4 + (coord(3)+1)*16 )
    !! Backward transformation form character to 3 integer:
    !! coord(1) = mod(ichar(offset_bit),4) - 1
    !! coord(2) = mod(ichar(offset_bit),16)/4 - 1
    !! coord(3) = ichar(offset_bit)/16 - 1
    !!
    !! If not present default is to center i.e offset_bit = achar(1+4+16)
    character, optional, intent(in) :: offset_bit(:)

    !> Level to which input points belong to
    integer, intent(in) :: iLevel

    !> global treelm mesh info
    type(treelmesh_type), intent(in) :: tree

    !> Number of points to add in method_data of this variable
    integer, intent(in) :: nPnts

    !> Index of points in the growing array and variable val array.
    !! Size: nPoints
    !!
    !! This must be stored in boundary or source depends on who
    !! calls this routine.
    !! This index is required to return a value using getValOfIndex.
    integer, intent(out) :: idx(:)
    ! ---------------------------------------------------------------------- !
    type(tem_varSys_op_data_type), pointer :: opData
    integer :: iPnt, iDep
    type(grw_intArray_type), allocatable :: inputIndex_loc(:)
    integer, allocatable :: idxPerPnt(:)
    ! ---------------------------------------------------------------------- !

    call C_F_POINTER( fun%method_Data, opData )

    ! allcoate the index array for all inpits
    if (.not. allocated(opData%input_pntIndex)) then
      allocate( opData%input_pntIndex(fun%nInputs) )
    end if
    ! allocate temporary inputIndex with size of nInputs and initialize
    ! growing array with length nPnts
    allocate(inputIndex_loc(fun%nInputs))

    ! Now fill in the index arrays for the inputs
    call tem_opVar_fill_inputIndex( fun        = fun,                &
      &                             varSys     = varSys,             &
      &                             point      = point,              &
      &                             offset_bit = offset_bit,         &
      &                             iLevel     = iLevel,             &
      &                             tree       = tree,               &
      &                             nPnts      = nPnts,              &
      &                             inputIndex = inputIndex_loc      )

    ! KM: Workaround for intel compiler in SuperMUC, the pointer
    ! gets corrupted after recursive call to depend variable
    call C_F_POINTER( fun%method_Data, opData )
    ! initialize index with zero to identify points which does not
    ! belong to subTree
    allocate(idxPerPnt(fun%nInputs))
    idx = 0
    do iPnt = 1, nPnts
      do iDep = 1, fun%nInputs
        idxPerPnt(iDep) = inputIndex_loc(iDep)%val(iPnt)
      end do
      ! set index only when any of dependent variable has valid index
      if (any(idxPerPnt > 0)) then
        do iDep = 1, fun%nInputs
          call append(me  = opData%input_pntIndex(iDep)%indexLvl(iLevel), &
            &         val = inputIndex_loc(iDep)%val(iPnt)              )
        end do
        ! set index to last position in input_pntIndex of dep var 1 of
        ! indexLvl of iLevel
        idx(iPnt) = opData%input_pntIndex(1)%indexLvl(iLevel)%nVals
      end if
    end do

    do iDep = 1, fun%nInputs
      call truncate (opData%input_pntIndex(iDep)%indexLvl(iLevel) )
    end do

  end subroutine tem_opVar_setupIndices
  ! ************************************************************************** !


  ! ************************************************************************** !
  !> subroutine to fill index for the setuo Index routine called for operation
  !! variables, it is also used by the solver
  recursive subroutine tem_opVar_fill_inputIndex( fun, varSys, point,       &
    &                                             offset_bit, iLevel, tree, &
    &                                             nPnts, inputIndex         )
    !---------------------------`----------------------------------------------!
    !> 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

    !> List of space coordinate points to store as growing array in
    !! method_data
    real(kind=rk), intent(in) :: point(:,:)

    !> Offset bit encoded as character for every point.
    !!
    !! Offset integer coord(3) is converted into a character with
    !! offset_bit = achar( (coord(1)+1) + (coord(2)+1)*4 + (coord(3)+1)*16 )
    !! Backward transformation form character to 3 integer:
    !! coord(1) = mod(ichar(offset_bit),4) - 1
    !! coord(2) = mod(ichar(offset_bit),16)/4 - 1
    !! coord(3) = ichar(offset_bit)/16 - 1
    !!
    !! If not present default is to center i.e offset_bit = achar(1+4+16)
    character, optional, intent(in) :: offset_bit(:)

    !> Level to which input points belong to
    integer, intent(in) :: iLevel

    !> global treelm mesh info
    type(treelmesh_type), intent(in) :: tree

    !> Number of points to add in method_data of this variable
    integer, intent(in) :: nPnts

    !> input index for dependent variables
    !! size: fun%nInputs
    type(grw_intArray_type), intent(out) :: inputIndex(:)
    !--------------------------------------------------------------------------!
    integer, allocatable :: idx_loc(:)
    integer :: iDep, posDepVar
    !--------------------------------------------------------------------------!
    ! allocate local index array, needed for setup_indices call
    allocate( idx_loc(nPnts) )

    do iDep = 1, fun%nInputs
      idx_loc = 0

      ! get position of dependent var in varSys
      posDepVar = fun%input_varPos(iDep)

      ! setup indices in dependent variable.
      ! idx output from any variables will be same so it just overwrites and
      ! returns the idx of last dependent variable
      call varSys%method%val(posDepVar)%setup_indices( &
        & varSys     = varSys,                         &
        & point      = point,                          &
        & offset_bit = offset_bit,                     &
        & iLevel     = iLevel,                         &
        & tree       = tree,                           &
        & nPnts      = nPnts,                          &
        & idx        = idx_loc                         )

      call append( me  = inputIndex(iDep), &
        &          val = idx_loc           )

      call truncate (inputIndex(iDep))
      write(logUnit(9),*) 'nIndex on level for input variable ',  &
        &                 trim(varSys%varname%val(posDepVar)),    &
        &                 ' on Lvl ', iLevel,                     &
        &                 ' are = ', inputIndex(iDep)%nVals
    end do
    deallocate(idx_loc)

  end subroutine tem_opVar_fill_inputIndex

end module tem_operation_var_module
! **************************************************************************** !
