tem_varSys_module.f90 Source File


This file depends on

sourcefile~~tem_varsys_module.f90~~EfferentGraph sourcefile~tem_varsys_module.f90 tem_varSys_module.f90 sourcefile~tem_dyn_array.f90 tem_dyn_array.f90 sourcefile~tem_varsys_module.f90->sourcefile~tem_dyn_array.f90 sourcefile~tem_time_module.f90 tem_time_module.f90 sourcefile~tem_varsys_module.f90->sourcefile~tem_time_module.f90 sourcefile~tem_logging_module.f90 tem_logging_module.f90 sourcefile~tem_varsys_module.f90->sourcefile~tem_logging_module.f90 sourcefile~tem_aux_module.f90 tem_aux_module.f90 sourcefile~tem_varsys_module.f90->sourcefile~tem_aux_module.f90 sourcefile~env_module.f90 env_module.f90 sourcefile~tem_varsys_module.f90->sourcefile~env_module.f90 sourcefile~treelmesh_module.f90 treelmesh_module.f90 sourcefile~tem_varsys_module.f90->sourcefile~treelmesh_module.f90 sourcefile~tem_dyn_array.f90->sourcefile~env_module.f90 sourcefile~tem_time_module.f90->sourcefile~env_module.f90 sourcefile~tem_logging_module.f90->sourcefile~env_module.f90 sourcefile~tem_aux_module.f90->sourcefile~tem_logging_module.f90 sourcefile~tem_aux_module.f90->sourcefile~env_module.f90 sourcefile~tem_lua_requires_module.f90 tem_lua_requires_module.f90 sourcefile~tem_aux_module.f90->sourcefile~tem_lua_requires_module.f90 sourcefile~tem_tools_module.f90 tem_tools_module.f90 sourcefile~tem_aux_module.f90->sourcefile~tem_tools_module.f90 sourcefile~tem_comm_env_module.f90 tem_comm_env_module.f90 sourcefile~tem_aux_module.f90->sourcefile~tem_comm_env_module.f90 sourcefile~tem_revision_module.f90 tem_revision_module.f90 sourcefile~tem_aux_module.f90->sourcefile~tem_revision_module.f90 sourcefile~treelmesh_module.f90->sourcefile~tem_logging_module.f90 sourcefile~treelmesh_module.f90->sourcefile~tem_aux_module.f90 sourcefile~treelmesh_module.f90->sourcefile~env_module.f90 sourcefile~tem_property_module.f90 tem_property_module.f90 sourcefile~treelmesh_module.f90->sourcefile~tem_property_module.f90 sourcefile~treelmesh_module.f90->sourcefile~tem_tools_module.f90 sourcefile~tem_global_module.f90 tem_global_module.f90 sourcefile~treelmesh_module.f90->sourcefile~tem_global_module.f90 sourcefile~tem_sparta_module.f90 tem_sparta_module.f90 sourcefile~treelmesh_module.f90->sourcefile~tem_sparta_module.f90 sourcefile~tem_topology_module.f90 tem_topology_module.f90 sourcefile~treelmesh_module.f90->sourcefile~tem_topology_module.f90 sourcefile~tem_property_module.f90->sourcefile~env_module.f90 sourcefile~tem_prophead_module.f90 tem_prophead_module.f90 sourcefile~tem_property_module.f90->sourcefile~tem_prophead_module.f90 sourcefile~tem_lua_requires_module.f90->sourcefile~env_module.f90 sourcefile~tem_tools_module.f90->sourcefile~env_module.f90 sourcefile~tem_global_module.f90->sourcefile~tem_logging_module.f90 sourcefile~tem_global_module.f90->sourcefile~tem_aux_module.f90 sourcefile~tem_global_module.f90->sourcefile~env_module.f90 sourcefile~tem_global_module.f90->sourcefile~tem_prophead_module.f90 sourcefile~tem_sparta_module.f90->sourcefile~tem_logging_module.f90 sourcefile~tem_sparta_module.f90->sourcefile~tem_aux_module.f90 sourcefile~tem_sparta_module.f90->sourcefile~env_module.f90 sourcefile~tem_float_module.f90 tem_float_module.f90 sourcefile~tem_sparta_module.f90->sourcefile~tem_float_module.f90 sourcefile~tem_topology_module.f90->sourcefile~env_module.f90 sourcefile~tem_float_module.f90->sourcefile~env_module.f90 sourcefile~tem_prophead_module.f90->sourcefile~env_module.f90

Files dependent on this one

sourcefile~~tem_varsys_module.f90~~AfferentGraph sourcefile~tem_varsys_module.f90 tem_varSys_module.f90 sourcefile~tem_logical_operator_test.f90 tem_logical_operator_test.f90 sourcefile~tem_logical_operator_test.f90->sourcefile~tem_varsys_module.f90 sourcefile~tem_logical_operation_var_module.f90 tem_logical_operation_var_module.f90 sourcefile~tem_logical_operation_var_module.f90->sourcefile~tem_varsys_module.f90 sourcefile~tem_convergence_module.f90 tem_convergence_module.f90 sourcefile~tem_convergence_module.f90->sourcefile~tem_varsys_module.f90 sourcefile~hvs_vtk_module.f90 hvs_vtk_module.f90 sourcefile~hvs_vtk_module.f90->sourcefile~tem_varsys_module.f90 sourcefile~tem_varsys_opvar_test.f90 tem_varSys_opVar_test.f90 sourcefile~tem_varsys_opvar_test.f90->sourcefile~tem_varsys_module.f90 sourcefile~tem_coordinate_module.f90 tem_coordinate_module.f90 sourcefile~tem_coordinate_module.f90->sourcefile~tem_varsys_module.f90 sourcefile~tem_variable_evaltype_test.f90 tem_variable_evaltype_test.f90 sourcefile~tem_variable_evaltype_test.f90->sourcefile~tem_varsys_module.f90 sourcefile~tem_variable_combine_test.f90 tem_variable_combine_test.f90 sourcefile~tem_variable_combine_test.f90->sourcefile~tem_varsys_module.f90 sourcefile~tem_varsys_statevar_test.f90 tem_varSys_stateVar_test.f90 sourcefile~tem_varsys_statevar_test.f90->sourcefile~tem_varsys_module.f90 sourcefile~tem_variable_extract_test.f90 tem_variable_extract_test.f90 sourcefile~tem_variable_extract_test.f90->sourcefile~tem_varsys_module.f90 sourcefile~tem_varmap_module.f90 tem_varMap_module.f90 sourcefile~tem_varmap_module.f90->sourcefile~tem_varsys_module.f90 sourcefile~tem_tracking_module.f90 tem_tracking_module.f90 sourcefile~tem_tracking_module.f90->sourcefile~tem_varsys_module.f90 sourcefile~tem_derived_module.f90 tem_derived_module.f90 sourcefile~tem_derived_module.f90->sourcefile~tem_varsys_module.f90 sourcefile~hvs_vtk_dummy.f90 hvs_vtk_dummy.f90 sourcefile~hvs_vtk_dummy.f90->sourcefile~tem_varsys_module.f90 sourcefile~tem_bc_module.f90 tem_bc_module.f90 sourcefile~tem_bc_module.f90->sourcefile~tem_varsys_module.f90 sourcefile~tem_operation_var_module.f90 tem_operation_var_module.f90 sourcefile~tem_operation_var_module.f90->sourcefile~tem_varsys_module.f90 sourcefile~tem_meshinfo_module.f90 tem_meshInfo_module.f90 sourcefile~tem_meshinfo_module.f90->sourcefile~tem_varsys_module.f90 sourcefile~hvs_ascii_module.f90 hvs_ascii_module.f90 sourcefile~hvs_ascii_module.f90->sourcefile~tem_varsys_module.f90 sourcefile~tem_varsys_derivevar_test.f90 tem_varSys_deriveVar_test.f90 sourcefile~tem_varsys_derivevar_test.f90->sourcefile~tem_varsys_module.f90 sourcefile~tem_depend_module.f90 tem_depend_module.f90 sourcefile~tem_depend_module.f90->sourcefile~tem_varsys_module.f90 sourcefile~tem_restart_module.f90 tem_restart_module.f90 sourcefile~tem_restart_module.f90->sourcefile~tem_varsys_module.f90 sourcefile~tem_varsys_stfunvar_test.f90 tem_varSys_stfunVar_test.f90 sourcefile~tem_varsys_stfunvar_test.f90->sourcefile~tem_varsys_module.f90 sourcefile~tem_varsys_test.f90 tem_varSys_test.f90 sourcefile~tem_varsys_test.f90->sourcefile~tem_varsys_module.f90 sourcefile~tem_reduction_spatial_module.f90 tem_reduction_spatial_module.f90 sourcefile~tem_reduction_spatial_module.f90->sourcefile~tem_varsys_module.f90 sourcefile~hvs_output_module.f90 hvs_output_module.f90 sourcefile~hvs_output_module.f90->sourcefile~tem_varsys_module.f90 sourcefile~tem_spacetime_var_module.f90 tem_spacetime_var_module.f90 sourcefile~tem_spacetime_var_module.f90->sourcefile~tem_varsys_module.f90 sourcefile~tem_variable_module.f90 tem_variable_module.f90 sourcefile~tem_variable_module.f90->sourcefile~tem_varsys_module.f90

Contents

Source Code


Source Code

! Copyright (c) 2014-2016, 2019 Harald Klimach <harald.klimach@uni-siegen.de>
! Copyright (c) 2014-2019 Kannan Masilamani <kannan.masilamani@uni-siegen.de>
! Copyright (c) 2014-2016, 2018-2019 Peter Vitt <peter.vitt2@uni-siegen.de>
! Copyright (c) 2014 Nikhil Anand <nikhil.anand@uni-siegen.de>
! Copyright (c) 2014-2016 Jiaxing Qi <jiaxing.qi@uni-siegen.de>
! Copyright (c) 2016-2017 Tobias Schneider <tobias1.schneider@student.uni-siegen.de>
! Copyright (c) 2016 Verena Krupp <verena.krupp@uni-siegen.de>
! Copyright (c) 2016 Daniel Fleischer <daniel.fleischer@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.
!> The variable system allows the description of arbitrary quantities and their
!! relations.
!!
!! Its meant to enable the extraction of most of the data in the solver, but
!! also to describe additional quantities to be used by the solvers.
module tem_varSys_module
  use iso_c_binding,            only: c_ptr
  use env_module,               only: labelLen, rk, minLength, zeroLength
  use tem_time_module,          only: tem_time_type
  use tem_dyn_array_module,     only: dyn_labelArray_type, init, append, &
    &                                 empty, PositionOfVal
  use treelmesh_module,         only: treelmesh_type
  use tem_aux_module,           only: tem_abort
  use tem_logging_module,       only: logUnit, llerror, lldebug

  use aotus_module,             only: flu_state, aot_get_val
  use aot_table_module,         only: aot_table_open, aot_table_close, &
    &                                 aot_table_length, aot_get_val
  use aot_out_module,           only: aot_out_type, aot_out_open, &
    &                                 aot_out_close, aot_out_val, &
    &                                 aot_out_open_table, aot_out_close_table

  implicit none

! Copyright (c) 2012-2013 Simon Zimny <s.zimny@grs-sim.de>
! Copyright (c) 2012 Manuel Hasert <m.hasert@grs-sim.de>
! Copyright (c) 2012-2016 Harald Klimach <harald.klimach@uni-siegen.de>
! Copyright (c) 2012, 2015-2016 Kannan Masilamani <kannan.masilamani@uni-siegen.de>
! Copyright (c) 2013 Daniel Harlacher <d.harlacher@grs-sim.de>
! Copyright (c) 2014 Kartik Jain <kartik.jain@uni-siegen.de>
! Copyright (c) 2014 Verena Krupp <verena.krupp@uni-siegen.de>
! Copyright (c) 2015-2017 Jiaxing Qi <jiaxing.qi@uni-siegen.de>
! Copyright (c) 2015-2016 Peter Vitt <peter.vitt2@uni-siegen.de>
! Copyright (c) 2016 Daniel Fleischer <daniel.fleischer@student.uni-siegen.de>
! Copyright (c) 2016 Tobias Schneider <tobias1.schneider@student.uni-siegen.de>
! Copyright (c) 2017 Daniel PetrĂ³ <daniel.petro@student.uni-siegen.de>
!
! Parts of this file were written by Harald Klimach, Simon Zimny and Manuel
! Hasert for German Research School for Simulation Sciences GmbH.
!
! Parts of this file were written by Harald Klimach, Kannan Masilamani,
! Daniel Harlacher, Kartik Jain, Verena Krupp, Jiaxing Qi, Peter Vitt,
! Daniel Fleischer, Tobias Girresser and Daniel PetrĂ³ for University Siegen.
!
! 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 file contains the source code for growing and dynamic arrays.
! This is used for arrays of primitives (int, long_int, real, ...) as well as
! for arrays of derived datatypes (tem_variable_type,...).
!
! To use these macros include the following to your source file.
!
! Smart growing array (GA) for ?tstring?
! Growing Arrays:
!
! declaration
!
!
! implementation
!

! -----------------------------------------------------------------
! 2d Array, which can grow in second dimension only (GA2d)
! tname ... indicates type of dynamic array (long, int, real, ...)


!
!------------------------------------------------------------------------------
!
! dynamic Arrays:
!
! declaration
!
!
! implementation
!

  ! **************************************************************************** !
  !> Description of the method how to obtain a variable
  type tem_varSys_op_type

    !> Position of this variable in the variable system.
    integer :: mypos

    !> Position of state variable in the state array
    !! Currently used only in MUSUBI to access one dimensional state array
    integer, allocatable :: state_varPos(:)

    !> Position of auxiliary variable in the auxilied field array
    !! Currently used only in MUSUBI to access one dimensional auxiliary array
    !! In Musubi, auxField vars are conserved macroscopic quantities computed
    !! from PDF state
    integer, allocatable :: auxField_varPos(:)

    !> Number of components for this variable.
    integer :: nComponents

    !> Number of variables, that are needed as input for the operation
    !! to obtain the variable.
    integer :: nInputs

    !> Position of the input variables in the variable system.
    !!
    !! There are as many entries as nInputs.
    integer, allocatable :: input_varPos(:)

    !> Component index of the input variable in the variable system.
    !! It is used only when there is only one input variable.
    !! Index values must not be zero and > nComponents of input variable
    integer, allocatable :: input_varIndex(:)

    !> Data that is required by the get method.
    type(c_ptr) :: method_data

    !> Operation type
    character(len=labelLen) :: operType

    !> Function to actually obtain the variable at a given point.
    !!
    !! This is either a function accessing a state variable directly,
    !! a function returning a space time function evaluation or a
    !! derived quantity, that computes a new variable out of others.
    procedure(tem_varSys_proc_point), pointer :: get_point => null()

    !> Function to actually obtain the variable in a given element.
    procedure(tem_varSys_proc_element), pointer :: get_element => null()

    !> Function to set parameter in the data_type stored in method_data.
    procedure(tem_varSys_proc_setParams), pointer :: set_params => null()

    !> Function to get parameter in the data_type stored in method_data
    procedure(tem_varSys_proc_getParams), pointer :: get_params => null()

    !> Function to setup points set for boundaries and sources.
    !! Pointe set are stored in method_data level wise 1D growing array for
    !! dimension X,Y and Z.
    !!  * For solver variables, points are stored in solver container.
    !!  * For spacetime variables, points are stored in spacetime function.
    !!  * For operation variables, points are passed down to its
    !!    input_variable.
    procedure(tem_varSys_proc_setupIndices), pointer :: &
      &                                      setup_indices => null()

    !> Function to get value for point set stored in method_data
    !! for requested index in point set.
    !! This function either returns a pre-stored value or compute value
    !! depends on variable type and spacetime function.
    !! For time-independent spacetime function, values are computed
    !! in setupIndices and growing array of points are deleted
    procedure(tem_varSys_proc_getValOfIndex), pointer :: &
      &                                       get_valOfIndex => null()

  end type tem_varSys_op_type
  ! **************************************************************************** !


  !> growing array type for type(tem_varsys_op_type)
  type grw_varoparray_type
    integer :: nvals = 0
    integer :: containersize = 0
    type(tem_varsys_op_type), allocatable :: val(:)
  end type

  !> initialize the dynamic array
  interface init
    module procedure init_ga_varop
  end interface

  !> truncate the array, meaning
  !! cut off the trailing empty entries
  interface truncate
    module procedure truncate_ga_varop
  end interface

  !> empty the entries  without changing arrays
  interface empty
    module procedure empty_ga_varop
  end interface

  !> destroy the dynamic array
  interface destroy
    module procedure destroy_ga_varop
  end interface

  !> insert an element at a given position
  interface placeat
    module procedure placeat_ga_varop
    module procedure placeat_ga_varop_vec
  end interface

  !> append a value to the dynamic array
  !! and return its position.
  interface append
    module procedure append_ga_varop
    module procedure append_ga_varop_vec
  end interface

  !> increase the size of the container
  !! for the array.
  interface expand
    module procedure expand_ga_varop
  end interface



  ! **************************************************************************** !
  !> Description of the variable system.
  type tem_varSys_type
    !> A descriptive name for this system of variables.
    character(len=LabelLen) :: SystemName

    !> Number of variables in the state.
    integer :: nStateVars = 0

    !> Number of scalars in the state.
    !!
    !! This keeps track of the length of the state array.
    integer :: nScalars = 0

    !> Number of auxField variables
    integer :: nAuxVars = 0

    !> Number of scalars in the auxField
    !! This keeps track of the length of the auxField array
    integer :: nAuxScalars

    !> Definition of how to obtain a variable.
    type(grw_varOpArray_type) :: method

    !> List of variables in the system.
    type(dyn_labelArray_type) :: varname

  end type tem_varSys_type
  ! **************************************************************************** !


  ! **************************************************************************** !
  !> A supporting data type to define a solver specific element evaluation for
  !! stfuns.
  type tem_varSys_solverData_evalElem_type
    !> Data from the solver for the evaluation.
    type(c_ptr) :: solver_bundle
    !> Callback function to set the appropriate function for stFun vars.
    procedure(tem_varsys_set_evalelem), pointer :: stFun_setter => null()
    !> Callback function to set the appropriate function for solver specific
    !! operation vars.
    procedure(tem_varsys_set_evalelem), pointer :: opVar_setter => null()
  end type tem_varSys_solverData_evalElem_type
  ! **************************************************************************** !


  ! **************************************************************************** !
  abstract interface
    !> Interface description for a variable access method (single point).
    !!
    !! To obtain values of a given variable, it is necessary to state the
    !! space and time at which the variable should be evaluated.
    !! The interface is vectorized and provides n values at n different space
    !! locations at the same time.
    !! Of course the variable system itself also needs to be passed in, to
    !! allow the computation of other derived quantities as needed.
    !! The method description itself is passed in automatically, and has not
    !! to be provided explicitly.
    subroutine tem_varSys_proc_point(fun, varsys, point, time, tree, nPnts, &
      &                              res)
      import :: tem_varSys_op_type, tem_varSys_type, tem_time_type, rk,        &
        &       treelmesh_type

      !> 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(:)
    end subroutine tem_varSys_proc_point


    !> Interface description for a variable access method (complete element).
    !!
    !! To obtain values of a given variable, it is necessary to state the
    !! treeID and time at which the variable should be evaluated.
    !! The interface is nDofs values to cover the all degrees of freedoms
    !! in the element.
    !! Of course the variable system itself also needs to be passed in, to
    !! allow the computation of other derived quantities as needed.
    !! The method description itself is passed in automatically, and has not
    !! to be provided explicitly.
    subroutine tem_varSys_proc_element(fun, varsys, elempos, time, tree, &
      &                                nElems, nDofs, res)
      import :: tem_varSys_op_type, tem_varSys_type, tem_time_type, rk, &
        &       treelmesh_type

      !> 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 element in tree%treeID to get the variable for.
      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 elements 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:
      !! (nComponents of resulting variable) x (nDegrees of freedom) x (nElems)
      !! Access: (iElem-1)*fun%nComponents*nDofs +
      !!         (iDof-1)*fun%nComponents + iComp
      real(kind=rk), intent(out) :: res(:)
    end subroutine tem_varSys_proc_element


    !> Interface description for a variable to set parameter in data type
    !! stored in method_data.
    !!
    !! To set the parameter, provide a string (recommended in lua format).
    !! Depends on variable procedure, input string will be processed and
    !! information will be store in method_data.
    !! For operation variable, pass the information down to its
    !! input_variable.
    subroutine tem_varSys_proc_setParams(fun, varSys, instring)
      import :: tem_varSys_op_type, tem_varSys_type

      !> 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
    end subroutine tem_varSys_proc_setParams

    !> Interface description for a variable to get parameter in data type
    !! stored in method_data.
    !!
    !! To get the parameter, provide requested parameter via
    !! instring (recommended in lua format) and routine returns the parameter
    !! via outstring in lua format.
    !! For operation variable, pass the information down to its
    !! input_variable.
    subroutine tem_varSys_proc_getParams(fun, varSys, instring, outstring)
      import :: tem_varSys_op_type, tem_varSys_type

      !> 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
    end subroutine tem_varSys_proc_getParams

    !> Interface description for a variable to setup point sets and return
    !! indices of points in the 1D growing array.
    !!
    !! Before setupIndices, use setParams to set information about boundary
    !! or sources in the spacetime function type.
    !! Points are stored in the 1D growing array for each dimension in the
    !! level wise list in the method_data.
    !! offset_bit is required for each point on boundaries so the remote domain
    !! can use this offset to identify rank for the point.
    !! If offset_bit is not present then default is set to element center
    !! i.e offset_bit = achar(1+4+16).
    !!
    !! For state variable, this routine will also compute elemPos in the
    !! given level and local_coord for every point. Elempos and local_coord are
    !! stored in same place as growing array of points.
    !!
    !! For spacetime function variable, if function is time independent then
    !! values are computed and stored level wise after all points are
    !! added to spacetime function.
    subroutine tem_varSys_proc_setupIndices(fun, varSys, point, offset_bit,    &
      &                                     iLevel, tree, nPnts, idx)
      import :: tem_varSys_op_type, tem_varSys_type, rk, treelmesh_type

      !> 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: n
      !!
      !! 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(:)
    end subroutine tem_varSys_proc_setupIndices


    !> Interface description for a variable to return a value at the given
    !! index position in the growing array points set stored in method_data.
    !!
    !! For spacetime function, if value is pre-stored, it will return a value
    !! at given index else if will evaluate a variable at index of a point
    !! for a given time and return a value.
    !!
    !! If index is not present and first is present then index is computed
    !! from first and number of return value (n).
    subroutine tem_varSys_proc_getValOfIndex(fun, varSys, time, iLevel, &
      &                                      idx, idxLen, nVals, res)
      import :: tem_varSys_op_type, tem_varSys_type, tem_time_type, rk

      !> Description of the method to obtain the variables, here some preset
      !! values might be stored, like the space time function to use or the
      !! required variables.
      class(tem_varSys_op_type), intent(in) :: fun

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

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

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

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

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

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

      !> Resulting values for the requested variable.
      !!
      !! Dimension: n requested entries x nComponents of this variable
      !! Access: (iElem-1)*fun%nComponents + iComp
      real(kind=rk), intent(out) :: res(:)
    end subroutine tem_varSys_proc_getValOfIndex


    subroutine tem_varsys_set_evalelem(set_elem_eval, fun)
      import :: tem_varSys_op_type, tem_varSys_solverData_evalElem_type
      !> Description on how to set the element retrieval function for stfuns.
      !! and solver specific operation variables
      class(tem_varSys_solverData_evalElem_type), intent(in) :: set_elem_eval

      !> 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.
      type(tem_varSys_op_type), intent(inout) :: fun
    end subroutine tem_varsys_set_evalelem

  end interface
  ! **************************************************************************** !


  ! **************************************************************************** !
  interface tem_varSys_load
    module procedure tem_varSys_load_vector
    module procedure tem_varSys_load_single
  end interface tem_varSys_load
  ! **************************************************************************** !

  ! **************************************************************************** !
  interface tem_varSys_dump
    module procedure tem_varSys_dump_vector
    module procedure tem_varSys_dump_single
  end interface tem_varSys_dump
  ! **************************************************************************** !

  ! **************************************************************************** !
  interface tem_varSys_out
    module procedure tem_varSys_out_vector
    module procedure tem_varSys_out_single
  end interface tem_varSys_out
  ! **************************************************************************** !

  interface assignment(=)
    module procedure copy_varOp
  end interface


contains


  subroutine init_ga_varop(me, length)
    type(grw_varoparray_type), intent(out) :: me !< dynamic array to init
    integer, intent(in), optional :: length !< initial length of the container

    if (present(length)) then
      me%containersize = length
    else
      me%containersize = zerolength
    end if
    ! deallocate ...
    if( allocated( me%val ))     &
      deallocate(me%val)
    ! ... and reallocate
    allocate(me%val(me%containersize))
    me%nvals = 0

  end subroutine init_ga_varop

  subroutine destroy_ga_varop(me)
    type(grw_varoparray_type), intent(inout) :: me !< dynamic array to destroy

    me%containersize = 0
    me%nvals = 0
    if( allocated( me%val ) ) deallocate(me%val)

  end subroutine destroy_ga_varop


  subroutine truncate_ga_varop(me)
    !------------------------------------------------------------------------
    type(grw_varoparray_type) :: me !< array to truncate
    !------------------------------------------------------------------------
    type(tem_varsys_op_type), allocatable :: tarray(:)
    !------------------------------------------------------------------------
    integer :: ii
    !------------------------------------------------------------------------

    ! nothing to do if container size is not larger than the number of values
    ! in the array.
    if (me%containersize > me%nvals) then
      allocate(tarray(me%nvals))
      do ii = 1, me%nvals
        tarray(ii) = me%val(ii)
      end do
      call move_alloc(tarray, me%val)
      me%containersize = me%nvals
    end if

  end subroutine truncate_ga_varop


  subroutine empty_ga_varop(me)
    !------------------------------------------------------------------------
    type(grw_varoparray_type) :: me !< array to sorttruncate
    !------------------------------------------------------------------------

    me%nvals = 0

  end subroutine empty_ga_varop

  !> adds the value to a given position inside the growing array.
  !!
  !! if the requested position is outside the current array bounds, the array
  !! will be resized accordingly. if it is inside the current array bounds, the
  !! element at the requested position will be replaced.
  subroutine placeat_ga_varop(me, val, pos, length)
    type(grw_varoparray_type) :: me !< array to place the value into
    type(tem_varsys_op_type), intent(in) :: val !< value to place at the given position
    integer, intent(in) :: pos !< predefined position
    !> optional length to expand the array
    integer, intent(in), optional :: length


    ! value to append is larger than all existing ones,
    ! just put it to the end of the list, this captures
    ! also the case of empty lists.
    ! in this case foundpos = me%nvals + 1 holds.
    if (pos > me%containersize) then
      ! expand the array, if its boundary is reached
      call expand(me = me, pos = pos, length = length)
    end if

    me%nvals = max( pos, me%nvals )
    me%val(pos) = val

  end subroutine placeat_ga_varop


  !> adds the values starting from a given position inside the growing array.
  !!
  !! if the requested position is outside the current array bounds, the array
  !! will be resized accordingly. if it is inside the current array bounds, the
  !! elements starting from the requested position will be replaced up to
  !! the element at position `pos + size(val) - 1`.
  subroutine placeat_ga_varop_vec(me, val, pos, length)
    type(grw_varoparray_type) :: me !< array to append the value to
    type(tem_varsys_op_type), intent(in) :: val(:) !< values to append
    integer, intent(in) :: pos !< predefined position
    !> optional length to expand the array
    integer, intent(in), optional :: length

    integer :: ub, ii

    if (me%nvals == huge(me%nvals)) then
      write(*,*) "reached end of integer range for growing array!"
      write(*,*) "aborting!!"
      stop
    end if

    ub = pos + size(val) - 1

    if (ub > me%containersize) then
      ! expand the array, if its boundary is reached
      call expand(me = me, pos = ub, length = length)
    end if

    me%nvals = max( ub, me%nvals )
    do ii = pos, ub
      me%val(ii) = val(1+ii-pos)
    end do

  end subroutine placeat_ga_varop_vec


  subroutine append_ga_varop(me, val, length)
    type(grw_varoparray_type) :: me !< array to append the value to
    type(tem_varsys_op_type), intent(in) :: val !< value to append
    !> optional length to expand the array
    integer, intent(in), optional :: length

    ! value to append is larger than all existing ones,
    ! just put it to the end of the list, this captures
    ! also the case of empty lists.
    ! in this case foundpos = me%nvals + 1 holds.
    if (me%nvals+1 > me%containersize) then
      ! expand the array, if its boundary is reached
      call expand(me = me, length = length)
    end if

    me%nvals = me%nvals+1
    me%val(me%nvals) = val

  end subroutine append_ga_varop

  subroutine append_ga_varop_vec(me, val, length)
    type(grw_varoparray_type) :: me !< array to append the value to
    type(tem_varsys_op_type), intent(in) :: val(:) !< values to append
    !> optional length to expand the array
    integer, intent(in), optional :: length

    integer :: lb, ub, ii

    if (me%nvals == huge(me%nvals)) then
      write(*,*) "reached end of integer range for growing array!"
      write(*,*) "aborting!!"
      stop
    end if

    lb = me%nvals + 1
    ub = lb + size(val) - 1

    if (ub > me%containersize) then
      ! expand the array, if its boundary is reached
      call expand(me = me, pos = ub, length = length)
    end if

    me%nvals = max( ub, me%nvals )
    do ii = lb, ub
      me%val(ii) = val(1+ii-lb)
    end do

  end subroutine append_ga_varop_vec


  subroutine expand_ga_varop(me, pos, length)
    type(grw_varoparray_type) :: me !< array to resize
    integer, intent(in), optional :: pos !< optional predefined position
    !> optional length to expand the array
    integer, intent(in), optional :: length

    type(tem_varsys_op_type), allocatable :: swpval(:)
    integer :: explen, ii

    explen = 0
    ! increase the container by the requested length of double it
    if( present(length) ) then
      explen = max( length, minlength )
    else
      ! set the global minimum length, if doubling would be smaller than that
      explen = max(me%containersize, minlength)
    end if

    ! if a position is given, increase the container to at least the size to
    ! fit the position.
    if( present(pos) ) explen = max(explen, pos-me%containersize)

    ! if the current size plus explen exceeds the max container size,
    ! reduce the size to the max container size.
    if( (huge(me%containersize) - explen) <= me%containersize) then
      ! set max container size
      me%containersize = huge(me%containersize)
    else
      ! set the new container size
      me%containersize = me%containersize + explen
    end if

    if ( me%nvals > 0 ) then
      allocate(swpval(me%containersize))
      do ii = 1, me%nvals
        swpval(ii) = me%val(ii)
      end do
      call move_alloc( swpval, me%val )
    else ! me%nvals == 0
      if ( allocated(me%val) ) deallocate( me%val )
      allocate( me%val(me%containersize) )
    end if

  end subroutine expand_ga_varop


  ! **************************************************************************** !
  !> Initialize a variable system.
  subroutine tem_varSys_init(me, systemName, length)
    ! --------------------------------------------------------------------------!
    !> Variable system to initialize.
    type(tem_varSys_type),           intent(out) :: me

    !> A descriptive name for this system of variables.
    character(len=*),                intent(in)  :: systemName

    !> Initial length for the arrays used in the variable system.
    integer,               optional, intent(in)  :: length
    ! --------------------------------------------------------------------------!

    me%systemName = trim(systemName)
    me%nStateVars = 0
    me%nScalars = 0
    me%nAuxVars = 0
    me%nAuxScalars = 0

    call init( me     = me%varname, &
      &        length = length      )

    call init( me     = me%method, &
      &        length = length     )

  end subroutine tem_varSys_init
  ! **************************************************************************** !


  ! **************************************************************************** !
  !> Append a new state variable to the variable system.
  !!
  !! Note that the actual data on how to access the state has to be encoded
  !! in the method_data pointer.
  !! The callers therefore have to ensure the consistency of those access
  !! definitions.
  !! nComponents has to be at least 1.
  subroutine tem_varSys_append_stateVar(me, varname, nComponents, method_data, &
    &                                   get_point, get_element, set_params,    &
    &                                   get_params, setup_indices,             &
    &                                   get_valOfIndex, pos, wasAdded)
    ! --------------------------------------------------------------------------!
    !> Variable system to append the state variable to.
    type(tem_varSys_type),             intent(inout) :: me

    !> Variable to append to the state.
    character(len=*),                  intent(in)    :: varname

    !> Number of components in this variable.
    integer, intent(in) :: nComponents

    !> Data that is required by the methods to obtain the variable.
    type(c_ptr), intent(in) :: method_data

    !> Procedure which allows the retrieval of the variable at given points.
    procedure(tem_varSys_proc_point), pointer :: get_point

    !> Procedure which allows the retrieval of the variable in an element.
    procedure(tem_varSys_proc_element), pointer :: get_element

    !> Procedure which allows to set parameter in method_data
    procedure(tem_varSys_proc_setParams), optional, pointer :: set_params

    !> Procedure which allows to get parameter in method_data
    procedure(tem_varSys_proc_getParams), optional, pointer :: get_params

    !> Procedure to setup growing array of points, variable value
    !! in method_data and return index of points set
    procedure(tem_varSys_proc_setupIndices), pointer :: setup_indices

    !> Procedure which allows to retrieval of the variable at point or val
    !! array index
    procedure(tem_varSys_proc_getValOfIndex), pointer :: get_valOfIndex

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

    !> Indicator, if the variable was actually added to the system.
    logical,                 optional, intent(out)   :: wasAdded
    ! --------------------------------------------------------------------------!
    logical :: newVar
    type(tem_varSys_op_type) :: state_access
    integer :: iComp, varPos
    ! --------------------------------------------------------------------------!

    call append( me       = me%varname, &
      &          val      = varname,    &
      &          pos      = varPos,     &
      &          wasAdded = newVar      )

    if (newVar) then
      state_access%mypos = varPos
      state_access%operType = 'state'
      state_access%nInputs = 0
      state_access%nComponents = nComponents
      allocate(state_access%state_varPos(nComponents))
      do iComp = 1, nComponents
        state_access%state_varPos(iComp) = me%nScalars + iComp
      end do
      allocate(state_access%input_varPos(0))

      state_access%method_data = method_data
      state_access%get_point => get_point
      state_access%get_element => get_element
      if (present(set_params)) then
        state_access%set_params => set_params
      else
        state_access%set_params => tem_varSys_setParams_dummy
      end if
      if (present(get_params)) then
        state_access%get_params => get_params
      else
        state_access%get_params => tem_varSys_getParams_dummy
      end if
      state_access%setup_indices => setup_indices
      state_access%get_valOfIndex => get_valOfIndex
      call append( me  = me%method,   &
        &          val = state_access )
      me%nStateVars = me%nStateVars+1
      me%nScalars = me%nScalars + nComponents
    end if

    if (present(pos)) pos = varPos
    if (present(wasAdded)) wasAdded = newVar

  end subroutine tem_varSys_append_stateVar
  ! **************************************************************************** !


  ! **************************************************************************** !
  !> Append a new auxiliaryField variable to the variable system.
  !!
  !! Note that the actual data on how to access the state has to be encoded
  !! in the method_data pointer.
  !! The callers therefore have to ensure the consistency of those access
  !! definitions.
  !! nComponents has to be at least 1.
  subroutine tem_varSys_append_auxFieldVar(me, varname, nComponents, &
    & method_data, get_point, get_element, set_params, get_params,   &
    & setup_indices, get_valOfIndex, pos, wasAdded)
    ! --------------------------------------------------------------------------!
    !> Variable system to append the state variable to.
    type(tem_varSys_type),             intent(inout) :: me

    !> Variable to append to the state.
    character(len=*),                  intent(in)    :: varname

    !> Number of components in this variable.
    integer, intent(in) :: nComponents

    !> Data that is required by the methods to obtain the variable.
    type(c_ptr), intent(in) :: method_data

    !> Procedure which allows the retrieval of the variable at given points.
    procedure(tem_varSys_proc_point), pointer :: get_point

    !> Procedure which allows the retrieval of the variable in an element.
    procedure(tem_varSys_proc_element), pointer :: get_element

    !> Procedure which allows to set parameter in method_data
    procedure(tem_varSys_proc_setParams), optional, pointer :: set_params

    !> Procedure which allows to get parameter in method_data
    procedure(tem_varSys_proc_getParams), optional, pointer :: get_params

    !> Procedure to setup growing array of points, variable value
    !! in method_data and return index of points set
    procedure(tem_varSys_proc_setupIndices), pointer :: setup_indices

    !> Procedure which allows to retrieval of the variable at point or val
    !! array index
    procedure(tem_varSys_proc_getValOfIndex), pointer :: get_valOfIndex

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

    !> Indicator, if the variable was actually added to the system.
    logical,                 optional, intent(out)   :: wasAdded
    ! --------------------------------------------------------------------------!
    logical :: newVar
    type(tem_varSys_op_type) :: auxField_access
    integer :: iComp, varPos
    ! --------------------------------------------------------------------------!

    call append( me       = me%varname, &
      &          val      = varname,    &
      &          pos      = varPos,     &
      &          wasAdded = newVar      )

    if (newVar) then
      auxField_access%mypos = varPos
      auxField_access%operType = 'auxfield'
      auxField_access%nInputs = 0
      auxField_access%nComponents = nComponents
      allocate(auxField_access%auxField_varPos(nComponents))
      do iComp = 1, nComponents
        auxField_access%auxField_varPos(iComp) = me%nAuxScalars + iComp
      end do
      allocate(auxField_access%input_varPos(0))

      auxField_access%method_data = method_data
      auxField_access%get_point => get_point
      auxField_access%get_element => get_element
      if (present(set_params)) then
        auxField_access%set_params => set_params
      else
        auxField_access%set_params => tem_varSys_setParams_dummy
      end if
      if (present(get_params)) then
        auxField_access%get_params => get_params
      else
        auxField_access%get_params => tem_varSys_getParams_dummy
      end if
      auxField_access%setup_indices => setup_indices
      auxField_access%get_valOfIndex => get_valOfIndex
      call append( me  = me%method,   &
        &          val = auxField_access )
      me%nAuxVars = me%nAuxVars+1
      me%nAuxScalars = me%nAuxScalars + nComponents
    end if

    if (present(pos)) pos = varPos
    if (present(wasAdded)) wasAdded = newVar

  end subroutine tem_varSys_append_auxFieldVar
  ! **************************************************************************** !

  ! **************************************************************************** !
  !> Append a new (non-state) variable to the variable system.
  !!
  !! Non-state variables might depend on other variables, and do not modify the
  !! internal state variable counter.
  subroutine tem_varSys_append_derVar(me, varname, operType, nComponents,    &
    &                                 input_varname, input_varindex,         &
    &                                 method_data, get_point, get_element,   &
    &                                 set_params, get_params, setup_indices, &
    &                                 get_valOfIndex, pos, wasAdded          )
    ! --------------------------------------------------------------------------!
    !> Variable system to append the state variable to.
    type(tem_varSys_type), intent(inout) :: me

    !> Variable to append to the state.
    character(len=*), intent(in) :: varname

    !> Operation type
    character(len=*), optional, intent(in) :: operType

    !> Number of components in this variable.
    integer, intent(in) :: nComponents

    !> List of variable names, this variable depends on.
    !!
    !! (Might be empty). The variable will only be appended, if all invars are
    !! already found in the variable system.
    character(len=*), optional, intent(in) :: input_varname(:)

    !> Component index to access from input variable nComponents
    !!
    !! (Might be empty). The variable will only be appended, if all index are
    !! available in input_varname
    integer, optional, intent(in) :: input_varIndex(:)

    !> Data that is required by the methods to obtain the variable.
    type(c_ptr), intent(in) :: method_data

    !> Procedure which allows the retrieval of the variable at given points.
    procedure(tem_varSys_proc_point), pointer :: get_point

    !> Procedure which allows the retrieval of the variable in an element.
    procedure(tem_varSys_proc_element), pointer :: get_element

    !> Procedure which allows to set parameter in method_data
    procedure(tem_varSys_proc_setParams), optional, pointer :: set_params

    !> Procedure which allows to get parameter in method_data
    procedure(tem_varSys_proc_getParams), optional, pointer :: get_params

    !> Procedure to setup growing array of points, variable value
    !! in method_data and return index of points set
    procedure(tem_varSys_proc_setupIndices), pointer :: setup_indices

    !> Procedure which allows to retrieval of the variable at point or val
    !! array index
    procedure(tem_varSys_proc_getValOfIndex), pointer :: get_valOfIndex

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

    !> Indicator, if the variable was actually added to the system.
    logical,  optional, intent(out)   :: wasAdded
    ! --------------------------------------------------------------------------!
    logical :: newVar, nonexisting_input
    type(tem_varSys_op_type) :: var_access
    integer :: iIn, nInputs, varPos, in_nComp, nVarIndex
    integer, allocatable :: inpos(:), input_varIndex_loc(:)
    ! --------------------------------------------------------------------------!

    if (present(input_varname)) then
      nInputs = size(input_varname)
    else
      nInputs = 0
    end if

    if (present(input_varindex)) then
      nVarIndex = size(input_varindex)
    else
      nVarIndex = 0
    end if

    allocate(inpos(nInputs))
    allocate(input_varIndex_loc(nVarIndex))

    nonexisting_input = .false.

    ! Check whether all input variables are present. If not, the variable can't
    ! be added, because variables it depends on are not present.
    do iIn=1,nInputs
      inpos(iIn) = PositionOfVal(me%varname, input_varname(iIn))
      nonexisting_input = (inpos(iIn) <= 0)
      if (nonexisting_input) then
        write(logUnit(llerror),*) 'Variable ' // trim(varname) &
          & // ' cant be added due to missing input variable '   &
          & // trim(input_varname(iIn))
        EXIT
      end if

      ! allocate and set input_varIndex only if nVarIndex > 0
      ! input_varIndex is used only for one nInputs
      if ( nVarIndex > 0 .and. nInputs == 1 ) then
        ! Check whether the component index requested from input_varname
        ! is available
        in_nComp = me%method%val(inpos(iIn))%nComponents
        if ( (maxval(input_varIndex) > in_nComp) .or. &
          &  (nVarIndex > in_nComp) ) then
          write(logUnit(llerror),*) 'WARNING: Component index requested ' &
            & // 'is not extractable from variable: '  &
            & // trim(input_varname(iIn))
          write(logUnit(llerror),*) ' Input variable name: ' &
            & // trim(input_varname(iIn))
          write(logUnit(llerror),*) ' nComponents of input_varname: ', &
            & in_nComp
          write(logUnit(llerror),*) ' Requested component index of ' &
            & // 'input var: ', input_varIndex
          EXIT
        end if
        ! Copy input_varIndex
        input_varIndex_loc = input_varIndex
      end if

    end do

    if (nonexisting_input) then

      varPos = 0
      newVar = .false.

    else

      ! All inputs are found in the variable system, thus this derived
      ! quantity can be added now.
      call append( me       = me%varname, &
        &          val      = varname,    &
        &          pos      = varPos,     &
        &          wasAdded = newVar      )

      if (newVar) then
        write(logUnit(lldebug),*) 'Variable ' // trim(varname) &
          & // ' was newly added'
        var_access%mypos = varPos
        if (present(operType)) then
          var_access%operType = trim(operType)
        else
          var_access%operType = 'derived'
        end if
        var_access%nInputs = nInputs
        var_access%nComponents = nComponents
        allocate(var_access%state_varPos(0))
        allocate(var_access%auxField_varPos(0))
        allocate(var_access%input_varPos(nInputs))
        var_access%input_varPos = inpos
        allocate(var_access%input_varIndex(nVarIndex))
        var_access%input_varIndex = input_varIndex_loc

        var_access%method_data = method_data
        var_access%get_point => get_point
        var_access%get_element => get_element
        if (present(set_params)) then
          var_access%set_params => set_params
        else
          var_access%set_params => tem_varSys_setParams_dummy
        end if
        if (present(get_params)) then
          var_access%get_params => get_params
        else
          var_access%get_params => tem_varSys_getParams_dummy
        end if
        var_access%setup_indices => setup_indices
        var_access%get_valOfIndex => get_valOfIndex
        call append( me  = me%method,   &
          &          val = var_access )
      else
        write(logUnit(lldebug),*) 'Variable ' // trim(varname) &
          & // ' is already present'
      end if

    end if

    if (present(pos)) pos = varPos
    if (present(wasAdded)) wasAdded = newVar

  end subroutine tem_varSys_append_derVar
  ! **************************************************************************** !


  ! **************************************************************************** !
  !> Load variable system(s) from a lua file
  !!
  subroutine tem_varSys_load_vector( me, conf, parent, key )
    ! ---------------------------------------------------------------------------
    !> The variable system to read in
    type(tem_varSys_type), intent(out), allocatable :: me(:)
    !> Lua handle connected to the script to read the table from
    type(flu_State) :: conf
    !> A parent table handle in which to look the current variables up
    integer, intent(in), optional :: parent
    !> load from this key
    character(len=*), optional, intent(in) :: key
    ! ---------------------------------------------------------------------------
    integer :: thandle, subhandle
    integer :: nSys, iSys
    character(len=LabelLen) :: local_key
    ! ---------------------------------------------------------------------------

    if( present( key )) then
      local_key = key
    else
      local_key = 'varsys'
    endif

    call aot_table_open( L       = conf,           &
      &                  parent  = parent,         &
      &                  thandle = thandle,        &
      &                  key     = trim(local_key) )
    !varsys table is present, load variable system
    if( thandle > 0 ) then
      ! Now the 'varsys' table is opened. Let's first check, if
      ! there is a single variable system directly inside, or if there are
      ! several, so we have to open a table first and check whether it has
      ! another table
      call aot_table_open( L       = conf,      &
        &                  parent  = thandle,   &
        &                  thandle = subhandle, &
        &                  pos     = 1          )
      ! more than one varSys
      if ( subhandle > 0 ) then
        call aot_table_close( L = conf, thandle = subhandle )
        nSys = aot_table_length( L = conf, thandle = thandle )
        allocate(me(nSys))
        ! load each varSys
        do iSys = 1, nSys
          call aot_table_open( L       = conf,      &
            &                  parent  = thandle,   &
            &                  thandle = subhandle, &
            &                  pos     = iSys       )
          call tem_varSys_load_single( me        = me(iSys),  &
            &                          conf      = conf,      &
            &                          parent    = subhandle, &
            &                          openTable = .false.    )
          call aot_table_close( L = conf, thandle = subhandle )
        end do
      else !single varSys
        nSys = 1
        allocate(me(nSys))
        call tem_varSys_load_single( me        = me(1),   &
          &                          conf      = conf,    &
          &                          parent    = thandle, &
          &                          openTable = .false.  )
      end if
    else
      write(logUnit(1),*) 'Error: Failed to load varsys table'
      call tem_abort()
    end if

    call aot_table_close( L = conf, thandle = thandle )

  end subroutine tem_varSys_load_vector
  ! **************************************************************************** !


  ! **************************************************************************** !
  !> load varSys from lua file.
  !! Required for harvester to load varSys from tracking or restart header
  !! file.
  subroutine tem_varSys_load_single(me, conf, parent, key, openTable)
    ! --------------------------------------------------------------------------!
    !> varSys to read from the Lua script(conf) and fill
    type(tem_varSys_type), intent(out) :: me

    !> Lua handle connected to the script to read the table from
    type(flu_state) :: conf

    !> A parent table handle in which to look the current variable up
    integer, optional, intent(in) :: parent

    !> load varsys from this key
    character(len=*), optional, intent(in) :: key

    ! if varsys table is already opened then openTable = .false.
    logical, optional, intent(in) :: openTable
    ! --------------------------------------------------------------------------!
    integer :: syshandle, varhandle, varsubhandle, iVar, nVars, pos
    type(tem_varSys_op_type) :: method
    character(len=labelLen) :: varname
    integer :: iError
    integer, allocatable :: vError(:)
    character(len=LabelLen) :: local_key
    logical :: openTable_loc
    ! --------------------------------------------------------------------------!
    ! if varsys table is already opened then openTable = .false.
    if( present(openTable) ) then
      openTable_loc = openTable
    else
      openTable_loc = .true.
    end if

    if( present( key )) then
      local_key = key
    else
      local_key = 'varsys'
    endif

    if (openTable_loc) then
      write(logUnit(1),*) 'Loading varsys table '
      ! Try to open the varsys table
      call aot_table_open( L       = conf,           &
        &                  parent  = parent,         &
        &                  thandle = syshandle,      &
        &                  key     = trim(local_key) )
      write(logUnit(1),*) 'Syshandle ', syshandle
    else
      ! if table is open before just use the parent as syshandle to load
      ! variable info
      syshandle = parent
    end if

    !varsys table is present, load variable system
    if (syshandle > 0) then
      ! Get the variable system name
      call aot_get_val( L       = conf,          &
        &               thandle = syshandle,     &
        &               val     = me%systemname, &
        &               ErrCode = iError,        &
        &               key     = 'systemname'   )

      ! get nStateVars
      call aot_get_val( L       = conf,          &
        &               thandle = syshandle,     &
        &               val     = me%nStateVars, &
        &               ErrCode = iError,        &
        &               key     = 'nStateVars'   )

      call aot_get_val( L       = conf,        &
        &               thandle = syshandle,   &
        &               val     = me%nScalars, &
        &               ErrCode = iError,      &
        &               key     = 'nScalars'   )

      call aot_table_open( L       = conf,      &
        &                  parent  = syshandle, &
        &                  thandle = varhandle, &
        &                  key     = 'variable' )

      if (varhandle > 0) then
        nVars = aot_table_length( L = conf, thandle = varhandle )
        call init( me     = me%varname, &
          &        length = nVars       )

        call init( me     = me%method, &
          &        length = nVars      )

        do iVar = 1, nVars
          method%myPos = iVar
          call aot_table_open( L       = conf,         &
            &                  parent  = varhandle,    &
            &                  thandle = varsubhandle, &
            &                  pos     = iVar          )

          call aot_get_val( L       = conf,         &
            &               thandle = varsubhandle, &
            &               val     = varname,      &
            &               ErrCode = iError,       &
            &               key     = 'name'        )

          call append( me       = me%varname, &
            &          val      = varname,    &
            &          pos      = pos         )

          call aot_get_val( L       = conf,               &
            &               thandle = varsubhandle,       &
            &               val     = method%nComponents, &
            &               ErrCode = iError,             &
            &               key     = 'ncomponents'       )

          call aot_get_val( L         = conf,                &
            &               thandle   = varsubhandle,        &
            &               val       = method%state_varPos, &
            &               maxLength = method%nComponents,  &
            &               ErrCode   = vError,              &
            &               key       = 'state_varpos'       )

          call append( me = me%method, val = method )
          call aot_table_close( L = conf, thandle = varsubhandle)
        end do
        call aot_table_close( L = conf, thandle = varhandle)
      else
        write(logUnit(1),*) 'Error: Failed to load variable table within ' &
          &                 //'varsys table'
        call tem_abort()
      end if ! variable table
    else
      write(logUnit(1),*) 'Error: Failed to load varsys table single'
      call tem_abort()
    end if !varsys table
    call aot_table_close( L = conf, thandle = syshandle)

  end subroutine tem_varSys_load_single
  ! **************************************************************************** !


  ! **************************************************************************** !
  !> Dumps array of varSys to given unit
  subroutine tem_varSys_dump_vector(me, outUnit, dumpVarPos)
    ! ---------------------------------------------------------------------------
    !> variable to write into the lua file
    type(tem_varSys_type), intent(in) :: me(:)

    !> unit to write to
    integer, intent(inout) :: outUnit

    !> Position of variables to dump
    integer, optional, intent(in) :: dumpVarPos(:)
    ! ---------------------------------------------------------------------------
    ! aotus type handling the output to the file in lua format
    type(aot_out_type) :: conf
    ! ---------------------------------------------------------------------------
    call aot_out_open( put_conf = conf, outUnit = outUnit )
    call tem_varSys_out_vector( me, conf, dumpVarPos )
    call aot_out_close( put_conf = conf )

  end subroutine tem_varSys_dump_vector
  ! **************************************************************************** !


  ! **************************************************************************** !
  !> Dump single varSys to given unit
  subroutine tem_varSys_dump_single(me, outUnit, dumpVarPos)
    ! ---------------------------------------------------------------------------
    !> variable to write into the lua file
    type(tem_varSys_type), intent(in) :: me

    !> unit to write to
    integer, intent(inout) :: outUnit

    !> Position of variables to dump
    integer, optional, intent(in) :: dumpVarPos(:)
    ! ---------------------------------------------------------------------------
    ! aotus type handling the output to the file in lua format
    type(aot_out_type) :: conf
    ! ---------------------------------------------------------------------------
    call aot_out_open( put_conf = conf, outUnit = outUnit )
    call tem_varSys_out_single( me, conf, dumpVarPos )
    call aot_out_close( put_conf = conf )

  end subroutine tem_varSys_dump_single
  ! **************************************************************************** !


  ! **************************************************************************** !
  !> Allows the output of array of varSys to lua out
  subroutine tem_varSys_out_vector(me, conf, dumpVarPos)
    ! ---------------------------------------------------------------------------
    !> variable to write into the lua file
    type(tem_varSys_type), intent(in) :: me(:)

    !> aotus type handling the output to the file in lua format
    type(aot_out_type), intent(inout) :: conf

    !> Position of variables to dump
    integer, optional, intent(in) :: dumpVarPos(:)
    ! ---------------------------------------------------------------------------
    integer :: iSys
    ! ---------------------------------------------------------------------------
    call aot_out_open_table( put_conf = conf, tname='varsys' )
    do iSys = 1, size(me)
      call tem_varSys_out_single( me         = me(iSys),   &
        &                         conf       = conf,       &
        &                         dumpVarPos = dumpVarPos, &
        &                         level      = 1           )
    end do
    call aot_out_close_table( put_conf = conf )

  end subroutine tem_varSys_out_vector
  ! **************************************************************************** !


  ! **************************************************************************** !
  !> Write the system of variables description into a Lua file.
  !! use the aotus_out functions for doing so, in order to obtain a neatly
  !! formatted lua file
  !!
  subroutine tem_varSys_out_single(me, conf, dumpVarPos, level)
    ! ---------------------------------------------------------------------------
    !> Variable system to write out
    type(tem_varSys_type), intent(in) :: me

    !> Aotus type handling the output to the file in lua format
    type(aot_out_type), intent(inout) :: conf

    !> Position of variables to dump
    integer, optional, intent(in) :: dumpVarPos(:)

    !> to dump varSys with key or without key
    integer, optional, intent(in) :: level
    ! ---------------------------------------------------------------------------
    integer :: iVar
    ! number of variables to be dumped (nVals or nVars)
    integer :: nVarDump, pos, level_loc, nScalars, nStateVars
    integer, allocatable :: dumpVarPos_loc(:)
    ! --------------------------------------------------------------------------
    if (present(level)) then
      level_loc = level
    else
      level_loc = 0
    end if

    ! dumpVarPos is present, then dump only variable at this position in varSys
    if (present(dumpVarPos)) then
      nVarDump = size(dumpVarPos)
      allocate(dumpVarPos_loc(nVarDump))
      dumpVarPos_loc = dumpVarPos
      nStateVars = nVarDump
      nScalars = sum(me%method%val(dumpVarPos_loc(:))%nComponents)
    else
      nVarDump = me%varname%nVals
      allocate(dumpVarPos_loc(nVarDump))
      dumpVarPos_loc = me%method%val(1:nVarDump)%myPos
      nStateVars = me%nStateVars
      nScalars = me%nScalars
    end if


    if( level_loc == 0) then
      call aot_out_open_table( put_conf = conf, tname = 'varsys' )
    else
      call aot_out_open_table( put_conf = conf )
    end if

    call aot_out_val( put_conf = conf, val = trim(me%SystemName),              &
      &               vname = 'systemname')

    call aot_out_open_table( put_conf = conf, tname = 'variable' )
    do iVar = 1, nVarDump
      call aot_out_open_table( put_conf = conf )

      pos = dumpVarPos_loc(iVar)

      call aot_out_val( put_conf = conf,                                       &
        &               val      = trim(me%varname%val(pos)),                  &
        &               vname    = 'name' )
      call aot_out_val( put_conf = conf,                                       &
        &               val      = me%method%val(pos)%nComponents,             &
        &               vname    = 'ncomponents' )

      !call aot_out_val( put_conf = conf,                                       &
      !  &               val      = me%method%val(pos)%myPos,                   &
      !  &               vname    = 'mypos' )

      if (allocated(me%method%val(pos)%state_varPos)) then
        if (size(me%method%val(pos)%state_varPos)>0) then
          call aot_out_val( put_conf = conf,                                   &
            &               val      = me%method%val(pos)%state_varPos,        &
            &               vname    = 'state_varpos' )
        end if
      end if

      !if (me%method%val(pos)%nInputs > 0) then
      !  call aot_out_val( put_conf = conf,                                    &
      !    &               val      = me%method%val(pos)%input_varPos,         &
      !    &               vname    = 'input_varPos' )
      !end if

      call aot_out_close_table( put_conf = conf )
    end do

    call aot_out_close_table( put_conf = conf )
    call aot_out_val( put_conf = conf, val = nScalars, vname = 'nScalars' )
    call aot_out_val( put_conf = conf, val = nStateVars, vname = 'nStateVars' )
    call aot_out_val( put_conf = conf, val = me%nAuxScalars, &
      & vname = 'nAuxScalars' )
    call aot_out_val( put_conf = conf, val = me%nAuxVars, vname = 'nAuxVars' )
    call aot_out_close_table( put_conf = conf )

  end subroutine tem_varSys_out_single
  ! **************************************************************************** !


  ! **************************************************************************** !
  !> A routine to evaluate chunk of elements for given list of variables
  !!
  !! If subTree is present, it will use map2Global from subTree else
  !! map2Global is created for current chunk of global mesh
  subroutine tem_get_element_chunk(varSys, varPos, elemPos, time, tree, nElems,&
    &                              nDofs, res)
    ! --------------------------------------------------------------------------!
    !> Variable system describing available data.
    type(tem_varsys_type), intent(in) :: varsys

    !> Position of variables to evaluate in varSys
    integer, intent(in) :: varPos(:)

    !> Mesh definition of the input data.
    type(treelmesh_type), intent(in) :: tree

    !> Time information for the current data.
    type(tem_time_type), intent(in) :: time

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

    !> Position of treeID of the element to get the variable for.
    integer, intent(in) :: elemPos(:)

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

    !> Output data
    !! size: io_buffer_size
    real(kind=rk), intent(out) :: res(:)
    ! --------------------------------------------------------------------------!
    integer :: maxComponents, nScalars, nComponents, elemsize, compOff
    integer :: iElem, iVar, iDof, nVars, res_size
    integer :: e_start, t_start, d_start
    real(kind=rk), allocatable :: tmpdat(:)
    ! --------------------------------------------------------------------------!
    ! number of variables to evaluate
    nVars = size(varPos)

    ! Number of scalars in current output
    nScalars = sum(varSys%method%val(varPos(:))%nComponents)

    ! Size of a single element
    elemsize = nScalars*nDofs

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

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

    compOff = 0
    vars: do iVar=1, nVars

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

      ! get the size of the needed part of the res array
      res_size = nElems * nDofs * nComponents

      ! derive the quantities for all the elements in the current chunk
      call varSys%method%val(varpos(iVar))%get_element(                     &
        &                                varSys  = varSys,                  &
        &                                elemPos = elemPos,                 &
        &                                time    = time,                    &
        &                                tree    = tree,                    &
        &                                nElems  = nElems,                  &
        &                                nDofs   = nDofs,                   &
        &                                res     = tmpdat(:res_size)        )

      ! copy the information to the right positions in the result array
      ! res contains results for all variables,
      ! tmpdat is only for one variable
      do iElem=1,nElems
        e_start = (iElem-1)*elemsize
        t_start = (iElem-1)*nComponents*nDofs
        do iDof=1,nDofs
          d_start = (iDof-1)*nScalars + compOff
          res( (e_start+d_start+1) : (e_start+d_start+nComponents) ) &
            &  = tmpdat( t_start + (iDof-1)*nComponents + 1 &
            &            :t_start + iDof*nComponents        )
        end do
      end do
      ! Increase the component offset for the next variables.
      compOff = compOff + nComponents
    end do vars

  end subroutine tem_get_element_chunk
  ! **************************************************************************** !


  ! **************************************************************************** !
  !> A routine to evaluate chunk of points for given list of variables
  !!
  !! If subTree is present, it will use map2Global from subTree else
  !! map2Global is created for current chunk of global mesh
  subroutine tem_get_point_chunk(varSys, varPos, point, time, tree, nPnts, &
    &                            res)
    ! --------------------------------------------------------------------------!
    !> Variable system describing available data.
    type(tem_varsys_type), intent(in) :: varsys

    !> Position of variables to evaluate in varSys
    integer, intent(in) :: varPos(:)

    !> Mesh definition of the input data.
    type(treelmesh_type), intent(in) :: tree

    !> Time information for the current data.
    type(tem_time_type), intent(in) :: time

    !> 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(:,:)

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

    !> Output data
    !! size: io_buffer_size
    real(kind=rk), intent(out) :: res(:)
    ! --------------------------------------------------------------------------!
    integer :: maxComponents, nScalars, nComponents, compOff
    integer :: iPnt, iVar, nVars, res_size
    integer :: e_start, t_start
    real(kind=rk), allocatable :: tmpdat(:)
    ! --------------------------------------------------------------------------!
    ! number of variables to evaluate
    nVars = size(varPos)

    ! Number of scalars in current output
    nScalars = sum(varSys%method%val(varPos(:))%nComponents)

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

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

    compOff = 0
    vars: do iVar=1, nVars

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

      ! get the size of the needed part of the res array
      res_size = nPnts * nComponents

      ! derive the quantities for all the elements in the current chunk
      call varSys%method%val(varpos(iVar))%get_point(               &
        &                                varSys = varSys,           &
        &                                point  = point,            &
        &                                time   = time,             &
        &                                tree   = tree,             &
        &                                nPnts  = nPnts,            &
        &                                res    = tmpdat(:res_size) )

      ! copy the information to the right positions in the result array
      ! res contains results for all variables,
      ! tmpdat is only for one variable
      do iPnt=1,nPnts
        e_start = (iPnt-1)*nScalars + compOff
        t_start = (iPnt-1)*nComponents
        res( (e_start+1) : (e_start+nComponents) )         &
          &  = tmpdat( t_start + 1 : t_start + nComponents )
      end do
      ! Increase the component offset for the next variables.
      compOff = compOff + nComponents
    end do vars

  end subroutine tem_get_point_chunk
  ! **************************************************************************** !


  ! **************************************************************************** !
  !> This function provides the assignment operator of tem_varSys_op_type
  subroutine copy_varOp(left, right)
    ! ---------------------------------------------------------------------------
    !> varSys op to copy to
    type(tem_varSys_op_type), intent(out) :: left
    !> varSys op to copy from
    type(tem_varSys_op_type), intent(in) :: right
    ! ---------------------------------------------------------------------------
    integer :: nVals
    ! ---------------------------------------------------------------------------
    left%myPos = right%myPos

    if (allocated(right%state_varPos)) then
      nVals = size(right%state_varPos)
      allocate(left%state_varPos(nVals))
      left%state_varPos = right%state_varPos
    end if
    if (allocated(right%auxField_varPos)) then
      nVals = size(right%auxField_varPos)
      allocate(left%auxField_varPos(nVals))
      left%auxField_varPos = right%auxField_varPos
    end if
    left%nComponents = right%nComponents
    left%nInputs = right%nInputs
    left%operType = right%operType

    if ( allocated(right%input_varPos) ) then
      allocate(left%input_varPos(size(right%input_varPos)))
      left%input_varPos = right%input_varPos
    end if

    if ( allocated(right%input_varIndex) ) then
      allocate(left%input_varIndex(size(right%input_varIndex)))
      left%input_varIndex = right%input_varIndex
    end if

    left%method_data = right%method_data
    left%get_element => right%get_element
    left%get_point => right%get_point
    left%set_params => right%set_params
    left%get_params => right%get_params
    left%setup_indices => right%setup_indices
    left%get_valOfIndex => right%get_valOfIndex

  end subroutine copy_varOp
  ! ************************************************************************ !
  ! ************************************************************************ !


  ! ************************************************************************ !
  !> Free a variable description.
  !!
  !! Deallocates arrays and nullifies pointers.
  !! Does not touch the method data. Thus you should free the method data
  !! before calling this routine.
  subroutine tem_free_varOp(fun)
    !> Variable method to free.
    type(tem_varSys_op_type), intent(inout) :: fun

    if (allocated(fun%state_varPos)) deallocate(fun%state_varPos)
    if (allocated(fun%auxField_varPos)) deallocate(fun%auxField_varPos)
    fun%nComponents = 0
    fun%nInputs = 0
    if (allocated(fun%input_varpos)) deallocate(fun%input_varpos)
    if (allocated(fun%input_varIndex)) deallocate(fun%input_varIndex)
    nullify(fun%get_point)
    nullify(fun%get_element)
    nullify(fun%set_params)
    nullify(fun%get_params)
    nullify(fun%setup_indices)
    nullify(fun%get_valOfIndex)

  end subroutine tem_free_varOp
  ! ************************************************************************ !
  ! ************************************************************************ !


  ! ************************************************************************ !
  !> Empty a variable system.
  !!
  !! This frees all variables and empties the method and varname arrays.
  !! Note, this does not free the memory used for method data.
  !! Any method data that is to be freed has to be dealt with beforehand.
  subroutine tem_empty_varsys(varsys)
    !> Variable system to empty.
    type(tem_varSys_type), intent(inout) :: varsys

    integer :: iMethod

    do iMethod=1,varsys%method%nVals
      call tem_free_varOp(varsys%method%val(iMethod))
    end do
    call empty(varsys%method)
    call empty(varsys%varname)
    varsys%nStateVars = 0
    varsys%nAuxVars = 0
    varsys%nScalars = 0
    varsys%nAuxScalars = 0
    varsys%systemName = ''

  end subroutine tem_empty_varsys
  ! ************************************************************************ !
  ! ************************************************************************ !


  ! ************************************************************************ !
  subroutine tem_varSys_getParams_dummy(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

    ! Nothing to do here ...
    outstring = trim(instring)
    call tem_abort('Not implemented for '//trim(varSys%varName%val(fun%myPos)) &
      &         // ':is dummy getParams')
  end subroutine tem_varSys_getParams_dummy
  ! ************************************************************************ !

  ! ************************************************************************ !
  subroutine tem_varSys_setParams_dummy(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

    ! Nothing to do here ...
    write(*,*) trim(instring)
    call tem_abort('Not implemented for '//trim(varSys%varName%val(fun%myPos)) &
      &         // ':is dummy setParams')
  end subroutine tem_varSys_setParams_dummy
  ! ************************************************************************ !

  ! ************************************************************************ !
  subroutine tem_varSys_setupIndices_dummy(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: n
    !!
    !! 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(:)

    write(*,*) 'point(1,1):', point(1,1)
    if (present(offset_bit)) write(*,*) 'offset_bit(1):', offset_bit(1)
    write(*,*) 'iLevel:', iLevel
    write(*,*) 'tree%nElems:', tree%nElems
    write(*,*) 'nPnts:', nPnts
    idx = 0
    call tem_abort('Not implemented for '//trim(varSys%varName%val(fun%myPos)) &
      &           //':is dummy setupIndices')
  end subroutine tem_varSys_setupIndices_dummy
  ! ************************************************************************ !


  ! ************************************************************************ !
  subroutine tem_varSys_getValOfIndex_dummy(fun, varSys, time, iLevel, &
    &                                       idx, idxLen, nVals, res)

    !> Description of the method to obtain the variables, here some preset
    !! values might be stored, like the space time function to use or the
    !! required variables.
    class(tem_varSys_op_type), intent(in) :: fun

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

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

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

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

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

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

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

    write(*,*) 'time%sim:', time%sim
    write(*,*) 'iLevel:', iLevel
    write(*,*) 'idx(1):', idx(1)
    if (present(idxLen)) write(*,*) 'offset_bit(1):', idxLen(1)
    write(*,*) 'nVals:', nVals
    res = 0.0_rk
    call tem_abort('Not implemented for '//trim(varSys%varName%val(fun%myPos)) &
      &           //':is dummy getValOfIndex')
  end subroutine tem_varSys_getValOfIndex_dummy
  ! ************************************************************************ !

  ! ************************************************************************ !
  subroutine tem_varSys_getPoint_dummy(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(:)

    write(*,*) 'point(1,1): ', point(1,1)
    write(*,*) 'time%sim:', time%sim
    write(*,*) 'tree%nElems: ', tree%nElems
    write(*,*) 'nPnts:', nPnts
    res = 0.0_rk
    call tem_abort('Not implemented for '//trim(varSys%varName%val(fun%myPos)) &
      &           //':is dummy getPoint')
  end subroutine tem_varSys_getPoint_dummy
  ! ************************************************************************ !


  ! ************************************************************************ !
  subroutine tem_varSys_getElement_dummy(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 element in tree%treeID to get the variable for.
    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 elements 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:
    !! (nComponents of resulting variable) x (nDegrees of freedom) x (nElems)
    !! Access: (iElem-1)*fun%nComponents*nDofs +
    !!         (iDof-1)*fun%nComponents + iComp
    real(kind=rk), intent(out) :: res(:)

    write(*,*) 'elempos(1): ', elempos(1)
    write(*,*) 'time%sim:', time%sim
    write(*,*) 'tree%nElems: ', tree%nElems
    write(*,*) 'nElems: ', nElems
    write(*,*) 'nDofs: ', nDofs
    res = 0.0_rk
    call tem_abort('Not implemented for '//trim(varSys%varName%val(fun%myPos)) &
      &           //':is dummy get_element')
  end subroutine tem_varSys_getElement_dummy
  ! ************************************************************************ !


  ! ************************************************************************** !
  subroutine tem_varSys_check_inArgs( fun, varSys, time, iLevel, idx, idxLen, &
    &                                 nVals, label                            )
    !------------------------------------------------------------------------!
    !> 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(:)

    integer, intent(in) :: nVals
    character(len=*), intent(in) :: label
    !------------------------------------------------------------------------!

    if (size(idx) /= nVals) then
      write(logunit(10),*) 'in ', trim(label), ' idx length /= nVals !'
      write(logunit(10),*) '  varname:', trim(varsys%varname%val(fun%mypos))
      write(logunit(10),*) '  time%sim:', time%sim
      write(logunit(10),*) '  iLevel:', iLevel
      write(logunit(10),*) '  nVals:', nVals
      write(logunit(10),*) '  size(idx):', size(idx)
      if (present(idxLen)) write(logunit(10),*) 'idxLen provided:', idxLen(1)
    end if
  end subroutine tem_varsys_check_inArgs
  ! ************************************************************************** !

end module tem_varSys_module