! See copyright notice at the end of this file.
! **************************************************************************** !
!> author: Harald Klimach
!! author: Kartik Jain
!!
!! TREe based ELemental Mesh module, providing the definition and methods to
!! act on a octree based mesh.
!!
!! The three letter abbrivation TEM is also used to refer to this mesh type.
!! All leaf elements are stored in the treelmesh_type.
!! The definition is designed to be used for parallel executions.
!! It includes a reference to tem_global_module::tem_global_type
!! and the number of elements,
!! which were loaded on the local partition, as well as the first and last
!! tree ID of each partition in the current communicator.
!! This enables the determination of any tree ID's position among the processes,
!! independent of their actual existence.
!! This is important to identify the location of neighbor elements when creating
!! connectivity lists for the linearized element lists locally on each process
!! within the solvers.
!!
!! ![apes_schematic](|media|/apes_schematic.png)
!!
!! See also
!! ---
!!
!! - [[tem_balance_module]] "Load Balancing"
!! - [[tem_dyn_array_module]] "Dynamic Data Structures"
!!
module treelmesh_module

  ! include treelm modules
  use mpi
  use env_module,                only: rk, long_k, long_k_mpi, rk_mpi, &
    &                                  newUnit, PathLen,               &
    &                                  tem_create_EndianSuffix,        &
    &                                  globalMaxLevels, io_buffer_size
  use tem_aux_module,            only: tem_abort, check_mpi_error
  use tem_global_module,         only: tem_global_type, load_tem_global,   &
    &                                  dump_tem_global, tem_global_mesh_free
  use tem_property_module,       only: tem_property_type, gather_property, &
    &                                  prp_hasbnd, prp_sendHalo, prp_hasQVal
  use tem_topology_module,       only: tem_path_type, tem_firstIdAtLevel, &
    &                                  tem_lastIdAtLevel, tem_levelOf,    &
    &                                  tem_IDofCoord
  use tem_logging_module,        only: logUnit
  use tem_tools_module,          only: tem_horizontalSpacer
  use tem_sparta_module,  only: tem_balance_sparta, &
    &                           tem_sparta_type,    &
    &                           tem_init_sparta,    &
    &                           tem_destroy_sparta, &
    &                           tem_exchange_sparta

  ! include aotus modules
  use aotus_module,     only: aot_get_val, flu_state, aoterr_NonExistent, &
    &                         aoterr_WrongType, aoterr_Fatal
  use aot_table_module, only: aot_table_open, aot_table_close, aot_table_length

  implicit none

  private

  public :: load_tem
  public :: unload_treelmesh
  public :: treelmesh_type, load_treelmesh, dump_treelmesh
  public :: free_treelmesh
  public :: generate_treelm_cube
  public :: treelmesh_fromList
  public :: tem_load_internal
  public :: dump_MeshHeader
  public :: exchange_elements
  public :: tem_dump_weights
  public :: tem_load_weights

  !> declaration of the overall mesh on the local partition
  type treelmesh_type
    !> This entry provides some global informations on the mesh, which is needed
    !! across all partitions.
    type(tem_global_type) :: global

    ! @todo: need update after mesh changes
    integer :: nElems !< Total number of Elements on this partition
    integer(kind=long_k) :: elemOffset !< Offset of this partition

    !> List of treeIDs of the first element on each partition, has a length of
    !! global%nParts
    integer(kind=long_k), allocatable :: Part_First(:)

    !> List of treeIDs of the last element on each partition, has a length of
    !! global%nParts
    integer(kind=long_k), allocatable :: Part_Last(:)

    !> The treeIDs of each element (the array has a length of nElems)
    !! these IDs identify the element in the complete mesh, with a breadth-first
    !! numbering scheme. Allowing easy computation of parents and children.
    !! [treeid] \ref treeID
    !!
    !! # TreeID {#treeID}
    !! The coordinate description of an element has three spatial entries and a
    !! fourth entry indicating the refinement level in the octree, resulting in
    !! a tuple of four integers: \f$(x, y, z, L)\f$.
    !!
    !! This coordinate fully describes the spatial shape and position of an
    !! element and is important to determine spatial relations between elements.
    !! It can be obtained from a given treeID by the following procedure.
    !! 1. First, the \ref tem_topology_module::tem_levelof "refinement level" is
    !!    found.
    !! 2. Then a \ref tem_topology_module::tem_firstidatlevel "level-offset" is
    !!    calculated.
    !! 3. After the offset is known, the Morton index of the node in question is
    !!    obtained by subtracting the offset \f$s\f$ from the treeID of the
    !!    node.
    !! 4. Finally a simple bit interleaving rule allows the
    !!    \ref tem_topology_module::tem_coordofid "transformation from the
    !!    Morton index to the three dimensional coordinates".
    !!
    !! The spatial indices are limited by the refinement level:
    !! \f$x, y, z \in \mathbf{Z}_{\ge 0}: x, y, z < 2^L\f$.
    !! To avoid undefined coordinates, movements through the mesh by additions
    !! to indices are done in a modulo(\f$2^L\f$) safeguard resulting in a
    !! periodic universe cube.
    !! This enables the movement in the mesh in horizontal direction (e.g.
    !! for stencil neighbor identification) by index alteration and translation
    !! into a treeID again.
    !! @todo: change to dynamic array
    integer(kind=long_k), allocatable :: treeID(:)

    !> Bit field of element properties for each element (the array has a length
    !! of nElems), each bit in the 8 byte integer indicates the presence or
    !! absence of a given property, which are further described in the
    !! Property component.
    !! @todo: change to dynamic array
    integer(kind=long_k), allocatable :: ElemPropertyBits(:)

    !> declaration of additional elemental properties, the array has a length of
    !! global%nProperties, each property provides a list of local element
    !! indices of elements with this property.
    type(tem_property_type), pointer :: Property(:) => null()

    !> Path for all elements in the treeID array
    type(tem_path_type), allocatable :: pathList(:)

    !> Name of the file that contains weight information for each element to
    !! allow load balancing.
    !! If empty string, balancing will be deactivated.
    character(len=pathLen) :: weights_file = ''

    !> File to write weights to after simulation.
    !! Defaults to the same name as weights_file, if you do not want to have
    !! the input file overwritten, you can provide a different file to write
    !! to by setting write_weights accordingly in the configuration.
    character(len=pathLen) :: write_weights

    !> Level weights as provided by the solver.
    real(kind=rk) :: levelWeight(globalMaxLevels)
  end type


contains


  ! ************************************************************************ !
  !> Load the treelmesh
  !!
  !! Depending on the configuration given in conf, this will either load a
  !! a \ref tem_distributed_octree "tree mesh" from disk with the given prefix
  !! or \ref tem_load_internal "create the mesh internally".
  !! @todo: to calculate levelWeight, global has to be loaded fisrt,
  !!        maybe load_global shoule be moved out of this routine
  subroutine load_tem( me, conf, myPart, nParts, comm, levelWeight, meshDir, &
    &                  parent)
    ! -------------------------------------------------------------------- !
    !> Structure to load the mesh to
    type(treelmesh_type), intent(out) :: me
    !> Directory containing the mesh informations
    type(flu_State) :: conf
    !> Partition to use on the calling process (= MPI Rank in comm)
    integer, intent(in) :: myPart
    !> Number of partitions, the mesh is partitioned into (= Number of MPI
    !! processes in comm).
    integer, intent(in) :: nParts
    !> MPI Communicator to use
    !! now defaults to the one given in the tree%global%comm
    integer, intent(in) :: comm
    !> Balancing weight for elements on different levels.
    !! If these weights are present, the mesh is partitioned initially according
    !! to the given weights.
    !! This can not be combined with given offsets. If offset is also provided
    !! they will overwrite the levelWeight distribution.
    real(kind=rk), intent(in), optional :: levelWeight(globalMaxLevels)
    !> output Mesh directory name
    character(len=*), intent(out), optional :: meshDir
    !> optional parent handle
    integer, intent(in), optional :: parent
    ! -------------------------------------------------------------------- !
    integer :: tem_handle
    character(len=pathLen) :: dirname
    character(len=4) :: EndianSuffix
    integer :: mesh_error, weights_error
    ! The communicator which is used
    integer :: commLocal
    real(kind=rk), allocatable :: weights(:)
    integer(kind=long_k) :: chunksize, remainder
    logical :: found_weights
    type( tem_sparta_type ) :: sparta
    ! -------------------------------------------------------------------- !

    ! Use the incoming communicator
    commLocal = comm

    call tem_horizontalSpacer(fUnit = logUnit(1))

    EndianSuffix = tem_create_EndianSuffix()

    ! Make sure the dirname has some sane value.
    dirname = ''

    call aot_table_open(L       = conf,       &
      &                 parent  = parent,     &
      &                 thandle = tem_handle, &
      &                 key     = 'mesh'      )


    if (tem_handle /= 0) then ! key mesh is a table

      ! The mesh is actually given as a table, parse it and
      ! generate a mesh internally.
      write(logUnit(1),*) 'Generating an internally defined mesh'
      call tem_load_internal(me, conf, tem_handle, myPart, nParts, commLocal)

    else ! tem_handle == 0, mesh is not a table
      ! The mesh is not a table, try to interpret it as a string.
      call aot_get_val( L       = conf,       &
        &               thandle = parent,     &
        &               key     = 'mesh',     &
        &               val     = dirname,    &
        &               ErrCode = mesh_error, &
        &               default = 'mesh/'     )

      if (btest(mesh_error, aoterr_Fatal)) then
        write(logUnit(0),*) 'FATAL Error occured, while retrieving mesh folder!'
        if (btest(mesh_error, aoterr_NonExistent))then
          write(logUnit(0),*) 'Variable not existent!'
          write(logUnit(0),*) 'Using the default value mesh folder: mesh/'
        else if (btest(mesh_error, aoterr_WrongType))then
          write(logUnit(0),*) 'Variable has wrong type!'
          write(logUnit(0),*) 'STOPPING'
          call tem_abort()
        end if
      end if

      ! Load global information
      ! @todo: this routine may already be called outside
      if (associated(me%global%property)) deallocate(me%global%property)
      call load_tem_global(me%global, dirname, myPart, nParts, commLocal)

      ! Homogeneous distribution
      chunksize = me%global%nElems / nParts
      remainder = mod(me%global%nElems, int(nParts, kind=long_k))
      me%nElems = int(chunksize)
      if (myPart < int(remainder)) then
        me%nElems = me%nElems + 1
      end if
      me%elemOffset = chunksize * int(myPart, kind=long_k)     &
        &             + min(int(myPart, kind=long_k), remainder)

      call load_treelmesh( me     = me,    &
        &                  nParts = nParts )

    end if

    ! If requested, provide the directory name to the caller.
    if ( present(meshDir) ) then
      meshDir = trim( dirname )
    end if

    ! Make sure each process has at least one element!
    if (me%global%nElems/me%global%nParts == 0) then
      write(logunit(1),*) ''
      write(logunit(1),*) 'FATAL Error occured:'
      write(logunit(1),*) 'Less elements then processes, this is not'
      write(logunit(1),*) 'supported! And you should reconfigure your'
      write(logunit(1),*) 'setup.'
      write(logunit(1),*) ''
      write(logunit(1),*) 'Restart this simulation with at most ', &
        &                 me%global%nElems, ' processes!'
      call tem_abort()
    end if

    me%levelWeight = 0.0_rk
    if (present(levelWeight)) then
      me%levelWeight = levelWeight
    end if

    ! Get the name for the weights file.
    ! This defaults to the mesh prefix with weights appended.
    call aot_get_val( L       = conf,            &
      &               key     = 'weights',       &
      &               thandle = parent,          &
      &               val     = me%weights_file, &
      &               ErrCode = weights_error,   &
      &               default = ''               )
    if (.not. btest(weights_error, aoterr_NonExistent)) then
      if (myPart == 0 .and. trim(me%weights_file) /= '') then
        inquire(file=trim(me%weights_file)//EndianSuffix, exist=found_weights)
        if (.not. found_weights) then
          write(logUnit(1),*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
          write(logUnit(1),*) 'Weight file '                         &
            &                 // trim(me%weights_file)//EndianSuffix &
            &                 // ' can NOT be found!'
          write(logUnit(1),*) 'May proceed WITHOUT BALANCING!'
          write(logUnit(1),*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
        end if
      end if
    end if

    call aot_get_val( L       = conf,             &
      &               key     = 'write_weights',  &
      &               thandle = parent,           &
      &               val     = me%write_weights, &
      &               ErrCode = weights_error,    &
      &               default = me%weights_file   )

    if (trim(me%write_weights) /= '') then
      write(logunit(1),*) 'The computational weight of each element will be'
      write(logunit(1),*) 'written at the end of the simulation into file'
      write(logunit(1),*) trim(me%write_weights)//EndianSuffix
      write(logunit(1),*) 'If this is not desired, set write_weights explicitly'
      write(logunit(1),*) 'to an empty string.'
    end if

    ! Do sparta balancing, but only if there are actually multiple processes.
    if (me%global%nParts > 1) then
      allocate( weights(me%nElems) )
      call tem_load_weights( me      = me,           &
        &                    weights = weights,      &
        &                    success = found_weights )

      if ( found_weights ) then
        call tem_init_sparta( sparta, nParts )
        call tem_balance_Sparta(weight  = weights,      &
          &                     myPart  = myPart,       &
          &                     nParts  = nParts,       &
          &                     comm    = commLocal,    &
          &                     myElems = me%nElems,    &
          &                     sparta  = sparta,       &
          &                     offset  = me%elemOffset )
        call exchange_elements( me, sparta )
        call tem_destroy_sparta( sparta )
      end if
      deallocate( weights )
    end if

    ! Output basic mesh infos to the screen
    call dump_meshHeader( me )

    call aot_table_close(L=conf, thandle=tem_handle)

    call tem_horizontalSpacer(fUnit = logUnit(1))

  end subroutine load_tem
  ! ************************************************************************ !


  ! ************************************************************************ !
  !> Unload the treelmesh
  !!
  !! Deallocate the data within tree, while keep the global information
  !!
  subroutine unload_treelmesh( me )
    ! -------------------------------------------------------------------- !
    !> Structure to load the mesh to
    type(treelmesh_type), intent(inout) :: me
    ! -------------------------------------------------------------------- !

    write(logUnit(6),*) 'Deallocate the tree data structure'

    if (allocated(me%Part_first)) deallocate( me%Part_First )
    if (allocated(me%Part_last)) deallocate( me%Part_Last  )
    if (allocated(me%treeID)) deallocate( me%treeID  )
    if (allocated(me%ElemPropertyBits)) deallocate( me%ElemPropertyBits  )
    if (associated(me%property)) deallocate(me%property)
    if (allocated(me%pathList)) deallocate( me%pathList  )

  end subroutine unload_treelmesh
  ! ************************************************************************ !


  ! ************************************************************************ !
  !> Completely free the memory of the treelmesh again.
  subroutine free_treelmesh(me)
    !> Structure to free
    type(treelmesh_type), intent(inout) :: me

    call unload_treelmesh(me)
    call tem_global_mesh_free(me%global)
  end subroutine free_treelmesh
  ! ************************************************************************ !


  ! ************************************************************************ !
  !> Load an internally generated mesh
  !!
  !! Currently, only full cubes and slices on a defined level can be created.
  !! The generated meshes are fully periodic.
  !!```lua
  !! -- use the predefined full cube
  !! mesh = { predefined = 'cube',
  !!          -- for a single element in Z-direction use slice:
  !!          -- predefined = 'slice',
  !!
  !!          -- for just a line (single element in Y and Z) use line:
  !!          -- predefined = 'line',
  !!          -- If you set the refinementLevel to 0, a single element will
  !!          -- be created with treeID=1 on level 1, with the bounding cube
  !!          -- redefined such, that the element fills the original
  !!          -- specification for the bounding cube.
  !!          -- The line can also be generated with an arbitrary number of
  !!          -- elements instead of a level. In this case the length refers
  !!          -- to the overall length of the given elements not the extent of
  !!          -- the cube. Periodicity is assumed in X direction as well.
  !!          -- To create the line with a certain number of elements set the
  !!          -- element_count instead of the refinementLevel. If both are
  !!          -- given, the element_count overrules the refinementLevel.
  !!
  !!          -- for a line with boundary conditions line_bounded:
  !!          -- predefined = 'line_bounded'
  !!          -- Like line, but needs to have elementcount given, not
  !!          -- refinementLevel. This mesh defines boundary conditions in the
  !!          -- X direction as west (on the left) and east (on the right), that
  !!          -- then need to be defined in the solver configuration.
  !! -- origin of the cube
  !! origin = {0.,0.,0.},
  !! -- length of the cube
  !! length = 10.,
  !! -- refinement level to resolve the cube
  !! refinementLevel = 4 }
  !!```
  !!
  subroutine tem_load_internal( me, conf, thandle, myPart, nParts, comm )
    ! -------------------------------------------------------------------- !
    !> Structure to load the mesh to
    type(treelmesh_type), intent(out) :: me
    !> Directory containing the mesh informations
    type(flu_State) :: conf
    !> Handle for the table to read the description
    !! of the mesh from.
    integer, intent(in) :: thandle
    !> Partition to use on the calling process (= MPI Rank in comm)
    integer, intent(in) :: myPart
    !> Number of partitions, the mesh is partitioned into (= Number of MPI
    !! processes in comm).
    integer, intent(in) :: nParts
    !> MPI Communicator to use
    integer, intent(in) :: comm
    ! -------------------------------------------------------------------- !
    character(len=40) :: meshtype
    real(kind=rk) :: origin(3)
    real(kind=rk) :: length
    integer :: elementcount
    integer :: level
    integer :: sub_handle
    integer :: iError, i
    integer :: orig_err(3)
    ! -------------------------------------------------------------------- !

    call aot_get_val( L       = conf,         &
      &               thandle = thandle,      &
      &               val     = meshtype,     &
      &               ErrCode = iError,       &
      &               key     = 'predefined', &
      &               default = 'cube'        )

    select case(trim(meshtype))
    case('cube')
      write(logUnit(1),*)'Creating a full cubical mesh without boundaries'

      ! Get the origin of the cube:
      call aot_table_open( L       = conf,       &
        &                  parent  = thandle,    &
        &                  thandle = sub_handle, &
        &                  key     = 'origin'    )
      if (aot_table_length(L = conf, thandle = sub_handle) == 3) then
        do i=1,3
          call aot_get_val( L       = conf,       &
            &               thandle = sub_handle, &
            &               val     = origin(i),  &
            &               ErrCode = iError,     &
            &               pos     = i           )
        end do
      else
        write(logUnit(1),*)'Origin specified for the mesh to generate has to ' &
          &            //'have 3 entries!'
        write(logUnit(1),*)'Default to {0., 0., 0.} '
        origin = 0.0_rk
      end if
      call aot_table_close( L=conf, thandle=sub_handle )

      ! Get the length of the cube:
      call aot_get_val( L       = conf,     &
        &               thandle = thandle,  &
        &               val     = length,   &
        &               ErrCode = iError,   &
        &               key     = 'length', &
        &               default = 1.0_rk    )

      ! Get the refinement level:
      call aot_get_val( L       = conf,             &
        &               thandle = thandle,          &
        &               val     = level,            &
        &               ErrCode = iError,           &
        &               key     = 'refinementLevel' )

      ! Now generate a mesh with the gathered configuration parameters
      if (level > 0) then
        ! Only if the level is greater than 0, actually create mesh.
        call generate_treelm_cube( me     = me,     &
          &                        origin = origin, &
          &                        length = length, &
          &                        level  = level,  &
          &                        myPart = myPart, &
          &                        nParts = nParts, &
          &                        comm   = comm    )
      else
        ! Meshes on level 0 have only one element, but we need treeids > 0,
        ! thus we create single elements on level 1 with treeid=1, and
        ! redefine the bounding cube accordingly.
        call generate_treelm_single( me     = me,       &
          &                          origin = origin,   &
          &                          length = 2*length, &
          &                          myPart = myPart,   &
          &                          nParts = nParts,   &
          &                          comm   = comm      )
      end if


    case('slice')
      write(logUnit(1),*)'Creating a slice (single element in Z) mesh' &
        &          //' without boundaries'

      ! Get the origin of the cube:
      origin = 0.0_rk
      call aot_table_open( L       = conf,       &
        &                  parent  = thandle,    &
        &                  thandle = sub_handle, &
        &                  key     = 'origin'    )
      if (aot_table_length(L = conf, thandle = sub_handle) >= 2) then
        do i=1,2
          call aot_get_val( L       = conf,       &
            &               thandle = sub_handle, &
            &               val     = origin(i),  &
            &               ErrCode = iError,     &
            &               pos     = i           )
        end do
      else
        write(logUnit(1),*)'Origin specified for the mesh to generate has to ' &
          &            // 'have 2 entries!'
        write(logUnit(1),*)'Default to {0., 0.} '
      end if
      call aot_get_val( L       = conf,       &
        &               thandle = sub_handle, &
        &               val     = origin(3),  &
        &               ErrCode = iError,     &
        &               Default = 0.0_rk,     &
        &               pos     = 3           )
      call aot_table_close( L=conf, thandle=sub_handle )

      ! Get the length of the cube:
      call aot_get_val( L       = conf,     &
        &               thandle = thandle,  &
        &               val     = length,   &
        &               ErrCode = iError,   &
        &               key     = 'length', &
        &               default = 1.0_rk    )

      ! Get the refinement level:
      call aot_get_val( L       = conf,             &
        &               thandle = thandle,          &
        &               val     = level,            &
        &               ErrCode = iError,           &
        &               key     = 'refinementLevel' )

      ! Now generate a mesh with the gathered configuration parameters
      if (level > 0) then
        ! Only if the level is greater than 0, actually create mesh.
        call generate_treelm_slice( me     = me,     &
          &                         origin = origin, &
          &                         length = length, &
          &                         level  = level,  &
          &                         myPart = myPart, &
          &                         nParts = nParts, &
          &                         comm   = comm    )
      else
        ! Meshes on level 0 have only one element, but we need treeids > 0,
        ! thus we create single elements on level 1 with treeid=1, and
        ! redefine the bounding cube accordingly.
        call generate_treelm_single( me     = me,       &
          &                          origin = origin,   &
          &                          length = 2*length, &
          &                          myPart = myPart,   &
          &                          nParts = nParts,   &
          &                          comm   = comm      )
      end if


    case('line')
      write(logUnit(1),*) 'Creating a line (single element in Y and Z) mesh' &
        &                 // ' without boundaries'

      ! Get the origin of the cube:
      origin = 0.0_rk
      call aot_table_open( L       = conf,       &
        &                  parent  = thandle,    &
        &                  thandle = sub_handle, &
        &                  key     = 'origin'    )
      if (aot_table_length(L = conf, thandle = sub_handle) >= 1) then
        call aot_get_val( L       = conf,       &
          &               thandle = sub_handle, &
          &               val     = origin(1),  &
          &               ErrCode = iError,     &
          &               pos     = 1           )
      else
        write(logUnit(1),*) 'Origin specified for the mesh to generate has' &
          &                 // ' to have 1 entry!'
        write(logUnit(1),*)'Default to {0.} '
      end if
      call aot_get_val( L       = conf,       &
        &               thandle = sub_handle, &
        &               val     = origin(2),  &
        &               ErrCode = iError,     &
        &               Default = 0.0_rk,     &
        &               pos     = 2           )
      call aot_get_val( L       = conf,       &
        &               thandle = sub_handle, &
        &               val     = origin(3),  &
        &               ErrCode = iError,     &
        &               Default = 0.0_rk,     &
        &               pos     = 3           )
      call aot_table_close( L=conf, thandle=sub_handle )

      ! Get the length of the cube:
      call aot_get_val( L       = conf,     &
        &               thandle = thandle,  &
        &               val     = length,   &
        &               ErrCode = iError,   &
        &               key     = 'length', &
        &               default = 1.0_rk    )

      ! Get the refinement level:
      call aot_get_val( L       = conf,              &
        &               thandle = thandle,           &
        &               val     = level,             &
        &               ErrCode = iError,            &
        &               key     = 'refinementLevel', &
        &               default = -1                 )

      ! Get the element count:
      call aot_get_val( L       = conf,            &
        &               thandle = thandle,         &
        &               val     = elementcount,    &
        &               ErrCode = iError,          &
        &               key     = 'element_count', &
        &               default = -1               )

      if ( (elementcount < 1) .and. (level < 0) ) then
        write(logUnit(1),*) 'For the line mesh you need to specify either'
        write(logUnit(1),*) 'the element_count or the refinementLevel!'
        write(logUnit(1),*) 'STOPPING...'
        call tem_abort()
      end if

      ! Now generate a mesh with the gathered configuration parameters
      if (elementcount > 1) then
        call generate_treelm_elements( me           = me,              &
          &                            origin       = origin,          &
          &                            length       = length,          &
          &                            elementcount = elementcount,    &
          &                            myPart       = myPart,          &
          &                            nParts       = nParts,          &
          &                            comm         = comm,            &
          &                            predefined   = trim(meshtype),  &
          &                            bclabel      = 'internal 1D BC' )
      else
        if ((level > 0) .and. (elementcount < 1)) then
          ! Only if the level is greater than 0, actually create mesh.
          call generate_treelm_line( me     = me,     &
            &                        origin = origin, &
            &                        length = length, &
            &                        level  = level,  &
            &                        myPart = myPart, &
            &                        nParts = nParts, &
            &                        comm   = comm    )
        else
          ! Meshes on level 0 have only one element, but we need treeids > 0,
          ! thus we create single elements on level 1 with treeid=1, and
          ! redefine the bounding cube accordingly.
          call generate_treelm_single( me     = me,       &
            &                          origin = origin,   &
            &                          length = 2*length, &
            &                          myPart = myPart,   &
            &                          nParts = nParts,   &
            &                          comm   = comm      )
        end if
      end if

    case('line_bounded')
      write(logUnit(1),*) 'Creating a line (single element in Y and Z) mesh' &
        &                 // ' WITH boundaries'
      write(logUnit(1),*) 'Boundaries west and east have to be provided!'

      ! Get the origin of the cube:
      origin = 0.0_rk
      call aot_table_open( L       = conf,       &
        &                  parent  = thandle,    &
        &                  thandle = sub_handle, &
        &                  key     = 'origin'    )
      if (aot_table_length(L = conf, thandle = sub_handle) >= 1) then
        call aot_get_val( L       = conf,       &
          &               thandle = sub_handle, &
          &               val     = origin(1),  &
          &               ErrCode = iError,     &
          &               pos     = 1           )
      else
        write(logUnit(1),*) 'Origin specified for the mesh to generate has' &
          &                 // ' to have 1 entry!'
        write(logUnit(1),*)'Default to {0.} '
      end if
      call aot_get_val( L       = conf,       &
        &               thandle = sub_handle, &
        &               val     = origin(2),  &
        &               ErrCode = iError,     &
        &               Default = 0.0_rk,     &
        &               pos     = 2           )
      call aot_get_val( L       = conf,       &
        &               thandle = sub_handle, &
        &               val     = origin(3),  &
        &               ErrCode = iError,     &
        &               Default = 0.0_rk,     &
        &               pos     = 3           )
      call aot_table_close( L=conf, thandle=sub_handle )

      ! Get the length of the cube:
      call aot_get_val( L       = conf,     &
        &               thandle = thandle,  &
        &               val     = length,   &
        &               ErrCode = iError,   &
        &               key     = 'length', &
        &               default = 1.0_rk    )

      ! Get the element count:
      call aot_get_val( L       = conf,            &
        &               thandle = thandle,         &
        &               val     = elementcount,    &
        &               ErrCode = iError,          &
        &               key     = 'element_count', &
        &               default = -1               )

      if ( (elementcount < 1) ) then
        write(logUnit(1),*) 'For the bounded line mesh you need to specify'
        write(logUnit(1),*) 'the element_count!'
        write(logUnit(1),*) 'STOPPING...'
        call tem_abort()
      end if

      ! Now generate a mesh with the gathered configuration parameters
      call generate_treelm_elements( me           = me,                      &
        &                            origin       = origin,                  &
        &                            length       = length,                  &
        &                            elementcount = elementcount,            &
        &                            myPart       = myPart,                  &
        &                            nParts       = nParts,                  &
        &                            comm         = comm,                    &
        &                            predefined   = trim(meshtype),          &
        &                            bclabel      = 'bounded internal 1D BC' )

    case('single')
      write(logUnit(1),*) 'Creating a single element mesh without boundaries'

      ! Get the origin of the cube:
      call aot_get_val( L       = conf,                    &
        &               thandle = thandle,                 &
        &               val     = origin,                  &
        &               ErrCode = orig_err,                &
        &               key     = 'origin',                &
        &               default = [0.0_rk, 0.0_rk, 0.0_rk] )

      ! Get the length of the cube:
      call aot_get_val( L       = conf,     &
        &               thandle = thandle,  &
        &               val     = length,   &
        &               ErrCode = iError,   &
        &               key     = 'length', &
        &               default = 1.0_rk    )

      ! Now generate a mesh with the gathered configuration parameters
      call generate_treelm_single( me     = me,     &
        &                          origin = origin, &
        &                          length = length, &
        &                          myPart = myPart, &
        &                          nParts = nParts, &
        &                          comm   = comm    )


    case('default')
      write(logUnit(1),*) 'Do not know how to generate mesh ' &
        &                 // trim(meshtype)

    end select

  end subroutine tem_load_internal
  ! ************************************************************************ !


  ! ************************************************************************ !
  !> This subroutine reads a mesh in treelm format from disk at the
  !! specified directory name.
  !!
  subroutine load_treelmesh(me, nParts )
    ! -------------------------------------------------------------------- !
    !> Structure to load the mesh to
    type(treelmesh_type), intent(inout) :: me
    !> Directory containing the mesh informations
    ! character(len=*), intent(in) :: dirname
    !> Partition to use on the calling process (= MPI Rank in comm)
    ! integer, intent(in) :: myPart
    !> Number of partitions, the mesh is partitioned into (= Number of MPI
    !! processes in comm).
    integer, intent(in) :: nParts
    ! -------------------------------------------------------------------- !
    logical :: ex
    integer :: iElem
    integer :: iProp
    integer :: fh, etype, ftype
    integer :: iostatus( MPI_STATUS_SIZE )
    integer :: iError
    integer :: file_rec_len
    integer :: typesize
    character(len=300)   :: ElemFileName
    character(len=4)     :: EndianSuffix
    integer(kind=long_k), allocatable :: buffer(:)
    integer(kind=long_k), allocatable :: globbuffer(:)
    integer(kind=MPI_OFFSET_KIND)     :: displacement
    ! -------------------------------------------------------------------- !

    EndianSuffix = tem_create_EndianSuffix()
    ElemFileName = trim(me%global%dirname)//'elemlist'//EndianSuffix
    write(logUnit(1), *) 'Load mesh from file: '//trim(ElemFileName)

    if (associated(me%property)) deallocate(me%property)
    allocate(me%Property(me%global%nProperties))
    allocate(me%Part_First(nParts))
    allocate(me%Part_Last(nParts))

    allocate(me%treeID(me%nElems))

    if (me%global%myPart == 0) then
      ! Only one process should check the file existence
      inquire(file=trim(ElemFileName), exist=ex)
      if (.not. ex) then
        write(logUnit(1),*) 'ERROR: File ' // trim(ElemFileName) &
          &                 // ' not found!'
        write(logUnit(1),*) 'Aborting.'
        call tem_abort()
      end if
    end if

    allocate(me%ElemPropertyBits(me%nElems))
    allocate(buffer(me%nElems*2))
    if (me%global%nElems*2 > io_buffer_size) then
      write(logUnit(1),*) 'using MPI to read'

      ! Create a contiguous type to describe the vector per element
      call MPI_Type_contiguous( 2, long_k_mpi, etype, iError )
      call check_mpi_error(iError,'type etype in load_treelmesh')
      call MPI_Type_commit(     etype, iError )
      call check_mpi_error(iError,'commit etype in load_treelmesh')
      call MPI_Type_size(etype, typesize, iError )
      call check_mpi_error(iError,'typesize in load_treelmesh')

      ! Create a MPI Subarray  as ftype for file view
      call MPI_Type_contiguous( me%nElems, etype, ftype, iError )
      call check_mpi_error(iError,'type ftype in load_treelmesh')
      call MPI_Type_commit(     ftype, iError )
      call check_mpi_error(iError,'commit ftype in load_treelmesh')

      ! Calculate displacement for file view
      displacement = me%elemOffset * typesize * 1_MPI_OFFSET_KIND

      ! Set the view for each process on the file above
      ! Open the binary file for MPI I/O (Write)
      call MPI_File_open( me%global%comm, trim(ElemFileName),        &
        &                 MPI_MODE_RDONLY, MPI_INFO_NULL, fh, iError )
      call check_mpi_error(iError,'file_open in load_treelmesh')

      call MPI_File_set_view( fh, displacement, etype, ftype, "native", &
        &                     MPI_INFO_NULL, iError                     )
      call check_mpi_error(iError,'set_view in load_treelmesh')
      ! Read data from the file
      call MPI_File_read_all( fh, buffer, me%nElems, etype, iostatus, iError )
      call check_mpi_error(iError,'read_all in load_treelmesh')

      !Free the MPI_Datatypes which were created and close the file
      call MPI_Type_free (etype, iError)
      call check_mpi_error(iError,'free etype in load_treelmesh')
      call MPI_Type_free (ftype, iError)
      call check_mpi_error(iError,'free ftype in load_treelmesh')
      call MPI_File_close(fh,    iError)
      call check_mpi_error(iError,'close file in load_treelmesh')

      ! END IO-part

    else

      !> The mesh is so small, it probably is better to read it on one process
      !! and distribute the resulting data via the network.
      write(logUnit(1),*) 'using one process to read and broadcast'
      allocate( globbuffer(me%global%nElems*2) )
      if (me%global%myPart == 0) then !! single proc read

        inquire(iolength = file_rec_len) globbuffer
        fh = newunit()
        open( unit = fh,                 &
          &   file = trim(ElemFileName), &
          &   recl = file_rec_len,       &
          &   action = 'read',           &
          &   access = 'direct',         &
          &   form   = 'unformatted'     )

        read(fh, rec=1) globbuffer

        close(fh)

      end if !! single proc read

      call MPI_Bcast( globbuffer, int(2*me%global%nElems), long_k_mpi, 0, &
        &             me%global%comm, iError                              )

      buffer = globbuffer(me%elemOffset*2+1:(me%ElemOffset+me%nElems)*2)
      deallocate(globbuffer)

    end if ! END IO-part

    ! Fill the tree with what was read in buffer from disk
    do iElem = 1, me%nElems
      me%treeID(iElem)           = buffer( (iElem-1)*2 + 1 )
      me%ElemPropertyBits(iElem) = buffer( (iElem-1)*2 + 2 )
    end do

    deallocate( buffer )

    ! clear the prp_sendHalo property since the distribution of elements
    ! might have changed
    do iElem=1,me%nElems
      me%ElemPropertyBits(iElem)  = ibclr( me%ElemPropertyBits(iElem), &
        &                                  prp_sendHalo                )
    end do

    call MPI_allgather( me%treeID(1), 1, long_k_mpi,  &
      &                 me%Part_First, 1, long_k_mpi, &
      &                 me%global%comm, iError        )
    call MPI_allgather( me%treeID(me%nElems), 1, long_k_mpi, &
      &                 me%Part_Last, 1, long_k_mpi,         &
      &                 me%global%comm, iError               )

    do iProp=1,me%global%nProperties
      call gather_Property( Property = me%Property(iProp),        &
        &                   Header   = me%global%Property(iProp), &
        &                   BitField = me%ElemPropertyBits,       &
        &                   comm     = me%global%comm             )
    end do

    write(logUnit(1),*) 'Done, reading the mesh!'

  end subroutine load_treelmesh
  ! ************************************************************************ !


  ! ************************************************************************ !
  !> Write a given mesh to disk. it is stored to the directory given in the
  !! [[tem_global_module:tem_global_type]].
  !!
  !! Dump treelmesh_type%global to header.lua
  !! Dump treeID and propertyBits to elemlist.lsb (little endian environments)
  !! or elemlist.msb (big endian environments)
  !!
  subroutine dump_treelmesh( me, root_only )
    ! -------------------------------------------------------------------- !
    !> Mesh to dump to disk
    type(treelmesh_type), intent(inout) :: me
    !> root dump global mesh when true and
    !! all process dump its own mesh when false
    logical, intent(in), optional       :: root_only
    ! -------------------------------------------------------------------- !
    integer(kind=long_k), allocatable   :: buffer(:)
    integer                             :: fh, etype, ftype
    integer                             :: iostatus( MPI_STATUS_SIZE )
    integer                             :: iElem
    integer                             :: iError
    integer                             :: typesize
    logical                             :: root_not_participating
    character(len=300)                  :: ElemFileName
    character(len=4)                    :: EndianSuffix
    integer(kind=long_k)                :: file_offset
    integer(kind=MPI_OFFSET_KIND)       :: displacement
    ! -------------------------------------------------------------------- !

    if(present(root_only)) then
      root_not_participating = root_only
    else
      root_not_participating = .true.
    endif

    EndianSuffix = tem_create_EndianSuffix()

    ! Dump global info to header.lua
    ! if root_not_participating = .false. that each process dumps its own mesh
    ! set global%nElems to local%nElems to dump local nElems
    ! in header file of each process
    if (.not. root_not_participating) me%global%nElems = me%nElems

    call dump_tem_global(me%global, root_only)
    ElemFileName = trim(me%global%dirname) // 'elemlist' // EndianSuffix

    file_offset = 0
    allocate( buffer(me%nElems * 2) )

    ! Get the offset for MPI File view
    if(root_not_participating) then
      call MPI_Exscan( int(me%nElems,long_k), file_offset, 1, &
        &              long_k_mpi, MPI_SUM, me%global%comm, iError)
    else
      file_offset = 0
    endif

    ! Fill the buffer with data to be written on disk
    do iElem = 1, me%nElems
      buffer((iElem-1)*2+1) = me%treeID(iElem)
      buffer((iElem-1)*2+2) = me%ElemPropertyBits(iElem)
    end do

    write(logUnit(1),*) 'using MPI to write'

    ! Open the binary file for MPI I/O (Write)
    call MPI_File_open( me%global%comm,                       &
      &                 trim(ElemFileName),                   &
      &                 ior(MPI_MODE_WRONLY,MPI_MODE_CREATE), &
      &                 MPI_INFO_NULL,                        &
      &                 fh, iError                            )

    ! Catch if there was an exception by MPI
    call check_mpi_error(iError,'file_open in dump_treelmesh')

    ! Create a contiguous type to describe the vector per element
    call MPI_Type_contiguous( 2, long_k_mpi, etype, iError )
    call check_mpi_error(iError,'type etype in dump_treelmesh')
    call MPI_Type_commit(     etype, iError )

    ! Catch MPI Exception(s)
    call check_mpi_error(iError,'commit etype in dump_treelmesh')

    !get size of etype
    call MPI_Type_size(etype, typesize, iError )
    call check_mpi_error(iError,'typesize in dump_treelmesh')

    ! Create a MPI Contiguous as ftype for file view
    call MPI_Type_contiguous( me%nElems, etype, ftype, iError )
    call check_mpi_error(iError,'type ftype in dump_treelmesh')
    call MPI_Type_commit(     ftype, iError )


    ! Catch MPI Exception(s)
    call check_mpi_error(iError,'commit ftype in dump_treelmesh')

    ! calculate dsplacement for file view
    displacement = file_offset * typesize * 1_MPI_OFFSET_KIND

    ! Set the view for each process on the file above
    call MPI_File_set_view( fh, displacement, etype, ftype, "native", &
      &                     MPI_INFO_NULL, iError                     )

    ! Catch MPI Exception(s)
    call check_mpi_error(iError,'set_view in dump_treelmesh')

    ! Write the Data to File
    call MPI_File_write_all( fh, buffer, me%nElems, etype, iostatus, &
      &                      iError                                  )

    ! Catch MPI Exception(s)
    call check_mpi_error(iError,'write_all in dump_treelmesh')

    !Free the MPI_Datatypes which were created and close the file
    call MPI_Type_free(etype, iError)
    call check_mpi_error(iError,'free etype in dump_treelmesh')
    call MPI_Type_free(ftype, iError)
    call check_mpi_error(iError,'free ftype in dump_treelmesh')
    call MPI_File_close(fh,    iError)
    call check_mpi_error(iError,'close file in dump_treelmesh')

    deallocate( buffer )

  end subroutine dump_treelmesh
  ! ************************************************************************ !

  ! ************************************************************************ !
  !> summary: Generate a predefined cube
  !! This serves as an simple grid generation for performance or scaling
  !! analysis without being obliged to use Seeder.
  !! You have to specify the generic grid parameters in the lua file insted of
  !! the mesh folder
  !!```lua
  !! mesh = { predefined='cube',
  !!          origin = {0.,0.,0.},
  !!          length = 10.,
  !!          refinementLevel = 6 }
  !!```
  !! You have to specify the shape 'cube', a bounding box origin, its length
  !! and also the refinement level, which results in different amount elements
  !! in the grid.
  !! The result of this routine is mainly the treeID list with the additional
  !! lists for saving the properties. They are all set to zero here, however.
  !! As we only have a simple cube which includes all the elements on this
  !! level, the treeID list just contains contiguously increasing integers
  !!
  subroutine generate_treelm_cube( me, origin, length, level,  myPart, nParts, &
    &                              comm )
    ! -------------------------------------------------------------------- !
    !> Mesh to generate
    type(treelmesh_type), intent(out) :: me
    !> Corner of the cube
    real(kind=rk), intent(in) :: origin(3)
    !> Length of cube
    real(kind=rk), intent(in) :: length
    !> Resolution level
    integer, intent(in) :: level
    !> Partition of the caller (starts with 0)
    integer, intent(in) :: myPart
    !> Number of partitions
    integer, intent(in) :: nParts
    !> communicator to be used
    integer, intent(in) :: comm
    ! -------------------------------------------------------------------- !
    integer(kind=long_k) :: firstID, lastID
    integer(kind=long_k) :: share
    integer :: remainder
    integer :: iPart, iElem
    ! -------------------------------------------------------------------- !

    me%global%nParts = nParts
    me%global%myPart = myPart
    me%global%comm = comm

    me%global%origin = origin
    me%global%BoundingCubeLength = length
    me%global%minLevel = level
    me%global%maxLevel = level
    me%global%label = 'Generic_Cube'
    me%global%predefined = 'cube'
    write(me%global%comment,'(a15,i7,a16,i2,a1)')                   &
      &     'Generated with ', nParts, ' Parts on Level ', level, '.'
    me%global%dirname = './'

    ! No properties in this mesh
    me%global%nProperties = 0
    if (associated(me%global%Property)) deallocate(me%global%Property)
    allocate(me%global%Property(me%global%nProperties))

    allocate(me%Part_First(nParts))
    allocate(me%Part_Last(nParts))

    ! Compute the treeIDs of the mesh:
    firstID = tem_firstIdAtLevel(level)
    lastID = tem_lastIdAtLevel(level)

    ! Total number of elements in this mesh
    me%global%nElems = lastID - firstID + 1_long_k

    share = me%global%nElems / int(nParts, kind=long_k)
    remainder = int(mod(me%global%nElems, int(nParts, kind=long_k)))

    ! The first partition starts always with the firstID
    me%Part_First(1) = firstID

    ! Up to remainder partitions have share + 1 elements
    do iPart=2,remainder+1
      me%Part_Last(iPart-1) = me%Part_First(iPart-1) + share
      me%Part_First(iPart) = me%Part_Last(iPart-1) + 1
    end do

    ! The remaining elements get exactly the share elements:
    do iPart=remainder+2,nParts
      me%Part_Last(iPart-1) = me%Part_First(iPart-1) + share - 1
      me%Part_First(iPart) = me%Part_Last(iPart-1) + 1
    end do

    ! The last partition ends always with the lastID
    me%Part_Last(nParts) = lastID

    ! Local data:
    me%nElems = int(me%Part_Last(myPart+1) - me%Part_First(myPart+1)) + 1
    me%elemOffset = me%Part_first(myPart+1) - firstID

    allocate(me%treeID(me%nElems))
    allocate(me%ElemPropertyBits(me%nElems))
    if (associated(me%property)) deallocate(me%property)
    allocate(me%Property(me%global%nProperties))

    ! No Properties:
    me%ElemPropertyBits = 0_long_k

    ! Filling the treeIDs:
    do iElem=1,me%nElems
      me%treeID(iElem) = me%Part_First(myPart+1) + int(iElem - 1, kind=long_k)
    end do

  end subroutine generate_treelm_cube
  ! ************************************************************************ !


  ! ************************************************************************ !
  function tem_CoordOf_2d_Id(id2d, level) result(coord)
    ! -------------------------------------------------------------------- !
    !> input element ID
    integer(kind=long_k), intent(in) :: Id2d
    integer, intent(in) :: level
    !> coordinate of element return value
    integer :: coord(4)
    ! -------------------------------------------------------------------- !
    integer(kind=long_k) :: fak(2)
    integer :: bitlevel
    integer :: i
    integer(kind=long_k) :: telem
    ! -------------------------------------------------------------------- !

    coord(:) = 0
    coord(4) = level
    fak(1) = 1
    fak(2) = 2
    bitlevel = 1

    tElem = id2d

    ! get x coordinate from
    !
    ! x = sum(iLevel=0, inf) { 2**iLevel * mod(tElem / (4**iLevel) ),2)   }
    ! y = sum(iLevel=0, inf) { 2**iLevel * mod(tElem /(2 * (4**iLevel) ),2)   }
    !
    do while ((tElem / fak(1)) > 0)
      do i=1,2
        coord(i) = coord(i) + bitlevel * int(mod(tElem / fak(i), 2_long_k))
      end do
      bitlevel = bitlevel * 2
      fak = fak * 4
    end do
  end function tem_CoordOf_2d_Id
  ! ************************************************************************ !


  ! ************************************************************************ !
  !> summary: Generate a predefined slice
  !! This serves as an simple grid generation for performance or scaling
  !! analysis without being obliged to use Seeder.
  !! You have to specify the generic grid parameters in the lua file instead of
  !! the mesh folder
  !!```lua
  !! mesh = { predefined='slice',
  !!          origin = {0.,0.,0.},
  !!          length = 10.,
  !!          refinementLevel = 6 }
  !!```
  !! You have to specify the shape 'slice', a bounding box origin, its length
  !! and also the refinement level, which results in different amount elements
  !! in the grid.
  !! The result of this routine is mainly the treeID list with the additional
  !! lists for saving the properties.
  !! The generated slice will be all elements in the first layer (XY plane) in
  !! Z-direction, with periodicity in Z direction.
  !!
  subroutine generate_treelm_slice( me, origin, length, level,  myPart, &
    &                               nParts, comm )
    ! -------------------------------------------------------------------- !
    !> Mesh to generate
    type(treelmesh_type), intent(out) :: me
    !> Corner of the cube
    real(kind=rk), intent(in) :: origin(3)
    !> Length of cube
    real(kind=rk), intent(in) :: length
    !> Resolution level
    integer, intent(in) :: level
    !> Partition of the caller (starts with 0)
    integer, intent(in) :: myPart
    !> Number of partitions
    integer, intent(in) :: nParts
    !> communicator to be used
    integer, intent(in) :: comm
    ! -------------------------------------------------------------------- !
    integer(kind=long_k) :: firstID, lastID
    integer(kind=long_k) :: share
    integer(kind=long_k) :: previous
    integer :: remainder
    integer :: iPart, iElem
    integer :: coord(4)
    integer :: lastcoord
    ! -------------------------------------------------------------------- !

    me%global%nParts = nParts
    me%global%myPart = myPart
    me%global%comm = comm

    me%global%origin = origin
    me%global%BoundingCubeLength = length
    me%global%minLevel = level
    me%global%maxLevel = level
    me%global%label = 'Generic_Slice'
    me%global%predefined = 'slice'
    write(me%global%comment,'(a15,i7,a16,i2,a1)')                   &
      &     'Generated with ', nParts, ' Parts on Level ', level, '.'
    me%global%dirname = './'

    ! Boundary property to define periodic boundary
    me%global%nProperties = 1
    if (associated(me%global%Property)) deallocate(me%global%property)
    if (associated(me%Property)) deallocate(me%property)
    allocate(me%global%Property(me%global%nProperties))
    allocate(me%Property(me%global%nProperties))

    allocate(me%Part_First(nParts))
    allocate(me%Part_Last(nParts))

    ! Compute the treeIDs of the mesh:
    firstID = tem_firstIdAtLevel(level)
    lastcoord = 2**level - 1
    lastID = tem_IdOfCoord( coord  = [lastcoord, lastcoord, 0, level], &
      &                     offset = firstID                           )

    ! Total number of elements in this mesh
    me%global%nElems = 4_long_k**level

    share = me%global%nElems / int(nParts, kind=long_k)
    remainder = int(mod(me%global%nElems, int(nParts, kind=long_k)))

    ! The first partition starts always with the firstID
    me%Part_First(1) = firstID

    ! Up to remainder partitions have share + 1 elements
    previous = 0_long_k
    do iPart=2,remainder+1
      previous = previous + int(share)
      coord = tem_coordof_2d_id(previous, level)
      me%Part_Last(iPart-1) = tem_idofcoord(coord, offset = firstID)
      previous = previous + 1
      coord = tem_coordof_2d_id(previous, level)
      me%Part_First(iPart) = tem_idofcoord(coord, offset = firstID)
    end do

    ! The remaining elements get exactly the share elements:
    do iPart=remainder+2,nParts
      previous = previous + int(share) - 1
      coord = tem_coordof_2d_id(previous, level)
      me%Part_Last(iPart-1) = tem_idofcoord(coord, offset = firstID)
      previous = previous + 1
      coord = tem_coordof_2d_id(previous, level)
      me%Part_First(iPart) = tem_idofcoord(coord, offset = firstID)
    end do

    ! The last partition ends always with the lastID
    me%Part_Last(nParts) = lastID

    ! Local data:
    if (myPart < remainder) then
      me%nElems = int(share+1)
      me%elemOffset = int(myPart, kind=long_k) * (share+1_long_k)
    else
      me%nElems = int(share)
      me%elemOffset = (int(myPart, kind=long_k) * share) &
        &           + int(remainder, kind=long_k)
    end if
    ! All elements have (periodic) boundaries
    me%Property(1)%nElems = me%nElems
    me%Property(1)%offset = me%elemOffset
    allocate(me%Property(1)%ElemID(me%Property(1)%nElems))
    ! Please note, that the tem_bc_prop_module will set boundary conditions
    ! based on this label accordingly!
    me%global%Property(1)%label = 'internal 2D BC'
    me%global%Property(1)%bitpos = prp_hasBnd
    me%global%Property(1)%nElems = me%Property(1)%nElems

    allocate(me%treeID(me%nElems))
    allocate(me%ElemPropertyBits(me%nElems))

    ! No Properties:
    me%ElemPropertyBits = ibset(0_long_k, prp_hasBnd)

    ! Filling the treeIDs:
    do iElem=1,me%nElems
      me%Property(1)%ElemID(iElem) = iElem
      coord = tem_coordof_2d_id(me%elemOffset + iElem - 1, level)
      me%treeID(iElem) = tem_idofcoord(coord, offset=firstID)
    end do

  end subroutine generate_treelm_slice
  ! ************************************************************************ !


  ! ************************************************************************ !
  !> summary: Generate a predefined line
  !> This serves as an simple grid generation for performance or scaling
  !! analysis without being obliged to use Seeder.
  !! You have to specify the generic grid parameters in the lua file instead of
  !! the mesh folder
  !! ~~~~~~~~~~~~~~~~~~~~~{.lua}
  !! mesh = { predefined='line',
  !!          origin = {0.,0.,0.},
  !!          length = 10.,
  !!          refinementLevel = 6 }
  !! ~~~~~~~~~~~~~~~~~~~~~
  !! \n
  !! You have to specify the shape 'line', a bounding box origin, its length
  !! and also the refinement level, which results in different amount elements
  !! in the grid.\n
  !! The result of this routine is mainly the treeID list with the additional
  !! lists for saving the properties.
  !! The generated line will be all elements along the X-Axis
  !! with periodicity in Y and Z direction.
  !!
  subroutine generate_treelm_line( me, origin, length, level,  myPart, &
    &                              nParts, comm )
    ! -------------------------------------------------------------------- !
    !> Mesh to generate
    type(treelmesh_type), intent(out) :: me
    !> Corner of the cube
    real(kind=rk), intent(in) :: origin(3)
    !> Length of cube
    real(kind=rk), intent(in) :: length
    !> Resolution level
    integer, intent(in) :: level
    !> Partition of the caller (starts with 0)
    integer, intent(in) :: myPart
    !> Number of partitions
    integer, intent(in) :: nParts
    !> communicator to be used
    integer, intent(in) :: comm
    ! -------------------------------------------------------------------- !
    integer(kind=long_k) :: firstID, lastID
    integer(kind=long_k) :: share
    integer :: remainder
    integer :: iPart, iElem
    integer :: coord(4)
    integer :: lastcoord
    ! -------------------------------------------------------------------- !

    me%global%nParts = nParts
    me%global%myPart = myPart
    me%global%comm = comm

    me%global%origin = origin
    me%global%BoundingCubeLength = length
    me%global%minLevel = level
    me%global%maxLevel = level
    me%global%label = 'Generic_Line'
    me%global%predefined = 'line'
    write(me%global%comment,'(a15,i7,a16,i2,a1)')                   &
      &     'Generated with ', nParts, ' Parts on Level ', level, '.'
    me%global%dirname = './'

    ! Boundary property to define periodic boundary
    me%global%nProperties = 1
    if (associated(me%global%property)) deallocate(me%global%property)
    if (associated(me%property)) deallocate(me%property)
    allocate(me%global%Property(me%global%nProperties))
    allocate(me%Property(me%global%nProperties))

    allocate(me%Part_First(nParts))
    allocate(me%Part_Last(nParts))

    ! Compute the treeIDs of the mesh:
    firstID = tem_firstIdAtLevel(level)
    lastcoord = 2**level - 1
    lastID = tem_IdOfCoord( coord = [lastcoord, 0, 0, level], offset = firstID )

    ! Total number of elements in this mesh
    me%global%nElems = 2_long_k**level

    share = me%global%nElems / int(nParts, kind=long_k)
    remainder = int(mod(me%global%nElems, int(nParts, kind=long_k)))

    ! The first partition starts always with the firstID
    me%Part_First(1) = firstID

    ! Up to remainder partitions have share + 1 elements
    coord = 0
    coord(4) = level
    do iPart=2,remainder+1
      coord(1) =  coord(1) + int(share)
      me%Part_Last(iPart-1) = tem_idofcoord(coord, offset = firstID)
      coord(1) = coord(1) + 1
      me%Part_First(iPart) = tem_idofcoord(coord, offset = firstID)
    end do

    ! The remaining elements get exactly the share elements:
    do iPart=remainder+2,nParts
      coord(1) = coord(1) + int(share) - 1
      me%Part_Last(iPart-1) = tem_idofcoord(coord, offset = firstID)
      coord(1) = coord(1) + 1
      me%Part_First(iPart) = tem_idofcoord(coord, offset = firstID)
    end do

    ! The last partition ends always with the lastID
    me%Part_Last(nParts) = lastID

    ! Local data:
    if (myPart < remainder) then
      me%nElems = int(share+1)
      me%elemOffset = int(myPart, kind=long_k) * (share+1_long_k)
    else
      me%nElems = int(share)
      me%elemOffset = (int(myPart, kind=long_k) * share) &
        &           + int(remainder, kind=long_k)
    end if
    ! All elements have (periodic) boundaries
    me%Property(1)%nElems = me%nElems
    me%Property(1)%offset = me%elemOffset
    allocate(me%Property(1)%ElemID(me%Property(1)%nElems))
    ! Please note, that the tem_bc_prop_module will set boundary conditions
    ! based on this label accordingly!
    me%global%Property(1)%label = 'internal 1D BC'
    me%global%Property(1)%bitpos = prp_hasBnd
    me%global%Property(1)%nElems = me%Property(1)%nElems

    allocate(me%treeID(me%nElems))
    allocate(me%ElemPropertyBits(me%nElems))

    ! Only has boundary property:
    me%ElemPropertyBits = ibset(0_long_k, prp_hasBnd)

    ! Filling the treeIDs:
    do iElem=1,me%nElems
      me%Property(1)%ElemID(iElem) = iElem
      ! We can only have 2**20 elements per dimension, so in this case where
      ! we create a line, we won't exceed the integer limit. Thus we safely can
      ! cast the long integer elemOffset to a normal integer.
      coord(1) = int(me%elemOffset) + iElem - 1
      me%treeID(iElem) = tem_idofcoord(coord, offset=firstID)
    end do

  end subroutine generate_treelm_line
  ! ************************************************************************ !


  ! ************************************************************************ !
  !> Generate a predefined line with a given number of elements
  !!
  !! This serves as an simple grid generation for performance or scaling
  !! analysis without being obliged to use Seeder.
  !! You have to specify the generic grid parameters in the lua file instead of
  !! the mesh folder
  !! ~~~~~~~~~~~~~~~~~~~~~{.lua}
  !! mesh = { predefined='line', -- or: 'line_bounded'
  !!          origin = {0.,0.,0.},
  !!          length = 10.,
  !!          elementcount = 6 }
  !! ~~~~~~~~~~~~~~~~~~~~~
  !!
  !! You have to specify the shape 'line', a bounding box origin, its length
  !! and the number of elements, which results in different amount elements
  !! in the grid.\n
  !! The result of this routine is mainly the treeID list with the additional
  !! lists for saving the properties.
  !! The generated line will be a line of elementcount elements along the
  !! X-Axis with periodicity in all directions.
  !!
  !! It is also possible to generate a mesh with boundary conditions in the X
  !! direction (west and east), by using the predefined 'line_bounded'.
  !! In this case these two boundary conditions need to be provided in the
  !! solver configuration.
  subroutine generate_treelm_elements( me, origin, length, elementcount, &
    &                                  myPart, nParts, comm, predefined, &
    &                                  bclabel )
    ! -------------------------------------------------------------------- !
    !> Mesh to generate
    type(treelmesh_type), intent(out) :: me
    !> Corner of the cube
    real(kind=rk), intent(in) :: origin(3)
    !> Length of cube
    real(kind=rk), intent(in) :: length
    !> Number of elements in the line
    integer, intent(in) :: elementcount
    !> Partition of the caller (starts with 0)
    integer, intent(in) :: myPart
    !> Number of partitions
    integer, intent(in) :: nParts
    !> communicator to be used
    integer, intent(in) :: comm
    !> Label describing the internal mesh.
    character(len=*), intent(in) :: predefined
    !> Label describing the boundary conditions to set for this mesh.
    character(len=*), intent(in) :: bclabel
    ! -------------------------------------------------------------------- !
    integer :: level
    integer(kind=long_k) :: firstID, lastID
    integer(kind=long_k) :: share
    integer :: remainder
    integer :: iPart, iElem
    integer :: coord(4)
    integer :: lastcoord
    integer :: xbound_pad
    ! -------------------------------------------------------------------- !

    me%global%nParts = nParts
    me%global%myPart = myPart
    me%global%comm = comm

    ! Face definitions make the first and last face coincide.
    ! To allow boundary definitions on both sides, we need to make sure, there
    ! is at least one additional element.
    if (predefined == 'line_bounded') then
      xbound_pad = 1
    else
      xbound_pad = 0
    end if

    ! Find an appropriate level of at least 1.
    level = max( ceiling(log(real(elementcount+xbound_pad,kind=rk)) &
      &                  / log(2.0_rk)),            1               )

    me%global%origin = origin
    me%global%BoundingCubeLength = length*(real(2**level, kind=rk) &
      &                                    /real(elementcount, kind=rk))
    me%global%minLevel = level
    me%global%maxLevel = level
    me%global%label = 'Generic_Line'
    me%global%predefined = predefined
    write(me%global%comment,'(a15,i7,a16,i2,a1)')                             &
      &               'Generated with ', nParts, ' Parts and ', elementcount, &
      &               ' elements.'
    me%global%dirname = './'

    ! Boundary property to define periodic boundary
    me%global%nProperties = 1
    if (associated(me%global%property)) deallocate(me%global%property)
    if (associated(me%property)) deallocate(me%property)
    allocate(me%global%Property(me%global%nProperties))
    allocate(me%Property(me%global%nProperties))

    allocate(me%Part_First(nParts))
    allocate(me%Part_Last(nParts))

    ! Compute the treeIDs of the mesh:
    firstID = tem_firstIdAtLevel(level)
    lastcoord = elementcount - 1
    lastID = tem_IdOfCoord(coord = [lastcoord, 0, 0, level], offset = firstID )

    ! Total number of elements in this mesh
    me%global%nElems = elementcount

    share = me%global%nElems / int(nParts, kind=long_k)
    remainder = int(mod(me%global%nElems, int(nParts, kind=long_k)))

    ! The first partition starts always with the firstID
    me%Part_First(1) = firstID

    ! Up to remainder partitions have share + 1 elements
    coord = 0
    coord(4) = level
    do iPart=2,remainder+1
      coord(1) =  coord(1) + int(share)
      me%Part_Last(iPart-1) = tem_idofcoord(coord, offset = firstID)
      coord(1) = coord(1) + 1
      me%Part_First(iPart) = tem_idofcoord(coord, offset = firstID)
    end do

    ! The remaining elements get exactly the share elements:
    do iPart=remainder+2,nParts
      coord(1) = coord(1) + int(share) - 1
      me%Part_Last(iPart-1) = tem_idofcoord(coord, offset = firstID)
      coord(1) = coord(1) + 1
      me%Part_First(iPart) = tem_idofcoord(coord, offset = firstID)
    end do

    ! The last partition ends always with the lastID
    me%Part_Last(nParts) = lastID

    ! Local data:
    if (myPart < remainder) then
      me%nElems = int(share+1)
      me%elemOffset = int(myPart, kind=long_k) * (share+1_long_k)
    else
      me%nElems = int(share)
      me%elemOffset = (int(myPart, kind=long_k) * share) &
        &           + int(remainder, kind=long_k)
    end if
    ! All elements have (periodic) boundaries
    me%Property(1)%nElems = me%nElems
    me%Property(1)%offset = me%elemOffset
    allocate(me%Property(1)%ElemID(me%Property(1)%nElems))
    ! Please note, that the tem_bc_prop_module will set boundary conditions
    ! based on this label accordingly!
    me%global%Property(1)%label = trim(bclabel)
    me%global%Property(1)%bitpos = prp_hasBnd
    me%global%Property(1)%nElems = me%Property(1)%nElems

    allocate(me%treeID(me%nElems))
    allocate(me%ElemPropertyBits(me%nElems))

    ! Only has Boundary Property:
    me%ElemPropertyBits = ibset(0_long_k, prp_hasBnd)

    ! Filling the treeIDs:
    do iElem = 1, me%nElems
      me%Property(1)%ElemID(iElem) = iElem
      ! We can only have 2**20 elements per dimension, so in this case where
      ! we create a line, we won't exceed the integer limit. Thus we safely can
      ! cast the long integer elemOffset to a normal integer.
      coord(1) = int(me%elemOffset) + iElem - 1
      me%treeID(iElem) = tem_idofcoord(coord, offset=firstID)
    end do

  end subroutine generate_treelm_elements
  ! ************************************************************************ !


  ! ************************************************************************ !
  !> summary: Generate a predefined single element
  !> This serves as an simple grid generation for performance or scaling
  !! analysis without being obliged to use Seeder.
  !! You have to specify the generic grid parameters in the lua file instead of
  !! the mesh folder
  !!```lua
  !! mesh = { predefined='single',
  !!          origin = {0.,0.,0.},
  !!          length = 10.}
  !!```
  !! You have to specify the shape 'single', a bounding box origin and its
  !! length. The generated element will have treeID 1 and is periodic in
  !! all directions.
  !! The given length will be the length of the element, and the bounding
  !! cube will be adapted accordingly, to map this to treeID 1.
  !!
  subroutine generate_treelm_single( me, origin, length, myPart, &
    &                                nParts, comm )
    ! -------------------------------------------------------------------- !
    !> Mesh to generate
    type(treelmesh_type), intent(out) :: me
    !> Corner of the cube
    real(kind=rk), intent(in) :: origin(3)
    !> Length of cube
    real(kind=rk), intent(in) :: length
    !> Partition of the caller (starts with 0)
    integer, intent(in) :: myPart
    !> Number of partitions
    integer, intent(in) :: nParts
    !> communicator to be used
    integer, intent(in) :: comm
    ! -------------------------------------------------------------------- !
    integer(kind=long_k) :: firstID, lastID
    integer(kind=long_k) :: share
    integer :: remainder
    integer :: iPart
    ! -------------------------------------------------------------------- !

    me%global%nParts = nParts
    me%global%myPart = myPart
    me%global%comm = comm

    me%global%origin = origin
    me%global%BoundingCubeLength = length
    me%global%minLevel = 1
    me%global%maxLevel = 1
    me%global%label = 'Generic_Single'
    me%global%predefined = 'single'
    write(me%global%comment,'(a15,i7,a16,i2,a1)')         &
      &     'Generated with ', nParts, ' Parts on Level 1.'
    me%global%dirname = './'

    ! Boundary property to define periodic boundary
    me%global%nProperties = 1
    if (associated(me%global%property)) deallocate(me%global%property)
    if (associated(me%property)) deallocate(me%property)
    allocate(me%global%Property(me%global%nProperties))
    allocate(me%Property(me%global%nProperties))

    allocate(me%Part_First(nParts))
    allocate(me%Part_Last(nParts))

    ! Compute the treeIDs of the mesh:
    firstID = 1_long_k
    lastID = firstID

    ! Total number of elements in this mesh
    me%global%nElems = 1_long_k

    share = me%global%nElems / int(nParts, kind=long_k)
    remainder = int(mod(me%global%nElems, int(nParts, kind=long_k)))

    ! The first partition starts always with the firstID
    me%Part_First(1) = firstID
    me%Part_Last(1) = lastID

    ! There is only one element, all partitions but the first one are empty.
    do iPart=2,nParts
      me%Part_First(iPart) = lastID
      me%Part_Last(iPart) = lastID - 1
    end do

    ! Local data:
    if (myPart == 0) then
      me%nElems = 1
      me%elemOffset = 0_long_k
    else
      me%nElems = 0
      me%elemOffset = 1_long_k
    end if

    ! All elements have (periodic) boundaries
    me%Property(1)%nElems = me%nElems
    me%Property(1)%offset = me%elemOffset
    allocate(me%Property(1)%ElemID(me%Property(1)%nElems))
    ! Please note, that the tem_bc_prop_module will set boundary conditions
    ! based on this label accordingly!
    me%global%Property(1)%label = 'internal 0D BC'
    me%global%Property(1)%bitpos = prp_hasBnd
    me%global%Property(1)%nElems = me%Property(1)%nElems

    allocate(me%treeID(me%nElems))
    allocate(me%ElemPropertyBits(me%nElems))

    ! No Properties:
    me%ElemPropertyBits = ibset(0_long_k, prp_hasBnd)

    ! There is only one element:
    me%treeID = firstID

  end subroutine generate_treelm_single
  ! ************************************************************************ !


  ! ************************************************************************ !
  !> This subroutine creates a simple mesh based on a treeID list
  !! requires the setting of me%global before hand
  !! Works on the communicator specified in me%global%comm
  !!
  subroutine treelmesh_fromList(me, tIDlist, comm, dirname)
    ! -------------------------------------------------------------------- !
    !> Structure to load the mesh to
    type(treelmesh_type), intent(inout) :: me
    !> treeID list to generate the global list in the mesh from
    integer(kind=long_k), intent(in) :: tIDlist(:)
    !> MPI Communicator to use, defaults to the one in me%global%comm if not
    !! specified
    integer, intent(in), optional :: comm
    !> Directory to store the mesh in. Is taken to be me%global%dirname if not
    !! specified
    character(len=*), intent(in), optional :: dirname
    ! -------------------------------------------------------------------- !
    integer :: commLoc
    integer(kind=long_k) :: offset, nElems
    integer :: iError
    integer :: nElemsList
    ! mpi variables
    integer :: commsize, rank
    ! -------------------------------------------------------------------- !

    if (present(dirname)) then
      me%global%dirname = dirname
    end if

    if (present(comm)) then
      commLoc = comm
    else
      commLoc = me%global%comm
    end if

    ! determine number of treeIDs in list
    nElemsList  = size(tIDlist)
    ! Get size of communicator and my rank in it
    call mpi_comm_size(commLoc, commsize, iError)
    call mpi_comm_rank(commLoc, rank, iError)
    me%global%nParts = commsize
    me%global%myPart = rank

    if (associated(me%property)) deallocate(me%property)
    allocate(me%Property(me%global%nProperties))
    allocate(me%treeID(nElemsList))
    allocate(me%ElemPropertyBits(nElemsList))
    allocate(me%Part_First(commsize))
    allocate(me%Part_Last(commsize))
    ! Assign treeID list to the new mesh. Sort first

    offset = 0
    nElems  = int(nElemsList, kind=long_k)
    ! Determine offset for each process (get nElems of ranks before)
    call MPI_Exscan( nElems, offset, 1, MPI_INTEGER8, MPI_SUM, commLoc, iError )
    me%elemOffset = offset
    me%nElems = nElemsList
    ! Last process offset + its nElems is total number of elements
    nElems = offset + nElems
    ! Distribute among new communicator
    call MPI_Bcast( nElems, 1, MPI_INTEGER8, commsize-1, commLoc, iError )
    ! Store in global number of elements
    me%global%nElems = nElems
    ! Copy treeID list to treeID in the mesh
    me%treeID = tIDlist

    ! Reset properties. Not needed
    me%ElemPropertyBits = 0
    ! Exchange Part_First and Part_Last
    call mpi_allgather( me%treeID(1),                   &
      &                 1, MPI_INTEGER8,                &
      &                 me%Part_First, 1, MPI_INTEGER8, &
      &                 commLoc, iError                 )

    call mpi_allgather( me%treeID(nElemsList),         &
      &                 1, MPI_INTEGER8,               &
      &                 me%Part_Last, 1, MPI_INTEGER8, &
      &                 commLoc, iError                )

  end subroutine treelmesh_fromList
  ! ************************************************************************ !


  ! ************************************************************************ !
  !> write mesh information
  subroutine dump_meshHeader( me )
    ! -------------------------------------------------------------------- !
    !> fluid tree from mesh
    type( treelmesh_type ), intent(in) :: me
    ! -------------------------------------------------------------------- !
    integer :: iProp
    ! -------------------------------------------------------------------- !
    write(logunit(1),*)  'Got a mesh with following properties:'
    write(logUnit(1), "('                 Mesh name: ',A)") &
      &               trim(me%global%label)
    if (trim(me%global%predefined) /= '') then
      write(logunit(1), "(' This is a predefined mesh: ', A)") &
        &               trim(me%global%predefined)
    end if
    write(logunit(1), "('  Total number of elements: ', I0)") me%global%nElems
    write(logunit(1), "('  Local number of elements: ', I0)") me%nElems
    write(logUnit(1), "('      Number of properties: ', I0)") &
      &               me%global%nProperties
    write(logunit(1), "('      Number of partitions: ', I0)") me%global%nParts
    write(logunit(1), "('             Minimum level: ', I0)") me%global%minlevel
    write(logunit(1), "('             Maximum level: ', I0)") me%global%maxlevel
    write(logunit(1), "('       Bounding Box Origin: ', 3F10.4 )") &
      &               me%global%origin
    write(logunit(1), "('       Bounding Box Length: ', F10.4 )") &
      &               real(me%global%BoundingCubeLength)
    write(logUnit(1), "('     Coarsest element size: ', F10.4 )") &
      &               real( me%global%BoundingCubeLength /      &
      &               real( 2**me%global%minlevel, kind=rk ) )
    write(logUnit(1), "('       Finest element size: ', F10.4 )") &
      &               real( me%global%BoundingCubeLength /      &
      &               real( 2**me%global%maxlevel, kind=rk ) )
    do iProp=1,me%global%nProperties
      if ( me%global%Property(iProp)%BitPos == prp_hasBnd ) then
        write(logUnit(1), "('  Number of boundary elements: ', I0 )") &
          &               me%Property(iProp)%nElems
      end if
      if ( me%global%Property(iProp)%BitPos == prp_hasQVal ) then
        write(logUnit(1), "('      Number of qVal elements: ', I0 )") &
          &               me%Property(iProp)%nElems
      end if
    end do
  end subroutine dump_meshHeader
  ! ************************************************************************ !


  ! ************************************************************************ !
  subroutine tem_load_weights( me, weights, success )
    ! -------------------------------------------------------------------- !
    type(treelmesh_type), intent(in) :: me
    real(kind=rk), intent(out) :: weights(me%nElems)
    logical, intent(out) :: success
    ! -------------------------------------------------------------------- !
    integer(kind=MPI_OFFSET_KIND) :: displacement
    integer :: fh, ftype, iError, typesize
    integer :: iostatus( MPI_STATUS_SIZE )
    integer :: level
    integer :: iElem
    logical :: ex
    character(len=PathLen) :: filename
    character(len=4) :: EndianSuffix
    ! -------------------------------------------------------------------- !

    EndianSuffix = tem_create_endianSuffix()
    if (trim(me%weights_file) /= '') then
      filename = trim(me%weights_file)//EndianSuffix
      ! check if a corresponding weight file exists
      if (me%global%myPart == 0) then
        inquire(file=trim(filename), exist=ex)
      end if

      call MPI_Bcast(ex, 1, MPI_LOGICAL, 0, me%global%comm, iError)
    else
      ex = .false.
    end if

    if (ex) then
      ! Found a weights file, which is used to read a weight for each
      ! element.
      write(logUnit(3),*) 'Loading Weights from file: '//trim(filename)
      ! Open the binary file for MPI I/O (Write)
      call MPI_File_open( me%global%comm, trim(filename),            &
        &                 MPI_MODE_RDONLY, MPI_INFO_NULL, fh, iError )
      call check_mpi_error(iError,'file_open in load_weights')

      ! Create a MPI Subarray  as ftype for file view
      call MPI_Type_contiguous( me%nElems, rk_mpi , ftype, iError )
      call check_mpi_error(iError,'type ftype in load_weights')
      call MPI_Type_commit( ftype, iError )
      call check_mpi_error(iError,'commit ftype in load_weights')

      !get size of etype
      call MPI_Type_size(rk_mpi, typesize, iError )
      call check_mpi_error(iError,'typesize in load_weights')

      ! calculate displacement
      displacement = me%elemoffset * typesize *1_MPI_OFFSET_KIND

      ! Set the view for each process on the file above
      call MPI_File_set_view( fh, displacement, rk_mpi, ftype, &
        &                     "native", MPI_INFO_NULL, iError  )
      call check_mpi_error(iError,'set_view in load_weights')

      ! Read data from the file
      call MPI_File_read_all( fh, weights, me%nElems, rk_mpi, iostatus, iError )
      call check_mpi_error(iError,'read_all in load_weights')

      !Free the MPI_Datatypes which were created and close the file
      call MPI_Type_free(ftype, iError)
      call check_mpi_error(iError,'free ftype in load_weights')
      call MPI_File_close(fh,    iError)
      call check_mpi_error(iError,'close file in load_weights')
      success = .true.
    else if (maxval(me%levelweight) > 0.0_rk) then
      write(logunit(3),*) 'Using levelwise weights.'
      do iElem=1,me%nElems
        level = tem_levelOf( me%treeID(iElem) )
        weights(iElem) = me%levelWeight(level)
      end do
      success = .true.
    else
      write(logunit(3),*) 'NO BALANCING!'
      success = .false.
    end if
    write(logUnit(3),*) 'Done loading weights.'

  end subroutine tem_load_weights
  ! ************************************************************************ !

  ! ************************************************************************ !
  subroutine tem_dump_weights( me, filename, weights )
    ! -------------------------------------------------------------------- !
    type(treelmesh_type), intent(in) :: me
    !> Weights file name
    character(len=*), intent(in) :: filename
    real(kind=rk), intent(in)    :: weights(me%nElems)
    ! -------------------------------------------------------------------- !
    integer                       :: nElems
    integer                       :: comm
    integer(kind=long_k)          :: offset, globElems
    integer(kind=long_k)          :: filesize
    integer(kind=MPI_OFFSET_KIND) :: displacement
    integer                       :: fh, ftype, iError
    integer                       :: iostatus( MPI_STATUS_SIZE )
    integer                       :: typesize
    ! -------------------------------------------------------------------- !

    nElems = me%nElems
    comm = me%global%comm
    offset = me%elemOffset
    globElems = me%global%nElems

    ! Found a weights file, which is used to read a weight for each
    ! element.
    write(logUnit(1),*) 'Dumping Weights to file: '//trim(filename)
    ! Open the binary file for MPI I/O (Write)
    call MPI_File_open( comm, trim(filename),                 &
      &                 ior(MPI_MODE_WRONLY,MPI_MODE_CREATE), &
      &                 MPI_INFO_NULL, fh, iError             )
    call check_mpi_error(iError,'file_open in dump_weights')

    ! Create a MPI Subarray  as ftype for file view
    call MPI_Type_contiguous( nElems, rk_mpi , ftype, iError )
    call check_mpi_error(iError,'type ftype in dump_weights')
    call MPI_Type_commit( ftype, iError )
    call check_mpi_error(iError,'commit ftype in dump_weights')

    !get size of etype
    call MPI_Type_size(rk_mpi, typesize, iError )
    call check_mpi_error(iError,'typesize in dump_weights')

    ! calculate displacement
    displacement= offset * typesize * 1_MPI_OFFSET_KIND

    ! set filesize
    filesize = globElems * typesize
    call MPI_File_set_size( fh, filesize, iError )

    ! Set the view for each process on the file above
    call MPI_File_set_view( fh, displacement, rk_mpi, ftype, &
      &                     "native", MPI_INFO_NULL, iError  )
    call check_mpi_error(iError,'set_view in dump_weights')

    ! Read data from the file
    call MPI_File_write_all( fh, weights, nElems, rk_mpi, iostatus, iError )
    call check_mpi_error(iError,'read_all in dump_weights')

    !Free the MPI_Datatypes which were created and close the file
    call MPI_Type_free(ftype, iError)
    call check_mpi_error(iError,'free ftype in dump_weights')
    call MPI_File_close(fh,    iError)
    call check_mpi_error(iError,'close file in dump_weights')
    write(logUnit(1),*) 'Done dumping weights.'

  end subroutine tem_dump_weights
  ! ************************************************************************ !


  ! ************************************************************************ !
  !> Exchange mesh elements with other process
  !! The following data within treelmesh_type is updated in this routine:
  !!    treeID, ElemPropertyBits, Part_First, Part_Last
  !! Data pathList is deallocated.
  subroutine exchange_elements( me, sparta )
    ! -------------------------------------------------------------------- !
    !> Tree Mesh type
    !! nElems and offset should be ready before this routine
    type( treelmesh_type ), intent(inout) :: me
    !> How many elements to exchange with other processes
    type( tem_sparta_type ), intent(in) :: sparta
    ! -------------------------------------------------------------------- !
    integer :: iError, iProp
    ! -------------------------------------------------------------------- !

    ! Exchange treeID and propertyBits ---------------------------------
    call tem_exchange_sparta( me = sparta, &
      &                       val = me%treeID, &
      &                       nComponents = 1, &
      &                       comm = me%global%comm )
    call tem_exchange_sparta( me = sparta, &
      &                       val = me%ElemPropertyBits, &
      &                       nComponents = 1, &
      &                       comm = me%global%comm )
    ! Exchange treeID and propertyBits ---------------------------------

    ! Update Part_first and Part_Last ------------------------
    call MPI_allgather( me%treeID(1), 1, long_k_mpi,  &
      &                 me%Part_First, 1, long_k_mpi, &
      &                 me%global%comm, iError        )
    call MPI_allgather( me%treeID(me%nElems), 1, long_k_mpi, &
      &                 me%Part_Last, 1, long_k_mpi,         &
      &                 me%global%comm, iError               )
    ! Update Part_first and Part_Last ------------------------

    do iProp=1,me%global%nProperties
      call gather_Property( Property = me%Property(iProp),                     &
        &                   Header   = me%global%Property(iProp),              &
        &                   BitField = me%ElemPropertyBits,                    &
        &                   comm     = me%global%comm )
    end do

    if (allocated(me%pathList)) deallocate( me%pathList  )

  end subroutine exchange_elements
  ! ************************************************************************ !

end module treelmesh_module
! **************************************************************************** !

!> \page motivation Motivation for the TreElm Library
!!
!! Mesh-based algorithms build one of the most important discretization schemes
!! for the numerical modeling of field problems.
!! They are also well suited for parallel computations, as the meshes can be
!! partitioned and the smaller problems solved concurrently.
!! Their basic idea is the discretization of the computational domain by
!! covering it fully with smaller, well defined elements.
!! For complex geometries, it is in general necessary to use an unstructured
!! mesh to discretize the computational domain.
!! In unstructured meshes each element has typically a different size and form.
!! Therefore, all information of element geometry and neighbor relations have to
!! be described explicitly.
!! This fully explicit description however, poses a severe hurdle for parallel
!! processing, as potentially any part of the mesh might be required in the
!! neighborhood of any process.
!! Due to the limited memory, providing the complete mesh information locally to
!! each partition is not an option, as the ability to process larger meshes is
!! often the motivation to switch from one to many computing units in the first
!! place.
!! Although, it is possible to move this limitation out of the solver itself,
!! Tu et al. pointed out the importance to avoid such bottlenecks throughout
!! the complete chain of tools in the simulation.
!! Thus, not only the main computational loop, but the entire simulation process
!! has to be considered, including mesh generation and post-processing.
!! The usage of hierarchical subdivisions in a fixed space resulting in and
!! octree representation offers an option to overcome the hurdle for parallel
!! processing of complex geometries with mesh-based methods.
!!
!! TreElM provides a basis for massively parallel mesh-based simulations, using
!! an octree discretization in combination with a space-filling curve.
!! The usage of octree meshes for flow simulations has long been proposed, for
!! example by Flaherty et al. in 1997.
!! By using octree meshes the main bottleneck of unstructured
!! meshes is avoided, while the flexibility of local refinement and resolution
!! of arbitrary complex geometries is maintained.
!! With the known topology of the octree, the needed information about other
!! partitions can be kept at a minimum and the connectivity can be computed
!! locally.
!! Direct usage of cubic octree elements has the advantage that they can be
!! computed efficiently in most schemes, or are even required as for example
!! in the lattice Boltzmann method.
!! \image html clumps_sm.png
!! The figure above shows an example for meshes, that are obtained by the
!! octree discretization.
!! It shows the filling of a channel with many small cylindric obstacles,
!! indicated by transparent objects.
!! The complete space covered by the octree is outlined by the blue grid lines
!! and shows the refinement towards the channel with obstacles.
!! Inside the channel, there is the actual computational domain, shown by the
!! green cubical elements.
!!
!! The basic concept in TreElM is to exploit the known hierarchical
!! topology, provided by the octree.
!! Another concept that is supported by the octree and beneficial for the
!! parallel operation is the strictly elementwise view on the mesh data.
!! As each element will be uniquely associated with a certain partition during
!! the computation, such a view enables the parallel treatment straight from
!! reading or creating the mesh onwards.
!! TreElM therefore maintains an elemental view on the mesh data in the data layout
!! and the operations on the data.
!!
!! \section toolchain The Simulation Tool Chain
!!
!! The tools in the framework range from the mesh generator *Seeder*
!! over a finite volume solver for compressible and a
!! lattice Boltzmann solver for incompressible flow to the post-processing tool
!! *Harvester*, that produces visualization files from simulation data
!! attached to the mesh.
!! This interaction is sketched in the following Figure.
!! \image html apes_schematic.png
!! Since neither input nor output (I/O) can be avoided to the largest part, even
!! with an integrated workflow, data structures are designed to be suitable for
!! efficient parallel reading and writing.
!!
!! \subsection see_also See Also
!! - \ref tem_distributed_octree "The distributed Octree"
!! - \ref tem_node_identification "Node identification in the tree"
!!
!! [treeid]: @ref treelmesh_module#treelmesh_type#treeid "treeID"
!!
!! \page tem_distributed_octree Distributed Octree
!! An [octree](http://en.wikipedia.org/wiki/Octree)
!! is obtained by recursive bisection in each dimension of a
!! three-dimensional cube.
!! As such it is the generalization of a binary tree into three dimensions, the
!! two dimensional case is usually referred to as quadtree.
!! We use quadtrees and two dimensional meshes for illustrations.
!! A three-dimensional picture of the actual octree meshes, obtained by this
!! method is provided in the Figure.
!! \image html clumps_sm.png
!! This illustration shows the subdivision of coarser elements into eight
!! smaller elements and the recursive repetition of this refinement where
!! necessary.
!! The outermost cube is also addressed as the universe, as it contains the
!! complete representable space of the octree and all elements have to live
!! inside of it.
!! The basic idea in TreElM for a scalable parallel deployment of the octree
!! discretization is a simple serialized tree representation and sketched in
!! the following figure.
!! \image html tree_structure_sm.png
!! The serialization turns the hierarchical structure of the tree, where each
!! node has several links to other nodes into a serial ordered list of one node
!! after the other.
!! This is achieved by a [space-filling curve](http://en.wikipedia.org/wiki/Space-filling_curve).
!! Space filling curves are designed to fill the three dimensional space
!! completely with a one-dimensional one.
!! They thus map a one-dimensional unit interval onto a higher dimensional unit
!! interval.
!! This feature can be exploited to define a serialization of the three
!! dimensional space.
!! They are also defined recursively and complete coverage of the continuous
!! space is achieved after infinite iterations.
!!
!! Finite iterations of the space filling curve definition can be interpreted as
!! discrete curves and can be used to describe the ordering for finite intervals
!! in the higher dimensional space.
!! This recursive definition of the space filling curve ordering nicely fits the
!! recursive definition of the spatial subdivisions by the octree.
!! We can thus interpret the space filling curve as a rule how the children of a
!! node in the octree should be numbered.
!! Consecutively numbering all nodes on each refinement level according to the
!! space filling curve results in a breadth first numbering of the complete
!! tree.
!! As these numbers identify the elements in the full tree, we refer to them as
!! [treeid]s.
!! Their properties are described in the following and the relation
!! between mesh and tree is illustrated with the help of a quadtree in the Figure
!! above.
!!
!! \section partitioning Partitioning of the Tree
!! With this numbering scheme for all nodes in the full tree it is possible to
!! deduce the geometrical position of each node purely by its assigned number.
!! Just the [treeid]s of the leaf
!! nodes, that are actually part of the
!! computational domain are stored with the described ordering to represent the
!! mesh.
!! Thus, we obtain a sparse octree, that contains only the relevant leaf nodes
!! for the computational nodes.
!! In the Figure of the 3D mesh, these are highlighted as green cubes.
!! None of the other tree nodes, including the parents of the elements in the
!! computational domain, are actually stored.
!! Due to the ordering, the serialized list of elements is efficiently
!! searchable.
!! The one dimensional list of elements can be easily partitioned by splitting
!! it into individual chunks.
!! This design is sketched in the Figure above for a two
!! dimensional mesh as shown on the left side of the figure.
!! There is an obstacle in the middle of the mesh, covering the four central
!! fine elements.
!! The mesh is refined towards this obstacle.
!! A line is used to indicate the connection between consecutive elements
!! according to the ordering with the space filling curve.
!! On the right side the according quadtree is shown and the partitioning into
!! 6 parts is indicated by regions with different colors.
!! Thus, we obtain a simple ordered linear representation of the sparse octree.
!! Each actually existing element in the mesh is numbered by position in
!! addition to its [treeid] and the resulting serialized tree representation is
!! given at the bottom of the Figure in the form of a simple [treeid] list.
!!
!! \subsection see_also See also
!! - \ref treelm_fileformat Treelm File Format
!! - \ref tem_balance_module "Load Balancing"
!!
!! [treeid]: @ref treelmesh_module#treelmesh_type#treeid "treeID"
!!
!! \page treelm_fileformat Treelm File Format
!! As already introduced in the previous sections, the
!! \ref tem_distributed_octree "octree mesh" is represented by the sparse
!! storage of the leaf nodes in the computational domain.
!! The storage of the mesh on disk is basically the sequence of [treeid]s
!! building the computational domain.
!! However, this is in general insufficient to fully describe the mesh and an
!! additional elemental information is introduced to attach further data to the
!! elements.
!! TreElM uses 64-bit signed integers to encode the [treeid]s for
!! greatest portability.
!! In addition, a second 64-bit signed integer is used to build a bit-mask for
!! properties like boundary conditions, that might be attached to each element.
!! This data is written in native binary format to disk, resulting in a file
!! with 16 bytes per element.
!! Due to this simple format, the data might be accessed on an elemental basis
!! with Fortran's direct IO, it can be easily converted between big and little
!! endian representations and allows fast reading and writing.
!!
!! This binary data is accompanied by a header file with descriptive information
!! about the mesh, like the total number of elements, the physical extent and
!! origin of the universal cube and descriptions for the attached properties.
!! The data in this header is provided in the form of a [Lua](http://www.lua.org/)
!! script, to allow for a flexible handling of additional data and future evolution
!! of the format.
!! With the help of the \ref aotus_module "Aotus library" data can be easily
!! retrieved from and written to such Lua scripts.
!!
!! [treeid]: @ref treelmesh_module#treelmesh_type#treeid "treeID"
!!
!!
!! \page implementation Implementation in the APES Framework
!! In distributed computations it is important to avoid any kind of bottleneck,
!! as the local computing resources in each computing unit are
!! usually quite limited in comparison to the total computing power.
!! For large scale computations that aim to use the full system for their
!! simulation, it therefore is essential that the complete processing chain is
!! distributable.
!! Otherwise any step from the mesh generation over the simulation to the
!! post-processing might prohibit a successful simulation.
!! The *TreElM* library builds the basis for a completely
!! integrated framework that allows the distributed processing of all necessary
!! simulation steps.
!! To enable a flexible and convenient user configuration, the \ref aotus_module
!! "Aotus library"
!! is not only used in *TreElM* header files, but also exposed to the
!! calling applications.
!! The layout of the overall framework with the central *TreElM* library is
!! shown by the schematic figure in the \ref motivation "Motivation".
!!
!! The solvers read the mesh data, using *TreElM* functions and process the
!! mesh into a solver specific data structure with attached simulation data.
!! Results are written in the form of restart files, that follow the elemental
!! design of the *TreElM* mesh by writing the elemental solution following the
!! same ordering to disk.
!! With the help of tracking objects it is possible to write subsets of the mesh
!! with attached simulation data in arbitrary intervals to disk.
!! The *TreElM* library not only provides the means to identify the halos to
!! be exchanged between partitions but also provides various communication
!! pattern that can be used for the actual exchange.
!! Due to this encapsulation of the communication, the communication layout can
!! be easily exchanged.
!! It is even possible to replace the complete parallel paradigm in the
!! deployment with only few changes to the code.
!!
!! Aside from these basic functionalities the library provides several auxiliary
!! routines, that provide for example logging and debugging facilities.
!! Finally the post processing tool *Harvester* is used to visualize the
!! results of the simulation.
!! Since *Harvester* is a stand alone tool, it can be deployed on specialized
!! machines, that are more suited for visualization.
!! Additionally, the separated design ensures the independence of the solvers
!! from any third party libraries, that are typically required for visualization
!! outputs.
!!
!! [treeid]: @ref treelmesh_module#treelmesh_type#treeid "treeID"
!!
!!
!! \page tem_requirements Requirements
!!
!! Below are the most important requirements for getting, compiling and running
!! solvers developed with the TreElm framework.
!!
!! Compilation
!! -----------
!! - [Mercurial](mercurial.selenic.com): Getting sources
!! - [Python](www.python.org): the build tool waf relies on python and most of
!!           the scripts inside the package are written with Python
!! - Compiler: the vast source files are written in Fortran 95/2003.
!!            We recommend [GNU Gfortran 4.6+](http://gcc.gnu.org/wiki/GFortran)
!!            and [Intel 12.1+](http://software.intel.com/en-us/articles/
!!                              intel-compilers/)
!! - MPI 2.0 enabled environment
!!           We recommend [OpenMPI 1.5+](http://www.open-mpi.org/)
!!
!! Post-Processing and Visualization
!! ---------------------------------
!! - [Gnuplot](http://www.gnuplot.info/)
!! - [Paraview](http://www.paraview.org/)
!! - [Tecplot](http://www.tecplot.com)
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Copyright 2012 German Research School for Simulation Sciences GmbH, Aachen.
! Copyright 2013 University of Siegen, Siegen.
! All rights reserved.
!
! 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 GERMAN RESEARCH SCHOOL FOR SIMULATION SCIENCES
! GMBH, AACHEN “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 GERMAN RESEARCH SCHOOL
! FOR SIMULATION SCIENCES GMBH, AACHEN 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 views and conclusions contained in the software and documentation are
! those of the authors and should not be interpreted as representing official
! policies, either expressed or implied, of German Research School for
! Simulation Sciences GmbH, Aachen.
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
