tem_spacetime_var_module.f90 Source File


This file depends on

sourcefile~~tem_spacetime_var_module.f90~~EfferentGraph sourcefile~tem_spacetime_var_module.f90 tem_spacetime_var_module.f90 sourcefile~env_module.f90 env_module.f90 sourcefile~tem_spacetime_var_module.f90->sourcefile~env_module.f90 sourcefile~tem_spacetime_fun_module.f90 tem_spacetime_fun_module.f90 sourcefile~tem_spacetime_var_module.f90->sourcefile~tem_spacetime_fun_module.f90 sourcefile~tem_varsys_module.f90 tem_varSys_module.f90 sourcefile~tem_spacetime_var_module.f90->sourcefile~tem_varsys_module.f90 sourcefile~treelmesh_module.f90 treelmesh_module.f90 sourcefile~tem_spacetime_var_module.f90->sourcefile~treelmesh_module.f90 sourcefile~tem_spatial_module.f90 tem_spatial_module.f90 sourcefile~tem_spacetime_var_module.f90->sourcefile~tem_spatial_module.f90 sourcefile~tem_logging_module.f90 tem_logging_module.f90 sourcefile~tem_spacetime_var_module.f90->sourcefile~tem_logging_module.f90 sourcefile~tem_dyn_array.f90 tem_dyn_array.f90 sourcefile~tem_spacetime_var_module.f90->sourcefile~tem_dyn_array.f90 sourcefile~tem_pointdata_module.f90 tem_pointData_module.f90 sourcefile~tem_spacetime_var_module.f90->sourcefile~tem_pointdata_module.f90 sourcefile~tem_tools_module.f90 tem_tools_module.f90 sourcefile~tem_spacetime_var_module.f90->sourcefile~tem_tools_module.f90 sourcefile~tem_param_module.f90 tem_param_module.f90 sourcefile~tem_spacetime_var_module.f90->sourcefile~tem_param_module.f90 sourcefile~tem_geometry_module.f90 tem_geometry_module.f90 sourcefile~tem_spacetime_var_module.f90->sourcefile~tem_geometry_module.f90 sourcefile~tem_topology_module.f90 tem_topology_module.f90 sourcefile~tem_spacetime_var_module.f90->sourcefile~tem_topology_module.f90 sourcefile~tem_variable_module.f90 tem_variable_module.f90 sourcefile~tem_spacetime_var_module.f90->sourcefile~tem_variable_module.f90 sourcefile~tem_time_module.f90 tem_time_module.f90 sourcefile~tem_spacetime_var_module.f90->sourcefile~tem_time_module.f90 sourcefile~tem_aux_module.f90 tem_aux_module.f90 sourcefile~tem_spacetime_var_module.f90->sourcefile~tem_aux_module.f90 sourcefile~tem_grow_array.f90 tem_grow_array.f90 sourcefile~tem_spacetime_var_module.f90->sourcefile~tem_grow_array.f90

Files dependent on this one

sourcefile~~tem_spacetime_var_module.f90~~AfferentGraph sourcefile~tem_spacetime_var_module.f90 tem_spacetime_var_module.f90 sourcefile~tem_derived_module.f90 tem_derived_module.f90 sourcefile~tem_derived_module.f90->sourcefile~tem_spacetime_var_module.f90 sourcefile~tem_operation_var_module.f90 tem_operation_var_module.f90 sourcefile~tem_derived_module.f90->sourcefile~tem_operation_var_module.f90 sourcefile~tem_varmap_module.f90 tem_varMap_module.f90 sourcefile~tem_varmap_module.f90->sourcefile~tem_spacetime_var_module.f90 sourcefile~tem_operation_var_module.f90->sourcefile~tem_varmap_module.f90 sourcefile~tem_convergence_module.f90 tem_convergence_module.f90 sourcefile~tem_convergence_module.f90->sourcefile~tem_varmap_module.f90 sourcefile~tem_restart_module.f90 tem_restart_module.f90 sourcefile~tem_restart_module.f90->sourcefile~tem_varmap_module.f90 sourcefile~tem_depend_module.f90 tem_depend_module.f90 sourcefile~tem_depend_module.f90->sourcefile~tem_varmap_module.f90 sourcefile~tem_bc_module.f90 tem_bc_module.f90 sourcefile~tem_bc_module.f90->sourcefile~tem_varmap_module.f90 sourcefile~tem_tracking_module.f90 tem_tracking_module.f90 sourcefile~tem_tracking_module.f90->sourcefile~tem_varmap_module.f90 sourcefile~hvs_output_module.f90 hvs_output_module.f90 sourcefile~tem_tracking_module.f90->sourcefile~hvs_output_module.f90 sourcefile~tem_simcontrol_module.f90 tem_simControl_module.f90 sourcefile~tem_tracking_module.f90->sourcefile~tem_simcontrol_module.f90 sourcefile~hvs_output_module.f90->sourcefile~tem_varmap_module.f90 sourcefile~hvs_output_module.f90->sourcefile~tem_restart_module.f90 sourcefile~tem_abortcriteria_module.f90 tem_abortCriteria_module.f90 sourcefile~tem_abortcriteria_module.f90->sourcefile~tem_convergence_module.f90 sourcefile~tem_ini_condition_module.f90 tem_ini_condition_module.f90 sourcefile~tem_ini_condition_module.f90->sourcefile~tem_depend_module.f90 sourcefile~tem_general_module.f90 tem_general_module.f90 sourcefile~tem_general_module.f90->sourcefile~tem_restart_module.f90 sourcefile~tem_general_module.f90->sourcefile~tem_abortcriteria_module.f90 sourcefile~tem_general_module.f90->sourcefile~tem_simcontrol_module.f90 sourcefile~tem_simcontrol_module.f90->sourcefile~tem_convergence_module.f90 sourcefile~tem_simcontrol_module.f90->sourcefile~tem_abortcriteria_module.f90

Contents


Source Code

! Copyright (c) 2016-2017, 2019 Kannan Masilamani <kannan.masilamani@uni-siegen.de>
! Copyright (c) 2016, 2018 Harald Klimach <harald.klimach@uni-siegen.de>
! Copyright (c) 2016-2018 Peter Vitt <peter.vitt2@uni-siegen.de>
! Copyright (c) 2016 Verena Krupp <verena.krupp@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) 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.
! ***************************************************************************** !
!> This module provides variable definitions for space time functions.
!!
module tem_spacetime_var_module
  use, intrinsic :: iso_c_binding,  only: c_ptr,       &
    &                                     c_loc,       &
    &                                     c_f_pointer

  ! include treelm modules
  use env_module,               only: rk, long_k, labelLen
  use tem_param_module,         only: qOffset_inChar, q000
  use tem_varSys_module,        only: tem_varSys_type,               &
    &                                 tem_varSys_op_type,            &
    &                                 tem_varSys_check_inArgs,       &
    &                                 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_tools_module,         only: tem_PositionInSorted
  use tem_geometry_module,      only: tem_CoordOfReal, &
    &                                 tem_posOfId
  use tem_topology_module,      only: tem_IdOfCoord
  use tem_variable_module,      only: tem_variable_type
  use tem_time_module,          only: tem_time_type
  use tem_spacetime_fun_module, only: tem_spacetime_for,          &
    &                                 tem_st_fun_linkedList_type, &
    &                                 append,                     &
    &                                 tem_st_fun_listElem_type,   &
    &                                 tem_spacetime_fun_type
  use tem_spatial_module,       only: tem_spatial_storeVal
  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: truncate, PositionOfVal
  use tem_grow_array_module,    only: init, append, truncate
  use tem_pointData_module,     only: init, append, truncate

  ! aotus modules
  use aotus_module,     only: aot_get_val, aot_top_get_val,                    &
    &                         aoterr_Fatal, aoterr_WrongType,                  &
    &                         aoterr_NonExistent, flu_State, open_config_chunk

  implicit none

  private

  public :: tem_varSys_append_stFun
  public :: evaluate_add_spacetime_scalarByTreeID
  public :: evaluate_add_spacetime_vectorByTreeID
  public :: evaluate_add_spacetime_scalarByCoordinate
  public :: evaluate_add_spacetime_vectorByCoordinate
  public :: evaluate_first_spacetime_scalarByCoordinate
  public :: evaluate_first_spacetime_scalarByTreeID
  public :: evaluate_first_spacetime_vectorByCoordinate
  public :: evaluate_first_spacetime_vectorByTreeID
  public :: get_valOfIndex_add_scalar_spacetime
  public :: get_valOfIndex_add_vector_spacetime
  public :: get_valOfIndex_first_scalar_spacetime
  public :: get_valOfIndex_first_vector_spacetime
  public :: set_params_spacetime
  public :: get_params_spacetime
  public :: setup_indices_spacetime
  public :: evaluate_FOAG_spacetime_scalarByTreeID
  public :: evaluate_FOAG_spacetime_vectorByTreeID
  public :: evaluate_FOAG_spacetime_scalarByCoordinate
  public :: evaluate_FOAG_spacetime_vectorByCoordinate
  public :: get_valOfIndex_FOAG_scalar_spacetime
  public :: get_valOfIndex_FOAG_vector_spacetime

  interface tem_varSys_append_stfun
    module procedure tem_varSys_append_stFunVar
    module procedure tem_varSys_append_stFun_raw
  end interface tem_varSys_append_stfun


contains


  ! *************************************************************************** !
  !> Returns the get_point and get_element pointer according to the requested
  !! evaluation type.
  subroutine tem_varSys_assignEvalType(evaltype, nComp, get_point, &
    &                                  get_element, get_valOfIndex )
    ! -------------------------------------------------------------------------
    character(len=*), intent(in) :: evaltype
    integer, intent(in) :: nComp
    !> The function pointer to the get_point subroutine for the given
    !! operation.
    procedure(tem_varSys_proc_point), pointer, intent(out) :: get_point
    !> The function pointer to the get_element subroutine for the given
    !! operation.
    procedure(tem_varSys_proc_element), pointer, intent(out) :: get_element
    !> The function pointer to the get_valOfIndex subroutine for the given
    !! operation.
    procedure(tem_varSys_proc_getValOfIndex), pointer, intent(out) :: &
      &                                       get_valOfIndex
    ! -------------------------------------------------------------------------

    select case(evalType)

    case ('add')
      ! if nComp = 1, use scalar function else use vector function
      if ( nComp == 1 ) then
        get_point => evaluate_add_spacetime_scalarByCoordinate
        get_element => evaluate_add_spacetime_scalarByTreeID
        get_valOfIndex => get_valOfIndex_add_scalar_spacetime
      else
        get_point => evaluate_add_spacetime_vectorByCoordinate
        get_element => evaluate_add_spacetime_vectorByTreeID
        get_valOfIndex => get_valOfIndex_add_vector_spacetime
      end if

    case ('first')
      ! if nComp = 1, use scalar function else use vector function
      if ( nComp == 1 ) then
        get_point => evaluate_first_spacetime_scalarByCoordinate
        get_element => evaluate_first_spacetime_scalarByTreeID
        get_valOfIndex => get_valOfIndex_first_scalar_spacetime
      else
        get_point => evaluate_first_spacetime_vectorByCoordinate
        get_element => evaluate_first_spacetime_vectorByTreeID
        get_valOfIndex => get_valOfIndex_first_vector_spacetime
      end if
    case ('firstonly_asglobal')
      ! Anonymous spacetime function variable created for boundary
      ! or source variable assumes global shape and single st-fun
      !
      ! if nComp = 1, use scalar function else use vector function
      if ( nComp == 1 ) then
        get_point => evaluate_FOAG_spacetime_scalarByCoordinate
        get_element => evaluate_FOAG_spacetime_scalarByTreeID
        get_valOfIndex => get_valOfIndex_FOAG_scalar_spacetime
      else
        get_point => evaluate_FOAG_spacetime_vectorByCoordinate
        get_element => evaluate_FOAG_spacetime_vectorByTreeID
        get_valOfIndex => get_valOfIndex_FOAG_vector_spacetime
      end if

    end select

  end subroutine tem_varSys_assignEvalType
  ! *************************************************************************** !


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

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

    !> contains spacetime functions of all variables
    type(tem_st_fun_linkedList_type), intent(inout) :: st_funList

    !> 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 :: nComp
    logical :: wasAdded
    type(c_ptr) :: method_data
    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(tem_st_fun_listElem_type), pointer  :: newElem
    ! -------------------------------------------------------------------------
    nullify(get_point, get_element, set_params, get_params, setup_indices, &
      &     get_valOfIndex)
    nComp = stFunVar%nComponents

    ! append space time function to linked list of spacetime functions
    call append( st_funList, stFunVar%st_fun, newElem )

    ! c pointer to list of spacetime functions of current variable
    !method_data = c_loc(st_funList%current)
    method_data = c_loc(newElem)

    ! assign function pointer depends on evaluation type
    call tem_varSys_assignEvalType( evaltype       = stfunVar%evaltype,    &
      &                             nComp          = stfunVar%nComponents, &
      &                             get_point      = get_point,            &
      &                             get_element    = get_element,          &
      &                             get_valOfIndex = get_valOfIndex        )

    set_params => set_params_spacetime
    get_params => get_params_spacetime
    setup_indices => setup_indices_spacetime

    if (.not. associated(get_point)) then
      call tem_abort( 'Error: No evaluation is defined for variable ' &
        & // trim(stfunvar%label) )
    end if

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

  end subroutine tem_varSys_append_stFunVar
  ! *************************************************************************** !


  ! *************************************************************************** !
  !> subroutine to add the variables from the input lua script to the varsys
  subroutine tem_varSys_append_stFun_raw( varSys, stFun, varname, nComp, &
    &                                     evaltype, st_funList,          &
    &                                     solverData_evalElem            )
    ! -------------------------------------------------------------------------
    !> global variable system to which stFunVar to be appended
    type(tem_varSys_type), intent(inout) :: varSys

    !> variables defined in the lua file
    type(tem_spacetime_fun_type), pointer, intent(in) :: stFun(:)

    character(len=*), intent(in) :: varname

    integer, intent(in), optional :: nComp

    character(len=*), intent(in), optional :: evaltype

    !> contains spacetime functions of all variables
    type(tem_st_fun_linkedList_type), intent(inout), optional :: st_funList

    !> 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
    ! -------------------------------------------------------------------------
    type(tem_st_fun_listElem_type), pointer :: stfun_listelem
    integer :: addedPos
    logical :: wasAdded
    type(c_ptr) :: method_data
    integer :: ncomp_loc
    character(len=labelLen) :: evaltype_loc
    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(tem_st_fun_listElem_type),pointer :: newElem
    ! -------------------------------------------------------------------------
    nullify(get_point, get_element, set_params, get_params, setup_indices, &
      &     get_valOfIndex)

    nComp_loc = 1
    if (present(nComp)) nComp_loc = nComp

    evaltype_loc = 'add'
    if (present(evaltype)) evaltype_loc = evaltype

    if (present(st_funList)) then
      ! append space time function to linked list of spacetime functions
      call append( st_funList, stfun, newElem )

      ! c pointer to list of spacetime functions of current variable
      method_data = c_loc(newElem)
    else
      allocate(stfun_listelem)
      stfun_listelem%val => stfun
      stfun_listelem%nvals = size(stfun)
      method_data = c_loc(stfun_Listelem)
    end if


    ! assign function pointer depends on evaluation type
    call tem_varSys_assignEvalType( evaltype       = evaltype_loc,  &
      &                             nComp          = nComp_loc,     &
      &                             get_point      = get_point,     &
      &                             get_element    = get_element,   &
      &                             get_valOfIndex = get_valOfIndex )

    set_params => set_params_spacetime
    get_params => get_params_spacetime
    setup_indices => setup_indices_spacetime

    if (.not. associated(get_point)) then
      write(logUnit(1),*) 'Error: No evaluation is defined for variable '//&
        &                 trim(varname)
      call tem_abort()
    end if

    ! append variable to varSys
    call tem_varSys_append_derVar( me             = varSys,         &
      &                            varName        = varName,        &
      &                            operType       = 'st_fun',       &
      &                            nComponents    = nComp_loc,      &
      &                            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 function is provided, override
        ! the get_element pointer and use the provided setter
        ! solverData_evalElem instead to define the get_element routine.
        call solverData_evalElem%stFun_setter(varSys%method%val(addedPos))
      end if
      write(logUnit(9),*) 'Successfully appended variable "' &
        & // trim(varname) // '" to the variable system'

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

  end subroutine tem_varSys_append_stFun_raw
  ! *************************************************************************** !


  ! *************************************************************************** !
  !> This subroutine directly evaluate spacetime function for given coordinate
  !! points on 1st st_fun assuming global shape.
  recursive subroutine evaluate_FOAG_spacetime_scalarByCoordinate( 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(:)
    !--------------------------------------------------------------------------!
    ! --------------------------------------------------------------------------
    type(tem_st_fun_listElem_type), pointer :: fPtr
    ! --------------------------------------------------------------------------
    !call C_F_POINTER( fun%method_Data, fPtr )
    ! Avoid compiler warning about unused varsys
    call C_F_POINTER( varsys%method%val(fun%mypos)%method_Data, fPtr )
    res = tem_spacetime_for( me    = fPtr%val(1), &
      &                      coord = point,       &
      &                      time  = time,        &
      &                      n     = nPnts        )

    if (tree%nElems < 0) write(logunit(10),*) 'Avoid unused dummy argument'

  end subroutine evaluate_FOAG_spacetime_scalarByCoordinate
  ! *************************************************************************** !


  ! *************************************************************************** !
  !> This subroutine add multiple spacetime function for given coordinate
  !! points. Spacetime function are evaluated only on the coordinate which
  !! belong to st_fun shape.
  recursive subroutine evaluate_add_spacetime_scalarByCoordinate( 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 :: iStFun, iPoint, elemPos, pos
    type(tem_st_fun_listElem_type), pointer :: fPtr
    integer :: coord(4)
    integer(kind=long_k):: treeID
    real(kind=rk) :: st_res(1)
    ! --------------------------------------------------------------------------
    !call C_F_POINTER( fun%method_Data, fPtr )
    ! Avoid compiler warning about unused varsys
    call C_F_POINTER( varsys%method%val(fun%mypos)%method_Data, fPtr )
    res = 0.0_rk

    ! The space time function has a custom shape, thus we have to check
    ! whether the coordinates are covered by the space time function.
    ! As we are dealing with vectorial access here, we have to loop over all
    ! requested points.
    stLoop: do iStFun = 1, fPtr%nVals
      if(fPtr%val(iStFun)%subTree%useGlobalMesh) then
        res = res + tem_spacetime_for( me      = fPtr%val(iStFun), &
          &                            coord   = point,            &
          &                            time    = time,             &
          &                            n       = nPnts             )
      else
        do iPoint = 1, nPnts
          ! The space time function has a custom shape, thus we have to check
          ! whether the coordinates are covered by the space time function.
          coord =  tem_CoordOfReal(tree, point(iPoint,:), &
            &                      tree%global%maxLevel)
          treeId = tem_IdOfCoord(coord)
          elemPos = tem_PosofId(treeId, tree%treeID)
          pos = tem_PositionInSorted(                      &
            & me    = fPtr%val(iStFun)%subTree%map2global, &
            & val   = elemPos                              )
          ! When the elemenet is covered by the space time function's shape,
          ! we evaluate it at the given coordinates.
          if (pos > 0) then
            st_res = tem_spacetime_for(                            &
              & me      = fPtr%val(iStFun),                        &
              & coord   = reshape( source = (/ point(iPoint,1),    &
              &                                point(iPoint,2),    &
              &                                point(iPoint,3) /), &
              &                    shape  = (/ 1, 3 /) ),          &
              & time    = time,                                    &
              & n       = 1                                        )

            res(iPoint) = res(iPoint) + st_res(1)
          end if ! point part of subtree
        end do !iPoint
      end if ! global mesh
    end do stLoop

  end subroutine evaluate_add_spacetime_scalarByCoordinate
  ! *************************************************************************** !


  ! *************************************************************************** !
  !> This subroutine returns first spacetime function on overlapping
  !! regions fo given coordinate points.
  !! Spacetime function are evaluated only on the coordinate which
  !! belong to st_fun shape.
  recursive subroutine evaluate_first_spacetime_scalarByCoordinate( 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 :: iStFun, iPnt, elemPos, pos
    logical :: found
    type(tem_st_fun_listElem_type), pointer :: fPtr
    integer :: coord(4)
    integer(kind=long_k):: treeID
    real(kind=rk) :: st_res(1)
    ! --------------------------------------------------------------------------
    !call C_F_POINTER( fun%method_Data, fPtr )
    ! avoid compiler warning about unused varsys
    call C_F_POINTER( varsys%method%val(fun%mypos)%method_Data, fPtr )
    res = 0.0_rk

    do iPnt = 1, nPnts
      found = .false.
      ! Loop over all the space-time functions, if value is set by one then go
      ! next point
      stLoop: do iStFun = 1, fPtr%nVals
        if(fPtr%val(iStFun)%subTree%useGlobalMesh) then
          found = .true.
        else
          ! The space time function has a custom shape, thus we have to check
          ! whether the coordinates are covered by the space time function.
          coord =  tem_CoordOfReal(tree, point(iPnt,:), tree%global%maxLevel )
          treeId = tem_IdOfCoord(coord)
          elemPos = tem_PosofId(treeId, tree%treeID)
          pos = tem_PositionInSorted(                      &
            & me    = fPtr%val(iStFun)%subTree%map2global, &
            & val   = elemPos                              )
          ! When the elemenet is covered by the space time function's shape, we
          ! evaluate it at the given coordinates.
          if (pos > 0) then
            found = .true.
          end if ! point part of subtree
        end if ! global mesh

        ! When the first space time function was able to fulfill the request, we
        ! stop looping over them, because we only want one value.
        if (found) then
          st_res = tem_spacetime_for(                          &
            & me      = fPtr%val(iStFun),                      &
            & coord   = reshape( source = (/ point(iPnt,1),    &
            &                                point(iPnt,2),    &
            &                                point(iPnt,3) /), &
            &                    shape  = (/ 1, 3 /) ),        &
            & time    = time,                                  &
            & n       = 1                                      )

          res(iPnt) = st_res(1)
          exit stLoop
        endif
      end do stLoop
    end do !iPnt

  end subroutine evaluate_first_spacetime_scalarByCoordinate
  ! *************************************************************************** !


  ! *************************************************************************** !
  !> This subroutine directly evaluate spacetime function for given coordinate
  !! points on 1st st_fun assuming global shape for vectorial variable
  recursive subroutine evaluate_FOAG_spacetime_vectorByCoordinate( 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(:)
    !--------------------------------------------------------------------------!
    ! --------------------------------------------------------------------------
    type(tem_st_fun_listElem_type), pointer :: fPtr
    integer :: iPoint
    real(kind=rk) :: st_res(nPnts, fun%nComponents)
    ! --------------------------------------------------------------------------
    !call C_F_POINTER( fun%method_Data, fPtr )
    ! avoid compiler warning about unused varsys
    call C_F_POINTER( varsys%method%val(fun%mypos)%method_Data, fPtr )

    st_res = tem_spacetime_for( me      = fPtr%val(1),    &
      &                         coord   = point,          &
      &                         time    = time,           &
      &                         n       = nPnts,          &
      &                         nComp   = fun%nComponents )

    do iPoint = 1, nPnts
      res( (iPoint-1)*fun%nComponents + 1 : iPoint*fun%nComponents ) &
        & = st_res(iPoint, :)
    end do
    if (tree%nElems < 0) write(logunit(10),*) 'Avoid unused dummy argument'

  end subroutine evaluate_FOAG_spacetime_vectorByCoordinate
  ! *************************************************************************** !


  ! *************************************************************************** !
  !> This subroutine add multiple spacetime function for given coordinate
  !! points. Spacetime function are evaluated only on the coordinate which
  !! belong to st_fun shape for vectorial variable.
  recursive subroutine evaluate_add_spacetime_vectorByCoordinate( 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 :: iStFun, iPoint, elemPos, pos
    type(tem_st_fun_listElem_type), pointer :: fPtr
    integer :: coord(4)
    integer(kind=long_k):: treeID
    real(kind=rk) :: st_res_pnt(1, fun%nComponents)
    real(kind=rk) :: st_res(nPnts, fun%nComponents)
    ! --------------------------------------------------------------------------
    !call C_F_POINTER( fun%method_Data, fPtr )
    ! avoid compiler warning about unused varsys
    call C_F_POINTER( varsys%method%val(fun%mypos)%method_Data, fPtr )

    st_res = 0.0_rk

    ! The space time function has a custom shape, thus we have to check
    ! whether the coordinates are covered by the space time function.
    ! As we are dealing with vectorial access here, we have to loop over all
    ! requested points.
    stLoop: do iStFun = 1, fPtr%nVals
      if (fPtr%val(iStFun)%subTree%useGlobalMesh) then
        st_res = st_res + tem_spacetime_for( me      = fPtr%val(iStFun), &
          &                                  coord   = point,            &
          &                                  time    = time,             &
          &                                  n       = nPnts,            &
          &                                  nComp   = fun%nComponents   )
      else
        ! The space time function has a custom shape, thus we have to check
        ! whether the coordinates are covered by the space time function.
        do iPoint = 1, nPnts
          coord =  tem_CoordOfReal(tree, point(iPoint,:), tree%global%maxLevel)
          treeId = tem_IdOfCoord(coord)
          elemPos = tem_PosofId(treeId, tree%treeID)
          pos = tem_PositionInSorted(                      &
            & me    = fPtr%val(iStFun)%subTree%map2global, &
            & val   = elemPos                              )
          ! When the element is covered by the space time function's shape, we
          ! evaluate it at the given coordinates.
          if (pos > 0) then
            st_res_pnt = tem_spacetime_for(                        &
              & me      = fPtr%val(iStFun),                        &
              & coord   = reshape( source = (/ point(iPoint,1),    &
              &                                point(iPoint,2),    &
              &                                point(iPoint,3) /), &
              &                    shape  = (/ 1, 3 /) ),          &
              & time    = time,                                    &
              & n       = 1,                                       &
              & nComp   = fun%nComponents                          )

            st_res(iPoint,:) = st_res(iPoint,:) +  st_res_pnt(1,:)
          end if ! point part of subtree
        end do !iPoint
      end if ! global mesh
    end do stLoop

    do iPoint = 1, nPnts
      res( (iPoint-1)*fun%nComponents + 1 : iPoint*fun%nComponents ) &
        & = st_res(iPoint, :)
    end do

  end subroutine evaluate_add_spacetime_vectorByCoordinate
  ! *************************************************************************** !


  ! *************************************************************************** !
  !> This subroutine returns first spacetime function on overlapping
  !! regions fo given coordinate points for vectorial variable.
  !! Spacetime function are evaluated only on the coordinate which
  !! belong to st_fun shape.
  recursive subroutine evaluate_first_spacetime_vectorByCoordinate( 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 :: iStFun, iPoint, elemPos, pos
    logical :: found
    type(tem_st_fun_listElem_type), pointer :: fPtr
    integer :: coord(4)
    integer(kind=long_k):: treeID
    real(kind=rk) :: st_res(1, fun%nComponents)
    ! --------------------------------------------------------------------------

    !call C_F_POINTER( fun%method_Data, fPtr )
    ! avoid compiler warning about unused varsys
    call C_F_POINTER( varsys%method%val(fun%mypos)%method_Data, fPtr )

    res = 0.0_rk

    ! We have to loop over all requested points. For each point we have to
    ! figure out which is the first space time function that can fulfill the
    ! request. Thus a vectorized access is not possible here.
    do iPoint = 1, nPnts
      found = .false.
      ! The space time function has a custom shape, thus we have to check
      ! whether the coordinates are covered by the space time function.
      ! As we are dealing with vectorial access here, we have to loop over all
      ! requested points.
      stLoop: do iStFun = 1, fPtr%nVals
        if(fPtr%val(iStFun)%subTree%useGlobalMesh) then
          found = .true.
        else
          ! The space time function has a custom shape, thus we have to check
          ! whether the coordinates are covered by the space time function.
          coord =  tem_CoordOfReal(tree, point(iPoint,:), tree%global%maxLevel)
          treeId = tem_IdOfCoord(coord)
          elemPos = tem_PosofId(treeId, tree%treeID)
          pos = tem_PositionInSorted(                      &
            & me    = fPtr%val(iStFun)%subTree%map2global, &
            & val   = elemPos                              )
          ! When the elemenet is covered by the space time function's shape, we
          ! evaluate it at the given coordinates.
          if (pos > 0) then
            found = .true.
          end if ! point part of subtree
        end if ! global mesh

        ! When the first space time function was able to fulfill the request, we
        ! stop looping over them, because we only want one value.
        if (found) then
          st_res = tem_spacetime_for(                            &
            & me      = fPtr%val(iStFun),                        &
            & coord   = reshape( source = (/ point(iPoint,1),    &
            &                                point(iPoint,2),    &
            &                                point(iPoint,3) /), &
            &                    shape  = (/ 1, 3 /) ),          &
            & time    = time,                                    &
            & n       = 1,                                       &
            & nComp   = fun%nComponents                          )

          res( (iPoint-1)*fun%nComponents + 1 : iPoint*fun%nComponents ) &
            & = st_res(1, :)
          exit stLoop
        end if ! found
      end do stLoop
    end do !iVal

  end subroutine evaluate_first_spacetime_vectorByCoordinate
  ! *************************************************************************** !


  ! *************************************************************************** !
  !> This subroutine directly evaluate spacetime function for given element
  !! position in global TreeID list on 1st st_fun assuming global shape.
  recursive subroutine evaluate_FOAG_spacetime_scalarByTreeID( 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_st_fun_listElem_type), pointer :: fPtr
    ! -------------------------------------------------------------------------!
    !call C_F_POINTER( fun%method_Data, fPtr )
    ! avoid compiler warning about unused varsys
    call C_F_POINTER( varsys%method%val(fun%mypos)%method_Data, fPtr )

    ! assuming nDofs = 1
    if (nDofs > 1) then
      write(logUnit(1),*) 'Spacetime function does not support nDofs>1'
      call tem_abort()
    end if

    res(:) = tem_spacetime_for( me      = fPtr%val(1),             &
         &                      treeIDs = tree%treeID(elemPos(:)), &
         &                      tree    = tree,                    &
         &                      time    = time,                    &
         &                      n       = nElems                   )

  end subroutine evaluate_FOAG_spacetime_scalarByTreeID
  ! *************************************************************************** !


  ! *************************************************************************** !
  !> Call the spacetime function for scalar variable,
  !! which are stored in method_data which intern
  !! points to tem_st_fun_listelem_type.
  !! this spacetime function can be used as a analytical solution for reference.
  !! note: ndofs is not used by this routine. so, to evaluate spacetime function
  !! for ndofs in an element, use tem_varsys_proc_point interface.
  !!
  !! \verbatim
  !! -- in lua file, one can define it as following:
  !! variable = {{
  !!   name = 'reference',
  !!   ncomponents = 1,
  !!   vartype = st_fun,
  !!   st_fun =  luafun
  !!   },
  !!   ...
  !!  }
  !! \endverbatim
  !!
  !! the interface has to comply to the abstract interface
  !! tem_varsys_module#tem_varsys_proc_element.
  !!
  recursive subroutine evaluate_add_spacetime_scalarByTreeID( 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 :: iStFun, iElem, pos, oldPos
    type(tem_st_fun_listElem_type), pointer :: fPtr
    ! element wise output
    real(kind=rk) :: st_res(1)
    integer(kind=long_k) :: treeID(1)
    ! -------------------------------------------------------------------------!
    !call C_F_POINTER( fun%method_Data, fPtr )

    ! avoid compiler warning about unused varsys
    call C_F_POINTER( varsys%method%val(fun%mypos)%method_Data, fPtr )

    res = 0.0_rk
    ! assuming nDofs = 1
    if (nDofs > 1) then
      write(logUnit(1),*) 'Spacetime function does not support nDofs>1'
      call tem_abort()
    end if

    do iStFun = 1, fPtr%nVals
      if (fPtr%val(iStFun)%subTree%useGlobalMesh) then
        res(:) = res(:)                                                &
         &     + tem_spacetime_for( me      = fPtr%val(iStFun),        &
         &                          treeIDs = tree%treeID(elemPos(:)), &
         &                          tree    = tree,                    &
         &                          time    = time,                    &
         &                          n       = nElems                   )
      else
        oldPos = 1
        do iElem = 1, nElems
          ! position of element matches with any of map2global then
          ! this element is part of this subTree
          pos = tem_PositionInSorted(  &
            &                 me    = fPtr%val(iStFun)%subTree%map2global, &
            &                 val   = elemPos(iElem),                      &
            &                 lower = oldPos                               )
          if (pos > 0) then
            oldPos = pos
            treeID = tree%treeID(elemPos(iElem))
            st_res = tem_spacetime_for( me      = fPtr%val(iStFun),     &
              &                         treeIDs = treeID,               &
              &                         tree    = tree,                 &
              &                         time    = time,                 &
              &                         n       = 1                     )
            res(iElem) = res(iElem) + st_res(1)
          end if ! elemPos(iElem) part of subtree
        end do ! iElem
      end if ! global mesh
    end do ! nr. st_funs

  end subroutine evaluate_add_spacetime_scalarByTreeID
  ! *************************************************************************** !

  ! *************************************************************************** !
  !> Call the spacetime function for scalar variable,
  !! which are stored in method_data which intern
  !! points to tem_st_fun_listelem_type.
  !! this spacetime function can be used as a analytical solution for reference.
  !! note: ndofs is not used by this routine. so, to evaluate spacetime function
  !! for ndofs in an element, use tem_varsys_proc_point interface.
  !!
  !! \verbatim
  !! -- in lua file, one can define it as following:
  !! variable = {{
  !!   name = 'reference',
  !!   ncomponents = 1,
  !!   vartype = st_fun,
  !!   st_fun =  luafun
  !!   },
  !!   ...
  !!  }
  !! \endverbatim
  !!
  !! the interface has to comply to the abstract interface
  !! tem_varsys_module#tem_varsys_proc_element.
  !!
  recursive subroutine evaluate_first_spacetime_scalarByTreeID( 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 :: iStFun, iElem, pos, oldPos
    logical :: found = .false.
    type(tem_st_fun_listElem_type), pointer :: fPtr
    ! element wise output
    real(kind=rk) :: st_res(1)
    integer(kind=long_k) :: treeID(1)
    ! -------------------------------------------------------------------------!
    call C_F_POINTER( fun%method_Data, fPtr )
    res = 0.0_rk
    ! assuming nDofs = 1
    if (nDofs > 1) then

    ! avoid compiler warning about unused varsys
      write(logunit(1),*) 'Variable:', trim(varsys%varname%val(fun%mypos))
      write(logunit(1),*) 'tree%nElems:', tree%nElems
      write(logUnit(1),*) 'Spacetime function does not support nDofs>1'
      call tem_abort()
    end if
    do iElem = 1, nElems
      found = .false.
      stLoop: do iStFun = 1, fPtr%nVals
        ! If the space time function has a global shape, we can leave out the
        ! search for the elemPos.
        if(fPtr%val(iStFun)%subTree%useGlobalMesh) then
          found = .true.
        else
          oldPos = 1
          ! position of element matches with any of map2global then
          ! this element is part of this subTree
          pos = tem_PositionInSorted(                                      &
            &                 me    = fPtr%val(iStFun)%subTree%map2global, &
            &                 val   = elemPos(iElem),                      &
            &                 lower = oldPos                               )
          if (pos > 0) then
            oldPos = pos
            found = .true.
          end if ! elemPos(iElem) part of subtree
        end if ! global mesh

        ! We found the one value we were looking for, so we can stop now
        if (found) then
          treeID = tree%treeID(elemPos(iElem))
          st_res = tem_spacetime_for( me      = fPtr%val(iStFun),     &
            &                         treeIDs = treeID,               &
            &                         tree    = tree,                 &
            &                         time    = time,                 &
            &                         n       = 1                     )
          res(iElem) = st_res(1)
          exit stLoop
        end if !found
      end do stLoop
    end do ! iElem
  end subroutine evaluate_first_spacetime_scalarByTreeID
  ! *************************************************************************** !


  ! *************************************************************************** !
  !> This subroutine directly evaluate spacetime function for given element
  !! position in global TreeID list on 1st st_fun assuming global shape for
  !! vectorial variable.
  recursive subroutine evaluate_FOAG_spacetime_vectorByTreeID( 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_st_fun_listElem_type), pointer :: fPtr
    integer :: iElem
    real(kind=rk) :: st_res(nElems, fun%nComponents)
    ! -------------------------------------------------------------------------!
    call C_F_POINTER( fun%method_Data, fPtr )
    ! assuming nDofs = 1
    if (nDofs > 1) then
      write(logunit(1),*) 'Variable:', trim(varsys%varname%val(fun%mypos))
      write(logunit(1),*) 'tree%nElems:', tree%nElems
      write(logUnit(1),*) 'Spacetime function does not support nDofs>1'
      call tem_abort()
    end if

    st_res = tem_spacetime_for( me      = fPtr%val(1),             &
      &                         treeIDs = tree%treeID(elemPos(:)), &
      &                         tree    = tree,                    &
      &                         time    = time,                    &
      &                         nComp   = fun%nComponents,         &
      &                         n       = nElems                   )


    do iElem = 1, nElems
      res( (iElem-1)*fun%nComponents + 1 : iElem*fun%nComponents )             &
        & = st_res(iElem, :)
    end do

  end subroutine evaluate_FOAG_spacetime_vectorByTreeID
  ! *************************************************************************** !


  ! *************************************************************************** !
  !> Call the spacetime function for vector variable,
  !! which are stored in method_data which intern
  !! points to tem_st_fun_listElem_type.
  !! This spacetime function can be used as a analytical solution for reference.
  !! NOTE: nDofs is not used by this routine. So, to evaluate spacetime function
  !! for nDofs in an element, use tem_varSys_proc_point interface.
  !!
  !! \verbatim
  !! -- in lua file, one can define it as following:
  !! variable = {{
  !!   name = 'reference',
  !!   ncomponents = 3,
  !!   vartype = st_fun,
  !!   st_fun =  luaFun
  !!   },
  !!   ...
  !!  }
  !! \endverbatim
  !!
  !! The interface has to comply to the abstract interface
  !! tem_varSys_module#tem_varSys_proc_element.
  !!
  recursive subroutine evaluate_add_spacetime_vectorByTreeID( 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 :: iStFun, iElem, pos, oldPos
    type(tem_st_fun_listElem_type), pointer :: fPtr
    ! element wise output
    real(kind=rk), allocatable :: st_res(:,:)
    real(kind=rk) :: st_res_elem(1, fun%nComponents)
    integer(kind=long_k) :: treeID(1)
    ! -------------------------------------------------------------------------!
    !write(*,*) 'derive variable label: ', trim(varSys%varname%val(fun%myPos))
    call C_F_POINTER( fun%method_Data, fPtr )
    res = 0.0_rk
    ! assuming nDofs = 1
    if (nDofs > 1) then
      write(logunit(1),*) 'Variable:', trim(varsys%varname%val(fun%mypos))
      write(logUnit(1),*) 'Spacetime function does not support nDofs>1'
      call tem_abort()
    end if

    allocate(st_res(nElems, fun%nComponents))
    st_res = 0.0_rk
    do iStFun = 1, fPtr%nVals
      if (fPtr%val(iStFun)%subTree%useGlobalMesh ) then
        !write(*,*) 'global mesh'
        st_res = st_res +                                         &
         &  tem_spacetime_for( me      = fPtr%val(iStFun),        &
         &                     treeIDs = tree%treeID(elemPos(:)), &
         &                     tree    = tree,                    &
         &                     time    = time,                    &
         &                     nComp   = fun%nComponents,         &
         &                     n       = nElems                   )
      else
        !write(*,*) 'subTree Mesh'
        oldPos = 1
        do iElem = 1, nElems
          ! position of element matches with any of map2global then
          ! this element is part of this subTree
          pos = tem_PositionInSorted(  &
            &                 me    = fPtr%val(iStFun)%subTree%map2global, &
            &                 val   = elemPos(iElem),                      &
            &                 lower = oldPos                               )
          ! pos>0 if elemPos(iElem) is part of subTree
          if (pos > 0) then
            oldPos = pos
            treeID = tree%treeID(elemPos(iElem))
            st_res_elem =                                           &
              &  tem_spacetime_for( me      = fPtr%val(iStFun),     &
              &                     treeIDs = treeID,               &
              &                     tree    = tree,                 &
              &                     time    = time,                 &
              &                     nComp   = fun%nComponents,      &
              &                     n       = 1                     )
            st_res(iElem,:) = st_res(iElem,:) +  st_res_elem(1,:)
          end if
        end do
      end if ! global mesh
    end do ! iStFun

    do iElem = 1, nElems
      res( (iElem-1)*fun%nComponents + 1 : iElem*fun%nComponents )             &
        & = st_res(iElem, :)
    end do

  end subroutine evaluate_add_spacetime_vectorByTreeID
  ! *************************************************************************** !


  ! *************************************************************************** !
  !> Call the spacetime function for vector variable,
  !! which are stored in method_data which intern
  !! points to tem_st_fun_listElem_type.
  !! This spacetime function can be used as a analytical solution for reference.
  !! NOTE: nDofs is not used by this routine. So, to evaluate spacetime function
  !! for nDofs in an element, use tem_varSys_proc_point interface.
  !!
  !! \verbatim
  !! -- in lua file, one can define it as following:
  !! variable = {{
  !!   name = 'reference',
  !!   ncomponents = 3,
  !!   vartype = st_fun,
  !!   st_fun =  luaFun
  !!   },
  !!   ...
  !!  }
  !! \endverbatim
  !!
  !! The interface has to comply to the abstract interface
  !! tem_varSys_module#tem_varSys_proc_element.
  !!
  recursive subroutine evaluate_first_spacetime_vectorByTreeID( 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 :: iStFun, iElem, pos, oldPos
    logical :: found
    type(tem_st_fun_listElem_type), pointer :: fPtr
    ! element wise output
    real(kind=rk) :: st_res(1, fun%nComponents)
    integer(kind=long_k) :: treeID(1)
    ! -------------------------------------------------------------------------!
    !write(*,*) 'derive variable label: ', trim(varSys%varname%val(fun%myPos))
    call C_F_POINTER( fun%method_Data, fPtr )
    res = 0.0_rk
    ! assuming nDofs = 1
    if (nDofs > 1) then
      write(logunit(1),*) 'Variable:', trim(varsys%varname%val(fun%mypos))
      write(logUnit(1),*) 'Spacetime function does not support nDofs>1'
      call tem_abort()
    end if

    do iElem = 1, nElems
      found = .false.
      stLoop: do iStFun = 1, fPtr%nVals
        if (fPtr%val(iStFun)%subTree%useGlobalMesh) then
          found = .true.
        else
          oldPos = 1
          ! position of element matches with any of map2global then
          ! this element is part of this subTree
          pos = tem_PositionInSorted(                      &
            & me    = fPtr%val(iStFun)%subTree%map2global, &
            & val   = elemPos(iElem),                      &
            & lower = oldPos                               )
          ! pos>0 if elemPos(iElem) is part of subTree
          if (pos > 0) then
            oldPos = pos
            found = .true.
          end if
        end if ! global mesh

        ! We found the one value we were looking for, so we can stop now
        if (found) then
          treeID = tree%treeID(elemPos(iElem))
          st_res = tem_spacetime_for( me      = fPtr%val(iStFun), &
            &                         treeIDs = treeID,           &
            &                         tree    = tree,             &
            &                         time    = time,             &
            &                         nComp   = fun%nComponents,  &
            &                         n       = 1                 )

          res( (iElem-1)*fun%nComponents + 1 : iElem*fun%nComponents ) &
            & = st_res(1,:)

          exit stLoop
        end if ! found

      end do stLoop
    end do !iElem

  end subroutine evaluate_first_spacetime_vectorByTreeID
  ! ************************************************************************** !


  ! ************************************************************************** !
  !> This routine process instring and store information in spacetime function.
  !!
  !! Abort if same spacetime function is used as boundaries and sources.
  recursive subroutine set_params_spacetime( 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
    ! -------------------------------------------------------------------- !
    type(flu_state) :: conf
    type(tem_st_fun_listElem_type), pointer :: fPtr
    logical :: isSurface_loc
    integer :: iStFun
    integer :: iError
    ! -------------------------------------------------------------------- !
    call C_F_POINTER( fun%method_Data, fPtr )

    call open_config_chunk(L = conf, chunk = instring)
    call aot_get_val( L       = conf,          &
      &               key     = 'isSurface',   &
      &               val     = isSurface_loc, &
      &               ErrCode = iError         )

    if (btest(iError,  aoterr_Fatal)) then
      write(logUnit(1),*) 'Error: In set_params to load "isSurface" '//&
        &                 'for spacetime variable ', &
        &                 trim(varSys%varname%val(fun%mypos))
      if (btest(iError, aoterr_WrongType)) then
        write(logUnit(1),*) '       "isSurface" is wrong type'
      end if
      if (btest(iError, aoterr_NonExistent)) then
        write(logUnit(1),*) '       "isSurface" does not exist'
      end if
      call tem_abort()
    end if

    ! KM Store bc_id and field id only for apesmate coupling
    do iStFun = 1, fPtr%nVals
      select case (trim(fPtr%val(iStFun)%fun_kind))
      case ('apesmate')
        if (fPtr%isSurface == -1) then
          if (isSurface_loc) then
            fPtr%isSurface = 0
          else
            fPtr%isSurface = 1
          end if
        else
          write(logUnit(10),*) 'Warning: ' &
            &                  // 'Set params stfun: isSurface is already set'
          if (isSurface_loc) then
            if (fPtr%isSurface == 1) then
              write(logUnit(1),*) 'Error: In set params for spacetime function'
              write(logUnit(1),*) '       This st_fun is already used for' &
                &                 // ' volume'
              write(logUnit(1),*) '       so cannot be used for surface.'
              call tem_abort()
            end if
          else
            if (fPtr%isSurface == 0) then
              write(logUnit(1),*) 'Error: In set params for spacetime function'
              write(logUnit(1),*) '       This st_fun is already used for' &
                &                 // ' surface'
              write(logUnit(1),*) '       so cannot be used for volume.'
              call tem_abort()
            end if
          end if
        end if

        ! store surface or volume information in coupling type
        fPtr%val(iStFun)%aps_coupling%isSurface = fPtr%isSurface

      end select
    end do

  end subroutine set_params_spacetime
  ! ************************************************************************** !


  ! ************************************************************************** !
  !> This routine process instring and return string with requested info
  !! from spacetime function
  recursive subroutine get_params_spacetime(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
    ! -------------------------------------------------------------------- !
    select case (trim(instring))
    case ('vartype')
      ! return vartype
      outstring = 'st_fun'
    case default
      outstring = trim(varsys%varname%val(fun%mypos)) // '_UNKNOWN'
    end select
  end subroutine get_params_spacetime
  ! ************************************************************************** !


  ! ************************************************************************** !
  !> This routine stores provided points in method_data of spacetime_listElem
  !! and return the indices of points or evaluated value in the growing array.
  !! If spacetime function is time-independent then pre-compute values
  !! and store in growing array of evalVal in tem_pointData_type.
  recursive subroutine setup_indices_spacetime( 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_st_fun_listElem_type), pointer :: fPtr
    integer :: iStFun, iPnt, posInTree, nUniquePnts, iVar
    character :: offset_bit_local
    logical, allocatable :: storePnt(:), storeOffsetBit(:), storeVal(:)
    integer :: elemPos
    integer(kind=long_k) :: treeID
    logical :: addPoint, wasAdded
    real(kind=rk) :: uniquePnts(nPnts,3)
    ! -------------------------------------------------------------------------!


    call C_F_POINTER( fun%method_Data, fPtr )

    allocate(storePnt(fPtr%nVals))
    allocate(storeVal(fPtr%nVals))
    allocate(storeOffsetBit(fPtr%nVals))

    ! Store points only for time dependent spacetime functions
    do iStFun = 1, fPtr%nVals
      select case (trim(fPtr%val(iStFun)%fun_kind))
      case ('none', 'const')
        ! time independent and no need to store points
        storePnt(iStFun) = .false.
        storeOffsetBit(iStFun) = .false.
        storeVal(iStFun) = .false.
      case ('combined')
        ! spatial is time independent so compute spatial value and store it.
        storePnt(iStFun) = .false.
        storeOffsetBit(iStFun) = .false.
        storeVal(iStFun) = .true.

      case ('lua_fun', 'miescatter_displacementfieldz',               &
        &   'miescatter_magneticfieldx', 'miescatter_magneticfieldy', &
        &   'cylindrical_wave')
        storePnt(iStFun) = .true.
        storeOffsetBit(iStFun) = .false.
        storeVal(iStFun) = .false.
      case ('apesmate')
        storePnt(iStFun) = .true.
        storeOffsetBit(iStFun) = (fPtr%val(iStFun)%aps_coupling%isSurface==0)
        storeVal(iStFun) = .false.
      case ('precice')
        storePnt(iStFun) = .true.
        storeOffsetBit(iStFun) = .true.
        storeVal(iStFun) = .true.
        ! for writing to precice, we need the position of the variables in
        ! the variable system
        do iVar = 1, fPtr%val(iStFun)%precice_coupling%writeVar%nVars
          fPtr%val(iStFun)%precice_coupling%writeVar%varPos(iVar) =     &
             PositionOfVal( me = varSys%varName,                        &
               &            val = fPtr%val(iStFun)%precice_coupling     &
               &                                  %writeVar%names(iVar) )
             if (fPtr%val(iStFun)%precice_coupling%writeVar%varPos(iVar) == 0) &
               & then
               write(*,*) 'position in Varsys for writing variable ',          &
                 & trim(fPtr%val(iStFun)%precice_coupling%writeVar%names(iVar)),&
                 & ' not found'
                 call tem_abort()
             end if
        end do

        ! same for reading from precice
        do iVar = 1, fPtr%val(iStFun)%precice_coupling%readVar%nVars
          fPtr%val(iStFun)%precice_coupling%readVar%varPos(iVar) =     &
             PositionOfVal( me = varSys%varName,                       &
               &            val = fPtr%val(iStFun)%precice_coupling    &
               &                                  %readVar%names(iVar) )
             if (fPtr%val(iStFun)%precice_coupling%readVar%varPos(iVar) == 0) &
               & then
               write(*,*) 'position in Varsys for reading ',                   &
                 & trim(fPtr%val(iStFun)%precice_coupling%readVar%names(iVar)),&
                 & 'not found'
                 call tem_abort()
             end if
        end do

      case default
        write(logUnit(1),*)'ERROR: Unknown spatial function in '// &
          &                'setup_indices.'
        call tem_abort()
      end select
    end do !iStFun

    ! initialize index with zero to identify points which does not
    ! belong to subTree
    idx = 0

    ! number of unique points added
    nUniquePnts = 0
    do iPnt = 1, nPnts
      addPoint = .false.

      ! get treeID from globalmaxLevel since points from ghost elements are
      ! also passed to this routine
      treeID = tem_IdOfCoord( tem_CoordOfReal(tree, point(iPnt, :)) )
      elemPos = tem_PosOfId( treeID, tree%treeID )
      !if (elemPos == 0)  then
       ! call tem_abort('Error: treeID not found in st-fun setup_indices')
     ! end if

      ! if any spacetime function has useGlobalMesh then store all points
      ! else store only points which belong to subTree of spacetime variable.
      if ( any(fPtr%val(:)%subTree%useGlobalMesh) ) then
        addPoint = .true.
      else
        stFunLoop: do iStFun = 1, fPtr%nVals
          posInTree = tem_PositionInSorted(                                &
            &                 me    = fPtr%val(iStFun)%subTree%map2global, &
            &                 val   = elemPos                              )
          if (posInTree > 0) then
            addPoint = .true.
            exit stFunLoop
          end if
        end do stFunLoop
      end if !use global mesh

      if (addPoint) then
        ! use center offset bit as default
        if (present(offset_bit)) then
          offset_bit_local = offset_bit(iPnt)
        else
          offset_bit_local = qOffset_inChar(q000)
        end if

        ! append point, offset_bit and elemPos to pointData type
        call append(me             = fPtr%pntData%pntLvl(iLevel), &
          &         point          = point(iPnt,:),               &
          &         storePnt       = any(storePnt),               &
          &         offset_bit     = offset_bit_local,            &
          &         storeOffsetBit = any(storeOffsetBit),         &
          &         elemPos        = elemPos,                     &
          &         tree           = tree,                        &
          &         pos            = idx(iPnt),                   &
          &         wasAdded       = wasAdded                     )

        if (wasAdded) then
          nUniquePnts = nUniquePnts + 1
          uniquePnts(nUniquePnts,:) = point(iPnt,:)
        end if

      end if ! add point
    end do !iPnt

    if (any(storePnt)) call truncate(fPtr%pntData%pntLvl(iLevel)%grwPnt)
    if (any(storeOffsetBit)) &
      & call truncate(fPtr%pntData%pntLvl(iLevel)%offset_bit)

    deallocate(storePnt)
    deallocate(storeOffsetBit)
    call truncate(fPtr%pntData%pntLvl(iLevel)%treeID)
    call truncate(fPtr%pntData%pntLvl(iLevel)%elemPos)

    !!! Store spatial value for unique points depends on stFun type
    if ( any(storeVal) .and. nUniquePnts > 0 ) then
      do iStFun = 1, fPtr%nVals
        select case (trim(fPtr%val(iStFun)%fun_kind))
        case ('combined')
          if (fun%nComponents == 1) then
            call tem_spatial_storeVal( me     = fPtr%val(iStFun)%spatial,    &
              &                        coord  = uniquePnts(1:nUniquePnts,:), &
              &                        nVals  = nUniquePnts,                 &
              &                        iLevel = iLevel                       )
          else
            call tem_spatial_storeVal( me     = fPtr%val(iStFun)%spatial,    &
              &                        coord  = uniquePnts(1:nUniquePnts,:), &
              &                        nVals  = nUniquePnts,                 &
              &                        iLevel = iLevel,                      &
              &                        nComps = fun%nComponents              )
          end if
        end select
      end do !iStFun
    end if !store spatial value
    deallocate(storeVal)

  end subroutine setup_indices_spacetime
  ! *************************************************************************** !


  ! ************************************************************************** !
  !> This routine returns value at the given index. If value is pre-computed
  !! and value at given index is returned else value is computing on the
  !! points for given index.
  !!
  !! If a variable has more than one spacetime function, then result of 1st
  !! spacetime function is returned with assumed global shape
  recursive subroutine get_valOfIndex_FOAG_scalar_spacetime( 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_st_fun_listElem_type), pointer :: fPtr
    ! -------------------------------------------------------------------------!
    call C_F_POINTER( fun%method_Data, fPtr )

    ! 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 = 'get_valOfIndex_FOAG_scalar_spacetime'          )

    res(:) = tem_spacetime_for( me     = fPtr%val(1),                &
      &                         grwPnt = fPtr%pntData%pntLvl(iLevel) &
      &                                      %grwPnt,                &
      &                         idx    = idx,                        &
      &                         nVals  = nVals,                      &
      &                         iLevel = ilevel,                     &
      &                         time   = time                        )
  end subroutine get_valOfIndex_FOAG_scalar_spacetime
  ! ************************************************************************ !


  ! *************************************************************************** !
  !> This routine returns value at the given index. If value is pre-computed
  !! and value at given index is returned else value is computing on the
  !! points for given index.
  !!
  !! If a variable has more than one spacetime function, then result of 1st
  !! spacetime function is returned with assumed global shape
  !!
  !! For vectorial variable.
  recursive subroutine get_valOfIndex_FOAG_vector_spacetime( 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_st_fun_listElem_type), pointer :: fPtr
    integer :: iVal
    real(kind=rk) :: st_res(nVals, fun%nComponents)
    ! -------------------------------------------------------------------------!
    call C_F_POINTER( fun%method_Data, fPtr )

    ! 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 = 'get_valOfIndex_FOAG_vector_spacetime'          )

    st_res = tem_spacetime_for( me     = fPtr%val(1),                &
      &                         grwPnt = fPtr%pntData%pntLvl(iLevel) &
      &                                      %grwPnt,                &
      &                         idx    = idx,                        &
      &                         nVals  = nVals,                      &
      &                         iLevel = ilevel,                     &
      &                         time   = time,                       &
      &                         nComps = fun%nComponents             )

    do iVal = 1, nVals
      res( (iVal-1)*fun%nComponents + 1 : iVal*fun%nComponents ) &
        & = st_res(iVal, :)
    end do

  end subroutine get_valOfIndex_FOAG_vector_spacetime
  ! ************************************************************************ !


  ! *************************************************************************** !
  !> This routine returns value at the given index. If value is pre-computed
  !! and value at given index is returned else value is computing on the
  !! points for given index.
  !!
  !! If a variable has more than one spacetime function, then result is sum
  !! of all spacetime functions
  recursive subroutine get_valOfIndex_add_scalar_spacetime( 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 :: iStFun, iVal, pos, oldPos
    type(tem_st_fun_listElem_type), pointer :: fPtr
    ! element wise output
    real(kind=rk) :: st_res(1)
    ! -------------------------------------------------------------------------!
    call C_F_POINTER( fun%method_Data, fPtr )

    ! 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 = 'get_valOfIndex_add_scalar_spacetime'           )

    res = 0.0_rk
    do iStFun = 1, fPtr%nVals
      if (fPtr%val(iStFun)%subTree%useGlobalMesh) then
        res(:) = res(:)                                                  &
          &    + tem_spacetime_for( me     = fPtr%val(iStFun),           &
          &                         grwPnt = fPtr%pntData%pntLvl(iLevel) &
          &                                      %grwPnt,                &
          &                         idx    = idx,                        &
          &                         nVals  = nVals,                      &
          &                         iLevel = ilevel,                     &
          &                         time   = time                        )
      else
        oldPos = 1
        do iVal = 1, nVals

          if (idx(iVal) > 0) then
            ! position of element matches with any of map2global then
            ! this element is part of this subTree
            pos = tem_PositionInSorted(  &
              &                 me    = fPtr%val(iStFun)%subTree%map2global, &
              &                 val   = fPtr%pntData%pntLvl(iLevel)          &
              &                             %elemPos%val(idx(iVal)),         &
              &                 lower = oldPos                               )
            if (pos > 0) then
              oldPos = pos
              st_res = tem_spacetime_for( me     = fPtr%val(iStFun),           &
                &                         grwPnt = fPtr%pntData%pntLvl(iLevel) &
                &                                      %grwPnt,                &
                &                         idx    = (/ idx(iVal) /),            &
                &                         nVals  = 1,                          &
                &                         iLevel = ilevel,                     &
                &                         time   = time                        )
              res(iVal) = res(iVal) + st_res(1)
            end if ! elemPos(iElem) part of subtree
          end if !idx > 0
        end do !iVal
      end if ! global mesh
    end do ! nr. st_funs

  end subroutine get_valOfIndex_add_scalar_spacetime
  ! ************************************************************************ !


  ! *************************************************************************** !
  !> This routine returns value at the given index. If value is pre-computed
  !! and value at given index is returned else value is computing on the
  !! points for given index.
  !!
  !! If a variable has more than one spacetime function, then result is sum
  !! of all spacetime functions
  !!
  !! For vectorial variable
  recursive subroutine get_valOfIndex_add_vector_spacetime( 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 :: iStFun, iVal, pos, oldPos
    type(tem_st_fun_listElem_type), pointer :: fPtr
    real(kind=rk) :: st_res(nVals, fun%nComponents)
    real(kind=rk) :: st_res_elem(1, fun%nComponents)
    ! -------------------------------------------------------------------------!
    call C_F_POINTER( fun%method_Data, fPtr )

    ! 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 = 'get_valOfIndex_add_vector_spacetime'           )

    st_res = 0.0_rk
    do iStFun = 1, fPtr%nVals
      if (fPtr%val(iStFun)%subTree%useGlobalMesh) then
        st_res = st_res                                                  &
          &    + tem_spacetime_for( me     = fPtr%val(iStFun),           &
          &                         grwPnt = fPtr%pntData%pntLvl(iLevel) &
          &                                      %grwPnt,                &
          &                         idx    = idx,                        &
          &                         nVals  = nVals,                      &
          &                         iLevel = ilevel,                     &
          &                         time   = time,                       &
          &                         nComps = fun%nComponents             )
      else
        oldPos = 1
        do iVal = 1, nVals
          if (idx(iVal) > 0) then
            ! position of element matches with any of map2global then
            ! this element is part of this subTree
            pos = tem_PositionInSorted(  &
              &                 me    = fPtr%val(iStFun)%subTree%map2global, &
              &                 val   = fPtr%pntData%pntLvl(iLevel)          &
              &                             %elemPos%val(idx(iVal)),         &
              &                 lower = oldPos                               )
            if (pos > 0) then
              oldPos = pos
              st_res_elem =                                               &
                & tem_spacetime_for( me     = fPtr%val(iStFun),           &
                &                    grwPnt = fPtr%pntData%pntLvl(iLevel) &
                &                                 %grwPnt,                &
                &                    idx    = (/ idx(iVal) /),            &
                &                    nVals  = 1,                          &
                &                    iLevel = ilevel,                     &
                &                    time   = time,                       &
                &                    nComps = fun%nComponents             )
              st_res(iVal,:) = st_res(iVal,:) + st_res_elem(1,:)
            end if ! elemPos(iElem) part of subtree
          end if
        end do !iVal
      end if ! global mesh
    end do ! nr. st_funs

    do iVal = 1, nVals
      res( (iVal-1)*fun%nComponents + 1 : iVal*fun%nComponents ) &
        & = st_res(iVal, :)
    end do

  end subroutine get_valOfIndex_add_vector_spacetime
  ! *************************************************************************** !


  ! *************************************************************************** !
  !> This routine returns value at the given index. If value is pre-computed
  !! and value at given index is returned else value is computing on the
  !! points for given index.
  !!
  !! If a variable has more than one spacetime function, then result is
  !! first spacetime function in the overlapping region
  recursive subroutine get_valOfIndex_first_scalar_spacetime( 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 :: iStFun, iVal, pos
    type(tem_st_fun_listElem_type), pointer :: fPtr
    ! element wise output
    real(kind=rk) :: st_res(1)
    logical :: found
    ! -------------------------------------------------------------------------!
    call C_F_POINTER( fun%method_Data, fPtr )
    res = 0.0_rk

    ! 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 = 'get_valOfIndex_first_scalar_spacetime'         )

    do iVal = 1, nVals
      found = .false.
      stLoop: do iStFun = 1, fPtr%nVals
        if (fPtr%val(iStFun)%subTree%useGlobalMesh) then
          found = .true.
        else
          if (idx(iVal) > 0) then
            ! position of element matches with any of map2global then
            ! this element is part of this subTree
            pos = tem_PositionInSorted(  &
              &                 me    = fPtr%val(iStFun)%subTree%map2global, &
              &                 val   = fPtr%pntData%pntLvl(iLevel)          &
              &                             %elemPos%val(idx(iVal))          )

            found = (pos > 0)
          end if !idx>0
        end if ! global mesh
        ! We found the one value we were looking for, so we can stop now
        if (found) then
          st_res = tem_spacetime_for( me     = fPtr%val(iStFun),           &
            &                         grwPnt = fPtr%pntData%pntLvl(iLevel) &
            &                                      %grwPnt,                &
            &                         idx    = (/ idx(iVal) /),            &
            &                         nVals  = 1,                          &
            &                         iLevel = ilevel,                     &
            &                         time   = time                        )
          res(iVal) = st_res(1)
          exit stLoop
        end if
      end do stLoop
    end do !iVal

  end subroutine get_valOfIndex_first_scalar_spacetime
  ! *************************************************************************** !


  ! *************************************************************************** !
  !> This routine returns value at the given index. If value is pre-computed
  !! and value at given index is returned else value is computing on the
  !! points for given index.
  !!
  !! If a variable has more than one spacetime function, then result is
  !! first spacetime function in the overlapping region.
  !!
  !! For vectorial variable
  recursive subroutine get_valOfIndex_first_vector_spacetime( 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 :: iStFun, iVal, pos
    type(tem_st_fun_listElem_type), pointer :: fPtr
    ! element wise output
    real(kind=rk) :: st_res(1, fun%nComponents)
    logical :: found
    ! -------------------------------------------------------------------------!
    call C_F_POINTER( fun%method_Data, fPtr )
    res = 0.0_rk

    ! 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 = 'get_valOfIndex_first_vector_spacetime'         )

    do iVal = 1, nVals
      found = .false.
      stLoop: do iStFun = 1, fPtr%nVals
        if (fPtr%val(iStFun)%subTree%useGlobalMesh) then
          found = .true.
        else
          if (idx(iVal) > 0) then
            ! position of element matches with any of map2global then
            ! this element is part of this subTree
            pos = tem_PositionInSorted(  &
              &                 me    = fPtr%val(iStFun)%subTree%map2global, &
              &                 val   = fPtr%pntData%pntLvl(iLevel)          &
              &                             %elemPos%val(idx(iVal))          )

            found = (pos > 0)
          end if
        end if ! global mesh

        ! We found the one value we were looking for, so we can stop now
        if (found) then
          st_res = tem_spacetime_for( me     = fPtr%val(iStFun),           &
            &                         grwPnt = fPtr%pntData%pntLvl(iLevel) &
            &                                      %grwPnt,                &
            &                         idx    = (/ idx(iVal) /),            &
            &                         nVals  = 1,                          &
            &                         iLevel = ilevel,                     &
            &                         time   = time,                       &
            &                         nComps = fun%nComponents             )

          res( (iVal-1)*fun%nComponents + 1 : iVal*fun%nComponents ) &
            & = st_res(1, :)

          exit stLoop
        end if ! found
      end do stLoop
    end do !iVal

  end subroutine get_valOfIndex_first_vector_spacetime
  ! *************************************************************************** !


end module tem_spacetime_var_module
! ****************************************************************************** !