! Copyright (c) 2012-2013, 2015 Kannan Masilamani <kannan.masilamani@uni-siegen.de>
! Copyright (c) 2012-2014, 2017, 2022 Harald Klimach <harald.klimach@dlr.de>
! Copyright (c) 2012, 2014, 2016 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
!
!> Module to describe attributes attached to geometrical objects.
!!
!! Attribute list will be a dynamic unique list. Thus, multiple definitions of
!! the same attribute in the configuration file does not result in multiple
!! storages.
module sdr_attribute_module
  use env_module,            only: minLength, labelLen, zeroLength, rk
  use tem_aux_module,        only: tem_abort
  use tem_dyn_array_module,  only: dyn_labelArray_type, init, append, &
    &                              PositionOfVal
  use tem_grow_array_module, only: grw_intArray_type, grw_logicalArray_type, &
    &                              init, append
  use tem_tools_module,      only: upper_to_lower
  use tem_logging_module,    only: logunit

  use flu_binding,        only: flu_State, flu_next
  use aotus_module,       only: aot_get_val, flu_State,                        &
    &                           aoterr_Fatal, aoterr_NonExistent,              &
    &                           aoterr_WrongType
  use aot_table_module,   only: aot_table_open, aot_table_close,               &
    &                           aot_table_length, aot_table_top, aot_table_first

  implicit none

  private

  public :: sdr_attribute_type
  public :: sdr_load_attribute
  public :: sdr_init_attribute
  public :: dyn_attributeArray_type
  public :: sdr_attrList_type
  public :: sdr_identify_bc_colors
  public :: sdr_identify_inv_colors
  public :: sdr_any_bc_subresolution
  public :: sdr_any_bc_distanceRefine
  public :: sdr_isPeriodicDefined

  public :: init, append, PositionOfVal, truncate, destroy, empty, placeAt
  public :: sortTruncate

  !> Maximal number of different object kinds to distinguish.
  !!
  !! This is used to allow allocations of arrays for all
  !! object kinds with direct access to each component
  !! by the sdr_*_object constants below.
  integer, parameter, public :: sdr_object_kinds_max = 5

  ! Some parameters to easily identify the different kinds
  ! objects, that are distinguished.
  ! They have to be in the range 1:sdr_object_max!

  !> Identifier for boundary objects
  integer, parameter, public :: sdr_Boundary_object = 1

  !> Identifier for seed objects
  integer, parameter, public :: sdr_Seed_object = 2

  !> Identifier for refinement objects
  integer, parameter, public :: sdr_Refinement_object = 3

  !> Identifier for solid objects, that can be fluidify in the solver
  integer, parameter, public :: sdr_Fluidifyable_object = 4

  !> Identifier for no Solidification objects
  integer, parameter, public :: sdr_noSolidification_object = 5

  !> This data type contains info for distance refinement defined
  !! in config file for each boundary attribute.
  type sdr_distanceRefine_type
    !> Defines the distance to which this boundary
    !! should be refined
    real(kind=rk) :: radius

    !> level to refine nodes with in given radius
    integer :: reach_level
  end type sdr_distanceRefine_type

  !> This data type describes type of attribute
  type sdr_attribute_type
    !> Kind of this attribute:
    !! - boundary
    !! - seed
    !! - refinement
    !! - fluidifyable
    !! - noSolidification
    integer :: kind

    !> Label to identify this attribute.
    character(len=labelLen) :: label

    !> Color to associate the object with.
    !!
    !! This is only relevant for seeds and boundaries.
    !! Refinements do not have colors.
    !! Color names are case-insensitive!
    !! Boundaries limit seeds of the same color, or all seeds, if they
    !! do not have a color ('none' or 'all').
    !! Only boundaries without a color will actually be written as boundary
    !! conditions all others are assumed to be covered internally.
    !! NOTICE, that this results in a somewhat special behavior, where seeds
    !! of color 'none' or 'all' are just usual colors (though none and all will
    !! be the same), while boundaries of this color will confine all colors,
    !! not just the seeds with 'none' color. This means, that 'none' colored
    !! seeds can not be confined on their own. Their boundaries will always
    !! bound all colors!
    !! Boundaries, with a color matching no seed will be ignored
    !! by the flooding but are still relevant for refinement.
    !! The default is 'none'.
    character(len=labelLen) :: color

    !> Refinement level to resolve this attribute with at least.
    integer :: level

    !> ID in the list of attributes of this kind.
    integer :: id = 0

    !> Uni-id of this attribute in the list of unique names for this kind.
    !! (labels for boundaries,
    !!  colors for seeds)
    integer :: uni_id = 0

    !> This object should be resolved further than the elements of the mesh.
    !!
    !! Only useful for boundary objects. And only effective, if this boundary
    !! has a specific color (not the none-color).
    logical :: subresolution = .false.

    !> Indicator of whether to calculate link-wise distance
    logical :: calc_dist = .false.

    !> Indicator of whether to flood diagonal node during
    !! flood_periphery
    logical :: flood_diagonal = .true.

    !> Uni-id to identify attribute with certain number of
    !! distance refinement
    integer :: distRefine_id = 0

    !> Start position of distance refine object of this attribute in
    !! growing array of distance refine
    integer :: distance_first = 0

    !> Last position of distance refine object of this attribute in
    !! growing array of distance refine
    integer :: distance_last = 0

  end type sdr_attribute_type


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

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

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

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

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

  !> insert an element at a given position
  interface placeat
    module procedure placeat_ga_distancerefine
    module procedure placeat_ga_distancerefine_vec
  end interface

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

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

! \brief smart dynamic array (da) for type(sdr_attribute_type)
!
! this datatype implements a dynamic array, which is capable of
! growing and adding of unique elements. it is available for
! various types, here we deal with $tstring$.
! sorted array contains the pointers to val array, instead of
! the actual values in val array. for example:
! val:     8, 6, 7, 5
! sorted:  4, 2, 3, 1
  !> dynamic array (da) type for type(sdr_attribute_type)
  type dyn_attributearray_type
    integer :: nvals = 0
    integer :: containersize = 0
    type(sdr_attribute_type), allocatable :: val(:)
    integer, allocatable :: sorted(:) !< pointers, not values
  end type

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

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

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

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

  !> empty the array, reset nvals to be 0
  interface empty
    module procedure empty_da_attribute
  end interface


  !> fix the dynamic array, meaning:
  !! store the array in the sorted order and cut
  !! off the trailing empty entries
  interface sorttruncate
    module procedure sorttruncate_da_attribute
  end interface

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

  !> return the position of a given value
  !! in the array val, which is what you usually want to know.
  !! it is the index of a given value
  interface positionofval
    module procedure posofval_attribute
  end interface

  !> return the position of a given value
  !! in the list 'sorted'. this is mainly for internal usage.
  !! the sorted list is only a pointer list to the actual values
  !! thus, in order to get the index of a given value, you
  !! need to look up the entry in the sorted list.
  !! this is done by the positionofval routine
  interface sortedposofval
    module procedure sortposofval_attribute
  end interface


  type sdr_attrList_type
    !> Dynamic array of all unique attributes.
    !!
    !! Attributes are distinguished by kind, label, color and level.
    type(dyn_attributeArray_type) :: dynArray

    !> Unique list of used names for the various attribute kinds.
    !!
    !! (Labels for boundaries, colors for seeds, not relevant otherwise)
    !! This is especially needed to find the actual number of different colors.
    type(dyn_labelArray_type) :: uni_name(sdr_object_kinds_max)

    !> uni_id of the 'none' color.
    !!
    !! The none color is special in that it will not be stored to disk.
    !! However, it will be used just like any other color during construction
    !! of the prototree. To filter the color later before writing the colors
    !! to disk, we store the index of the 'none' color here.
    !! If there is no 'none' color, this id is -1.
    integer :: none_color_id = -1

    !> ID of unique colors for each kindpos(sdr_Boundary_object).
    !!
    !! Required to identify the colors to which a certain boundary should
    !! apply.
    !! 'none' and 'all' boundaries: negative value
    !! matching seed labels: positive uni_id of that color
    !! non-matching: 0 (affecting no flooding color at all)
    integer, allocatable :: bc_color_id(:)

    !> Flag for each unique color, wether it should be inverted.
    !!
    !! Inverted colors will change reverse their flooded flags after flooding.
    !! Only elements that are flooded by the none color are considered, and
    !! the none color itself can not be inverted.
    logical, allocatable :: color_inverted(:)

    !> Growing array with flags on whether to calculate distances for
    !! each unique boundary.
    !!
    !! There might be conflicting calc_dist settings for attributes with
    !! the same boundary label. The `bc_uni_calcdist` consolidates this.
    !! If any of the attributes with the same label has the calc_dist set
    !! to true, this will be true and all attributes with that label will
    !! be changed accordingly.
    type(grw_logicalArray_type) :: bc_uni_calcdist

    !> A list of all the attribute positions for each kind.
    !! This provides the mapping of the kind id to the attribute position.
    type(grw_intArray_type) :: kindpos(sdr_object_kinds_max)

    !> Growing array of distance refinement info for boundary attribute.
    !! Use dynArray%distRefine_id to create attribute with distance refine
    type(grw_distanceRefineArray_type) :: distRefine

  end type sdr_attrList_type

!HK: To work around Cray compiler bug, we might need to define the assignment
!HK  ourselves...
!HK!  interface assignment(=)
!HK!    module procedure Copy_var
!HK!  end interface


  interface operator(==)
    module procedure isEqual
  end interface

  interface operator(/=)
    module procedure isUnequal
  end interface

  interface operator(<)
    module procedure isSmaller
  end interface

  interface operator(<=)
    module procedure isSmallerOrEqual
  end interface

  interface operator(>)
    module procedure isGreater
  end interface

!  interface operator(>=)
!    module procedure isGreaterOrEqual
!  end interface


contains


  !> Subroutine to initialize the list for all attributes
  subroutine sdr_init_attribute(me)
    type(sdr_attrList_type), intent(out) :: me

    me%none_color_id = -1

  end subroutine sdr_init_attribute



  ! ****************************************************************************!
  !> This routine loads the attribute information from the config file.
  subroutine sdr_load_attribute( attrList, conf, thandle, subres_colors, &
    &                            attr_pos )
    ! --------------------------------------------------------------------------!
    !> dynamic array of attribute type
    type(sdr_attrList_type), intent(inout) :: attrList
    type(flu_State) :: conf !< lua state
    integer, intent(in) :: thandle !<  handle to load attribute from
    !> List of colors which should by default use subelement resolution for
    !! their boundaries.
    type(dyn_labelArray_type), intent(in) :: subres_colors
    integer, intent(out) :: attr_pos !< position of attribute in attribute array
    ! --------------------------------------------------------------------------!
    integer :: iError
    character(len=labelLen) :: attr_kind
    integer :: attr_handle
    logical :: wasAdded
    logical :: newuni
    type(sdr_attribute_type) :: attribute ! save attibute info temporarily
    integer :: uni_pos
    logical :: merged_cd
    logical :: defsub
    integer :: iDistRefine
    integer :: iBC, iAttr
    type(grw_intArray_type) :: prev_bcs
    ! --------------------------------------------------------------------------!

    write(logunit(1),*) 'Loading attribute'
    ! Open the attribute table.
    call aot_table_open(L = conf, parent = thandle, thandle = attr_handle, &
      &                 key = 'attribute')

    ! Load attribute only if the table exists otherwise prompt an error and
    ! abort.
    if (attr_handle .gt. 0 ) then
      ! load KIND
      call aot_get_val(L = conf, thandle = attr_handle, key = 'kind', &
        &              val = attr_kind, ErrCode = iError )

      if (btest(iError, aoterr_Fatal)) then
        write(logUnit(0),*) &
          &  'FATAL Error occured, while retrieving attribute kind:'
        if (btest(iError, aoterr_NonExistent)) &
          &  write(logunit(0),*) 'Attribute kind not existent!'
        if (btest(iError, aoterr_WrongType)) &
          &  write(logunit(0),*)'Attribute kind has wrong type!'
        call tem_abort()
      end if

      select case(upper_to_lower(trim(attr_kind)))
        case('boundary','periodic')
          attribute%kind = sdr_Boundary_object

        case('seed')
          attribute%kind = sdr_Seed_object

        case('refinement')
          attribute%kind = sdr_Refinement_object

        case('fluidifyable')
          attribute%kind = sdr_Fluidifyable_object

        case('nosolidification')
          attribute%kind = sdr_noSolidification_object

        case default
          write(logunit(0),*) 'ERROR: Unknown attribute kind: '//trim(attr_kind)
          call tem_abort()

      end select

      ! Load LABEL of this attribute.
      call aot_get_val(L = conf, thandle = attr_handle, key = 'label', &
        &              val = attribute%label, ErrCode = iError, &
        &              default='sdr_nolabel' )

      if (btest(iError, aoterr_Fatal)) then
        write(logunit(0),*) &
          &  'FATAL Error occured, while retrieving attribute label:'
        if (btest(iError, aoterr_WrongType)) &
          &  write(logunit(0),*) 'Label has wrong type!'
        call tem_abort()
      end if

      ! Load COLOR of this attribute.
      call aot_get_val(L = conf, thandle = attr_handle, key = 'color', &
        &              val = attribute%color, ErrCode = iError, &
        &              default = 'none' )

      if (btest(iError, aoterr_Fatal)) then
        write(logunit(0),*) &
          &  'FATAL Error occured, while retrieving attribute label:'
        if (btest(iError, aoterr_WrongType)) &
          &  write(logunit(0),*) 'Label has wrong type!'
        call tem_abort()
      end if

      ! Color 'all' is to be the same as 'none'. Also, make this parameter case
      ! insensitive.
      attribute%color = upper_to_lower(trim(attribute%color))
      if (trim(attribute%color) == 'all') attribute%color = 'none'

      ! Load LEVEL associated with this attribute.
      ! The level describes the refinement level, that should be used at least
      ! to refine the objects with this attribute to.
      ! The default of 0 does not impose any refinement at all.
      call aot_get_val(L = conf, thandle = attr_handle, key = 'level', &
        &              val = attribute%level, ErrCode = iError, default = 0)
      if (btest(iError, aoterr_WrongType)) then
        write(logunit(0),*) &
          &  'FATAL Error occured, while retrieving attribute level:'
        write(logunit(0),*) 'Level has wrong type, should be an Integer!'
        call tem_abort()
      end if

      ! Loading boundary specific settings.
      bc_settings: if (attribute%kind == sdr_Boundary_object) then

        if (trim(attribute%color) /= 'none') then
          ! Set the default for the subresolution according to the color.
          defsub = (PositionOfVal(me=subres_colors, val=attribute%color) > 0)

          ! Should the objects with this attribute be resolved on a subelement
          ! level?
          !>\todo HK: subresolution should get its default from a color-specific
          !!          setting and only needs to be read, if
          !!          subelement_resolution > 0.
          call aot_get_val( L       = conf,                    &
            &               thandle = attr_handle,             &
            &               key     = 'subresolution',         &
            &               val     = attribute%subresolution, &
            &               ErrCode = iError,                  &
            &               default = defsub                   )

          if (btest(iError, aoterr_WrongType)) then
            write(logunit(0),*) &
              & 'FATAL Error occured, while retrieving attribute subresolution:'
            write(logunit(0),*) 'subresolution should be a logical!'
            call tem_abort()
          end if
        end if

        ! Load CALC_DIST indicator
        call aot_get_val( L       = conf,                &
          &               thandle = attr_handle,         &
          &               key     = 'calc_dist',         &
          &               val     = attribute%calc_dist, &
          &               ErrCode = iError,              &
          &               default = .false.              )

        if (btest(iError, aoterr_WrongType)) then
          write(logunit(0),*) &
            & 'FATAL Error occured, while retrieving attribute calc_dist:'
          write(logunit(0),*) 'calc_dist should be a logical!'
          call tem_abort()
        end if

        if (attribute%calc_dist) then
          ! Load flood_diagonal indicator (only relevant if calc_dist is active)
          call aot_get_val( L       = conf,                     &
            &               thandle = attr_handle,              &
            &               key     = 'flood_diagonal',         &
            &               val     = attribute%flood_diagonal, &
            &               ErrCode = iError,                   &
            &               default = .true.                    )

          if (btest(iError, aoterr_WrongType)) then
            write(logunit(0),*)  'FATAL Error occured, ' &
              &                  // 'while retrieving attribute flood_diagonal:'
            write(logunit(0),*) 'flood_diagonal should be a logical!'
            call tem_abort()
          end if
        end if

        ! load distance refine table from attribute table only for
        ! boundary attribute and append to growing
        ! array of distance refine in attrList
        call load_distanceRefine(attrList, conf, attr_handle, attribute)

      end if bc_settings

    else

      write(logunit(0),*) 'ERROR loading attribute for spatial object: '
      write(logunit(0),*) '  attribute table not defined'
      call tem_abort()

    endif

    call aot_table_close( L = conf, thandle = attr_handle )

    if (attribute%kind == sdr_Boundary_object) then
      ! Check whether this boundary is already present.
      ! If we already have boundaries we need to make sure that all of
      ! them have consistent calc_dist settings.
      ! (If any boundary with of the same label has calc_dist set, we
      !  need to set it for all boundaries of this label).
      uni_pos =  PositionOfVal( me = attrList%uni_name(sdr_Boundary_object), &
        &                       val = trim(attribute%label)                  )
      if (uni_pos > 0) then
        ! Found a boundary of the same name.
        merged_cd = attribute%calc_dist &
          &         .or. attrList%bc_uni_calcdist%val(uni_pos)
        if (merged_cd .neqv. attribute%calc_dist) then
          write(logunit(1),*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
          write(logunit(1),*) 'WARNING you set calc_dist=false here for'
          write(logunit(1),*) 'boundary '//trim(attribute%label)
          write(logunit(1),*) 'But it has already been defined as true by'
          write(logunit(1),*) 'another boundary object with this label.'
          write(logunit(1),*) 'This attribute will assume calc_dist=true!'
          write(logunit(1),*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
          attribute%calc_dist = merged_cd
        end if
        if (merged_cd .neqv. attrList%bc_uni_calcdist%val(uni_pos)) then
          write(logunit(1),*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
          write(logunit(1),*) 'WARNING you set calc_dist=true here for'
          write(logunit(1),*) 'boundary '//trim(attribute%label)
          write(logunit(1),*) 'But it has already been defined as false by'
          write(logunit(1),*) 'another boundary object with this label.'
          prev_bcs = sdr_attr_of_uni( attribute   = attrList,            &
            &                         object_kind = sdr_boundary_object, &
            &                         id          = uni_pos              )
          write(logunit(1),*) 'ATTENTION: setting calc_dist to TRUE for the' &
            &                 // ' following attributes:'
          write(logunit(1),*) prev_bcs%val(:prev_bcs%nVals)
          write(logunit(1),*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
          do iBC=1,prev_bcs%nVals
            iAttr = prev_bcs%val(iBC)
            attrList%dynArray%val(iAttr)%calc_dist = merged_cd
          end do
          attrList%bc_uni_calcdist%val(uni_pos) = merged_cd
        end if
      end if
    end if

    ! Append attribute to dynamic array
    call append( attrList%dynArray, attribute, &
      &          pos = attr_pos,               &
      &          wasAdded = wasAdded           )

    if (wasAdded) then

      write(logunit(1),*) 'This is a new attribute and was added at ', &
        &                 attr_pos
      write(logunit(1),"(A)") '           kind: '//trim(attr_kind)
      write(logunit(1),"(A)") '          label: '//trim(attribute%label)
      write(logunit(1),"(A)") '          color: '//trim(attribute%color)
      write(logunit(1),"(A,I0)") '          level: ', attribute%level
      if (attribute%kind == sdr_Boundary_object) then
        write(logunit(1),*) '  subresolution: ', attribute%subresolution
        write(logunit(1),*) '  calc distance: ', attribute%calc_dist
        if (attribute%calc_dist) then
          write(logunit(1),*) ' flood diagonal: ', attribute%flood_diagonal
        end if
        if (attribute%distRefine_id > 0) then
          write(logUnit(1),"(A,I0)") '  Nr. of distance_refine: ', &
            &   attribute%distance_last - attribute%distance_first + 1
          do iDistRefine = attribute%distance_first, attribute%distance_last
            write(logunit(1),"(A,I0)")    '    iDistRefine: ', iDistRefine
            write(logUnit(1),"(A,E13.6)") '         radius: ',      &
              &        attrList%distRefine%val(iDistRefine)%radius
            write(logUnit(1),"(A,I0)") '    reach_level: ', &
              &        attrList%distRefine%val(iDistRefine)%reach_level
          end do
        end if
      end if

      ! Update the counters for this specific kind of object.
      if (trim(attr_kind) /= 'periodic') then
        call append(attrList%kindpos(attribute%kind), attr_Pos)
        attrList%dynArray%val(attr_pos)%id &
          &  = attrList%kindpos(attribute%kind)%nVals

        select case(attribute%kind)
        case (sdr_Seed_object)
          call append(attrList%uni_name(attribute%kind), attribute%color, &
            &         pos = uni_pos)
          ! If this was the 'none' color, store its position.
          if (trim(attribute%color) == 'none') attrList%none_color_id = uni_pos

        case default
          call append(attrList%uni_name(attribute%kind), attribute%label, &
            &         pos = uni_pos, wasAdded = newuni)

          if (newuni) then
            write(logunit(1),*) 'Added a new unique label to the the list of' &
              &                 //' kind "'//trim(attr_kind)//'":'
            write(logunit(1),*) trim(attribute%label)

            if (attribute%kind == sdr_Boundary_object) then
              ! Add the unique calc_dist setting for this boundary.
              call append(attrList%bc_uni_calcdist, attribute%calc_dist)
            end if
          end if

        end select

        attrList%dynArray%val(attr_pos)%uni_id = uni_pos
        write(logunit(1),*) 'Unique identifier position of this attribute:', &
          &                 uni_pos

      else
        !do not count boundary Id counter in the attribute for periodic planes
        attrList%dynArray%val(attr_pos)%id = -1
      end if

    else
      write(logunit(1),*) 'This is an old attribute found at ', attr_pos
    end if

  end subroutine sdr_load_attribute
! ******************************************************************************!


  ! ****************************************************************************!
  !> This routine loads the distance refine table from attribute table and
  !! appends the loaded distance in radius and level to refine nodes within
  !! radius the given radius.
  !!
  !! to array of distance refine and store nVals in distance refine to
  !! distRefine_id in attribute to create unique list of attribute
  subroutine load_distanceRefine(attrList, conf, thandle, attribute)
    ! --------------------------------------------------------------------------!
    !> dynamic array of attribute type
    type(sdr_attrList_type), intent(inout) :: attrList
    type(flu_State) :: conf !< lua state
    integer, intent(in) :: thandle !<  current attribute handle
    !> save attibute info temporarily
    type(sdr_attribute_type), intent(inout) :: attribute
    ! --------------------------------------------------------------------------!
    integer :: dr_handle, dr_subhandle
    integer :: dr_length, iDistRefine
    ! --------------------------------------------------------------------------!
    write(logunit(1),*) 'Loading distance_refine'

    ! Open distance refine table
    call aot_table_open(L = conf, parent = thandle, thandle = dr_handle, &
      &                 key = 'distance_refine')

    ! load distance refine only if the table exists otherwise prompt
    ! warning message and set attribute%distRefine_id = 0
    if (dr_handle > 0) then
      ! set current attribute distance_refine start position in growing
      ! array of distance_refine
      attribute%distance_first = attrList%distRefine%nVals + 1

      ! Load distance refine sub handle to check between single table or multiple table
      call aot_table_open( L = conf, parent = dr_handle,   &
        &                  thandle = dr_subhandle, pos = 1 )
      if (dr_subhandle == 0) then
        ! distance refine is a single table
        call aot_table_close( L = conf, thandle = dr_subHandle )
        call load_distanceRefine_single(attrList  = attrList,  &
          &                             conf      = conf,      &
          &                             thandle   = dr_handle, &
          &                             attribute = attribute  )
      else
        ! distance refine is a multiple table
        call aot_table_close( L = conf, thandle = dr_subHandle )
        dr_length = aot_table_length(L = conf, thandle = dr_handle)
        do iDistRefine = 1, dr_length
          call aot_table_open( L = conf, parent = dr_handle,             &
            &                  thandle = dr_subhandle, pos = iDistRefine )
          call load_distanceRefine_single(attrList  = attrList,     &
            &                             conf      = conf,         &
            &                             thandle   = dr_subHandle, &
            &                             attribute = attribute     )
          call aot_table_close( L = conf, thandle = dr_subHandle )
        end do
      end if
      ! set nVals of distance_refine as distRefine_id to create
      ! unique list of attribute
      attribute%distRefine_id = attrList%distRefine%nVals

      ! set current attribute distance_refine last position in growing
      ! array of distance_refine
      attribute%distance_last = attrList%distRefine%nVals
    else
      write(logUnit(5),*) 'WARNING: No distance_refine defined for attribute'
      write(logUnit(5),*) '         boundary label: '//trim(attribute%label)

      ! no distance refine
      attribute%distRefine_id = 0
      attribute%distance_first = 0
      attribute%distance_last = 0
    end if

    call aot_table_close(L = conf, thandle = dr_handle)

  end subroutine load_distanceRefine
  ! ****************************************************************************!


  ! ****************************************************************************!
  !> load single distance refine table
  subroutine load_distanceRefine_single(attrList, conf, thandle, attribute)
    ! --------------------------------------------------------------------------!
    !> dynamic array of attribute type
    type(sdr_attrList_type), intent(inout) :: attrList
    type(flu_State) :: conf !< lua state
    integer, intent(in) :: thandle !<  handle to load distance refine from
    !> save attibute info temporarily
    type(sdr_attribute_type), intent(inout) :: attribute
    ! --------------------------------------------------------------------------!
    integer :: iError
    type(sdr_distanceRefine_type) :: distRefine
    integer :: level_offset
    ! --------------------------------------------------------------------------!

    ! load radius which defines the distance to apply the level
    call aot_get_val(L = conf, thandle = thandle, key = 'radius', &
      &              val = distRefine%radius, ErrCode = iError )

    if (btest(iError, aoterr_Fatal)) then
      write(logUnit(0),*) &
        &  'FATAL Error occured, while retrieving distance_refine radius:'
      if (btest(iError, aoterr_NonExistent)) &
        &  write(logunit(0),*) 'distance_refine radius not existent!'
      if (btest(iError, aoterr_WrongType)) &
        &  write(logunit(0),*)'distance_refine has wrong type!'
      call tem_abort()
    end if

    ! load level to refine node within given radius.
    ! if level is not defined, set default to current boundary level
    call aot_get_val(L = conf, thandle = thandle, key = 'level_offset', &
      &              val = level_offset, ErrCode = iError,              &
      &              default = 0)
    if (btest(iError, aoterr_WrongType)) then
      write(logunit(0),*) &
        &  'FATAL Error occured, while retrieving distance_refine &
        &  level_offset:'
      write(logunit(0),*) 'level_offset has wrong type, should be an Integer!'
      call tem_abort()
    end if

    if (level_offset > 0) then
      write(logUnit(0),*) 'ERROR: level_offset is greater than zero'
      write(logUnit(0),*) 'SOLUTION: level_offset must be <= 0'
      call tem_abort()
    end if

    ! compute reach level from level_offset
    distRefine%reach_level = max(attribute%level + level_offset,0)

    if (distRefine%reach_level==0) then
      write(logUnit(1),*) 'WARNING: Distance refine level is negative'
      write(logUnit(1),*) '         i.e level+offset_level < 0'
      write(logUnit(1),*) '         This distance refine will be omitted'
    end if

    ! append distance refine
    call append(me = attrList%distRefine, val = distRefine)

  end subroutine load_distanceRefine_single
  ! ****************************************************************************!


! ******************************************************************************!
  !> Identify all unique boundary colors, after all boundary attributes and
  !! seeds are known.
  !!
  !! This has to be done separately from reading the attributes, due to the
  !! interdependence of the boundaries and seed objects. It has to be ensured
  !! that all seed objects are known at this point.
  !!
  !! The bc_colors will correspond to the unique color name list, or are 0
  !! if the bc color does not exist in that list.
  !! 'none' bc colors are marked negatively as they are to be applied to all
  !! seed colors.
  subroutine sdr_identify_bc_colors(attribute)
    type(sdr_attrList_type), intent(inout) :: attribute

    integer :: nBCs
    integer :: colorpos
    integer :: iBC, iAttr
    character(len=labelLen) :: bc_color

    nBCs = attribute%kindpos(sdr_Boundary_Object)%nVals
    allocate(attribute%bc_color_id(nBCs))

    do iBC=1,nBCs
      iAttr = attribute%kindpos(sdr_Boundary_Object)%val(iBC)
      bc_color = attribute%dynArray%val(iAttr)%color
      if (trim(bc_color) == 'none') then
        attribute%bc_color_id(iBC) = -1
      else
        ! colorpos is 0 if the bc_color is not found in the uni_name list,
        ! otherwise it is the corresponding index in the list of unique color
        ! names.
        colorpos = PositionOfVal(me  = attribute%uni_name(sdr_Seed_Object), &
          &                      val = bc_color)
        attribute%bc_color_id(iBC) = colorpos
      end if
    end do

  end subroutine sdr_identify_bc_colors
  ! ****************************************************************************!


  ! ****************************************************************************!
  !> Look up colors that should be inverted and set their inversion flag
  !! accordingly.
  subroutine sdr_identify_inv_colors(attribute, invert_color)
    type(sdr_attrList_type), intent(inout) :: attribute
    type(dyn_labelArray_type), intent(in) :: invert_color

    integer :: nColors
    integer :: iColor
    character(len=labelLen) :: colname

    nColors = attribute%uni_name(sdr_seed_object)%nVals

    allocate(attribute%color_inverted(nColors))

    attribute%color_inverted = .false.

    do iColor=1,nColors
      colname = attribute%uni_name(sdr_seed_object)%val(iColor)
      if (trim(colname) /= 'none') then
        ! If the color is found in the list of colors to invert, mark it
        ! accordingly now.
        attribute%color_inverted(iColor) &
          &  = (PositionOfVal(invert_color, trim(colname)) > 0)
      end if
    end do

  end subroutine sdr_identify_inv_colors
  ! ****************************************************************************!


  ! ****************************************************************************!
  !> Returns if any bc attribute has the subresolution option set.
  !!
  !! The subresolution option only plays a role for boundary objects and only
  !! for those with a specific color (not none). This routine iterates over all
  !! boundary definitions and checks the subresolution flag. If any flag is
  !! true, the result of the function will be true.
  !! Only colors, which do have a corresponding seed are considered here.
  !! Note, that the subresolution flag will only be set for specific colors,
  !! thus, we do not need to check for 'none' colors here separately.
  function sdr_any_bc_subresolution(attribute) result(any_subres)
    !> The attribute list to check th boundaries in.
    type(sdr_attrList_type), intent(in) :: attribute

    !> Returned value, true if any boundary object has the subresolution flag
    !! set.
    logical :: any_subres

    integer :: nBCs
    integer :: iBC, iAttr

    nBCs = attribute%kindpos(sdr_Boundary_Object)%nVals
    any_subres = .false.

    do iBC=1,nBCs
      ! Only colors that actually do have an active seed (bc_color_id>0),
      ! need to be considered.
      if (attribute%bc_color_id(iBC) > 0) then
        iAttr = attribute%kindpos(sdr_Boundary_Object)%val(iBC)
        any_subres &
          & = (any_subres .or. attribute%dynArray%val(iAttr)%subresolution)
      end if
    end do

  end function sdr_any_bc_subresolution
  ! ****************************************************************************!


  ! ****************************************************************************!
  !> Returns if any bc attribute has the distance refine option with
  !! reach_level>0.
  !!
  !! The distance refine option only plays a role for boundary objects.
  !! This routine iterates over all
  !! boundary definitions and checks the distRefine_id>0 and
  !! distRefine%reach_level>0.
  function sdr_any_bc_distanceRefine(attribute) result(any_distRef)
    !> The attribute list to check th boundaries in.
    type(sdr_attrList_type), intent(in) :: attribute

    !> Returned value, true if any boundary object has the
    !! valid distance refine
    logical :: any_distRef

    integer :: nBCs
    integer :: iBC, iAttr, distRef_first, distRef_last

    nBCs = attribute%kindpos(sdr_Boundary_Object)%nVals
    any_distRef = .false.

    do iBC=1,nBCs
      iAttr = attribute%kindpos(sdr_Boundary_Object)%val(iBC)
      if (attribute%dynArray%val(iAttr)%distRefine_id > 0) then
        distRef_first = attribute%dynArray%val(iAttr)%distance_first
        distRef_last = attribute%dynArray%val(iAttr)%distance_last

        any_distRef = ( any_distRef .or.                             &
          & any(attribute%distRefine%val(distRef_first:distRef_last) &
          &                         %reach_level>0) )
      end if
    end do

  end function sdr_any_bc_distanceRefine
  ! ****************************************************************************!


  ! ****************************************************************************!
  !> Returns if periodic bc attribute is defined
  function sdr_isPeriodicDefined(attribute) result(isPeriodic)
    !> The attribute list to check th boundaries in.
    type(sdr_attrList_type), intent(in) :: attribute

    !> Returned value, true if periodic attribute is defined
    logical :: isPeriodic

    integer :: nAttrs

    nAttrs = attribute%dynArray%nVals
    isPeriodic = any(attribute%dynArray%val(1:nAttrs)%id == -1)
  end function sdr_isPeriodicDefined
  ! ****************************************************************************!

  ! ------------------------------------------------------------------------ !
  !> Get all attributes of the given object kind and (unique) id from the
  !! list of attributes in "attribute".
  function sdr_attr_of_uni(attribute, object_kind, id) result(allattr)
    ! -------------------------------------------------------------------- !
    !> The list of attributes to look up the unique entries in.
    type(sdr_attrList_type), intent(in) :: attribute

    !> The kind of objects to look for (boundaries, seeds, refinements...)
    integer, intent(in) :: object_kind

    !> The unique id, to look for (all attributes with this id will be
    !! returned).
    integer, intent(in) :: id

    !> Resulting list of attributes with the same unique id.
    type(grw_IntArray_type) :: allattr
    ! -------------------------------------------------------------------- !
    integer :: iAttr
    ! -------------------------------------------------------------------- !

    call init(allattr)

    do iAttr=1,attribute%dynArray%nVals
      if (attribute%dynArray%val(iAttr)%kind == object_kind) then
        if (attribute%dynArray%val(iAttr)%uni_id == id) then
          call append(me = allattr, val = iAttr)
        end if
      end if
    end do

  end function sdr_attr_of_uni
  ! ------------------------------------------------------------------------ !
  ! ------------------------------------------------------------------------ !


  subroutine init_ga_distancerefine(me, length)
    type(grw_distancerefinearray_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_distancerefine

  subroutine destroy_ga_distancerefine(me)
    type(grw_distancerefinearray_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_distancerefine


  subroutine truncate_ga_distancerefine(me)
    !------------------------------------------------------------------------
    type(grw_distancerefinearray_type) :: me !< array to truncate
    !------------------------------------------------------------------------
    type(sdr_distancerefine_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_distancerefine


  subroutine empty_ga_distancerefine(me)
    !------------------------------------------------------------------------
    type(grw_distancerefinearray_type) :: me !< array to sorttruncate
    !------------------------------------------------------------------------

    me%nvals = 0

  end subroutine empty_ga_distancerefine

  !> 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_distancerefine(me, val, pos, length)
    type(grw_distancerefinearray_type) :: me !< array to place the value into
    type(sdr_distancerefine_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_distancerefine


  !> 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_distancerefine_vec(me, val, pos, length)
    type(grw_distancerefinearray_type) :: me !< array to append the value to
    type(sdr_distancerefine_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_distancerefine_vec


  subroutine append_ga_distancerefine(me, val, length)
    type(grw_distancerefinearray_type) :: me !< array to append the value to
    type(sdr_distancerefine_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_distancerefine

  subroutine append_ga_distancerefine_vec(me, val, length)
    type(grw_distancerefinearray_type) :: me !< array to append the value to
    type(sdr_distancerefine_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_distancerefine_vec


  subroutine expand_ga_distancerefine(me, pos, length)
    type(grw_distancerefinearray_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_distancerefine_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_distancerefine

! ***************************************************************************** !
  !> initialization of a dynamic array
  !!
  !! before a dynamic array can be used, it has to be initialized
  !! with this routine. the initial length provided here, can
  !! avoid reallocations and memory copying, if approximated
  !! correctly enough. if none is specified, the provided container
  !! initially will be of size 0.
  subroutine init_da_attribute(me, length)
    !-----------------------------------------------------------------
    type(dyn_attributearray_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)
    if( allocated( me%sorted ) ) deallocate(me%sorted)
    ! ... and reallocate
    allocate(me%val(me%containersize))
    allocate(me%sorted(me%containersize))
    me%nvals = 0

  end subroutine init_da_attribute

  !> destruction of a dynamic array
  !!
  !! this subroutine takes care of a proper destruction of a
  !! dynamic array, it frees the allocated memory and resets
  !! the internal counts to 0.
  subroutine destroy_da_attribute(me)
    type(dyn_attributearray_type), intent(inout) :: me !< dynamic array to init

    me%containersize = 0
    me%nvals         = 0

    if( allocated( me%val ) ) deallocate(me%val)
    if( allocated( me%sorted ) ) deallocate(me%sorted)
  end subroutine destroy_da_attribute
! ***************************************************************************** !


! ***************************************************************************** !
  !> appending a value to the dynamic array
  !!
  !! with this subroutine, a given value can be added to the
  !! dynamic array. the actual position of this value in the
  !! dynamic array will be returned, so it can be found again
  !! easily later. with the wasadded flag, it is indicated,\n
  !! wasadded = true,  if this entry had to be added,\n
  !! wasadded = false, if this was already found in the array.
  subroutine append_da_attribute(me, val, length, pos, wasadded )
    !------------------------------------------------------------------------
    type(dyn_attributearray_type) :: me   !< array to append the value to
    type(sdr_attribute_type), intent(in)     :: val  !< value to append
    !> optional length to expand the array
    integer, intent(in), optional :: length
    !> position in the array, if the value is found
    integer, intent(out), optional :: pos
    !> flag to indicate, if val was newly added
    logical, intent(out), optional :: wasadded
    !------------------------------------------------------------------------
    integer :: foundpos
    integer :: i
    !------------------------------------------------------------------------

    ! do a binary search on existing entries (returns closest entry next to
    ! it if not found).
    foundpos = sortedposofval(me, val, .true.)
    if( present( wasadded ) ) wasadded = .false.

    ! if it found the value, the position is smaller than nvals
    if (foundpos <= me%nvals) then

      ! the returned position might actually be the right entry already or
      ! not, check for this here.
      if ( me%val(me%sorted(foundpos)) == val ) then

        ! found the value in a list of unique values,
        ! nothing to do, just return its position.
        if( present( pos ) ) pos = me%sorted(foundpos)

      else

        ! need to append a new value!

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

        if( present( wasadded ) ) wasadded = .true.
        if (me%nvals == me%containersize) then

          ! container is full, need to expand it
          call expand(me = me, length = length)
        end if
        me%nvals = me%nvals + 1

        ! put the new value into the last position in the
        ! array.
        me%val(me%nvals) = val
        do while( foundpos < me%nvals )
          if(me%val(me%sorted(foundpos)) /= val) then
            exit
          end if
          ! in case of multiple entries with the same value
          ! move on to the first differing entry.
          foundpos = foundpos + 1
        end do
        ! shift the sorted list of indices, to create a
        ! whole for the value to be inserted, at position
        ! foundpos.
        do i=me%nvals-1,foundpos,-1
          me%sorted(i+1) = me%sorted(i)
        end do
        ! put the index of the new value into the
        ! sorted list at the now freed position.
        me%sorted(foundpos) = me%nvals
        if( present( pos ) ) pos = me%nvals

      end if

    else

      ! 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( present( wasadded ) ) wasadded = .true.
      if (foundpos > me%containersize) then
        ! expand the array, if its boundary is reached
        call expand(me = me, length = length)
      end if
      me%nvals = foundpos
      me%val(foundpos) = val
      me%sorted(foundpos) = foundpos
      if( present( pos ) ) pos = foundpos

    end if

  end subroutine append_da_attribute
! ***************************************************************************** !


! ***************************************************************************** !
  !> appending a sorted list of values to the dynamic array
  !!
  !! with this subroutine, a given list of sorted values can be added to the
  !! dynamic array. the actual positions of these values in the
  !! dynamic array will be returned, so it can be found again
  !! easily later. with the wasadded flag, it is indicated,\n
  !! wasadded = true,  if this entry had to be added,\n
  !! wasadded = false, if this was already found in the array.
  subroutine append_da_vecattribute(me, val, length, pos, wasadded )
    !------------------------------------------------------------------------
    type(dyn_attributearray_type) :: me   !< array to append the value to
    type(sdr_attribute_type), intent(in)     :: val(:)  !< values to append
    !> optional length to expand the array
    integer, intent(in), optional :: length
    !> position in the array, the values are found at.
    integer, intent(out), optional :: pos(:)
    !> flag to indicate, if val was newly added
    logical, intent(out), optional :: wasadded(:)
    !------------------------------------------------------------------------
    type(sdr_attribute_type) :: lastval
    logical :: addedval(size(val))
    integer :: i
    integer :: veclen
    integer :: maxlen
    integer :: nappend
    integer :: rem_app
    integer :: curval, ival, iold, iadd
    integer, allocatable :: newsorted(:)
    !------------------------------------------------------------------------

    if (size(val) == 0) return

    veclen = size(val)
    maxlen = veclen + me%nvals

    allocate(newsorted(maxlen))

    addedval = .false.

    iold = 1
    iadd = 1

    nappend = 0
    curval = 0

    ! select the first entry before the loop unconditionally without checks
    ! for uniqueness (nothing to check against yet).
    if ( me%val(me%sorted(iold)) <= val(iadd) ) then
      curval = curval + 1
      newsorted(curval) = me%sorted(iold)
      lastval = me%val(me%sorted(iold))
      iold = iold + 1
    else
      curval = curval + 1
      nappend = nappend + 1
      newsorted(curval) = me%nvals + nappend
      lastval = val(iadd)
      if (present(pos))  pos(iadd) = newsorted(curval)
      addedval(iadd) = .true.
      iadd = iadd + 1
    end if

    do ival=2,maxlen

      if ( (iadd <= veclen) .and. (iold <= me%nvals) ) then

        if ( me%val(me%sorted(iold)) <= val(iadd) ) then

          ! the original list's values are appended to newsorted before
          ! the additional list is appended.
          curval = curval + 1
          newsorted(curval) = me%sorted(iold)
          lastval = me%val(me%sorted(iold))
          iold = iold + 1

        else

          ! only append the value to unique lists, if it is not yet in the list.
          ! (if it is already in the list, it has to be the previous (curval-1)
          !  entry.)
          if ( lastval < val(iadd) ) then
            nappend = nappend + 1
            curval = curval + 1
            newsorted(curval) = me%nvals + nappend
            lastval = val(iadd)
            addedval(iadd) = .true.
          end if
          if (present(pos)) pos(iadd) = newsorted(curval)
          iadd = iadd + 1

        end if

      else

        ! reached the end of one or both of the sorted lists.
        exit

      end if

    end do

    if (iold <= me%nvals) then
      ! still some values from the original list left.
      newsorted(curval+1:me%nvals+nappend) = me%sorted(iold:me%nvals)
    end if

    if (iadd <= veclen) then
      ! still some values from the list to append left.
      rem_app = iadd
      do i = rem_app,veclen
        if ( lastval < val(iadd) ) then
          nappend = nappend + 1
          curval = curval + 1
          newsorted(curval) = me%nvals + nappend
          lastval = val(iadd)
          addedval(iadd) = .true.
        end if
        if (present(pos)) pos(iadd) = newsorted(curval)
        iadd = iadd + 1
      end do
    end if

    if (me%nvals > huge(me%nvals)-nappend) then
       write(*,*) "reached end of integer range for dynamic array!"
       write(*,*) "aborting!!"
       stop
    end if

    if (me%nvals + nappend > me%containersize) then
      call expand( me        = me,      &
        &          increment = nappend, &
        &          length    = length   )
    end if
    me%sorted(:me%nvals+nappend) = newsorted(:me%nvals+nappend)
    curval = me%nvals
    do iadd=1,veclen
      if (addedval(iadd)) then
        curval = curval + 1
        me%val(curval) = val(iadd)
      end if
    end do
    me%nvals = me%nvals + nappend

    if( present( wasadded ) ) wasadded = addedval

  end subroutine append_da_vecattribute
! ***************************************************************************** !


! ***************************************************************************** !
  !> truncate the array after the last valid entry and hence cut off the empty
  !! trailing empty entries
  !!
  subroutine truncate_da_attribute(me)
    !------------------------------------------------------------------------
    type(dyn_attributearray_type) :: me !< array to sorttruncate
    !------------------------------------------------------------------------
    type(sdr_attribute_type), allocatable :: swpval(:)
    integer, allocatable :: swpsort(:)
    !------------------------------------------------------------------------

    if (me%nvals < me%containersize) then
      allocate(swpval(me%nvals))
      allocate(swpsort(me%nvals))

      swpval = me%val(:me%nvals)
      swpsort = me%sorted(:me%nvals)

      call move_alloc(swpval, me%val)
      call move_alloc(swpsort, me%sorted)

      me%containersize = me%nvals
    end if

  end subroutine truncate_da_attribute
! ***************************************************************************** !


! ***************************************************************************** !
  !> empty all contents of the array without changing the size or status of any
  !! array
  !!
  subroutine empty_da_attribute(me)
    !------------------------------------------------------------------------
    type(dyn_attributearray_type) :: me !< array to sorttruncate
    !------------------------------------------------------------------------
    ! reset the number of entries
    me%nvals = 0
  end subroutine empty_da_attribute
! ***************************************************************************** !


! ***************************************************************************** !
  !> fixing the dynamic array
  !!
  !! truncate the array after the last valid entry and hence cut off the empty
  !! trailing empty entries
  !! store the array in the sorted order according to the sorted( ) array
  !!
  subroutine sorttruncate_da_attribute(me)
    !------------------------------------------------------------------------
    type(dyn_attributearray_type) :: me !< array to sorttruncate
    !------------------------------------------------------------------------
    type(dyn_attributearray_type) :: tarray !< temporary array
    integer :: ival
    integer :: dpos
    !------------------------------------------------------------------------
    ! allocate the temporary array
    call init( me = tarray, length = me%nvals )
    ! copy the entries in a sorted fashion into the temporary array
    do ival = 1, me%nvals
      call append( me = tarray, val = me%val( me%sorted( ival )), &
           &       pos = dpos)
    enddo
    call destroy( me = me )

    me = tarray
    call destroy( me = tarray )

  end subroutine sorttruncate_da_attribute
! ***************************************************************************** !


! ***************************************************************************** !
  !> expanding the dynamic array
  !!
  !! this is a helping subroutine, which doubles the container
  !! of the given dynamic array. as the container might be
  !! initially 0-sized, a module variable minlength has been introduced, which
  !! is used here, to at least create a container of this size.
  subroutine expand_da_attribute(me, increment, length)
    !------------------------------------------------------------------------
    type(dyn_attributearray_type) :: me !< array to resize
    integer, optional :: increment !< used for vector append
    !> optional length to expand the array
    integer, intent(in), optional :: length
    !------------------------------------------------------------------------
    type(sdr_attribute_type), allocatable :: swpval(:)
    integer, allocatable :: swpsort(:)
    !------------------------------------------------------------------------
    integer :: addvals, explen
    !------------------------------------------------------------------------

    addvals = 1
    if (present(increment)) addvals = increment

    if (addvals > 0) then

      ! if length is present, use that, otherwise double the size
      if( present( length ) ) then
        explen = length
      else
        ! set the global minimum length, if doubling would be smaller than that
        explen = max(me%containersize, minlength)
      end if

      ! check whether all elements will fit
      if( addvals > explen ) then
        explen = addvals
      end if

      ! check whether the new size will exceed the max container size.
      if( (huge(me%containersize) - explen) <= me%containersize ) then
        ! if so, expand to the maximum size
        me%containersize = huge(me%containersize)
      else
        ! if not, expand to the calculated size
        me%containersize = me%containersize + explen
      end if

      ! only need to do something, if there are actually values to append.
      if (me%nvals > 0) then

        allocate(swpval(me%containersize))
        swpval(1:me%nvals) = me%val(1:me%nvals)
        call move_alloc( swpval, me%val )

        allocate(swpsort(me%containersize))
        swpsort(1:me%nvals) = me%sorted(1:me%nvals)
        call move_alloc( swpsort, me%sorted )

      else ! me%nvals == 0

        if( allocated(me%val) ) &
          deallocate(me%val)
        allocate(me%val(me%containersize))
        if( allocated(me%sorted) ) &
          deallocate(me%sorted)
        allocate(me%sorted(me%containersize))

      end if

    end if

  end subroutine expand_da_attribute
! ***************************************************************************** !


! ***************************************************************************** !
  !> return the sorted position of a value in the given dynamic array
  !!
  !! if the value was not found,
  !!  - return 0 if nextifnotfound = .false.
  !!  - return position at the end if nextifnotfound = .true.
  function sortposofval_attribute(me, val, nextifnotfound, lower, upper) result(pos)
    !------------------------------------------------------------------------
    type(dyn_attributearray_type), intent(in) :: me !< dynamic array
    type(sdr_attribute_type), intent(in) :: val !< value to look for
    !> flag to indicate, if the next entry in the list should be returned,
    !! if the searched one is not found.
    logical, intent(in), optional :: nextifnotfound
    integer, intent(in), optional :: lower !< lower search limit
    integer, intent(in), optional :: upper !< upper search limit
    integer :: pos !< position of val in the sorted list, 0 if not found
    !------------------------------------------------------------------------
    logical :: retnext
    integer :: lb, ub
    integer :: mid
    type(sdr_attribute_type) :: lb_val, ub_val
    type(sdr_attribute_type) :: mid_val
    !------------------------------------------------------------------------

    retnext = .false.
    if (present(nextifnotfound)) retnext = nextifnotfound

    lb = 1
    ub = me%nvals

    if( present( lower ) ) lb = lower
    if( present( upper ) ) ub = upper

    pos = 0
    if (retnext) pos = lb

    !> binary search on sorted list
    do while(ub >= lb)
      lb_val = me%val(me%sorted(lb))

      if (val < lb_val) then
        if (retnext) pos = lb
        exit
      end if

      ub_val = me%val(me%sorted(ub))

      if (val > ub_val) then
        if (retnext) pos = ub+1
        exit
      end if

      ! safe guard against integer limit overflow
      mid = lb + (ub-lb) / 2
      mid_val = me%val(me%sorted(mid))
      if (val == mid_val) then
        pos = mid
        exit
      end if
      if (val > mid_val) then
        lb = mid + 1
      else
        ub = mid - 1
      end if
    end do
  end function sortposofval_attribute
! ***************************************************************************** !


! ***************************************************************************** !
  !> the actual position of a given value in the dynamic array
  !!
  !! most likely this is what you need in codes, using this
  !! data structure, it first does the binary search on the sorted
  !! values with sortposofval_attribute and then returns the looked
  !! up position in the original unsorted array, which corresponds
  !! to the position returned by the append routine.
  function posofval_attribute(me, val, nextifnotfound, lower, upper) result(pos)
    !------------------------------------------------------------------------
    type(dyn_attributearray_type), intent(in) :: me !< dynamic array
    type(sdr_attribute_type), intent(in) :: val !< value to search for
    !> flag to indicate, if the position of the next entry in the sorted
    !! list should be returned instead, if val is not found.
    logical, intent(in), optional :: nextifnotfound
    integer, intent(in), optional :: lower !< lower search limit
    integer, intent(in), optional :: upper !< upper search limit
    integer :: pos !< position in the array of the searche value, 0 if not found
    !------------------------------------------------------------------------
    integer :: sortpos
    integer :: lb, ub
    !------------------------------------------------------------------------

    lb = 1
    ub = me%nvals

    if( present( lower ) ) lb = lower
    if( present( upper ) ) ub = upper

    sortpos = sortedposofval(me, val, nextifnotfound, lb, ub)

    ! if result (sorted pos)
    if ((sortpos <= me%nvals) .and. (sortpos > 0)) then
      pos = me%sorted(sortpos)
    else
      pos = sortpos
    end if
  end function posofval_attribute
! ***************************************************************************** !

! ******************************************************************************!
  !> This function provides the test for equality of two attributes.
  !!
  !! Two attributes are considered to be equal if their kind, label, level,
  !! color, subresolution and distancerefine are all equal.
  function isEqual(left, right) result(equality)
    ! ---------------------------------------------------------------------------
    !> attribute to compare
    type(sdr_attribute_type), intent(in) :: left
    !> attribute to compare against
    type(sdr_attribute_type), intent(in) :: right
    !> is equal??
    logical :: equality
    ! ---------------------------------------------------------------------------

    equality = ( left%kind == right%kind ) &
      &      .and. ( trim(left%label) == trim(right%label) ) &
      &      .and. ( trim(left%color) == trim(right%color) ) &
      &      .and. ( left%level == right%level ) &
      &      .and. ( left%distRefine_id == right%distRefine_id ) &
      &      .and. ( left%subresolution .eqv. right%subresolution )

  end function isEqual
! ******************************************************************************!



! ******************************************************************************!
  !> This function provides the test for unequality of two attributes.
  !!
  !! Two attributes are considered unequal if any of their kind, label, level,
  !! color, subresolution or distancerefine are not equal.
  function isUnequal(left, right) result(unequality)
    ! ---------------------------------------------------------------------------
    !> attribute to compare
    type(sdr_attribute_type), intent(in) :: left
    !> attribute to compare against
    type(sdr_attribute_type), intent(in) :: right
    !> is unequal??
    logical :: unequality
    ! ---------------------------------------------------------------------------

    unequality = ( left%kind /= right%kind ) &
      &        .or. ( trim(left%label) /= trim(right%label) ) &
      &        .or. ( trim(left%color) /= trim(right%color) ) &
      &        .or. ( left%level /= right%level ) &
      &        .or. ( left%distRefine_id /= right%distRefine_id ) &
      &        .or. ( left%subresolution .neqv. right%subresolution )

  end function isUnequal
! ******************************************************************************!


! ******************************************************************************!
  !> This function provides a comparison of two attributes.
  !!
  !! Check if left attribute is smaller than right attribute
  !! by checking smaller for attribute kind, label color and level
  function isSmaller(left, right) result(small)
    ! ---------------------------------------------------------------------------
    !> attribute to compare
    type(sdr_attribute_type), intent(in) :: left
    !> attribute to compare against
    type(sdr_attribute_type), intent(in) :: right
    !> is smaller??
    logical :: small
    ! ---------------------------------------------------------------------------

    small = .false.
    if (left%kind < right%kind) then
      small = .true.
    else if (left%kind == right%kind) then
      if (trim(left%label) < trim(right%label)) then
        small = .true.
      else if (trim(left%label) == trim(right%label)) then
        if (trim(left%color) < trim(right%color)) then
          small = .true.
        else if (trim(left%color) == trim(right%color)) then
          if (left%level < right%level) then
            small = .true.
          else if (left%level == right%level) then
            if (left%distRefine_id < right%distRefine_id) then
              small = .true.
            else if (left%distRefine_id == right%distRefine_id) then
              small = ( (left%subresolution .neqv. right%subresolution) &
                &       .and. right%subresolution )
            end if
          end if
        end if
      end if
    end if

  end function isSmaller
! ******************************************************************************!


! ******************************************************************************!
  !> This function provides a comparison of two attributes.
  !!
  !! Check if left attribute is smaller or equal than right attribute
  !! by checking smaller or equal for attribute kind, label color and level
  function isSmallerOrEqual(left, right) result(small)
    ! ---------------------------------------------------------------------------
    !> attribute to compare
    type(sdr_attribute_type), intent(in) :: left
    !> attribute to compare against
    type(sdr_attribute_type), intent(in) :: right
    !> is smaller??
    logical :: small
    ! ---------------------------------------------------------------------------

    small = .false.
    if (left%kind < right%kind) then
      small = .true.
    else if (left%kind == right%kind) then
      if (trim(left%label) < trim(right%label)) then
        small = .true.
      else if (trim(left%label) == trim(right%label)) then
        if (trim(left%color) < trim(right%color)) then
          small = .true.
        else if (trim(left%color) == trim(right%color)) then
          if (left%level < right%level) then
            small = .true.
          else if (left%level == right%level) then
            if (left%distRefine_id < right%distRefine_id) then
              small = .true.
            else if (left%distRefine_id == right%distRefine_id) then
              small = ( (left%subresolution .eqv. right%subresolution) &
                &       .or. right%subresolution )
            end if
          end if
        end if
      end if
    end if

  end function isSmallerOrEqual
! ******************************************************************************!


! ******************************************************************************!
  !> This function provides a comparison of two attributes.
  !!
  !! Check if left attribute is greater than right attribute
  !! by checking greater for attribute kind, label, color and level
  function isGreater(left, right) result(great)
    ! ---------------------------------------------------------------------------
    !> attribute to compare
    type(sdr_attribute_type), intent(in) :: left
    !> attribute to compare against
    type(sdr_attribute_type), intent(in) :: right
    !> is greater??
    logical :: great
    ! ---------------------------------------------------------------------------

    great = .false.
    if (left%kind > right%kind) then
      great = .true.
    else if (left%kind == right%kind) then
      if (trim(left%label) > trim(right%label)) then
        great = .true.
      else if (trim(left%label) == trim(right%label)) then
        if (trim(left%color) > trim(right%color)) then
          great = .true.
        else if (trim(left%color) == trim(right%color)) then
          if (left%level > right%level) then
            great = .true.
          else if (left%level == right%level) then
            if (left%distRefine_id > right%distRefine_id) then
              great = .true.
            else if (left%distRefine_id == right%distRefine_id) then
              great = ( (left%subresolution .neqv. right%subresolution) &
                &       .and. left%subresolution )
            end if
          end if
        end if
      end if
    end if

  end function isGreater
! ******************************************************************************!


! ******************************************************************************!
  !> This function provides a comparison of two attributes.
  !!
  !! Check if left attribute is greater or equal than right attribute
  !! by checking greater for attribute kind, label, color and level
!  function isGreaterOrEqual(left, right) result(great)
!    ! ---------------------------------------------------------------------------
!    !> attribute to compare
!    type(sdr_attribute_type), intent(in) :: left
!    !> attribute to compare against
!    type(sdr_attribute_type), intent(in) :: right
!    !> is greater??
!    logical :: great
!    ! ---------------------------------------------------------------------------
!
!    great = .false.
!    if (left%kind > right%kind) then
!      great = .true.
!    else if (left%kind == right%kind) then
!      if (trim(left%label) > trim(right%label)) then
!        great = .true.
!      else if (trim(left%label) == trim(right%label)) then
!        if (trim(left%color) > trim(right%color)) then
!          great = .true.
!        else if (trim(left%color) == trim(right%color)) then
!          if (left%level > right%level) then
!            great = .true.
!          else if (left%level == right%level) then
!            if (left%distRefine_id > right%distRefine_id) then
!              great = .true.
!            else if (left%distRefine_id == right%distRefine_id) then
!              great = ( (left%subresolution .eqv. right%subresolution) &
!                &       .or. left%subresolution )
!            end if
!          end if
!        end if
!      end if
!    end if
!
!  end function isGreaterOrEqual
! ******************************************************************************!

end module sdr_attribute_module

!> \page attribute Attribute
!! Attribute defines the property for the geometries defined in the
!! spatial object. Following attribute kinds are supported,
!! - \ref boundary "Boundary"
!! - \ref seed "Seed"
!! - \ref refinement "Refinement".
!! \n
!! Attribute list is dynamic unique list. Thus multiple definitions of
!! the same attribute in the configuration file does not result in multiple
!! storages.
