text:
!##############################################################################
!# ****************************************************************************
!# <name> hadaptaux </name>
!# ****************************************************************************
!#
!# <purpose>
!#
!# WARNING: Do not USE this module in your applications unless you really
!#          know what you are doing. This module does no error checking!!!
!#
!# This module contains the basic data structures and auxiliary routines
!# which are required for performing h-adaptivity in.
!#
!# The following routines are available:
!#
!#  1.) hadapt_getNVE / hadapt_getNNVE
!#      -> return the (maximum ) number of vertices per element
!#
!#  2.) hadapt_setVertexCoords1D
!#      -> Set the coordinates of vertices to the adaptivity structure in 1D
!#
!#  3.) hadapt_getVertexCoords1D
!#      -> Get the coordinates of vertices from the adaptivity structure in 1D
!#
!#  4.) hadapt_setVertexCoords2D
!#      -> Set the coordinates of vertices to the adaptivity structure in 2D
!#
!#  5.) hadapt_getVertexCoords2D
!#      -> Get the coordinates of vertices from the adaptivity structure in 2D
!#
!#  6.) hadapt_setVertexCoords3D
!#      -> Set the coordinates of vertices to the adaptivity structure in 3D
!#
!#  7.) hadapt_getVertexCoords3D
!#      -> Get the coordinates of vertices from the adaptivity structure in 3D
!#
!#  8.) hadapt_setVerticesAtElement
!#      -> Set the "vertices-at-element" structure to the adaptivity structure
!#
!#  9.) hadapt_getVerticesAtElement
!#      -> Get the "vertices-at-element" structure from the adaptivity structure
!#
!# 10.) hadapt_setNeighboursAtElement
!#      -> Set the "neighbours-at-element" structure to the adaptivity structure
!#
!# 11.) hadapt_getNeighboursAtElement
!#      -> Get the "neighbours-at-element" structure from the adaptivity structure
!#
!# 12.) hadapt_setNelOfType
!#      -> Set the "InelOfType" array to the adaptivity structure
!#
!# 13.) hadapt_getNelOfType
!#      -> Get the "InelOfType" array from the adaptivity structure
!#
!# 14.) hadapt_setBoundary
!#      -> Set the boundary structure to the adaptivity structure
!#
!# 15.) hadapt_getBoundary
!#      -> Get the boundary structure form the adaptivity structure
!#
!# 16.) hadapt_setNodalProperty
!#      -> Set the "nodal property" list to the adaptivity structure
!#
!# 17.) hadapt_getNodalProperty
!#      -> Get the "nodal property" list from the adaptivity structure
!#
!# 18.) hadapt_genElementsAtVertex
!#      -> Generate the "elements-at-vertex" data structure
!#
!# </purpose>
!##############################################################################

module hadaptaux

!$ use omp_lib
  use arraylistInt
  use basicgeometry
  use fsystem
  use genoutput
  use linearalgebra
  use mapInt_DP
  use octreeDP
  use quadtreeDP
  use sort
  use storage
  use triangulation

  implicit none

  private
  public :: t_hadapt
  public :: hadapt_genElementsAtVertex
  public :: hadapt_getSubdivisonLevel
  public :: hadapt_getBoundary
  public :: hadapt_getNVE
  public :: hadapt_getNNVE
  public :: hadapt_getNeighboursAtElement
  public :: hadapt_getNelOfType
  public :: hadapt_getNodalProperty
  public :: hadapt_getVertexCoords1D
  public :: hadapt_getVertexCoords2D
  public :: hadapt_getVertexCoords3D
  public :: hadapt_getVerticesAtElement
  public :: hadapt_setBoundary
  public :: hadapt_setNeighboursAtElement
  public :: hadapt_setNelOfType
  public :: hadapt_setNodalProperty
  public :: hadapt_setVertexCoords1D
  public :: hadapt_setVertexCoords2D
  public :: hadapt_setVertexCoords3D
  public :: hadapt_setVerticesAtElement

!<constants>

!<constantblock description="Global flags for mesh refinement/coarsening">

  ! No refinement and coarsening
  integer, parameter, public :: HADAPT_NOADAPTATION = 0

  ! Red-Green refinement and coarsening strategy (R. Banks et al.)
  integer, parameter, public :: HADAPT_REDGREEN     = 1

  ! Longest edge bisection strategy (M. Rivara)
  integer, parameter, public :: HADAPT_LONGESTEDGE  = 2

!</constantblock>

!<constantblock description="Bitfield identifiers for state of adaptation">

  ! Adaptation is undefined
  integer(I32), parameter, public :: HADAPT_NONE             = 0_I32

  ! Parameters of adaptivity structure are initialised
  integer(I32), parameter, public :: HADAPT_HAS_PARAMETERS   = 2_I32**0

  ! Map/quadtree/octree for vertex coordinates is generated
  integer(I32), parameter, public :: HADAPT_HAS_COORDS       = 2_I32**1

  ! Array for IverticesAtElement is generated
  integer(I32), parameter, public :: HADAPT_HAS_VERTATELEM   = 2_I32**2

  ! Array for IneighboursAtElement is generated
  integer(I32), parameter, public :: HADAPT_HAS_NEIGHATELEM  = 2_I32**3

  ! Array for ImidneighboursAtElement is generated
  integer(I32), parameter, public :: HADAPT_HAS_MIDNEIGH     = 2_I32**4

  ! Boundary data is generated
  integer(I32), parameter, public :: HADAPT_HAS_BOUNDARY     = 2_I32**5

  ! Nodal property is generated
  integer(I32), parameter, public :: HADAPT_HAS_NODALPROP    = 2_I32**6

  ! Number of elements for predefined type
  integer(I32), parameter, public :: HADAPT_HAS_NELOFTYPE    = 2_I32**7

  ! Array for IelementsAtVertex is generated
  integer(I32), parameter, public :: HADAPT_HAS_ELEMATVERTEX = 2_I32**8

  ! Dynamic data structures in 1D are all generated
  integer(I32), parameter, public :: HADAPT_HAS_DYNAMICDATA1D  = HADAPT_HAS_PARAMETERS+&
                                                                 HADAPT_HAS_COORDS+&
                                                                 HADAPT_HAS_VERTATELEM+&
                                                                 HADAPT_HAS_NEIGHATELEM+&
                                                                 HADAPT_HAS_NODALPROP+&
                                                                 HADAPT_HAS_NELOFTYPE+&
                                                                 HADAPT_HAS_ELEMATVERTEX

  ! Dynamic data structures in 2D are all generated
  integer(I32), parameter, public :: HADAPT_HAS_DYNAMICDATA2D  = HADAPT_HAS_PARAMETERS+&
                                                                 HADAPT_HAS_COORDS+&
                                                                 HADAPT_HAS_VERTATELEM+&
                                                                 HADAPT_HAS_NEIGHATELEM+&
                                                                 HADAPT_HAS_MIDNEIGH+&
                                                                 HADAPT_HAS_NODALPROP+&
                                                                 HADAPT_HAS_NELOFTYPE+&
                                                                 HADAPT_HAS_ELEMATVERTEX

  ! Dynamic data structures in 3D are all generated
  integer(I32), parameter, public :: HADAPT_HAS_DYNAMICDATA3D  = HADAPT_HAS_PARAMETERS+&
                                                                 HADAPT_HAS_COORDS+&
                                                                 HADAPT_HAS_VERTATELEM+&
                                                                 HADAPT_HAS_NEIGHATELEM+&
                                                                 HADAPT_HAS_NODALPROP+&
                                                                 HADAPT_HAS_NELOFTYPE+&
                                                                 HADAPT_HAS_ELEMATVERTEX

  ! Cells are marked for refinement
  integer(I32), parameter, public :: HADAPT_MARKEDREFINE     = 2_I32**9

  ! Cells are marked for coarsening
  integer(I32), parameter, public :: HADAPT_MARKEDCOARSEN    = 2_I32**10

  ! Cells are marked
  integer(I32), parameter, public :: HADAPT_MARKED           = HADAPT_MARKEDREFINE+&
                                                               HADAPT_MARKEDCOARSEN

!</constantblock>

!<constantblock description="Global constants for grid modification operations">

  ! Operation identifier for initialization of callback function
  integer, parameter, public :: HADAPT_OPR_INITCALLBACK     = -1

  ! Operation identifier for finalization of callback function
  integer, parameter, public :: HADAPT_OPR_DONECALLBACK     = -2

  ! Operation identifier for adjustment of vertex dimension
  integer, parameter, public :: HADAPT_OPR_ADJUSTVERTEXDIM  = 1

  ! Operation identifier for vertex insertion at edge midpoint
  integer, parameter, public :: HADAPT_OPR_INSERTVERTEXEDGE = 2

  ! Operation identifier for vertex insertion at element centroid
  integer, parameter, public :: HADAPT_OPR_INSERTVERTEXCENTR= 3

  ! Operation identifier for vertex removal
  integer, parameter, public :: HADAPT_OPR_REMOVEVERTEX     = 4

  ! Operation identifier for refinment: 1-tria : 2-tria
  integer, parameter, public :: HADAPT_OPR_REF_TRIA2TRIA    = 5

  ! Operation identifier for refinment: 1-tria : 3-tria
  integer, parameter, public :: HADAPT_OPR_REF_TRIA3TRIA12  = 6

  ! Operation identifier for refinment: 1-tria : 3-tria
  integer, parameter, public :: HADAPT_OPR_REF_TRIA3TRIA23  = 7

  ! Operation identifier for refinment: 1-tria : 4-tria
  integer, parameter, public :: HADAPT_OPR_REF_TRIA4TRIA    = 8

  ! Operation identifier for refinment: 1-quad : 2-quad
  integer, parameter, public :: HADAPT_OPR_REF_QUAD2QUAD    = 9

  ! Operation identifier for refinment: 1-quad : 3-tria
  integer, parameter, public :: HADAPT_OPR_REF_QUAD3TRIA    = 10

  ! Operation identifier for refinment: 1-quad : 4-tria
  integer, parameter, public :: HADAPT_OPR_REF_QUAD4TRIA    = 11

  ! Operation identifier for refinment: 1-quad : 4-quad
  integer, parameter, public :: HADAPT_OPR_REF_QUAD4QUAD    = 12

  ! Operation identifier for conversion: 2-tria : 4-tria
  integer, parameter, public :: HADAPT_OPR_CVT_TRIA2TRIA    = 13

  ! Operation identifier for conversion: 2-quad : 4-tria
  integer, parameter, public :: HADAPT_OPR_CVT_QUAD2QUAD    = 14

  ! Operation identifier for conversion: 3-tria : 4-quad
  integer, parameter, public :: HADAPT_OPR_CVT_QUAD3TRIA    = 15

  ! Operation identifier for conversion: 4-tria : 4-quad
  integer, parameter, public :: HADAPT_OPR_CVT_QUAD4TRIA    = 16

  ! Operation identifier for coarsening: 2-tria : 1-tria
  integer, parameter, public :: HADAPT_OPR_CRS_2TRIA1TRIA   = 17

  ! Operation identifier for coarsening: 4-tria : 1-tria
  integer, parameter, public :: HADAPT_OPR_CRS_4TRIA1TRIA   = 18

  ! Operation identifier for coarsening: 4-tria : 2-tria
  integer, parameter, public :: HADAPT_OPR_CRS_4TRIA2TRIA1  = 19

  ! Operation identifier for coarsening: 4-tria : 2-tria
  integer, parameter, public :: HADAPT_OPR_CRS_4TRIA2TRIA2  = 20

  ! Operation identifier for coarsening: 4-tria : 2-tria
  integer, parameter, public :: HADAPT_OPR_CRS_4TRIA2TRIA3  = 21

  ! Operation identifier for coarsening: 4-quad : 1-quad
  integer, parameter, public :: HADAPT_OPR_CRS_4QUAD1QUAD   = 22

  ! Operation identifier for coarsening: 4-quad : 2-quad
  integer, parameter, public :: HADAPT_OPR_CRS_4QUAD2QUAD   = 23

  ! Operation identifier for coarsening: 4-quad : 3-tria
  integer, parameter, public :: HADAPT_OPR_CRS_4QUAD3TRIA   = 24

  ! Operation identifier for coarsening: 4-quad : 4-tria
  integer, parameter, public :: HADAPT_OPR_CRS_4QUAD4TRIA   = 25

  ! Operation identifier for coarsening: 2-quad : 1-quad
  integer, parameter, public :: HADAPT_OPR_CRS_2QUAD1QUAD   = 26

  ! Operation identifier for coarsening: 2-quad : 3-tria
  integer, parameter, public :: HADAPT_OPR_CRS_2QUAD3TRIA   = 27

  ! Operation identifier for coarsening: 3-tria : 1-quad
  integer, parameter, public :: HADAPT_OPR_CRS_3TRIA1QUAD   = 28

  ! Operation identifier for coarsening: 4-tria : 1-quad
  integer, parameter, public :: HADAPT_OPR_CRS_4TRIA1QUAD   = 29

  ! Operation identifier for coarsening: 4-tria : 3-tria
  integer, parameter, public :: HADAPT_OPR_CRS_4TRIA3TRIA2  = 30

  ! Operation identifier for coarsening: 4-tria : 3-tria
  integer, parameter, public :: HADAPT_OPR_CRS_4TRIA3TRIA3  = 31

  ! Operation identifier for refinement: 1-line : 2-line
  integer, parameter, public :: HADAPT_OPR_REF_LINE2LINE    = 32

  ! Operation identifier for coarsening: 2-line : 1-line
  integer, parameter, public :: HADAPT_OPR_CRS_2LINE1LINE   = 33

!</constantblock>

!<constantblock
! description="Duplication flags. Specifies which information is shared between adaptivity structures">

  integer(I32), parameter, public :: HADAPT_SHARE_IMARKER            = 2_I32**0
  integer(I32), parameter, public :: HADAPT_SHARE_IVERTEXAGE         = 2_I32**1
  integer(I32), parameter, public :: HADAPT_SHARE_INODALPROPERTY     = 2_I32**2
  integer(I32), parameter, public :: HADAPT_SHARE_IVERTICESATELEMENT = 2_I32**3
  integer(I32), parameter, public :: HADAPT_SHARE_INEIGHATELEMENT    = 2_I32**4
  integer(I32), parameter, public :: HADAPT_SHARE_IMIDNEIGHATELEMENT = 2_I32**5
  integer(I32), parameter, public :: HADAPT_SHARE_RVERTEXCOORDINATES = 2_I32**6
  integer(I32), parameter, public :: HADAPT_SHARE_RBOUNDARY          = 2_I32**7
  integer(I32), parameter, public :: HADAPT_SHARE_RELEMENTSATVERTEX  = 2_I32**8

!</constantblock>

!</constants>

  ! ***************************************************************************
  ! ***************************************************************************
  ! ***************************************************************************

!<types>

  !<typeblock>

  ! This type contains all data structures to handle
  ! adaptive grid refinement and grid coarsening.
  type :: t_hadapt
    ! Format Tag: Specifies the state of adaptation
    integer(I32) :: iSpec = HADAPT_NONE

    ! Duplication flag. Bitfield that indicates which information is
    ! shared with another adaptivity structure.
    ! When a bit is set to 1, the corresponding array is
    ! maintained by another adaptivity structure and must
    ! not be deleted by hadapt_releaseAdaptation.
    ! When the bit is 0, the array is a real copy of another array
    ! and must be deleted in hadapt_releaseAdaptation.
    integer(I32) :: iduplicationFlag = 0_I32

    ! Tag: Specifies the strategy for grid refinement and coarsening
    integer :: iadaptationStrategy = HADAPT_NOADAPTATION

    ! Maximum number of subdivisions from the original mesh
    integer :: nsubdividemax = 0

    ! Total number of grid refinement steps
    integer :: nRefinementSteps = 0

    ! Total number of grid coarsening steps
    integer :: nCoarseningSteps = 0

    ! Total number of grid smoothing steps
    integer :: nSmoothingSteps = 0

    ! Tolerance for refinement
    real(DP) :: drefinementTolerance = 0

    ! Tolerance for coarsening
    real(DP) :: dcoarseningTolerance = 0

    ! Dimension of the triangulation
    integer :: ndim = 0

    ! Total number of vertices (initially)
    integer :: NVT0 = 0

    ! Total number of vertices
    integer :: NVT = 0

    ! Increment of vertices
    integer :: increaseNVT = 0

    ! Total number of boundary vertives (initially)
    integer :: NVBD0 = 0

    ! Total number of boundary vertives
    integer :: NVBD = 0

    ! Total number of boundary components (should not change)
    integer :: NBCT = 0

    ! Total number of elements (initially)
    integer :: NEL0 = 0

    ! Total number of elements
    integer :: NEL = 0

    ! Maximum number of elements (before reallocation)
    integer :: NELMAX = 0

    ! Total number of green elements (required internally)
    integer :: nGreenElements = 0

    ! Nuber of elements with a defined number of vertices per element.
    ! InelOfType(TRIA_NVELINE1D) = number of lines in the mesh (1D).
    ! InelOfType(TRIA_NVETRI2D)  = number of triangles in the mesh (2D).
    ! InelOfType(TRIA_NVEQUAD2D) = number of quadrilaterals in the mesh (2D).
    integer, dimension(TRIA_MAXNVE) :: InelOfType = 0

    ! Same as InelOfType but this array stores the number of elements
    ! which are initially present in the mesh
    integer, dimension(TRIA_MAXNVE) :: InelOfType0 = 0

    ! Element marker array.
    ! Handle to
    !       p_Imarker = array [1..NEL] of integer.
    ! For each element there is one identifier [0,..,31] in the marker.
    ! BIT0: Is 0 for triangle and 1 for quadrilateral
    ! BIT1: Is 0 if first edge is not subdivided, 1 otherwise.
    ! BIT2: Is 0 if second edge is not subdivided, 1 otherwise.
    ! BIT3: Is 0 if third edge is not subdivided, 1 otherwise.
    ! BIT4: Is 0 if fourth edge is not subdivided, 1 otherwise.
    integer :: h_Imarker = ST_NOHANDLE

    ! Vertex age array.
    ! Handle to
    !       p_IvertexAge = array [1..NVT] of integer
    ! Each node has an individual age which is used to determine
    ! the type of the element (green/red triangle/quadrilateral).
    ! The nodes of the initial grid are assigned age 0 by default.
    ! Each node which is introduced at the midpoint of one edge (I,J)
    ! is given the age MAX(AGE(I),AGE(J))+1. For vertices which are
    ! inserted at the center of a quadrilateral, the maximum age of
    ! all four surrounding vertices is adopted and increased by one.
    ! A negative age means, the the vertex is lockes, that is, it
    ! cannot be removed. By definition, vertices of the initial
    ! triangulation cannot be removed and are always locked.
    integer :: h_IvertexAge = ST_NOHANDLE

    ! Pointer to h_IvertexAge.
    ! This array is introduced to increase performance and must
    ! not be modified by the user
    integer, dimension(:), pointer :: p_IvertexAge => null()

    ! Nodal property array.
    ! Handle to
    !       p_InodalProperty = array [1..NVT+NMT] of integer
    ! p_InodalProperty(i) defines for each vertex i=(1..NVT)
    ! and each edge i=(NVT+1..NVT+NMT) its function inside of the
    ! geometry. Generally said, the range of the p_InodalProperty-array
    ! characterises the type of the node (=vertex/edge):
    ! = 0    : The vertex/edge is an inner vertex/edge
    ! > 0    : The vertex/edge is a boundary vertex/edge on the real
    !           boundary. KNPR(.) defines the number of the boundary
    !           component.
    integer :: h_InodalProperty = ST_NOHANDLE

    ! Pointer to h_InodalProperty.
    ! This array is introduced to increase performance and must
    ! not be modified by the user
    integer, dimension(:), pointer :: p_InodalProperty => null()

    ! Vertices adjacent to an element.
    ! Handle to
    !       p_IverticesAtElement = array [1..TRIA_MAXNVE2D,1..NEL] of integer.
    ! For each element the node numbers of the corner-vertices
    ! in mathematically positive sense.
    integer :: h_IverticesAtElement = ST_NOHANDLE

    ! Pointer to h_IverticesAtElement.
    ! This array is introduced to increase performance and must
    ! not be modified by the user
    integer, dimension(:,:), pointer :: p_IverticesAtElement => null ()

    ! Neighbour elements adjacent to an element.
    ! Handle to
    !       p_IneighboursAtElement = array [1..TRIA_MAXNME2D,1..NEL] of integer
    ! For each element, the numbers of adjacent elements in mathematically
    ! positive sense, metting the element in an edge.
    integer :: h_IneighboursAtElement = ST_NOHANDLE

    ! Pointer to h_IneighboursAtElement.
    ! This array is introduced to increase performance and must
    ! not be modified by the user
    integer, dimension(:,:), pointer :: p_IneighboursAtElement => null ()

    ! Midneighbour elements adjacent to an element.
    ! Handle to
    !       p_ImidneighboursAtElement = array [1..TRIA_MAXNME2D,1..NEL] of integer
    ! H-adaptivity is performed element-by-element. If one element is
    ! refined, then its neighbors need to be informed that there
    ! are possibly two elements adjacent along one edge. This is a
    ! nonconforming state. When time comes to process the other
    ! element (with hanging node), the element knows which elements
    ! are adjacent along the first and the second half of the edge.
    integer :: h_ImidneighboursAtElement = ST_NOHANDLE

    ! Pointer to h_ImidneighboursAtElement.
    ! This array is introduced to increase performance and must
    ! not be modified by the user
    integer, dimension(:,:), pointer :: p_ImidneighboursAtElement => null ()

    ! Handle to h_DVertexCoords in 1D
    integer :: h_DvertexCoords1D

    ! Pointer to handle h_DvertexCoords in 1D
    real(DP), dimension(:,:), pointer :: p_DvertexCoords1D => null()

    ! Quadtree storing the nodal coordinates in 2D
    type(t_quadtreeDP) :: rVertexCoordinates2D

    ! Octree storing the nodal coordinates in 2D
    type(t_octreeDP) :: rVertexCoordinates3D

    ! Array of maps storing the boundary data
    ! p_IboundaryCpIdx and p_IverticesAtBoundary
    type(t_mapInt_DP), dimension(:), pointer :: rBoundary => null()

    ! Arraylist for elements-meeting-at-vertex structure
    type(t_arraylistInt) :: rElementsAtVertex
  end type t_hadapt

  !</typeblock>

!</types>

contains

  ! ***************************************************************************
  ! ***************************************************************************
  ! ***************************************************************************

!<function>

  pure function hadapt_getNVE(rhadapt, iel) result(nve)

!<description>
    ! This function returns the number of vertices present in the given element
!</description>

!<input>
    ! Adaptivity structure
    type(t_hadapt), intent(in) :: rhadapt

    ! Number of the element
    integer, intent(in) :: iel
!</input>

!<result>
    ! Number of vertices per element
    integer :: nve
!</result>
!</function>

    do nve = ubound(rhadapt%p_IverticesAtElement,1), 1, -1
      if (rhadapt%p_IverticesAtElement(nve,iel) .ne. 0) return
    end do
  end function hadapt_getNVE

  ! ***************************************************************************

!<function>

  pure function hadapt_getNNVE(rhadapt) result(nnve)

!<description>
    ! This function returns the maximum number of vertices
!</description>

!<input>
    ! Adaptivity structure
    type(t_hadapt), intent(in) :: rhadapt
!</input>

!<result>
    ! Maximum number of vertices per element
    integer :: nnve
!</result>
!</function>

    do nnve = ubound(rhadapt%InelOfType,1), 1, -1
      if (rhadapt%InelOfType(nnve) .ne. 0) return
    end do
  end function hadapt_getNNVE

  ! ***************************************************************************

!<subroutine>

  subroutine hadapt_setVertexCoords1D(rhadapt, h_DvertexCoords, nvt)

!<description>
    ! This subroutine sets the vertex coordinates given by the handle
    ! h_DvertexCoords that points to the two-dimensional array
    ! p_DvertexCoords and stores them in the map.
!</description>

!<input>
    ! Handle to the vertex coordinates
    integer, intent(in) :: h_DvertexCoords

    ! Number of vertices
    integer, intent(in) :: nvt
!</input>

!<inputoutput>
    ! Adaptivity structure
    type(t_hadapt), intent(inout) :: rhadapt
!</inputoutput>
!</subroutine>

    ! Check if handle is not empty
    if (h_DvertexCoords .eq. ST_NOHANDLE) then
      call output_line('Invalid handle!',&
                       OU_CLASS_ERROR,OU_MODE_STD,'hadapt_setVertexCoords1D')
      call sys_halt()
    end if

    ! Check if structure is already generated, then remove old structure first
    if (iand(rhadapt%iSpec, HADAPT_HAS_COORDS) .eq.&
                            HADAPT_HAS_COORDS) then
      if (rhadapt%h_DvertexCoords1D .ne. ST_NOHANDLE)&
          call storage_free(rhadapt%h_DvertexCoords1D)
    end if

    ! Copy handle and set pointer
    call storage_copy(h_DvertexCoords, rhadapt%h_DvertexCoords1D)
    call storage_getbase_double2D(rhadapt%h_DvertexCoords1D,&
                                  rhadapt%p_DvertexCoords1D)

    ! Set specifier for structure
    rhadapt%iSpec = ior(rhadapt%iSpec, HADAPT_HAS_COORDS)

    ! Set dimensions
    rhadapt%ndim = NDIM1D
    rhadapt%NVT  = nvt
  end subroutine hadapt_setVertexCoords1D

  ! ***************************************************************************

!<subroutine>

  subroutine hadapt_getVertexCoords1D(rhadapt, h_DvertexCoords, nvt, ndim)

!<description>
    ! This subroutine gets the vertex coordinates from the structure
    ! and stores them in the two-dimensional array associated to the
    ! handle h_DvertexCoords. Note that the allocated memory will
    ! be reallocated if is does not provide the correct dimensions.
!</description>

!<input>
    ! Adaptivity structure
    type(t_hadapt), intent(in) :: rhadapt
!</input>

!<inputoutput>
    ! Handlt to the vertex coordinate vector
    integer, intent(inout) :: h_DvertexCoords
!</inputoutput>

!<output>
    ! OPTIONAL: number of vertices
    integer, intent(out), optional :: nvt

    ! OPTIONAL: number of spatial dimensions
    integer, intent(out), optional :: ndim
!</output>
!</subroutine>

    ! local variables
    real(DP), dimension(:,:), pointer :: p_DvertexCoords
    real(DP), dimension(:,:), pointer :: p_Ddata
    integer, dimension(2) :: Isize

    ! Check if coordinates exists
    if (iand(rhadapt%iSpec, HADAPT_HAS_COORDS) .ne.&
                            HADAPT_HAS_COORDS) then
      call output_line('Coords do not exist!',&
                       OU_CLASS_ERROR,OU_MODE_STD,'hadapt_getVertexCoords1D')
      call sys_halt()
    end if

    ! Check if coordinates are given in 1D
    if (rhadapt%ndim .ne. NDIM1D) then
      call output_line('Invalid spatial dimension!',&
                       OU_CLASS_ERROR,OU_MODE_STD,'hadapt_getVertexCoords1D')
      call sys_halt()
    end if

    ! Check if given handle is valid
    if (h_DvertexCoords .eq. ST_NOHANDLE) then
      Isize=(/1,rhadapt%NVT/)
      call storage_new('hadapt_getVertexCoords1D', 'DCORVG', Isize,&
                       ST_DOUBLE, h_DvertexCoords, ST_NEWBLOCK_NOINIT)
    else
      call storage_getsize(h_DvertexCoords, Isize)
      if (Isize(1) .ne. 1) then
        call output_line('First component of DvertexCoords must be 1!',&
                         OU_CLASS_ERROR,OU_MODE_STD,'hadapt_getVertexCoords1D')
        call sys_halt()
      end if

      if (Isize(2) .lt. rhadapt%NVT) then
        call storage_realloc('hadapt_getVertexCoords1D', rhadapt%NVT,&
                             h_DvertexCoords, ST_NEWBLOCK_NOINIT, .false.)
      end if
    end if

    ! Set pointers
    call storage_getbase_double2D(rhadapt%h_DvertexCoords1D,&
                                  p_Ddata, rhadapt%NVT)
    call storage_getbase_double2D(h_DvertexCoords, p_DvertexCoords)

    ! Copy data
    call lalg_copyVector(p_Ddata, p_DvertexCoords, 1, rhadapt%NVT)

    ! Set dimension
    if (present(ndim)) ndim = rhadapt%ndim
    if (present(nvt))  nvt  = rhadapt%NVT
  end subroutine hadapt_getVertexCoords1D

  ! ***************************************************************************

!<subroutine>

  subroutine hadapt_setVertexCoords2D(rhadapt, h_DvertexCoords, nvt)

!<description>
    ! This subroutine sets the vertex coordinates given by the handle
    ! h_DvertexCoords that points to the two-dimensional array
    ! p_DvertexCoords and stores them in the quadtree.
!</description>

!<input>
    ! Handle to the vertex coordinates
    integer, intent(in) :: h_DvertexCoords

    ! Number of vertices
    integer, intent(in) :: nvt
!</input>

!<inputoutput>
    ! Adaptivity structure
    type(t_hadapt), intent(inout) :: rhadapt
!</inputoutput>
!</subroutine>

    ! local variables
    real(DP), dimension(:,:), pointer :: p_DvertexCoords
    real(DP) :: xmin,xmax,ymin,ymax
    integer  :: nnode

    ! Check if handle is not empty
    if (h_DvertexCoords .eq. ST_NOHANDLE) then
      call output_line('Invalid handle!',&
                       OU_CLASS_ERROR,OU_MODE_STD,'hadapt_setVertexCoords2D')
      call sys_halt()
    end if

    ! Check if quadtree is already generated, then remove old quadtree/octree first
    if (iand(rhadapt%iSpec, HADAPT_HAS_COORDS) .eq.&
                            HADAPT_HAS_COORDS) then
      call qtree_release(rhadapt%rVertexCoordinates2D)
    end if

    ! Set pointer
    call storage_getbase_double2D(h_DvertexCoords, p_DvertexCoords, nvt)

    ! Get outer bounding-box of vertex coordinates
    xmin = minval(p_DvertexCoords(1,:))
    xmax = maxval(p_DvertexCoords(1,:))
    ymin = minval(p_DvertexCoords(2,:))
    ymax = maxval(p_DvertexCoords(2,:))

    ! Estimate number of initial quadrilaterals
    nnode = int(0.5_DP*nvt)

    ! Create quadtree for vertices
    call qtree_create(rhadapt%rVertexCoordinates2D, nvt,&
                      nnode, xmin, ymin, xmax, ymax)

    ! Copy vertex coordinates to quadtree
    call qtree_copy(p_DvertexCoords, rhadapt%rVertexCoordinates2D)

    ! Set specifier for quadtree
    rhadapt%iSpec = ior(rhadapt%iSpec, HADAPT_HAS_COORDS)

    ! Set dimensions
    rhadapt%ndim = NDIM2D
    rhadapt%NVT  = nvt
  end subroutine hadapt_setVertexCoords2D

  ! ***************************************************************************

!<subroutine>

  subroutine hadapt_getVertexCoords2D(rhadapt, h_DvertexCoords, nvt, ndim)

!<description>
    ! This subroutine gets the vertex coordinates from the quadtree
    ! and stores them in the two-dimensional array associated to the
    ! handle h_DvertexCoords. Note that the allocated memory will
    ! be reallocated if is does not provide the correct dimensions.
!</description>

!<input>
    ! Adaptivity structure
    type(t_hadapt), intent(in) :: rhadapt
!</input>

!<inputoutput>
    ! Handlt to the vertex coordinate vector
    integer, intent(inout) :: h_DvertexCoords
!</inputoutput>

!<output>
    ! OPTIONAL: number of vertices
    integer, intent(out), optional :: nvt

    ! OPTIONAL: number of spatial dimensions
    integer, intent(out), optional :: ndim
!</output>
!</subroutine>

    ! Check if coordinates exists
    if (iand(rhadapt%iSpec, HADAPT_HAS_COORDS) .ne.&
                            HADAPT_HAS_COORDS) then
      call output_line('Quadtree does not exist!',&
                       OU_CLASS_ERROR,OU_MODE_STD,'hadapt_getVertexCoords2D')
      call sys_halt()
    end if

    ! Check if coordinates are given in 2D
    if (rhadapt%ndim .ne. NDIM2D) then
      call output_line('Invalid spatial dimension!',&
                       OU_CLASS_ERROR,OU_MODE_STD,'hadapt_getVertexCoords2D')
      call sys_halt()
    end if

    ! Copy quadtree to handle h_DvertexCoords
    call qtree_copy(rhadapt%rVertexCoordinates2D, h_DvertexCoords)

    ! Set dimension
    if (present(ndim)) ndim = rhadapt%ndim
    if (present(nvt))  nvt  = rhadapt%NVT
  end subroutine hadapt_getVertexCoords2D

  ! ***************************************************************************

!<subroutine>

  subroutine hadapt_setVertexCoords3D(rhadapt, h_DvertexCoords, nvt)

!<description>
    ! This subroutine sets the vertex coordinates given by the handle
    ! h_DvertexCoords that points to the three-dimensional array
    ! p_DvertexCoords and stores them in the octree.
!</description>

!<input>
    ! Handle to the vertex coordinates
    integer, intent(in) :: h_DvertexCoords

    ! Number of vertices
    integer, intent(in) :: nvt
!</input>

!<inputoutput>
    ! Adaptivity structure
    type(t_hadapt), intent(inout)    :: rhadapt
!</inputoutput>
!</subroutine>

    ! local variables
    real(DP), dimension(:,:), pointer :: p_DvertexCoords
    real(DP) :: xmin,xmax,ymin,ymax,zmin,zmax
    integer  :: nnode

    ! Check if handle is not empty
    if (h_DvertexCoords .eq. ST_NOHANDLE) then
      call output_line('Invalid handle!',&
                       OU_CLASS_ERROR,OU_MODE_STD,'hadapt_setVertexCoords3D')
      call sys_halt()
    end if

    ! Check if quadtree is already generated, then remove old quadtree/octree first
    if (iand(rhadapt%iSpec, HADAPT_HAS_COORDS) .eq.&
                            HADAPT_HAS_COORDS) then
      call otree_release(rhadapt%rVertexCoordinates3D)
    end if

    ! Set pointer
    call storage_getbase_double2D(h_DvertexCoords, p_DvertexCoords, nvt)

    ! Get outer bounding-box of vertex coordinates
    xmin = minval(p_DvertexCoords(1,:))
    xmax = maxval(p_DvertexCoords(1,:))
    ymin = minval(p_DvertexCoords(2,:))
    ymax = maxval(p_DvertexCoords(2,:))
    zmin = minval(p_DvertexCoords(3,:))
    zmax = maxval(p_DvertexCoords(3,:))

    ! Estimate number of initial quadrilaterals
    nnode = int(0.5_DP*nvt)

    ! Create octree for vertices
    call otree_create(rhadapt%rVertexCoordinates3D, nvt,&
                            nnode, xmin, ymin, zmin, xmax, ymax, zmax)

    ! Copy vertex coordinates to octree
    call otree_copy(p_DvertexCoords, rhadapt%rVertexCoordinates3D)

    ! Set specifier for quadtree
    rhadapt%iSpec=ior(rhadapt%iSpec, HADAPT_HAS_COORDS)

    ! Set dimensions
    rhadapt%ndim = NDIM3D
    rhadapt%NVT  = nvt
  end subroutine hadapt_setVertexCoords3D

  ! ***************************************************************************

!<subroutine>

  subroutine hadapt_getVertexCoords3D(rhadapt, h_DvertexCoords, nvt, ndim)

!<description>
    ! This subroutine gets the vertex coordinates from the octree
    ! and stores them in the three-dimensional array associated to the
    ! handle h_DvertexCoords. Note that the allocated memory will
    ! be reallocated if is does not provide the correct dimensions.
!</description>

!<input>
    ! Adaptivity structure
    type(t_hadapt), intent(in) :: rhadapt
!</input>

!<inputoutput>
    ! Handlt to the vertex coordinate vector
    integer, intent(inout) :: h_DvertexCoords
!</inputoutput>

!<output>
    ! OPTIONAL: number of vertices
    integer, intent(out), optional :: nvt

    ! OPTIONAL: number of spatial dimensions
    integer, intent(out), optional :: ndim
!</output>
!</subroutine>

    ! Check if coordinates exists
    if (iand(rhadapt%iSpec, HADAPT_HAS_COORDS) .ne.&
                            HADAPT_HAS_COORDS) then
      call output_line('Octree does not exist!',&
                       OU_CLASS_ERROR,OU_MODE_STD,'hadapt_getVertexCoords3D')
      call sys_halt()
    end if

    ! Check if coordinates are given in 2D
    if (rhadapt%ndim .ne. NDIM3D) then
      call output_line('Invalid spatial dimension!',&
                       OU_CLASS_ERROR,OU_MODE_STD,'hadapt_getVertexCoords3D')
      call sys_halt()
    end if

    ! Copy octree to handle h_DvertexCoords
    call otree_copy(rhadapt%rVertexCoordinates3D, h_DvertexCoords)

    ! Set dimension
    if (present(ndim)) ndim = rhadapt%ndim
    if (present(nvt))  nvt  = rhadapt%NVT
  end subroutine hadapt_getVertexCoords3D

  ! ***************************************************************************

!<subroutine>

  subroutine hadapt_setVerticesAtElement(rhadapt, h_IverticesAtElement, nel)

!<description>
    ! This routine assigns the handle to the "vertices-at-element" array
    ! to the adaptivity structure and sets up internal links.
!</description>

!<input>
    ! Handle to p_IverticesAtElement
    integer, intent(in) :: h_IverticesAtElement

    ! Total number of elements
    integer, intent(in) :: nel
!</input>

!<inputoutput>
    ! adaptivity structure
    type(t_hadapt), intent(inout) :: rhadapt
!</inputoutput>
!</subroutine>

    ! Check if handle is not empty
    if (h_IverticesAtElement .eq. ST_NOHANDLE) then
      call output_line('Invalid handle!',&
                       OU_CLASS_ERROR,OU_MODE_STD,'hadapt_setVerticesAtElement')
      call sys_halt()
    end if

    rhadapt%h_IverticesAtElement = h_IverticesAtElement
    call storage_getbase_int2D(rhadapt%h_IverticesAtElement,&
                               rhadapt%p_IverticesAtElement)

    ! Set specifier for IverticesAtElement
    rhadapt%iSpec = ior(rhadapt%iSpec, HADAPT_HAS_VERTATELEM)

    ! Set dimensions
    rhadapt%NEL = nel
  end subroutine hadapt_setVerticesAtElement

  ! ***************************************************************************

!<subroutine>

  subroutine hadapt_getVerticesAtElement(rhadapt, h_IverticesAtElement, nel)

!<description>
    ! This routine assigns the "vertices-at-element" array from the
    ! adaptivity structure to the given handle.
!</description>

!<input>
    ! Adaptivity structure
    type(t_hadapt), intent(in) :: rhadapt
!</input>

!<inputoutput>
    ! Hande to p_IverticesAtElement
    integer, intent(inout) :: h_IverticesAtElement
!</inputoutput>

!<output>
    ! OPTIONAL: number of elements
    integer, intent(out), optional :: nel
!</output>
!</subroutine>

    ! Check if "vertices-at-element" array exists
    if (iand(rhadapt%iSpec, HADAPT_HAS_VERTATELEM) .ne.&
                            HADAPT_HAS_VERTATELEM) then
      call output_line('Structure does not exist!',&
                       OU_CLASS_ERROR,OU_MODE_STD,'hadapt_getVerticesAtElement')
      call sys_halt()
    end if

    ! Check if handle needs to be freed first
    if ((h_IverticesAtElement .ne. ST_NOHANDLE) .and.&
        (h_IverticesAtElement .ne. rhadapt%h_IverticesAtElement))&
        call storage_free(h_IverticesAtElement)

    ! Assign handle
    h_IverticesAtElement = rhadapt%h_IverticesAtElement

    ! Set dimensions
    if(present(nel))   nel = rhadapt%NEL
  end subroutine hadapt_getVerticesAtElement

  ! ***************************************************************************

!<subroutine>

  subroutine hadapt_setNeighboursAtElement(rhadapt, h_IneighboursAtElement)

!<description>
    ! This routine assigns the handle to the "neighbours-at-element" array
    ! to the adaptivity structure and sets up internal links.
!</description>

!<input>
    ! Handle to p_IneighboursAtElement
    integer, intent(in) :: h_IneighboursAtElement
!</input>

!<inputoutput>
    ! adaptivity structure
    type(t_hadapt), intent(inout) :: rhadapt
!</inputoutput>
!</subroutine>

    ! Check if handle is not empty
    if (h_IneighboursAtElement .eq. ST_NOHANDLE) then
      call output_line('Invalid handle!',&
                       OU_CLASS_ERROR,OU_MODE_STD,'hadapt_setNeighboursAtElement')
      call sys_halt()
    end if

    rhadapt%h_IneighboursAtElement = h_IneighboursAtElement
    call storage_getbase_int2D(rhadapt%h_IneighboursAtElement,&
                               rhadapt%p_IneighboursAtElement)

    ! What spatial dimension are we?
    select case(rhadapt%ndim)
    case (NDIM1D,&
          NDIM3D)
      ! Do nothing

    case(NDIM2D)
      ! Create structure of "mid-adjacent" neighbours
      call storage_copy(rhadapt%h_IneighboursAtElement,&
                        rhadapt%h_ImidneighboursAtElement)
      call storage_getbase_int2D(rhadapt%h_ImidneighboursAtElement,&
                                 rhadapt%p_ImidneighboursAtElement)

      ! Set specifier for ImidneighboursAtElement
      rhadapt%iSpec = ior(rhadapt%iSpec, HADAPT_HAS_MIDNEIGH)

    case default
      call output_line('Unsupported spatial dimension!',&
                       OU_CLASS_ERROR,OU_MODE_STD,'hadapt_setNeighboursAtElement')
      call sys_halt()
    end select

    ! Set specifier for IverticesAtElement
    rhadapt%iSpec = ior(rhadapt%iSpec, HADAPT_HAS_NEIGHATELEM)
  end subroutine hadapt_setNeighboursAtElement

  ! ***************************************************************************

!<subroutine>

  subroutine hadapt_getNeighboursAtElement(rhadapt,&
                                           h_IneighboursAtElement, nel)

!<description>
    ! This routine assigns the "neighbours-at-element" array from the
    ! adaptivity structure to the given handle
!</description>

!<input>
    ! Adaptivity structure
    type(t_hadapt), intent(in) :: rhadapt
!</input>

!<inputoutput>
    ! Hande to p_IneighboursAtElement
    integer, intent(inout) :: h_IneighboursAtElement
!</inputoutput>

!<output>
    ! OPTIONAL: number of elements
    integer, intent(out), optional :: nel
!</output>
!</subroutine>

    ! Check if "neighbours-at-element" array exists
    if (iand(rhadapt%iSpec, HADAPT_HAS_NEIGHATELEM) .ne.&
                            HADAPT_HAS_NEIGHATELEM) then
      call output_line('Structure does not exist!',&
                       OU_CLASS_ERROR,OU_MODE_STD,'hadapt_getNeighboursAtElement')
      call sys_halt()
    end if

    ! Check if handle needs to be freed first
    if ((h_IneighboursAtElement .ne. ST_NOHANDLE) .and.&
        (h_IneighboursAtElement .ne. rhadapt%h_IneighboursAtElement))&
        call storage_free(h_IneighboursAtElement)

    ! Assign handle
    h_IneighboursAtElement = rhadapt%h_IneighboursAtElement

    ! Set dimensions
    if(present(nel))   nel = rhadapt%NEL
  end subroutine hadapt_getNeighboursAtElement

  ! ***************************************************************************

!<subroutine>

  subroutine hadapt_setNelOfType(rhadapt, InelOfType)

!<description>
    ! This subroutine sets the number of elements with a defined number
    ! of vertices per elements, that is, InelOfType(TRIA_MAXNVE)
!</description>

!<input>
    ! Number of elements with a defined number of vertices per element.
    integer, dimension(TRIA_MAXNVE), intent(in) :: InelOfType
!</input>

!<inputoutput>
    ! Adaptivity structure
    type(t_hadapt), intent(inout) :: rhadapt
!</inputoutput>
!</subroutine>

    rhadapt%InelOfType = InelOfType

    ! Set specifier for InelOfType
    rhadapt%iSpec = ior(rhadapt%iSpec, HADAPT_HAS_NELOFTYPE)
  end subroutine hadapt_setNelOfType

  ! ***************************************************************************

!<subroutine>

  subroutine hadapt_getNelOfType(rhadapt, InelOfType)

!<description>
    ! This subroutine returns the number of elements with a defined number
    ! of vertices per elements, that is, InelOfType(TRIA_MAXNVE)
!</description>

!<input>
    ! Adaptivity structure
    type(t_hadapt), intent(in) :: rhadapt
!</input>

!<output>
    ! Number of elements with a defined number of vertices per element.
    integer, dimension(TRIA_MAXNVE), intent(out) :: InelOfType
!</output>
!</subroutine>

    ! Check if "InelOfType" array exists
    if (iand(rhadapt%iSpec, HADAPT_HAS_NELOFTYPE) .ne.&
                            HADAPT_HAS_NELOFTYPE) then
      call output_line('Structure does not exist!',&
                       OU_CLASS_ERROR,OU_MODE_STD,'hadapt_getNelOfType')
      call sys_halt()
    end if

    InelOfType = rhadapt%InelOfType
  end subroutine hadapt_getNelOfType

  ! ***************************************************************************

!<subroutine>

  subroutine hadapt_setBoundary(rhadapt, h_IboundaryCpIdx, h_IverticesAtBoundary,&
                                h_DvertexParameterValue, nbct, nvbd)

!<description>
    ! This subroutine sets the boundary structure to the
    ! adaptivity structure and initialises internal data
!</description>

!<input>
    ! Handle to p_IboundaryCpIdx
    integer, intent(in) :: h_IboundaryCpIdx

    ! Handle to p_IverticesAtBoundary
    integer, intent(in) :: h_IverticesAtBoundary

    ! Handle to p_DvertexParameterValue
    integer, intent(in) :: h_DvertexParameterValue

    ! Number of boundary components
    integer, intent(in) :: nbct

    ! Number of vertices at the boundary
    integer, intent(in) :: nvbd
!</input>

!<inputoutput>
    ! Adaptivity structure
    type(t_hadapt), intent(inout) :: rhadapt
!</inputoutput>
!</subroutine>

    ! local variables
    type(it_mapInt_DP) :: rmapIter
    real(DP), dimension(:), pointer   :: p_DvertexParameterValue
    integer, dimension(:), pointer    :: p_IboundaryCpIdx
    integer, dimension(:), pointer    :: p_IverticesAtBoundary
    real(DP), dimension(3) :: Ddata
    integer :: ibct,ioff,ivbd,ivbdEnd,ivbdStart,lvbd

    ! Check if handle are not empty
    if ((h_IboundaryCpIdx .eq. ST_NOHANDLE) .or.&
        (h_IverticesAtBoundary .eq. ST_NOHANDLE) .or.&
        (h_DvertexParameterValue .eq. ST_NOHANDLE)) then
      call output_line('Invalid handle!',&
                       OU_CLASS_ERROR,OU_MODE_STD,'hadapt_setBoundary')
      call sys_halt()
    end if

    ! Check if boundary structure exists and remove it
    if (associated(rhadapt%rBoundary)) then
      do ibct = 1, size(rhadapt%rBoundary,1)
        call map_release(rhadapt%rBoundary(ibct))
      end do
      deallocate(rhadapt%rBoundary)
    end if

    ! Set pointers
    call storage_getbase_double(h_DvertexParameterValue,&
                                p_DvertexParameterValue, nvbd)
    call storage_getbase_int(h_IboundaryCpIdx,&
                             p_IboundaryCpIdx, nbct+1)
    call storage_getbase_int(h_IverticesAtBoundary,&
                             p_IverticesAtBoundary, nvbd)

    ! Allocate array of maps
    allocate(rhadapt%rBoundary(nbct))

    ! Initialization
    ioff = 0

    do ibct = 1, nbct

      ! Create a separate map for each boundary component
      lvbd = p_IboundaryCpIdx(ibct+1)-p_IboundaryCpIdx(ibct)
      call map_create(rhadapt%rBoundary(ibct), lvbd, 3)

      ! Set subdimensions
      ivbdStart = ioff+1
      ivbdEnd   = ioff+lvbd
      ioff      = ioff+lvbd

      ! Insert first element into map
      Ddata(1) =      p_DvertexParameterValue(ivbdStart)
      Ddata(2) = real(p_IverticesAtBoundary(ivbdEnd),DP)
      Ddata(3) = real(p_IverticesAtBoundary(ivbdStart+1),DP)

      rmapIter = map_insert(rhadapt%rBoundary(ibct),&
                            p_IverticesAtBoundary(ivbdStart), Ddata)

      ! Insert second to last but one element into map
      do ivbd = ivbdStart+1, ivbdEnd-1
        Ddata(1) =      p_DvertexParameterValue(ivbd)
        Ddata(2) = real(p_IverticesAtBoundary(ivbd-1),DP)
        Ddata(3) = real(p_IverticesAtBoundary(ivbd+1),DP)

        rmapIter = map_insert(rhadapt%rBoundary(ibct),&
                              p_IverticesAtBoundary(ivbd), Ddata)
      end do

      ! Insert last element into map
      Ddata(1) =      p_DvertexParameterValue(ivbdEnd)
      Ddata(2) = real(p_IverticesAtBoundary(ivbdEnd-1),DP)
      Ddata(3) = real(p_IverticesAtBoundary(ivbdStart),DP)

      rmapIter = map_insert(rhadapt%rBoundary(ibct),&
                            p_IverticesAtBoundary(ivbdEnd), Ddata)
    end do

    ! Set dimensions
    rhadapt%NBCT  = nbct
    rhadapt%NVBD  = nvbd
    rhadapt%NVBD0 = nvbd

    ! Set specifier for boundary
    rhadapt%iSpec = ior(rhadapt%iSpec, HADAPT_HAS_BOUNDARY)

  end subroutine hadapt_setBoundary

  ! ***************************************************************************

!<function>

  function hadapt_getSubdivisonLevel(rhadapt) result(isubdivisionLevel)

!<description>
    ! This function calculates the deepest subdivision level present
    ! in the adaptivity structure. If the structure is not
    ! initialised, then this routine returns -1.
!</description>

!<input>
    ! Adaptivity structure
    type(t_hadapt), intent(in) :: rhadapt
!</input>

!<result>
    ! Deepest subdivision level
    integer :: IsubdivisionLevel
!</result>
!</function>

    ! local variable
    integer, dimension(:), pointer :: p_IvertexAge

    if (rhadapt%h_IvertexAge .eq. ST_NOHANDLE) then
      isubdivisionLevel = -1
    else
      call storage_getbase_int(rhadapt%h_IvertexAge, p_IvertexAge)
      isubdivisionLevel = maxval(p_IvertexAge)
    end if

  end function hadapt_getSubdivisonLevel

  ! ***************************************************************************

!<subroutine>

  subroutine hadapt_getBoundary(rhadapt, h_IboundaryCpIdx, h_IverticesAtBoundary,&
                                h_DvertexParameterValue, nvbd, nbct)

!<description>
    ! This subroutine extracts the boundary data from the adaptivity structure
    ! and generates the the arrays p_IboundaryCpIdx, p_IverticesAtBoundary and
    ! p_DvertexParameterValue. If the optional parameter nvbd is given
    ! then the number of boundary vertices is assigned.
!</description>

!<input>
    ! Adaptivity structure
    type(t_hadapt), intent(in) :: rhadapt
!</input>

!<inputoutput>
    ! Handle to p_IboundaryCpIdx
    integer, intent(inout) :: h_IboundaryCpIdx

    ! Handle to p_IverticesAtBoundary
    integer, intent(inout) :: h_IverticesAtBoundary

    ! Handle to p_DvertexParameterValue
    integer, intent(inout) :: h_DvertexParameterValue

    ! OPTIONAL: number of vertices at the boundary
    integer, intent(out), optional :: nvbd

    ! OPTIONAL: number of boundary components
    integer, intent(out), optional :: nbct
!</inputoutput>
!</subroutine>

    ! local variables
    type(it_mapInt_DP) :: rmapIter
    real(DP), dimension(:), pointer :: p_DvertexParameterValue, p_Ddata
    integer, dimension(:), pointer  :: p_IboundaryCpIdx
    integer, dimension(:), pointer  :: p_IverticesAtBoundary
    integer :: ibct,ioff,isize,ivbd,ivbdEnd,ivbdStart,lvbd

    ! Check if boundary data exists
    if (iand(rhadapt%iSpec, HADAPT_HAS_BOUNDARY) .ne.&
                            HADAPT_HAS_BOUNDARY) then
      call output_line('Structure does not exist!',&
                       OU_CLASS_ERROR,OU_MODE_STD,'hadapt_getBoundary')
      call sys_halt()
    end if

    ! Due to the fact, that the boundary data may be stored in
    ! multiple maps (one for each component) the static arrays are
    ! pre-allocated and filled step-by-step from the dynamic data in
    ! the boundary map(s).

    ! Check if the arrays p_IboundaryCpIdx, p_IverticesAtBoundary and
    ! p_DvertexParameterValue have correct dimension
    if (h_IboundaryCpIdx .eq. ST_NOHANDLE) then
      call storage_new('hadapt_getBoundary', 'p_IboundaryCpIdx',&
                       rhadapt%NBCT+1, ST_INT, h_IboundaryCpIdx,&
                       ST_NEWBLOCK_NOINIT)
    else
      call storage_getsize(h_IboundaryCpIdx, isize)
      if (isize < rhadapt%NBCT+1) then
        call storage_realloc('hadapt_getBoundary', rhadapt%NBCT+1,&
                             h_IboundaryCpIdx, ST_NEWBLOCK_NOINIT, .false.)
      end if
    end if

    if (h_IverticesAtBoundary .eq. ST_NOHANDLE) then
      call storage_new('hadapt_getBoundary', 'p_IverticesAtBoundary',&
                       rhadapt%NVBD, ST_INT, h_IverticesAtBoundary,&
                       ST_NEWBLOCK_NOINIT)
    else
      call storage_getsize(h_IverticesAtBoundary, isize)
      if (isize < rhadapt%NVBD) then
        call storage_realloc('hadapt_getBoundary', rhadapt%NVBD,&
                             h_IverticesAtBoundary, ST_NEWBLOCK_NOINIT, .false.)
      end if
    end if

    if (h_DvertexParameterValue .eq. ST_NOHANDLE) then
      call storage_new('hadapt_getBoundary', 'p_DvertexParameterValue',&
                       rhadapt%NVBD, ST_DOUBLE, h_DvertexParameterValue,&
                       ST_NEWBLOCK_NOINIT)
    else
      call storage_getsize(h_DvertexParameterValue, isize)
      if (isize < rhadapt%NVBD) then
        call storage_realloc('hadapt_getBoundary', rhadapt%NVBD,&
                             h_DvertexParameterValue, ST_NEWBLOCK_NOINIT, .false.)
      end if
    end if

    ! Set pointers
    call storage_getbase_int(h_IboundaryCpIdx,&
                             p_IboundaryCpIdx, rhadapt%NBCT+1)
    call storage_getbase_int(h_IverticesAtBoundary,&
                             p_IverticesAtBoundary, rhadapt%NVBD)
    call storage_getbase_double(h_DvertexParameterValue,&
                                p_DvertexParameterValue, rhadapt%NVBD)

    ! Initialization
    p_IboundaryCpIdx(1) = 1
    ioff = 0
    ivbd = 0

    do ibct = 1, rhadapt%NBCT

      ! Set subdimensions
      lvbd      = map_size(rhadapt%rBoundary(ibct))
      ivbdStart = ioff+1
      ivbdEnd   = ioff+lvbd
      ioff      = ioff+lvbd

      ! Set index for next boundary component
      p_IboundaryCpIdx(ibct+1) = ioff+1

      ! Get static boundary data from map
      rmapIter = map_begin(rhadapt%rBoundary(ibct))

      do while(.not.map_isNull(rmapIter))

        ! Increase counter
        ivbd = ivbd+1

        ! Get key value
        p_IverticesAtBoundary(ivbd) = map_get(rhadapt%rBoundary(ibct), rmapIter)

        ! Get auxiliary data
        call map_getbase_data(rhadapt%rBoundary(ibct), rmapIter, p_Ddata)
        p_DvertexParameterValue(ivbd) = p_Ddata(1)

        ! Increase iterators
        call map_next(rmapIter)
      end do

      ! Consistency check
      if (ivbd .ne. ivbdEnd) then
        call output_line('An internal error occured!',&
            OU_CLASS_ERROR,OU_MODE_STD,'hadapt_getBoundary')
        call sys_halt()
      end if

      ! Sort array w.r.t. increasing parameter values
      call sort_dp(p_DvertexParameterValue(ivbdStart:ivbdEnd),&
                   SORT_QUICK, p_IverticesAtBoundary(ivbdStart:ivbdEnd))
    end do

    ! Set dimension
    if(present(nvbd)) nvbd = rhadapt%NVBD
    if(present(nbct)) nbct = rhadapt%NBCT
  end subroutine hadapt_getBoundary

  ! ***************************************************************************

!<subroutine>

  subroutine hadapt_setNodalProperty(rhadapt, h_InodalProperty)

!<description>
    ! This subroutine sets the nodal property list to the adaptivity structure
!</description>

!<input>
    ! Handle to p_InodalProperty
    integer, intent(in) :: h_InodalProperty
!</input>

!<inputoutput>
    ! Adaptivity structure
    type(t_hadapt), intent(inout) :: rhadapt
!</inputoutput>
!</subroutine>

    ! Check if handle is not empty
    if (h_InodalProperty .eq. ST_NOHANDLE) then
      call output_line('Invalid handle!',&
                       OU_CLASS_ERROR,OU_MODE_STD,'hadapt_setNodalProperty')
      call sys_halt()
    end if

    rhadapt%h_InodalProperty = h_InodalProperty
    call storage_getbase_int(rhadapt%h_InodalProperty,&
                             rhadapt%p_InodalProperty)

    ! Set specifier for InodalProperty
    rhadapt%iSpec = ior(rhadapt%iSpec, HADAPT_HAS_NODALPROP)
  end subroutine hadapt_setNodalProperty

  ! ***************************************************************************

!<subroutine>

  subroutine hadapt_getNodalProperty(rhadapt, h_InodalProperty)

!<description>
    ! This subroutine assignes the "nodal property" array from the
    ! adaptivity structure to the given handle
!</description>

!<input>
    ! Adaptivity structure
    type(t_hadapt), intent(in) :: rhadapt
!</input>

!<inputoutput>
    ! Hande to p_InodalProperty
    integer, intent(inout) :: h_InodalProperty
!</inputoutput>
!</subroutine>

    ! Check if "nodal property" list exits
    if (iand(rhadapt%iSpec, HADAPT_HAS_NODALPROP) .ne.&
                            HADAPT_HAS_NODALPROP) then
      call output_line('Structure does not exist!',&
                       OU_CLASS_ERROR,OU_MODE_STD,'hadapt_getNodalProperty')
      call sys_halt()
    end if

    ! Check if handle needs to be freed first
    if ((h_InodalProperty .ne. ST_NOHANDLE) .and.&
        (h_InodalProperty .ne. rhadapt%h_InodalProperty))&
        call storage_free(h_InodalProperty)

    ! Assign handle
    h_InodalProperty = rhadapt%h_InodalProperty
  end subroutine hadapt_getNodalProperty

  ! ***************************************************************************

!<subroutine>

  subroutine hadapt_genElementsAtVertex(rhadapt)

!<description>
    ! This subroutine generates a list of arrays for the "elements-at-vertex"
    ! structure. This is an internal structure of the adaptivity structure
    ! which is used, e.g., in the removel of vertices.
!</description>

!<inputoutput>
    ! Adaptivity structur
    type(t_hadapt), intent(inout) :: rhadapt
!</inputoutput>
!</subroutine>

    ! local variables
    integer :: iel,ive

    ! Check if "vertices-at-element" list exists
    if (iand(rhadapt%iSpec, HADAPT_HAS_VERTATELEM) .ne.&
                            HADAPT_HAS_VERTATELEM) then
      call output_line('Structure does not exist!',&
                       OU_CLASS_ERROR,OU_MODE_STD,'hadapt_genElementsAtVertex')
      call sys_halt()
    end if

    ! Create arraylist
    call alst_create(rhadapt%relementsAtVertex, 2*rhadapt%NVT, 8*rhadapt%NEL)

    ! Fill arraylist
    do iel = 1, rhadapt%NEL
      do ive = 1, hadapt_getNVE(rhadapt,iel)
        call alst_push_back(rhadapt%relementsAtVertex,&
                            rhadapt%p_IverticesAtElement(ive,iel), iel)
      end do
    end do

    ! Set specifier for relementsAtVertex
    rhadapt%iSpec = ior(rhadapt%iSpec, HADAPT_HAS_ELEMATVERTEX)
  end subroutine hadapt_genElementsAtVertex

end module hadaptaux

max_stars_repo_licenses: Intel Unlicense

max_issues_repo_licenses: Intel Unlicense

max_forks_repo_licenses: Intel Unlicense


text:
program t
  implicit none
! inquire stmt 'exist' specifier
  logical::a=.true.
  inquire (file='tmpfile', exist=a)
  print *,a
  open (95, status='new', file='tmpfile')
  inquire (file='tmpfile', exist=a)
  print *,a
  close (95,status='delete')
   endprogram t

max_stars_repo_name: maddenp/ppp

max_issues_repo_name: maddenp/ppp

max_forks_repo_name: maddenp/ppp


! { dg-do compile }
! { dg-options "-O2 -fdump-tree-original" }

subroutine foo(x)
  integer :: x(4)
  x(:) = (/ 3, 1, 4, 1 /)
end subroutine

subroutine bar(x)
  integer :: x(4)
  x = (/ 3, 1, 4, 1 /)
end subroutine

! { dg-final { scan-tree-dump-times "memcpy|ref-all\[^\\n\]*ref-all" 2 "original" } }


text:
! { dg-do compile }
! { dg-options "-O2" }
!
! Tests the fix for PR80965 in which the use of the name 'loc'
! for the dummy argument of 'xyz' caused an ICE. If the module
! was used, the error "DUMMY attribute conflicts with INTRINSIC
! attribute in ‘loc’ at (1)" was emitted. Note that although 'loc'
! is a GNU extension and so can be over-ridden, this is not very
! good practice.
!
! Contributed by David Sagan  <david.sagan@gmail.com>
!
module mode3_mod
contains
  subroutine xyz (loc)
    implicit none
    class(*) :: loc
    real x(6)
    integer ix_use
    select type (loc)
      type is (integer)
        x = 0
        print *, "integer"
      type is (real)
        ix_use = 0
        print *, "real"
    end select
  end subroutine xyz
end module mode3_mod

max_stars_repo_name: TUDSSL/TICS

max_issues_repo_name: TUDSSL/TICS

max_forks_repo_name: TUDSSL/TICS


text:
! NCLFORTSTART
SUBROUTINE DCOMPUTEWSPD(wspd, u, v, nx, ny)

    IMPLICIT NONE

    !f2py threadsafe
    !f2py intent(in,out) :: wspd

    INTEGER, INTENT(IN) :: nx, ny
    REAL(KIND=8), DIMENSION(nx,ny), INTENT(OUT) :: wspd
    REAL(KIND=8), DIMENSION(nx,ny), INTENT(IN) :: u, v
! NCLEND

    INTEGER i, j

    !$OMP PARALLEL DO COLLAPSE(2) SCHEDULE(runtime)
    DO j = 1,ny
        DO i = 1,nx
            wspd(i,j) = SQRT(u(i,j)*u(i,j) + v(i,j)*v(i,j))
        END DO
    END DO
    !$OMP END PARALLEL DO

END SUBROUTINE DCOMPUTEWSPD

! NCLFORTSTART
SUBROUTINE DCOMPUTEWDIR(wdir, u, v, nx, ny)
    USE wrf_constants, ONLY : DEG_PER_RAD

    IMPLICIT NONE

    !f2py threadsafe
    !f2py intent(in,out) :: wdir

    INTEGER, INTENT(IN) :: nx, ny
    REAL(KIND=8), DIMENSION(nx,ny), INTENT(OUT) :: wdir
    REAL(KIND=8), DIMENSION(nx,ny), INTENT(IN) :: u, v
! NCLEND

    INTEGER i, j

    !$OMP PARALLEL DO COLLAPSE(2) SCHEDULE(runtime)
    DO j = 1,ny
        DO i = 1,nx
            wdir(i,j) = MOD(270.0 - ATAN2(v(i,j), u(i,j)) * DEG_PER_RAD, 360.)
        END DO
    END DO
    !$OMP END PARALLEL DO

END SUBROUTINE DCOMPUTEWDIR

max_stars_repo_name: bbonenfant/wrf-python

max_issues_repo_name: mostamndi/wrf-python

max_forks_repo_name: mostamndi/wrf-python


text:
	SUBROUTINE GG_GAIRM ( dattim, icolor, ssize, iwidth, lwidth,
     +			      itbclr2, lvlfl, lvfil, iflags, iret )
C************************************************************************
C* GG_AIRM								*
C*									*
C* This subroutine plots the current g-airmets for the lower 48 states  *
C* and the adjacent coastal waters.  G-airmets are rendered at 3-hour   *
C* intervals to depict snapshots                                        *
C*									*
C* GG_GAIRM ( DATTIM, ICOLOR, SSIZE, IWIDTH, LWIDTH, ITBVCLR2, LVLFL,	*
C* 	     LVFIL, IFLAGS, IRET )					*
C*									*
C* Input parameters:							*
C*	DATTIM		CHAR*		Ending time for airmet          *
C*	ICOLOR (5)	INTEGER		Line and symbol colors		*
C*	SSIZE  (5)	REAL		Symbol sizes			*
C*	IWIDTH (5)	INTEGER		Symbol widths			*
C*	LWIDTH (5)	INTEGER		Line widths			*
C*	ITBCLR2		INTEGER		Color2 for turbulence		*
C*	LVLFL		INTEGER		Flight level			*
C*	LVFIL  (2)	INTEGER		Lower/upper filter level	*
C*	IFLAGS (4)	INTEGER		Flags for labels		*
C*					  0 = false			*
C*					  1 = true			*
C*									*
C* Output parameters:							*
C*	IRET		INTEGER		Return code			*
C*									*
C**									*
C* Log:                                                                 *
C* L. Hinson/AWC         3/10   Derived from GG_AIRM                    *
C* L. Hinson/AWC         11/10  Increase NW Parameter to 6000           *
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
C*
	PARAMETER	( NW = 6000 )
C*
	PARAMETER	( NR = 6 )
C
	PARAMETER	( JIFR = 1, JMTO = 2, JTRB = 3, JICE = 4,
     +                    JSSW = 5, JLWS = 6 )
	PARAMETER	( JSYM = 1, JTIM = 2, JFLT = 3 )
        PARAMETER       ( IOFFM2HR = 4, IOFFM1HR = 5 )
        PARAMETER       ( IOFF0HR = 6, IOFF1HR = 7, IOFF2HR = 8 )
C*
	CHARACTER*(*)	dattim
	INTEGER		icolor(*), iwidth(*), lwidth(*), iflags(*),
     +			lvfil(*)
	REAL		ssize(*)
C*
	CHARACTER	path*25, templ*(MXTMPL), cdttm*20, dattm2*20,
     +			buffer*100, tfile*128, label*4, dattm4*20,
     +			tmstp4*20, tmstr4*20, carr (11)*20, tid*4,
     +			ttype*2
	CHARACTER*(MXFLSZ)      filnam, files (MXNMFL)
C*
	CHARACTER	type (NW)*2, timstr (NW)*20, timstp (NW)*20,
     +			waid (NW)*4, flevel (2,NW)*4, strir*3, strmo*4,
     +                  strws*4, timiss (NW)*20
        CHARACTER       timstpwrk*20, timstrwrk*20
	INTEGER		npt (NW), itest (NW), iupdt (NW), level(2,NW)
	REAL		rlat (50,NW), rlon (50,NW), tlat (50), tlon (50)
C*
	CHARACTER	stime*20, flstrt*160
	INTEGER		itarr (5), jtarr (5), itype
	INTEGER         iflagsmin(8)
	LOGICAL		done, flvls, filter
C-----------------------------------------------------------------------
	iret = 0
C       These are the corresponding minutes offsets to the IFLAGS
C       parameter settings.
	iflagsmin(IOFFM2HR) = -120
	iflagsmin(IOFFM1HR) = -60
	iflagsmin(IOFF0HR) = 0
	iflagsmin(IOFF1HR) = 60
	iflagsmin(IOFF2HR) = 120
C
C*	Scan the directory for all of the airmet data files.
C
	filnam = 'GAIRM'
	path  = ' '
	templ = ' '
	CALL ST_NULL ( filnam, filnam, nf, ier )
	CALL CTB_DTGET ( filnam, path, templ,
     +			 ic, is, if, ir, ii, ion, ihb, mnb, iha, mna,
     +			 mstrct, idtmch,ier)
	CALL ST_RNUL ( path, path, lens, ier )
	CALL ST_RNUL ( templ, templ, lens, ier )
	CALL ST_LSTR ( path, lenp, ier )
        nexp   = MXNMFL
	iorder = 1
	CALL FL_SCND ( path, templ, iorder, nexp, files, nfile, ier )
C
C*	Check for the last file requested by the user.
C
	CALL ST_LCUC ( dattim, dattim, ier )
	itype = 1
	IF  ( ( dattim .eq. 'LAST' ) .or.
     +	      ( dattim .eq. 'ALL' ) )  THEN
	    CALL CSS_GTIM ( itype, dattm2, ier )
	    CALL FL_MNAM ( dattm2, templ, filnam, ier )
	  ELSE
	    CALL CSS_GTIM ( itype, cdttm, ier )
	    CALL TI_STAN ( dattim, cdttm, dattm2, ier )
            IF ( ier .ne. 0 ) THEN
		CALL ER_WMSG ( 'TI', ier, dattim, ierr )
		iret = ier
		return
	    END IF
	    CALL FL_MNAM ( dattm2, templ, filnam, ier )
	END IF
C
C*	Find the earliest file to start searching. For ALL times
C*	go back 10 days, for any other entry for dattim subtract
C*	1 day from the time given.
C
	IF  ( dattim .eq. 'ALL' )  THEN
	    minuts = 14400
	  ELSE
	    minuts = 1440
	END IF
	CALL TI_CTOI ( dattm2, itarr, ier )
	CALL TI_SUBM ( itarr, minuts, jtarr, ier )
	CALL TI_ITOC ( jtarr, stime, ier )
	CALL FL_MNAM ( stime, templ, flstrt, ier )
C
C*	Get 4-digit year to compare dates.
C
	CALL TI_DTM4 ( dattm2, dattm4, ier )
C
C*	Decode each file until the end time is reached.
C
	nairm = 0
	done  = .false.
	ifl   = 1
	DO WHILE  ( ( ifl .le. nfile ) .and. ( .not. done ) )
	    IF  ( files(ifl) .gt. filnam )  THEN
		done = .true.
	      ELSE IF  ( nairm .ge. NW )  THEN
		done = .true.
	      ELSE
		IF  ( files (ifl) .ge. flstrt )  THEN
		  tfile = path ( :lenp ) // '/' // files ( ifl )
		  CALL FL_SOPN ( tfile, lunf, ier )
C
		  iostat = 0
C
C*		  Loop on the records in the file.
C
		  DO WHILE  ( iostat .eq. 0 )
		    READ ( lunf, 2, IOSTAT = iostat ) buffer
2		    FORMAT ( A )
		    IF  ( iostat .eq. 0 )  THEN
		      IF  ( buffer ( 1:1 ) .eq. '|' )  THEN
			  CALL ST_CLST ( buffer, '|', ' ', 11, carr,
     +					 num, ier )
C
C
			  CALL ST_NUMB ( carr ( 11 ), jcancl, ier )
			  CALL ST_NUMB ( carr ( 10 ), jflag, ier )
			  jcorr = MOD ( jflag, 3 )
			  jtest = jflag / 3
			  ttype = carr ( 2 )
			  tid   = carr ( 6 )
			  CALL ST_NUMB ( carr ( 7 ), jupdt, ier )
			  IF ( ier .lt. 0 ) jupdt = 0
                          DO ii = 1, nairm
				  IF ( (carr (5) .eq. timstp (ii)) .and.
     +				       ( ttype .eq. type (ii) ) .and.
     +				       ( jupdt .ne. iupdt (ii) )) THEN
				      timstp ( ii ) = carr ( 3 )
				  END IF
			  END DO
C
C*			  Read the lat/lon coordinates from the file.
C
			  knt = 0
			  jostat = 0
			  DO WHILE  ( jostat .eq. 0 )
			    READ ( lunf, 2, IOSTAT = jostat ) buffer
			    IF  ( jostat .eq. 0 )  THEN
				IF  ( buffer ( 1:1 ) .eq. '|' )  THEN
				    CALL FL_BKSP ( lunf, ier )
				    jostat = -1
				  ELSE
				    knt = knt + 1
				    READ ( buffer, 1000 )
     +					   tlat ( knt ), tlon ( knt )
1000				    FORMAT ( 2F9.2 )
				END IF
			    END IF
			  END DO
C
C*			  Add this report if it is not a cancellation.
C
			  IF ( jcancl .eq. 0 ) THEN
			    nairm               = nairm + 1
			    type ( nairm )      = ttype
			    timiss ( nairm )    = carr ( 3 )
                            timstr ( nairm )    = carr ( 4 )
			    timstp ( nairm )    = carr ( 5 )
			    waid ( nairm )      = tid
			    flevel ( 1, nairm ) = carr ( 8 )
			    flevel ( 2, nairm ) = carr ( 9 )
			    CALL ST_NUMB (flevel(1, nairm),
     +			                  level(1, nairm), ier)
                            CALL ST_NUMB (flevel(2, nairm),
     +                                    level(2, nairm), ier)
			    itest ( nairm )     = jtest
			    iupdt ( nairm )     = jupdt
			    DO kk = 1, knt
				rlat ( kk, nairm ) = tlat ( kk )
				rlon ( kk, nairm ) = tlon ( kk )
			    END DO
			    npt ( nairm ) = knt
			  END IF
		      END IF
		    END IF
		  END DO
C
		  CALL FL_CLOS ( lunf, ier )
C
		END IF
C
	    END IF
	    ifl = ifl + 1
	END DO
C
C*	Query current settings.
C
	CALL GQCOLR ( jcolr, ier )
	CALL GQLINE ( jtype, jtyhw, jwidth, jwdhw, ier )
	CALL GQTEXT ( jtxfn, jtxhw, siztx, jtxwid, jbrdr,
     +		      jrrotn, jjust, ier )
	CALL GQTURB ( szturb, jtuwid, ier )
	CALL GQICNG ( szicng, jicwid, ier )
	CALL GQSPCL ( szsswn, jswwid, ier )
C
C*	Set attributes.
C
	CALL GSTEXT ( 21, 2, 1.0, 1, 111, 1, 1, ier )
C
	DO  i = 1, NR
	    IF  ( ( ssize (i) .le.  0.0 ) .or.
     +		  ( ssize (i) .gt. 10.0 ) ) THEN
		IF ( ( i .eq. JTRB ) .or. ( i .eq. JICE ) .or.
     +               ( i .eq. JSSW ) .or. ( i .eq. JLWS ) ) THEN
		    ssize (i) = 1.5
		  ELSE
		    ssize (i) = 1.0
		END IF
	    END IF
C
	    IF  ( ( iwidth (i) .le.  0 ) .or.
     +		  ( iwidth (i) .gt. 10 ) ) THEN
		IF ( ( i .eq. JTRB ) .or. ( i .eq. JICE ) .or.
     +               ( i .eq. JSSW ) .or. ( i .eq. JLWS ) ) THEN
		    iwidth (i) = 2
		  ELSE
		    iwidth (i) = 1
		END IF
	    END IF
C
	    IF  ( ( lwidth (i) .le.  0 ) .or.
     +		  ( lwidth (i) .gt. 10 ) )  lwidth (i) = 2
	END DO
C
	syttb = 4.
	sytic = 5.
	syssw = 50.
	strir = 'IFR'
	strmo = 'MTOS'
	strws = 'LLWS'
C
C*	Plot the graphic for each valid report.
C
	DO ip = 1, nairm
C
C*	  Filter the levels.
C
	  filter = .false.
	  IF ( lvfil(1) .gt. 0 ) THEN
	      IF ( level(2,ip) .gt. 0 .and. level(2,ip) .le. lvfil(1))
     +            filter = .true.
	      IF ( level(1,ip) .gt. 0 .and. level(2,ip) .le. 0 .and.
     +		  level(1,ip) .lt. lvfil(1) ) filter = .true.
	  END IF
C
	  IF ( lvfil(2) .gt. 0 ) THEN
	      IF ( level(1,ip) .ge. lvfil(2) ) filter = .true.
              IF ( level(2,ip) .gt. 0 .and. level(1,ip) .le. 0 .and.
     +             level(2,ip) .gt. lvfil(2) ) filter = .true.
          END IF
C*
	  IF ( .not. filter ) THEN
	    DO ihr =  IOFFM2HR, IOFF2HR
	      IF ( iflags(ihr) .NE. 0) THEN
                minuts = iflagsmin(ihr)
		CALL TI_CTOI ( timstp ( ip ), itarr, ier )
		IF (minuts .gt. 0) THEN
		  CALL TI_ADDM ( itarr, minuts, jtarr, ier )
		ELSE
		  CALL TI_SUBM ( itarr, ABS(minuts), jtarr, ier)
		ENDIF
		CALL TI_ITOC ( jtarr, timstpwrk, ier )
		CALL TI_CTOI ( timstr ( ip ), itarr, ier )
		IF (minuts .gt. 0) THEN
		  CALL TI_ADDM ( itarr, minuts, jtarr, ier )
		ELSE
		  CALL TI_SUBM ( itarr, ABS(minuts), jtarr, ier)
		ENDIF
		CALL TI_ITOC ( jtarr, timstrwrk, ierr )
	        CALL TI_DTM4 ( timstpwrk, tmstp4, ier )
	        CALL TI_DTM4 ( timstrwrk, tmstr4, ier )

		IF  ( ( dattim .eq. 'ALL' ) .or.
     +		  ( ( tmstp4 .gt. dattm4 ) .and.
     +		    ( tmstr4 .le. dattm4 ) ) )  THEN
C
C*		Set the line type for a test or an actual report.
C
		    IF ( itest ( ip ) .eq. 0 ) THEN
			iltyp = 1
		      ELSE
			iltyp = 2
		    END IF
C
C*		Set the color and symbol based on the phenomenon type.
C
		    IF  ( type ( ip ) .eq. 'IR' )  THEN
			ic = icolor (JIFR)
			CALL GSLINE ( iltyp, 0, lwidth (JIFR), 0, ier )
C
		      ELSE IF  ( type ( ip ) .eq. 'MO' )  THEN
			ic = icolor (JMTO)
			CALL GSLINE ( iltyp, 0, lwidth (JMTO), 0, ier )
C
		      ELSE IF ( type ( ip ) .eq. 'TB' ) THEN
			ic = icolor (JTRB)
                	IF (lvlfl .ge. 0 ) THEN
			   IF ( level(1,ip) .ge. lvlfl ) ic = itbclr2
			   IF ( level(1,ip) .le. 0 .and. level(2,ip)
     +			    .gt. lvlfl ) ic = itbclr2
                	END IF

			CALL GSTURB ( ssize (JTRB), iwidth (JTRB), ier )
			CALL GSLINE ( iltyp, 0, lwidth (JTRB), 0, ier )
C
		      ELSE IF ( type ( ip ) .eq. 'IC' ) THEN
			ic = icolor (JICE)
			CALL GSICNG ( ssize (JICE), iwidth (JICE), ier )
			CALL GSLINE ( iltyp, 0, lwidth (JICE), 0, ier )
		      ELSE IF ( type ( ip ) .eq. 'SW' ) THEN
			ic = icolor (JSSW)
			CALL GSSPCL ( ssize (JSSW), iwidth (JSSW), ier )
			CALL GSLINE ( iltyp, 0, lwidth (JSSW), 0, ier )
		      ELSE IF ( type ( ip ) .eq. 'WS' ) THEN
			ic = icolor (JLWS)
			CALL GSLINE ( iltyp, 0, lwidth (JLWS), 0, ier )
		    END IF
C
C*		Draw the shape for a bounded phenomenon and get the
C*		label location.
C
		    IF ( ( ic .ne. 0 ) .and. ( npt ( ip ) .ne. 0 ) ) THEN
			npts = npt ( ip )
			CALL GSCOLR ( ic, ier )
			CALL GLINE ( 'M', npts , rlat ( 1, ip ),
     +			         rlon ( 1, ip ), ier )
C
			IF ( ( iflags ( JSYM ) .ne. 0 ) .or.
     +			 ( iflags ( JFLT ) .ne. 0 ) .or.
     +			 ( iflags ( JTIM ) .ne. 0 ) ) THEN
		            iyoff = -2
			    flvls = .false.
		            CALL GG_WLBL ( npts, rlat ( 1, ip ),
     +				       rlon ( 1, ip ), alat, alon, ier )
			    IF ( ( type ( ip ) .eq. 'TB' ) .or.
     +			     ( type ( ip ) .eq. 'IC' ) ) flvls = .true.
C
		            IF ( iflags ( JSYM ) .ne. 0 ) THEN
C
C*		          Plot the symbol.
C
			      IF ( ( iflags ( JFLT ) .gt. 0 ) .and.
     +			       flvls )  iyoff = -3
			      ixoff = 0
		              IF ( type ( ip ) .eq. 'IR' )  THEN
	                	CALL GSTEXT ( 21, 2, ssize(JIFR),
     +					  iwidth(JIFR), 111, 1, 1, ier )
		        	CALL GTEXT ( 'M', alat, alon, strir, 0.,
     +				         ixoff, iyoff, ier )
	                	CALL GSTEXT ( 21, 2, 1.0, 1, 111, 1, 1,
     +					  ier )
		               ELSE IF ( type ( ip ) .eq. 'MO' )  THEN
	                	CALL GSTEXT ( 21, 2, ssize(JMTO),
     +					  iwidth(JMTO), 111, 1, 1, ier )
		        	CALL GTEXT ( 'M', alat, alon, strmo, 0.,
     +				         ixoff, iyoff, ier )
	                	CALL GSTEXT ( 21, 2, 1.0, 1, 111, 1, 1,
     +					  ier )
		               ELSE IF ( type ( ip ) .eq. 'TB' )  THEN
		        	CALL GTURB ( 'M', 1, syttb, alat, alon,
     +				         ixoff, iyoff, ier )
		               ELSE IF ( type ( ip ) .eq. 'IC' )  THEN
	                	CALL GICNG ( 'M', 1, sytic, alat, alon,
     +				         ixoff, iyoff, ier )
		               ELSE IF ( type ( ip ) .eq. 'SW' )  THEN
	                	CALL GSPCL ( 'M', 1, syssw, alat, alon,
     +				         ixoff, iyoff, ier )
				iyoff = iyoff - 2
		               ELSE IF ( type ( ip ) .eq. 'WS' )  THEN
				CALL GSTEXT ( 21, 2, ssize(JLWS),
     +                                    iwidth(JLWS), 111, 1, 1, ier )
                        	CALL GTEXT ( 'M', alat, alon, strws, 0.,
     +                                   ixoff, iyoff, ier )
                        	CALL GSTEXT ( 21, 2, 1.0, 1, 111, 1, 1,
     +                                    ier )
		              END IF
			      iyoff = iyoff - 2
		            END IF
C
		            IF ( iflags ( JFLT ) .ne. 0 ) THEN
			      IF ( flvls ) THEN
C
C*			    Plot the flight level(s).
C
				ixoff = 4
				iyoff = -5
				DO ii = 1, 2
		                    label  = flevel ( ii, ip )
		                    CALL GTEXT ( 'M', alat, alon, label, 0.,
     +				             ixoff, iyoff, ier )
				    iyoff = iyoff + 4
				END DO
				IF ( flevel ( 1, ip ) .ne. ' ' ) THEN
				    iyoff = -7
				  ELSE
				    iyoff = -5
				END IF
			      END IF
		            END IF
C
		            IF ( iflags ( JTIM ) .ne. 0 ) THEN
C
C*		          Plot the time.
C
		              label  = timstr ( ip ) ( 8:11 )
		              CALL GTEXT ( 'M', alat, alon, label, 0.0, 0,
     +					iyoff, ier )
			    END IF
			END IF
		    END IF
		END IF
	      END IF
	    END DO
C*********************************************************************
	  END IF
	END DO
C
C*	Reset the saved attributes.
C
	CALL GSCOLR ( jcolr, ier )
	CALL GSLINE ( jtype, jtyhw, jwidth, jwdhw, ier )
	CALL GSTEXT ( jtxfn, jtxhw, siztx, jtxwid, jbrdr,
     +		      jrrotn, jjust, ier )
	CALL GSTURB ( szturb, jtuwid, ier )
	CALL GSICNG ( szicng, jicwid, ier )
	CALL GSSPCL ( szsswn, jswwid, ier )
C*
	RETURN
	END

ext: f

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


!========================================================================
!
!                   S P E C F E M 2 D  Version 7 . 0
!                   --------------------------------
!
!     Main historical authors: Dimitri Komatitsch and Jeroen Tromp
!                        Princeton University, USA
!                and CNRS / University of Marseille, France
!                 (there are currently many more authors!)
! (c) Princeton University and CNRS / University of Marseille, April 2014
!
! This software is a computer program whose purpose is to solve
! the two-dimensional viscoelastic anisotropic or poroelastic wave equation
! using a spectral-element method (SEM).
!
! This software is governed by the CeCILL license under French law and
! abiding by the rules of distribution of free software. You can use,
! modify and/or redistribute the software under the terms of the CeCILL
! license as circulated by CEA, CNRS and Inria at the following URL
! "http://www.cecill.info".
!
! As a counterpart to the access to the source code and rights to copy,
! modify and redistribute granted by the license, users are provided only
! with a limited warranty and the software's author, the holder of the
! economic rights, and the successive licensors have only limited
! liability.
!
! In this respect, the user's attention is drawn to the risks associated
! with loading, using, modifying and/or developing or reproducing the
! software by the user in light of its specific status of free software,
! that may mean that it is complicated to manipulate, and that also
! therefore means that it is reserved for developers and experienced
! professionals having in-depth computer knowledge. Users are therefore
! encouraged to load and test the software's suitability as regards their
! requirements in conditions enabling the security of their systems and/or
! data to be ensured and, more generally, to use and operate it in the
! same conditions as regards security.
!
! The full text of the license is available in file "LICENSE".
!
!========================================================================

! From array 'surface' (element, type : node/edge, node(s) ) that describes the
! elastic fixed surface, determines the points (ixmin, ixmax, izmin and izmax) on the surface
! for each element.
! We chose to have ixmin <= ixmax and izmin <= izmax, so as to be able to have DO loops on it with
! an increment of +1.
!
subroutine construct_elastic_fixed_surface ()

  use specfem_par, only : ngnod, knods, nelem_elastic_fixed_surface, elastic_fixed_edges, elastic_fixed_surface

  implicit none

  integer  :: i, k
  integer  :: ixmin, ixmax
  integer  :: izmin, izmax
  integer, dimension(ngnod)  :: n
  integer  :: e1, e2
  integer  :: type_elastic

  do i = 1, nelem_elastic_fixed_surface
     elastic_fixed_surface(1,i) = elastic_fixed_edges(1,i) ! Here we do a copy
     type_elastic = elastic_fixed_edges(2,i)
     e1 = elastic_fixed_edges(3,i)
     e2 = elastic_fixed_edges(4,i)
     do k = 1, ngnod
        n(k) = knods(k,elastic_fixed_surface(1,i))
     enddo

     call get_elastic_fixed_edge ( ngnod, n, type_elastic, e1, e2, ixmin, ixmax, izmin, izmax )

     elastic_fixed_surface(2,i) = ixmin
     elastic_fixed_surface(3,i) = ixmax
     elastic_fixed_surface(4,i) = izmin
     elastic_fixed_surface(5,i) = izmax

  enddo

end subroutine construct_elastic_fixed_surface

!-----------------------------------------------
! Get the points (ixmin, ixmax, izmin and izmax) on an node/edge for one element.
!-----------------------------------------------
subroutine get_elastic_fixed_edge ( ngnod, n, type, e1, e2, ixmin, ixmax, izmin, izmax )

  implicit none
  include "constants.h"

  integer, intent(in)  :: ngnod
  integer, dimension(ngnod), intent(in)  :: n
  integer, intent(in)  :: type, e1, e2
  integer, intent(out)  :: ixmin, ixmax, izmin, izmax

  if ( type == 1 ) then
     if ( e1 == n(1) ) then
        ixmin = 1
        ixmax = 1
        izmin = 1
        izmax = 1
     endif
     if ( e1 == n(2) ) then
        ixmin = NGLLX
        ixmax = NGLLX
        izmin = 1
        izmax = 1
     endif
     if ( e1 == n(3) ) then
        ixmin = NGLLX
        ixmax = NGLLX
        izmin = NGLLZ
        izmax = NGLLZ
     endif
     if ( e1 == n(4) ) then
        ixmin = 1
        ixmax = 1
        izmin = NGLLZ
        izmax = NGLLZ
     endif

  else
     if ( e1 ==  n(1) ) then
        ixmin = 1
        izmin = 1
        if ( e2 == n(2) ) then
           ixmax = NGLLX
           izmax = 1

        endif
        if ( e2 == n(4) ) then
           ixmax = 1
           izmax = NGLLZ

        endif
     endif
     if ( e1 == n(2) ) then
        ixmin = NGLLX
        izmin = 1
        if ( e2 == n(3) ) then
           ixmax = NGLLX
           izmax = NGLLZ

        endif
        if ( e2 == n(1) ) then
           ixmax = ixmin
           ixmin = 1
           izmax = 1

        endif
     endif
     if ( e1 == n(3) ) then
        ixmin = NGLLX
        izmin = NGLLZ
        if ( e2 == n(4) ) then
           ixmax = ixmin
           ixmin = 1
           izmax = NGLLZ

        endif
        if ( e2 == n(2) ) then
           ixmax = NGLLX
           izmax = izmin
           izmin = 1

        endif
     endif
     if ( e1 == n(4) ) then
        ixmin = 1
        izmin = NGLLZ
        if ( e2 == n(1) ) then
           ixmax = 1
           izmax = izmin
           izmin = 1

        endif
        if ( e2 == n(3) ) then
           ixmax = NGLLX
           izmax = NGLLZ

        endif
     endif
  endif

end subroutine get_elastic_fixed_edge

!-----------------------------------------------
! set the accel_elastic, veloc_elastic, and displ_elastic on the points (ixmin, ixmax, izmin and izmax)
! that on an node/edge for one element to be zero
!-----------------------------------------------

subroutine enforce_elastic_fixed_surface(accel_elastic,veloc_elastic,displ_elastic)

  use specfem_par, only: elastic_fixed_surface,ibool,nelem_elastic_fixed_surface,nglob,&
                         this_ibool_is_a_periodic_edge

  implicit none

  include "constants.h"

  real(kind=CUSTOM_REAL), dimension(3,nglob) :: accel_elastic,veloc_elastic,displ_elastic

!---
!--- local variables
!---

  integer :: ispec_elastic_fixed_surface,ispec,i,j,iglob

  do ispec_elastic_fixed_surface = 1, nelem_elastic_fixed_surface

    ispec = elastic_fixed_surface(1,ispec_elastic_fixed_surface)

    do j = elastic_fixed_surface(4,ispec_elastic_fixed_surface), elastic_fixed_surface(5,ispec_elastic_fixed_surface)
      do i = elastic_fixed_surface(2,ispec_elastic_fixed_surface), elastic_fixed_surface(3,ispec_elastic_fixed_surface)
        iglob = ibool(i,j,ispec)
        ! make sure that an acoustic free surface is not enforced on periodic edges
        if(.not. this_ibool_is_a_periodic_edge(iglob)) then
          displ_elastic(:,iglob) = ZERO
          veloc_elastic(:,iglob) = ZERO
          accel_elastic(:,iglob) = ZERO
        endif
      enddo
    enddo

  enddo

  end subroutine enforce_elastic_fixed_surface