! Copyright (c) 2012-2014, 2022 Harald Klimach <harald.klimach@dlr.de>
! Copyright (c) 2012-2013, 2015, 2019 Kannan Masilamani <kannan.masilamani@uni-siegen.de>
! Copyright (c) 2012 Jiaxing Qi <jiaxing.qi@uni-siegen.de>
! Copyright (c) 2016 Tobias Girresser <tobias.girresser@student.uni-siegen.de>
! Copyright (c) 2016 Peter Vitt <peter.vitt2@uni-siegen.de>
! Copyright (c) 2018 Daniel Fleischer <daniel.fleischer@student.uni-siegen.de>
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions are met:
!
! 1. Redistributions of source code must retain the above copyright notice, this
! list of conditions and the following disclaimer.
!
! 2. Redistributions in binary form must reproduce the above copyright notice,
! this list of conditions and the following disclaimer in the documentation
! and/or other materials provided with the distribution.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
! AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
! IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
! DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
! FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
! SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
! OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
! OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
! Copyright (c) 2012-2013 Simon Zimny <s.zimny@grs-sim.de>
! Copyright (c) 2012 Manuel Hasert <m.hasert@grs-sim.de>
! Copyright (c) 2012-2016 Harald Klimach <harald.klimach@uni-siegen.de>
! Copyright (c) 2012, 2015-2016 Kannan Masilamani <kannan.masilamani@uni-siegen.de>
! Copyright (c) 2013 Daniel Harlacher <d.harlacher@grs-sim.de>
! Copyright (c) 2014 Kartik Jain <kartik.jain@uni-siegen.de>
! Copyright (c) 2014 Verena Krupp <verena.krupp@uni-siegen.de>
! Copyright (c) 2015-2017 Jiaxing Qi <jiaxing.qi@uni-siegen.de>
! Copyright (c) 2015-2016 Peter Vitt <peter.vitt2@uni-siegen.de>
! Copyright (c) 2016 Daniel Fleischer <daniel.fleischer@student.uni-siegen.de>
! Copyright (c) 2016 Tobias Schneider <tobias1.schneider@student.uni-siegen.de>
! Copyright (c) 2017 Daniel Petró <daniel.petro@student.uni-siegen.de>
!
! Parts of this file were written by Harald Klimach, Simon Zimny and Manuel
! Hasert for German Research School for Simulation Sciences GmbH.
!
! Parts of this file were written by Harald Klimach, Kannan Masilamani,
! Daniel Harlacher, Kartik Jain, Verena Krupp, Jiaxing Qi, Peter Vitt,
! Daniel Fleischer, Tobias Girresser and Daniel Petró for University Siegen.
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions are met:
!
! 1. Redistributions of source code must retain the above copyright notice, this
! list of conditions and the following disclaimer.
!
! 2. Redistributions in binary form must reproduce the above copyright notice,
! this list of conditions and the following disclaimer in the documentation
! and/or other materials provided with the distribution.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
! AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
! IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
! DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
! FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
! SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
! OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
! OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

! This file contains the source code for growing and dynamic arrays.
! This is used for arrays of primitives (int, long_int, real, ...) as well as
! for arrays of derived datatypes (tem_variable_type,...).
!
! To use these macros include the following to your source file.
!
! Smart growing array (GA) for ?tstring?
! Growing Arrays:
!
! declaration
!
!
! implementation
!

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


!
!------------------------------------------------------------------------------
!
! dynamic Arrays:
!
! declaration
!
!
! implementation
!
! ******************************************************************************!
!> This module implements the description of a node in the tree that is going to
!! be created by Seeder.
module sdr_node_module
  use env_module,            only: rk, long_k, minLength, zeroLength
  use tem_param_module,      only: qQQQ
  use tem_grow_array_module, only: grw_int2dArray_type, grw_longArray_type,    &
    &                              grw_intArray_type, &
    &                              init, append, placeAt, truncate
  use tem_dyn_array_module,  only: dyn_longArray_type, init, append, truncate
  use tem_color_prop_module, only: colors_per_char
  use tem_geometry_module,   only: tem_eligibleChildren
  use tem_mergesort_module,  only: mrgrnk

  use sdr_geometry_module,   only: sdr_geometry_type

  implicit none

  private

  public :: sdr_node_type

  public :: init, append, truncate, destroy, empty, placeAt
  public :: sdr_wetNeighborsFace
  public :: sdr_mark_floodNode
  public :: sdr_set_nodeProp_bit
  public :: sdr_clear_nodeProp_bit
  public :: sdr_nodeProp_btest
  public :: sdr_nodeColors
  public :: sdr_bitfieldColors
  public :: sdr_color_log2char
  public :: sdr_inHeritBnd_eligibleChildren
  public :: sdr_intersectObjPos_type
  public :: grw_intersectObjPosArray_type
  public :: sdr_append_childIntersectedObject

  ! Some parameters to encode node properties in the PropertyBits bitmask

  !> This bit is set, if the node is a leaf, that is there are no further
  !! descendants.
  integer, parameter, public :: isLeaf_bit = 1

  !> If this bit is set, the node intersects some computational domain boundary.
  integer, parameter, public :: IntersectsBoundary_bit = 2

  !> If any of 26 direct neighbor is intersected by boundary
  integer, parameter, public :: hasBoundary_bit = 3

  !> This bit is set, if the node is a fluidifyable that can be
  !! fluidify in solver
  integer, parameter, public :: isFluidifyable_bit = 4

  !> If the node is flooded by any color, this bit is set.
  !! It is useful to identify elements, that belong to the computational
  !! domain.
  integer, parameter, public :: isFlooded_bit = 5

  !> Indicating if the node is colored with an non-none color.
  !!
  !! This is used to filter none colors when transferring the information from
  !! the prototree to the treelm data structure.
  integer, parameter, public :: isColored_bit = 6

  !> Indicates if this node is one that is to be included in the final mesh.
  !!
  !! It is needed to distinguish such nodes, which are refined further for
  !! subelement resolution of boundaries, from ordinary virtual nodes.
  integer, parameter, public :: isTarget_bit = 7

  !> This bit is set, if solidification is not allowed for the node
  !!
  integer, parameter, public :: isNoSolidification_bit = 0

  !> This data type describes first and last intersected object
  !! of a node in growing array of intersected_objects.
  !!
  !! Childs intersected objects are filled in parents intersected
  !! objects position until childs objects fit.
  !! Children are filled with sorted fashsion from child with
  !! max nObjects to least
  type sdr_intersectObjPos_type
    !> First position of intersected object in 1st list
    integer :: first

    !> Last position of intersected object in 1st list
    integer :: last
  end type sdr_intersectObjPos_type


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

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

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

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

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

  !> insert an element at a given position
  interface placeat
    module procedure placeat_ga_intersectobjpos
    module procedure placeat_ga_intersectobjpos_vec
  end interface

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

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

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

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

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

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

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

  !> insert an element at a given position
  interface placeat
    module procedure placeat_ga_grwint
    module procedure placeat_ga_grwint_vec
  end interface

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

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



  !> This data type describes a node in the tree to be created.
  type sdr_node_type

    !> Number of nodes in this list.
    integer :: nNodes = 0

    !> Number of integers to encode the properties.
    !! This depends on the number of colors and defines the size of the
    !! the first dimension (width) in the grw_int2dArray propertyBits.
    integer :: propLength = 1

    !> Auxilary value to provide the number of bytes available per integer
    !! in the propertybits.
    integer :: bytes_per_int

    !> Indicator for the starting bit of the last byte in the propertybits.
    integer :: lastbyte_pos

    !> Counter for the number of colors in the mesh.
    integer :: nColors

    !> Id of the 'none' color, this color will not be considered in the final
    !! color information for treelm.
    !!
    !! If there is no 'none' color this will be -1, otherwise it is the id
    !! of the color with the label 'none' in the unique list of color labels.
    integer :: none_color_id

    !> Number of characters needed to hold all colors.
    integer :: nColorChars

    !> Number of levels, that need to be further refined for colors with
    !! subelement resolution.
    !!
    !! If there are no colors with subelement resolution at all, this setting
    !! is 0.
    integer :: subelement_resolution

    !> True if distance refine objects are created for
    !! boundary attribute.
    !!
    !! If there is no distance refine object then no need to append
    !! distance_first and distance_last for every node
    logical :: distanceRefine = .false.

    !> List of treeIDs for all the nodes in the built tree.
    !!
    !! We are exploiting here the fact, that new nodes are created in order
    !! level by level. Therefore, it is not necessary to use an
    !! additional index list to maintain the ranking of the entries, as it is
    !! done in the dyn_arrays.
    type(dyn_longArray_type) :: treeID

    !> The property bits can be used to attach various information to the
    !! element.
    !! They can span multiple integers (width of the 2d growing array) and all
    !! bytes except the very last one are used to encode colors (1 byte per
    !! color).
    !!
    !! In the last byte other properties are encoded (Bits 0-7) of last byte.
    !! (The last byte is given by
    !!  ibits(PropertyBits%val(propLength,iNode), pos=bit_size(int)-8, len=8))
    !! The following bits are used:
    !! Bit 0: NoSolidification is set, if solidification is not allowed for node
    !! Bit 1: isLeaf Bit is set, if there are no further descendants
    !! Bit 2: IntersectsBoundaries, is set, if the element intersects boundaries
    !! Bit 3: IntersectsAll Set, if all objects are seperate intersected_object
    !! Bit 4: isFluidifyable is set, if the node is fluidifyable.
    !! Bit 5: isFlooded is set, if node is flooded by any color.
    !! Bit 6: isColored is set, if there is another than the none color in the
    !!        node.
    !! Bit 7: isTarget is set, if this node should be in the final mesh, even
    !!        so it has further children.
    !!
    !! To encode a color, we need a representation for each of the 6 sides and
    !! one for the volume, resulting in a total of 7 bits for each color.
    !! Thus, a complete Byte is used for each color. The first color is stored
    !! in the first Byte of the first integer and the other follow
    !! consecutively.
    !! Ordering of the sides as given by the definition described for the
    !! linkPos component.
    !! We also need to store intersections on a per color basis, as we want to
    !! get the flooding stopped by differently colored boundaries.
    !! The layout of the bits for each color ar therefore:
    !! Bit   0: intersects a relevant geometry (same color or geom. color = -1)
    !! Bit   1: node is flooded with this color
    !! Bit 2-7: colored wetness of the 6 node sides
    type(grw_int2dArray_type) :: PropertyBits

    !> Store which of 26 directions of node has boundary.
    !! Uses first 26 bits of the integer with each bit
    !! corresponds to treelm direction given by qOffset in tem_param_module
    type(grw_intArray_type) :: hasBndBits

    !> Minimum boundary ID intersect by the node
    !!
    !! Used to set boundaryID based on thuming rule when a node
    !! is intersected by more than one boundary object
    type(grw_intArray_type) :: minBCID

    !> List of all objects intersecting this node.
    !!
    !! If children are created for this node, the number of children
    !! with intersected objects thats fits in to the memory of its parent
    !! starting from child with highest number of intersected object
    !! will be copied to the memory position of its parent to save memory.
    type(grw_intArray_type) :: intersected_object

    !> First and last position of this node's intersected objects in
    !! growing array of intersected_object
    type(grw_intersectObjPosArray_type) :: userObjPos

    !> List of all sphere_distance objects intersecting this node
    type(grw_intArray_type) :: intersected_distance

    !> First and last position of this node's intersected distance refine
    !! objects in growing array of intersected_distance
    type(grw_intersectObjPosArray_type) :: distObjPos

    !> Memory unused by children from parent intersected object list
    integer :: memLeft_userObj

    !> Memory unused by children from parent intersected distance list
    integer :: memLeft_distObj

    !> minLevel for each node.
    !! Used in refineleaf and inherit distance refine object.
    !! If intersected object level is <= minLevel then that object
    !! is not added to intersected object list for children.
    !! This should reduce memory to inHerit lot of intersected objects
    !! to children
    type(grw_intArray_type) :: minLevel

    !> This is used to store the location of the 6 neighbors for leaf nodes and
    !! the first child for virtual nodes at linkpos(1).
    !!
    !! 1: neighbor x-1 / first child node
    !! 2: neighbor x+1
    !! 3: neighbor y-1
    !! 4: neighbor y+1
    !! 5: neighbor z-1
    !! 6: neighbor z+1
    type(grw_intArray_type) :: linkpos(6)

    !> If subelement resolution is required, the sublevel keeps track of how
    !! deep the subelement resolution is below the target element.
    !!
    !! For nodes above isTarget nodes it is set to a negative value. In the
    !! target node it is set to the subelement_resolution.
    !! On each level this counter is decreased by 1, when passing it on to the
    !! children. Boundaries, that require subelement resolution will stop in
    !! refinement, when sublevel has reached 0, or if a total level of 20 is
    !! reached.
    !! If subelement_resolution is 0, this array is not used
    !! (and not allocated).
    type(grw_intArray_type) :: sublevel

  end type sdr_node_type


  interface sdr_wetNeighborsFace
    module procedure sdr_wetNeighborsFace_single
    module procedure sdr_wetNeighborsFace_all
  end interface

  interface init
    module procedure init_node
  end interface

  interface append
    module procedure append_newNode
  end interface

  interface truncate
    module procedure truncate_node
  end interface

  interface sdr_append_childIntersectedObject
    module procedure sdr_append_childIntersectedObjectAll
    module procedure sdr_append_childIntersectedObjectGTminLevel
  end interface sdr_append_childIntersectedObject

contains

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

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

  end subroutine init_ga_intersectobjpos

  subroutine destroy_ga_intersectobjpos(me)
    type(grw_intersectobjposarray_type), intent(inout) :: me !< dynamic array to destroy

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

  end subroutine destroy_ga_intersectobjpos


  subroutine truncate_ga_intersectobjpos(me)
    !------------------------------------------------------------------------
    type(grw_intersectobjposarray_type) :: me !< array to truncate
    !------------------------------------------------------------------------
    type(sdr_intersectobjpos_type), allocatable :: tarray(:)
    !------------------------------------------------------------------------
    integer :: ii
    !------------------------------------------------------------------------

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

  end subroutine truncate_ga_intersectobjpos


  subroutine empty_ga_intersectobjpos(me)
    !------------------------------------------------------------------------
    type(grw_intersectobjposarray_type) :: me !< array to sorttruncate
    !------------------------------------------------------------------------

    me%nvals = 0

  end subroutine empty_ga_intersectobjpos

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


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

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

  end subroutine placeat_ga_intersectobjpos


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

    integer :: ub, ii

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

    ub = pos + size(val) - 1

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

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

  end subroutine placeat_ga_intersectobjpos_vec


  subroutine append_ga_intersectobjpos(me, val, length)
    type(grw_intersectobjposarray_type) :: me !< array to append the value to
    type(sdr_intersectobjpos_type), intent(in) :: val !< value to append
    !> optional length to expand the array
    integer, intent(in), optional :: length

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

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

  end subroutine append_ga_intersectobjpos

  subroutine append_ga_intersectobjpos_vec(me, val, length)
    type(grw_intersectobjposarray_type) :: me !< array to append the value to
    type(sdr_intersectobjpos_type), intent(in) :: val(:) !< values to append
    !> optional length to expand the array
    integer, intent(in), optional :: length

    integer :: lb, ub, ii

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

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

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

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

  end subroutine append_ga_intersectobjpos_vec


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

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

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

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

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

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

  end subroutine expand_ga_intersectobjpos

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

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

  end subroutine init_ga_grwint

  subroutine destroy_ga_grwint(me)
    type(grw_grwintarray_type), intent(inout) :: me !< dynamic array to destroy

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

  end subroutine destroy_ga_grwint


  subroutine truncate_ga_grwint(me)
    !------------------------------------------------------------------------
    type(grw_grwintarray_type) :: me !< array to truncate
    !------------------------------------------------------------------------
    type(grw_intarray_type), allocatable :: tarray(:)
    !------------------------------------------------------------------------
    integer :: ii
    !------------------------------------------------------------------------

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

  end subroutine truncate_ga_grwint


  subroutine empty_ga_grwint(me)
    !------------------------------------------------------------------------
    type(grw_grwintarray_type) :: me !< array to sorttruncate
    !------------------------------------------------------------------------

    me%nvals = 0

  end subroutine empty_ga_grwint

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


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

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

  end subroutine placeat_ga_grwint


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

    integer :: ub, ii

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

    ub = pos + size(val) - 1

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

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

  end subroutine placeat_ga_grwint_vec


  subroutine append_ga_grwint(me, val, length)
    type(grw_grwintarray_type) :: me !< array to append the value to
    type(grw_intarray_type), intent(in) :: val !< value to append
    !> optional length to expand the array
    integer, intent(in), optional :: length

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

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

  end subroutine append_ga_grwint

  subroutine append_ga_grwint_vec(me, val, length)
    type(grw_grwintarray_type) :: me !< array to append the value to
    type(grw_intarray_type), intent(in) :: val(:) !< values to append
    !> optional length to expand the array
    integer, intent(in), optional :: length

    integer :: lb, ub, ii

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

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

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

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

  end subroutine append_ga_grwint_vec


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

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

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

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

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

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

  end subroutine expand_ga_grwint


  !> Initialize the node type.
  subroutine init_node(me, nColors, none_color_id, distanceRefine, nSublevels, &
    &                  length)
    !> The nodelist to initialize.
    type(sdr_node_type), intent(out) :: me

    !> Number of different colors that are to be used.
    integer, intent(in) :: nColors

    !> Id of the color with the 'none' label, this color will be ignored in the
    !! treelm color data. It has to be -1, if there is no none color.
    integer, intent(in) :: none_color_id

    !> Number of levels to refine beyond target element levels for boundaries,
    !! that require this. Has to be 0, if there are no such boundaries.
    integer, intent(in) :: nSublevels

    !> If any boundary attribute has distance refine
    logical, intent(in) :: distanceRefine

    !> Initial length for the list of nodes.
    integer, intent(in), optional :: length

    integer :: iSide
    integer :: nNonecolors

    me%nNodes = 0
    me%nColors = nColors
    me%none_color_id = none_color_id
    me%subelement_resolution = nSublevels
    me%distanceRefine = distanceRefine
    me%memLeft_userObj = 0
    me%memLeft_distObj = 0

    ! Number of none colors (0 or 1), can be found by evaluating none_color_id
    ! as it is -1 if there is none and positive otherwise:
    nNonecolors = min(none_color_id,0) + 1

    me%nColorChars = ceiling(real(nColors-nNonecolors)/real(colors_per_char))

    ! Number of bytes there are in a default integer
    me%bytes_per_int = bit_size(me%bytes_per_int)/8

    ! Start position of the last byte in default integers
    me%lastbyte_pos = (me%bytes_per_int-1)*8

    ! We need nColors + 1 bytes at least, thus the proplength is given by:
    me%propLength = ceiling(real(nColors+1,kind=rk)/me%bytes_per_int)

    call init(me%treeID, length=length)
    call init(me%PropertyBits, length = length, width = me%propLength)
    call init(me%hasBndBits, length = length)
    call init(me%minLevel, length)
    call init(me%intersected_object, length)
    call init(me%userObjPos, length)
    if (me%distanceRefine) then
      call init(me%intersected_distance, length)
      call init(me%distObjPos, length)
    end if
    do iSide=1,6
      call init(me%linkpos(iSide), length)
    end do
    call init(me%minBCID, length)
    if (me%subelement_resolution > 0) call init(me%sublevel, length)

  end subroutine init_node


  !> Truncate the growing arrays in the node list to their actual size.
  subroutine truncate_node(me)
    !> The nodelist to initialize.
    type(sdr_node_type), intent(inout) :: me

    integer :: iSide

    !!call truncate(me%treeID)
    !!call truncate(me%PropertyBits)
    call truncate(me%hasBndBits)
    call truncate(me%minLevel)
    call truncate(me%intersected_object)
    call truncate(me%userObjPos)
    if (me%distanceRefine) then
      call truncate(me%intersected_distance)
      call truncate(me%distObjPos)
    end if
    do iSide=1,6
      call truncate(me%linkpos(iSide))
    end do

    call truncate(me%minBCID)
    if (me%subelement_resolution > 0) call truncate(me%sublevel)

  end subroutine truncate_node


  ! ****************************************************************************
  !> Append a new node to the protoTree.
  !!
  !! Please note that the chosen data structures exploit the fact, that the
  !! algorithm to build the tree is doing this ordered by treeIDs. Therfore,
  !! this routine appends a new node at the end of the list of existing nodes
  !! without any checks.
  subroutine append_newNode(me, treeID, PropertyBits, sublevel, minlevel,      &
    &                       pos, grwTreeID)
    ! --------------------------------------------------------------------------!
    !> The nodelist to which the node is to be appended.
    type(sdr_node_type), intent(inout) :: me

    !> The treeID of the node to append.
    integer(kind=long_k), intent(in) :: treeID

    !> The Bitmask of initial properties to set.
    integer, intent(in) :: PropertyBits(:)

    !> Sublevel indication, counts down after the target node. Has to be
    !! negative if not a subelement. Only relevant if subelement_resolution > 0.
    !!
    !! Refinement will be stopped, when sublevel 0 is reached.
    integer, intent(in) :: sublevel

    !> Minlevel to refine this node
    integer, intent(in) :: minLevel

    !> Position in the list, where the node was appended.
    integer, intent(out) :: pos

    !> growing array of treeID, append treeID to this list if present
    type(grw_longArray_type), intent(inout), optional :: grwTreeID
    ! --------------------------------------------------------------------------!
    integer :: iSide
    type(sdr_intersectObjPos_type) :: objPos
    ! --------------------------------------------------------------------------!

    if (present(grwTreeID)) then
      call append(grwTreeID, treeID)
      me%nNodes = me%nNodes + 1
      pos = me%nNodes
    else
      call append(me%treeID, treeID, pos=pos)
      me%nNodes = me%treeID%nVals
    end if

    call append(me%PropertyBits, PropertyBits)
    call append(me%hasBndBits, 0)
    do iSide=1,6
      call append(me%linkpos(iSide), 0)
    end do

    call append(me%minLevel, minLevel)

    !KM: No need to append empty intersected_object
    !call append(me%intersected_object, 0)
    objPos%first = 1
    objPos%last = 0
    call append(me%userObjPos, objPos)

    if (me%distanceRefine) then
      call append(me%distObjPos, objPos)
    end if

    call append(me%minBCID, 0)

    if (me%subelement_resolution > 0) call append(me%sublevel, sublevel)

  end subroutine append_newNode
  ! ****************************************************************************


  ! ****************************************************************************
  !> Set a bit in the last byte of the node properties in a given node.
  !!
  !! This routine can be used to set a property bit in the last byte of a node,
  !! which is used for the flags that are independent of the colors.
  subroutine sdr_set_nodeProp_bit(node, iNode, bit)
    type(sdr_node_type), intent(inout) :: node
    integer, intent(in) :: iNode
    integer, intent(in) :: bit

    integer :: nodeprop

    nodeprop = node%PropertyBits%val(node%propLength, iNode)
    node%PropertyBits%val(node%proplength, iNode) &
      &  = ibSet(nodeprop, node%lastbyte_pos + bit)
  end subroutine sdr_set_nodeProp_bit
  ! ****************************************************************************


  ! ****************************************************************************
  !> Clear a bit in the last byte of the node properties in a given node.
  !!
  !! This routine can be used to clear a property bit in the last byte of a
  !! node, which is used for the flags that are independent of the colors.
  subroutine sdr_clear_nodeProp_bit(node, iNode, bit)
    type(sdr_node_type), intent(inout) :: node
    integer, intent(in) :: iNode
    integer, intent(in) :: bit

    integer :: nodeprop

    nodeprop = node%PropertyBits%val(node%propLength, iNode)
    node%PropertyBits%val(node%proplength, iNode) &
      &  = ibClr(nodeprop, node%lastbyte_pos + bit)
  end subroutine sdr_clear_nodeProp_bit
  ! ****************************************************************************


  ! ****************************************************************************
  !> Set a bit in the last byte of the node properties in a given node.
  !!
  !! This routine can be used to set a property bit in the last byte of a node,
  !! which is used for the flags that are independent of the colors.
  pure function sdr_nodeProp_btest(node, iNode, bit) result(isSet)
    type(sdr_node_type), intent(in) :: node
    integer, intent(in) :: iNode
    integer, intent(in) :: bit
    logical :: isSet

    integer :: nodeprop

    nodeprop = node%PropertyBits%val(node%propLength, iNode)
    isSet = btest(nodeprop, node%lastbyte_pos + bit)
  end function sdr_nodeProp_btest
  ! ****************************************************************************


  ! ****************************************************************************
  !> This routine will wet the neighbors sides of the neighbor node
  !!
  !! col_int is the integer in which the flags for this color are encoded.
  !! col_bit is the 0th bit of this color within the integer.
  !! Thus, for the first color it should be col_int = 1, col_bit = 0.
  !! iNode is the index of the node whose neighbors are to be marked as
  !! wet. Note, that this will be done unconditionally, that is regardless of
  !! the state in iNode itself.
  !! Be aware, that the indices of neighbor need to be available in the linkpos
  !! array when calling this routine!
  subroutine sdr_wetNeighborsFace_single(node, iNode, col_int, col_bit)
    ! --------------------------------------------------------------------------!
    !> List of all nodes in the tree
    type(sdr_node_type), intent(inout) :: node
    !> Current node
    integer, intent(in) :: iNode
    integer, intent(in) :: col_int
    integer, intent(in) :: col_bit
    ! --------------------------------------------------------------------------!
    integer :: iDir, iSide, dir_off, neighbor_pos, neighborside
    integer :: mySide
    integer :: neighborbit_pos
    ! --------------------------------------------------------------------------!

    dirLoop: do iDir=1,3
      dir_off = (iDir-1)*2
      sideLoop: do iSide=1,2
        mySide = dir_off + iSide
        neighbor_pos = node%linkpos(mySide)%val(iNode)

        ! Copy the color of the current node to the adjacent side of the
        ! neighboring element. This side can be found with the following
        ! modulo operation.
        ! CAUTION: neighborside is defined with 0,1 where mySide is 1,2.
        !          the modulo maps 1 to 0 and 2 to 1, resulting in the
        !          expected behavior.
        neighborside = dir_off + mod(iSide, 2)

        ! Bit 0 is used for the intersection status,
        ! Bit 1 is used for the flooding status.
        ! Wetness is stored from Bit 2 onwards,
        ! neighborside ranges from 0 to 5 and col_bit
        ! provides the index of Bit 0 for this color.
        ! The bit, that has to be set for this side is therefore given by:
        neighborbit_pos = col_bit + neighborside + 2

        node%PropertyBits%val(col_int, neighbor_pos) &
          &   = ibset(node%PropertyBits%val(col_int, neighbor_pos), &
          &           neighborbit_pos)

      end do sideLoop
    end do dirLoop

  end subroutine sdr_wetNeighborsFace_single
  ! ****************************************************************************


  ! ****************************************************************************
  !> This routine will wet the neighbors sides of the neighbor node in all
  !! colors, that are flooded.
  !!
  !! The neighbors have to be available in the linkpos array prior calling
  !! this routine!
  subroutine sdr_wetNeighborsFace_all(node, iNode)
    ! --------------------------------------------------------------------------!
    !> List of all nodes in the tree
    type(sdr_node_type), intent(inout) :: node
    !> Current node
    integer, intent(in) :: iNode
    ! --------------------------------------------------------------------------!
    integer :: iColor
    integer :: col_int
    integer :: col_bit
    logical :: isflooded
    ! --------------------------------------------------------------------------!

    ! Extract the color of the current node from the property bits.
    do iColor=1,node%nColors
      col_int = (iColor-1) / node%bytes_per_int + 1
      col_bit = mod(iColor-1, node%bytes_per_int)*8
      isFlooded = btest(node%PropertyBits%val(col_int,iNode), &
        &               col_bit+1)
      if (isFlooded) then
        call sdr_wetNeighborsFace_single(node, iNode, col_int, col_bit)
      end if
    end do

  end subroutine sdr_wetNeighborsFace_all
  ! ****************************************************************************


  ! ****************************************************************************
  !> This routine floods the node with the given color and increases
  !! nFloodedLeaves.
  subroutine sdr_mark_floodNode( node, iNode, nFloodedLeaves, color )
    ! --------------------------------------------------------------------------!
    !> Node object
    type(sdr_node_type), intent(inout) :: node
    !> Node to modify
    integer, intent(in) :: iNode
    !> Number of flooded leaves to increase
    integer, intent(inout) :: nFloodedLeaves
    !> Color to set in the given bitfield
    integer, intent(in) :: color
    ! --------------------------------------------------------------------------!
    integer :: col_int, col_bit
    integer :: nodeFlooded_bit

    nodeFlooded_bit = node%lastbyte_pos + isFlooded_bit
    ! If this node was not flooded before, increment the counter for flooded
    ! leaves, and mark the node as flooded now.
    if ( .not. btest(node%PropertyBits%val(node%proplength, iNode), &
      &              nodeFlooded_bit) ) then
      nFloodedLeaves = nFloodedLeaves + 1

      node%PropertyBits%val(node%proplength, iNode) &
        &  = ibset(node%PropertyBits%val(node%proplength, iNode), &
        &          nodeFlooded_bit)
    end if
    col_int = (Color-1) / node%bytes_per_int + 1
    col_bit = mod(Color-1, node%bytes_per_int)*8 + 1
    node%PropertyBits%val(col_int,iNode) &
      &   = ibset(node%PropertyBits%val(col_int,iNode), col_bit)
  end subroutine sdr_mark_floodNode
  ! ****************************************************************************


  ! ****************************************************************************
  !> If parent has hasBoundary_bit then this function will inherit
  !! this property to eligible childrens
  pure function sdr_inHeritBnd_eligibleChildren(node, iNode) &
    &  result(child_hasBnd)
    ! --------------------------------------------------------------------------!
    !> Description of the nodes in the prototree.
    type(sdr_node_type), intent(in) :: node
    !> Index of the node to evaluate.
    integer, intent(in) :: iNode
    !> result contains which child has boundary
    logical :: child_hasBnd(8)
    ! --------------------------------------------------------------------------!
    integer, allocatable :: eligible_childs(:)
    integer :: iDir, iChild
    ! --------------------------------------------------------------------------!
    ! set default
    child_hasBnd = .false.
    if ( sdr_nodeProp_btest(node  = node,             &
      &                     iNode = iNode,            &
      &                     bit   = hasBoundary_bit ) ) then

      ! Inherit hasBoundary bit only to eligible children
      do iDir = 1, qQQQ
        if ( btest(node%hasBndBits%val(iNode), iDir) ) then
          call tem_eligibleChildren(eligible_childs, iDir)
          do iChild = 1, size(eligible_childs)
            child_hasBnd( eligible_childs(iChild) ) = .true.
          end do
          deallocate(eligible_childs)
        end if
      end do
    end if

  end function sdr_inHeritBnd_eligibleChildren
  ! ****************************************************************************


  ! ****************************************************************************
  !> Return all colors of the given node encoded in an array of characters.
  !!
  !! Each character is supposed to hold an ASCII symbol and by setting the
  !! byte bits, we can store up to 7 colors in each character.
  function sdr_nodeColors(node, iNode) result(colchar)
    ! --------------------------------------------------------------------------!
    !> Description of the nodes in the prototree.
    type(sdr_node_type), intent(in) :: node
    !> Index of the node to evaluate.
    integer, intent(in) :: iNode
    !> Gathered color information in an array of characters.
    character :: colchar(node%nColorChars)
    ! --------------------------------------------------------------------------!

    colchar = sdr_bitfieldColors( node = node,                              &
      &                           bitfield =node%PropertyBits%val(:, iNode) )

  end function sdr_nodeColors
  ! ****************************************************************************


  ! ****************************************************************************
  !> Return all colors of the given bitfield encoded in an array of characters.
  !!
  !! The bitfield has to correspond with the color definition in nodes.
  !! Each character is supposed to hold an ASCII symbol and by setting the
  !! byte bits, we can store up to 7 colors in each character.
  function sdr_bitfieldColors(node, bitfield) result(colchar)
    ! --------------------------------------------------------------------------!
    !> Description of the nodes in the prototree.
    type(sdr_node_type), intent(in) :: node
    !> Index of the node to evaluate.
    integer, intent(in) :: bitfield(:)
    !> Gathered color information in an array of characters.
    character :: colchar(node%nColorChars)
    ! --------------------------------------------------------------------------!
    integer :: charcount, intcount
    integer :: charbit, intbit
    integer :: iColor, iNotNone
    integer :: tmpColor(node%nColorChars)
    ! --------------------------------------------------------------------------!

    tmpColor = 0

    iNotNone = 0

    do iColor=0,node%nColors-1

      ! Only process non-none colors.
      if (iColor+1 /= node%none_color_id) then

        ! Integer and bit in integer for this color (bit 1 decides flooding)
        intcount = iColor/node%bytes_per_int + 1
        intbit = mod(iColor, node%bytes_per_int)*8 + 1

        if ( btest(bitfield(intcount), intbit) ) then
          ! Character and bit in character for this color
          charcount = iNotNone/colors_per_char + 1
          charbit = mod(iNotNone, colors_per_char)

          tmpColor(charcount) = ibset(tmpColor(charcount), charbit)
        end if

        iNotNone = iNotNone + 1

      end if

    end do

    colchar = achar(tmpColor)

  end function sdr_bitfieldColors
  ! ****************************************************************************


  ! ****************************************************************************
  !> Return all colors encoded in an array of characters based on an array of
  !! logicals indicating for each color wether it should be set or not.
  !!
  !! Each character is supposed to hold an ASCII symbol and by setting the
  !! byte bits, we can store up to 7 colors in each character.
  function sdr_color_log2char(node, logicals) result(colchar)
    ! --------------------------------------------------------------------------!
    !> Description of the nodes in the prototree.
    type(sdr_node_type), intent(in) :: node

    !> Flag for each color wether it should be set or not.
    logical, intent(in) :: logicals(:)

    !> Gathered color information in an array of characters.
    character :: colchar(node%nColorChars)
    ! --------------------------------------------------------------------------!
    integer :: charcount
    integer :: charbit
    integer :: iColor, iNotNone
    integer :: tmpColor(node%nColorChars)
    ! --------------------------------------------------------------------------!

    tmpColor = 0

    iNotNone = 0

    do iColor=0,node%nColors-1

      ! Only process non-none colors.
      if (iColor+1 /= node%none_color_id) then

        if ( logicals(iColor+1) ) then
          ! Character and bit in character for this color
          charcount = iNotNone/colors_per_char + 1
          charbit = mod(iNotNone, colors_per_char)

          tmpColor(charcount) = ibset(tmpColor(charcount), charbit)
        end if

        iNotNone = iNotNone + 1

      end if

    end do

    colchar = achar(tmpColor)

  end function sdr_color_log2char
  ! ****************************************************************************


  ! ****************************************************************************
  !> This routine appends temporary child intersected object to actual growing
  !! array of intersected objects. To reduce memory usuage, the child with
  !! maximum number of intersected objects from parent is appended at same
  !! position as its parent. Rest of the childrens intersected objects are
  !! appended to the end of growing array
  subroutine sdr_append_childIntersectedObjectAll(node, parent, testAll,       &
    & intersected_object, grwObjPos, child_nodePos, child_nObjects,            &
    & child_intersected_object, child_objPos, memLeft)
    ! --------------------------------------------------------------------------!
    !> contains information about all nodes in protoTree
    type(sdr_node_type), intent(in) :: node
    !> Position of parent node
    integer, intent(in) :: parent
    !> true for level 1 nodes
    logical, intent(in) :: testAll
    !> Growing array of intersected objects.
    !! Could be user defined  or distance refine spatial objects
    type(grw_intArray_type), intent(inout) :: intersected_object
    !> First and last position of intersected object of all nodes in
    !! intersected_object list
    type(grw_intersectObjPosArray_type), intent(inout) :: grwObjPos
    !> 8 children node position in protoTree
    integer, intent(in) :: child_nodePos(8)
    !> number of intersected objets in 8 children
    integer, intent(in) :: child_nObjects(8)
    !> Temporary array of intersected objects for 8 children
    type(grw_intArray_type), intent(inout) :: child_intersected_object
    !> first and last index of intersected objects of 8 children in
    !! child_intersected_object list
    type(sdr_intersectObjPos_type), intent(in) :: child_objPos(8)
    !> memory of parent intersected object unused by children
    integer, intent(out) :: memLeft
    ! --------------------------------------------------------------------------!
    integer :: iChild, childIDX, child_sorted(8)
    integer :: child_intersected_first_inParent
    logical :: parent_isTarget
    ! --------------------------------------------------------------------------!
    ! Append child with maximum number of intersected object
    ! in the position of parent to reduce memory only if parent
    ! is not a target bit.

    ! if parent is target bit do not write intersected object
    ! of children in parent position
    parent_isTarget = sdr_nodeProp_btest(node  = node,       &
        &                                iNode = parent,     &
        &                                bit   = isTarget_bit)

    ! Find the sorted list
    call mrgrnk(child_nObjects, child_sorted)

    ! left space is number of objects we can still fit in parent position
    ! Initially, We could fill upto nObjects of parent
    memLeft = grwObjPos%val(parent)%last - grwObjPos%val(parent)%first + 1

    ! Start position of a child in parent list
    child_intersected_first_inParent = grwObjPos%val(parent)%first

    ! Find the sorted list, mrgrnk sorted min to max so use
    ! offset 8-iChild to append child with max nObjects first
    do iChild = 8, 1, -1
      ! append child according to sorted list
      childIDX = child_sorted(iChild)

      ! append child only when it has intersected objects
      ! if childs nObjects is <= left_space than fill this child
      ! in parent position and set intersected first and last accordingly.
      ! Also update left_space
      if (child_nObjects(childIDX) <= memLeft .and.&
        & .not. testAll .and. .not. parent_isTarget) then
        ! Set the first and last index of this childs intesected object
        ! in the array
        grwObjPos%val(child_nodePos(childIDX))%first =             &
          &                       child_intersected_first_inParent
        grwObjPos%val(child_nodePos(childIDX))%last =              &
          &                       child_intersected_first_inParent &
          &                     + child_nObjects(childIDX) - 1

        call placeAt( me  = intersected_object,               &
          &           val = child_intersected_object%val(     &
          &                    child_objPos(childIDX)%first : &
          &                    child_objPos(childIDX)%last ), &
          &           pos = child_intersected_first_inParent  )

        ! update left_space
        memLeft = memLeft - child_nObjects(childIDX)
        ! update next childs start position in parent list
        child_intersected_first_inParent =                &
          & grwObjPos%val(child_nodePos(childIDX))%last + 1

      else
        ! append to end of intersected object list

        ! Set the first and last index of this childs intesected object
        ! in the array
        grwObjPos%val(child_nodePos(childIDX))%first =       &
          &                       intersected_object%nVals + 1

        grwObjPos%val(child_nodePos(childIDX))%last =     &
          &                      intersected_object%nVals &
          &                    + child_nObjects(childIDX)
        call append( me  = intersected_object,               &
          &          val = child_intersected_object%val(     &
          &                   child_objPos(childIDX)%first : &
          &                   child_objPos(childIDX)%last )  )

      end if
    end do

  end subroutine sdr_append_childIntersectedObjectAll
  ! ****************************************************************************


  ! ****************************************************************************
  !> This routine appends temporary child intersected object to actual growing
  !! array of intersected objects. To reduce memory usuage, the child with
  !! maximum number of intersected objects from parent is appended at same
  !! position as its parent. Rest of the childrens intersected objects are
  !! appended to the end of growing array.
  !!
  !! We append only objects with level greater than children minLevel
  !! to reduce memory.
  subroutine sdr_append_childIntersectedObjectGTminLevel(geometry, node,       &
    & parent, testAll, intersected_object, grwObjPos, child_nodePos,           &
    & child_intersected_object, child_objPos, memLeft)
    ! --------------------------------------------------------------------------!
    !> type which contains all geometry object infos
    type(sdr_geometry_type), intent(in) :: geometry
    !> contains information about all nodes in protoTree
    type(sdr_node_type), intent(in) :: node
    !> Position of parent node
    integer, intent(in) :: parent
    !> true for level 1 nodes
    logical, intent(in) :: testAll
    !> Growing array of intersected objects.
    !! Could be user defined  or distance refine spatial objects
    type(grw_intArray_type), intent(inout) :: intersected_object
    !> First and last position of intersected object of all nodes in
    !! intersected_object list
    type(grw_intersectObjPosArray_type), intent(inout) :: grwObjPos
    !> 8 children node position in protoTree
    integer, intent(in) :: child_nodePos(8)
    !> Temporary array of intersected objects for 8 children
    type(grw_intArray_type), intent(inout) :: child_intersected_object
    !> first and last index of intersected objects of 8 children in
    !! child_intersected_object list
    type(sdr_intersectObjPos_type), intent(in) :: child_objPos(8)
    !> memory of parent intersected object unused by children
    integer, intent(out) :: memLeft
    ! --------------------------------------------------------------------------!
    integer :: iChild, childIDX, child_sorted(8)
    integer :: child_objPos_first_inParent, child_objPos_next_inParent
    integer :: child_nObjectsNew(8)
    logical :: parent_isTarget
    integer :: iObject, obj_pos, attr_pos, objLevel
    ! --------------------------------------------------------------------------!
    ! Append child with maximum number of intersected object
    ! in the position of parent to reduce memory only if parent
    ! is not a target bit.

    ! if parent is target bit do not write intersected object
    ! of children in parent position
    parent_isTarget = sdr_nodeProp_btest(node  = node,       &
        &                                iNode = parent,     &
        &                                bit   = isTarget_bit)

    ! left space is number of objects we can still fit in parent position
    ! Initially, We could fill upto nObjects of parent
    memLeft = grwObjPos%val(parent)%last - grwObjPos%val(parent)%first + 1

    ! Start position of a child in parent list
    child_objPos_first_inParent = grwObjPos%val(parent)%first
    child_objPos_next_inParent = child_objPos_first_inParent

    child_nObjectsNew = 0!child_nObjects
    do iChild = 1, 8
      do iObject = child_objPos(iChild)%first, child_objPos(iChild)%last
        obj_pos = child_intersected_object%val(iObject)
        attr_pos = geometry%spatialObj%val(obj_pos)%attribute_position
        objLevel = geometry%attribute%dynArray%val(attr_pos)%level

        if (objLevel > node%minLevel%val(child_nodePos(iChild))) then
          ! number of objects with level > minLevel
          child_nObjectsNew(iChild) = child_nObjectsNew(iChild) + 1
        end if !objlevel>minLevel
      end do ! iObject
    end do !iChild

    ! Find the sorted list
    call mrgrnk(child_nObjectsNew, child_sorted)

    ! Find the sorted list, mrgrnk sorted min to max so use
    ! offset 8-iChild to append child with max nObjects first
    do iChild = 8, 1, -1
      ! append child according to sorted list
      childIDX = child_sorted(iChild)

      ! append child only when it has intersected objects
      ! if childs nObjects is <= left_space than fill this child
      ! in parent position and set intersected first and last accordingly.
      ! Also update left_space
      if (child_nObjectsNew(childIDX) <= memLeft .and.&
        & .not. testAll .and. .not. parent_isTarget) then

        ! append objects with level > child minlevel
        ! Count new number of objects per child which has level
        ! greater than childs minLevel and
        do iObject = child_objPos(childIDX)%first, child_objPos(childIDX)%last
          obj_pos = child_intersected_object%val(iObject)
          attr_pos = geometry%spatialObj%val(obj_pos)%attribute_position
          objLevel = geometry%attribute%dynArray%val(attr_pos)%level

          if (objLevel > node%minLevel%val(child_nodePos(childIDX))) then
            call placeAt( me  = intersected_object,        &
              &           val = obj_pos,                   &
              &           pos = child_objPos_next_inParent )

            ! Position of next childObjPos in parent
            child_objPos_next_inParent = child_objPos_next_inParent + 1
          end if !objlevel>minLevel
        end do ! iObject

        ! Set the first and last index of this childs intesected object
        ! in the array
        grwObjPos%val(child_nodePos(childIDX))%first =             &
          &                       child_objPos_first_inParent
        grwObjPos%val(child_nodePos(childIDX))%last =           &
          &                       child_objPos_first_inParent   &
          &                     + child_nObjectsNew(childIDX) - 1

        ! update left_space
        memLeft = memLeft - child_nObjectsNew(childIDX)
        ! update next childs start position in parent list
        child_objPos_first_inParent =                     &
          & grwObjPos%val(child_nodePos(childIDX))%last + 1

      else
        ! append to end of intersected object list

        ! Set the first and last index of this childs intesected object
        ! in the array
        grwObjPos%val(child_nodePos(childIDX))%first =       &
          &                       intersected_object%nVals + 1
        ! append objects with level > child minlevel
        ! Count new number of objects per child which has level
        ! greater than childs minLevel and
        do iObject = child_objPos(childIDX)%first, child_objPos(childIDX)%last
          obj_pos = child_intersected_object%val(iObject)
          attr_pos = geometry%spatialObj%val(obj_pos)%attribute_position
          objLevel = geometry%attribute%dynArray%val(attr_pos)%level

          if (objLevel > node%minLevel%val(child_nodePos(childIDX))) then
            call append( me  = intersected_object, val = obj_pos )
          end if !objlevel>minLevel
        end do ! iObject

        ! Set the last index of this childs intesected object
        ! in the array
        grwObjPos%val(child_nodePos(childIDX))%last =      &
          &                      intersected_object%nVals

      end if
    end do

  end subroutine sdr_append_childIntersectedObjectGTminLevel
  ! ****************************************************************************


end module sdr_node_module
