text:
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!
!                          Futility Development Group                          !
!                             All rights reserved.                             !
!                                                                              !
! Futility is a jointly-maintained, open-source project between the University !
! of Michigan and Oak Ridge National Laboratory.  The copyright and license    !
! can be found in LICENSE.txt in the head directory of this repository.        !
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!
!> @brief Module defines objects for representing an XDMF file as a hierarchical
!> mesh
!>
!> This module reads an XDMF file and stores the information in a hierarchical
!> mesh type. It can also write the hierarchical mesh to XDMF.
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++!
MODULE XDMFMesh
#include "Futility_DBC.h"
USE ExceptionHandler
USE Futility_DBC
USE ISO_FORTRAN_ENV
USE IntrType
USE Strings
USE FileType_XML
USE FileType_HDF5
USE ParameterLists

IMPLICIT NONE
PRIVATE

! Public members
PUBLIC :: eXDMF
#ifdef FUTILITY_HAVE_HDF5
PUBLIC :: XDMFMeshType
PUBLIC :: XDMFTopologyList
PUBLIC :: ImportXDMFMesh
PUBLIC :: ExportXDMFMesh
PUBLIC :: ASSIGNMENT(=)
!> The module name
CHARACTER(LEN=*),PARAMETER :: modName='FILETYPE_XDMF'
#endif

!> Exception handler for the module
TYPE(ExceptionHandlerType),SAVE :: eXDMF

#ifdef FUTILITY_HAVE_HDF5
!> Parameter list that holds XDMF topology names, ids, etc.
TYPE(ParamType),SAVE :: XDMFTopologyList

!> Type to hold the vertices that make up a mesh cell
TYPE :: XDMFCell
  !> The cell type id followed by the vertex ids
  !> XDMF ID, v1, v2, ..., v_n
  INTEGER(SLK), ALLOCATABLE :: vertex_list(:)
ENDTYPE XDMFCell

!> Type to hold a list of cell IDs that make up a named set.
TYPE :: XDMFCellSet
  !> The name of the set
  TYPE(StringType) :: name
  !> The cell IDs
  INTEGER(SLK), ALLOCATABLE :: cell_list(:)
ENDTYPE XDMFCellSet

!> Type to hold the XDMF mesh data
TYPE :: XDMFMeshType
  !> The name of the mesh
  TYPE(StringType) :: name
  !> If the mesh cells are all the same topology
  LOGICAL(SBK) :: singleTopology=.FALSE.
  !> The bounding box for the mesh
  !> xmin, xmax, ymin, ymax
  REAL(SDK) :: boundingBox(4) = 0.0_SDK
  !> The vertices that compose the mesh
  !> Looks like:
  !> x1, x2, x3, ..., xn
  !> y1, y2, y3, ..., yn
  !> z1, z2, z3, ..., zn
  !> Therefore vertices will be of shape (3, N)
  REAL(SDK), ALLOCATABLE :: vertices(:, :)
  !> The mesh cells
  TYPE(XDMFCell), ALLOCATABLE :: cells(:)
  !> Material for each mesh cell
  INTEGER(SNK), ALLOCATABLE :: material_ids(:)
  !> Named sets within the mesh
  TYPE(XDMFCellSet), ALLOCATABLE :: cell_sets(:)
  !> Parent mesh
  TYPE(XDMFMeshType), POINTER :: parent => NULL()
  !> Child meshes
  TYPE(XDMFMeshType), POINTER :: children(:) => NULL()
  CONTAINS
    !> @copybrief XDMFMeshType::clear_XDMFMeshType
    !> @copydoc XDMFMeshType::clear_XDMFMeshType
    PROCEDURE,PASS :: clear => clear_XDMFMeshType
    !> @copybrief XDMFMeshType::distanceToLeaf_XDMFMeshType
    !> @copydoc XDMFMeshType::distanceToLeaf_XDMFMeshType
    PROCEDURE,PASS :: distanceToLeaf => distanceToLeaf_XDMFMeshType
    !> @copybrief XDMFMeshType::recomputeBoundingBox_XDMFMeshType
    !> @copydoc XDMFMeshType::recomputeBoundingBox_XDMFMeshType
    PROCEDURE,PASS :: recomputeBoundingBox => recomputeBoundingBox_XDMFMeshType
ENDTYPE XDMFMeshType

!> @brief Interface for assignment operator (=)
INTERFACE ASSIGNMENT(=)
  !> @copybrief FileType_XDMF::assign_XDMFMeshType
  !> @copydoc FileType_XDMF::assign_XDMFMeshType
  MODULE PROCEDURE assign_XDMFMeshType
ENDINTERFACE
!
!===============================================================================
CONTAINS
!
!-------------------------------------------------------------------------------
!> @brief Initializes the XDMFTopologyList
!>
SUBROUTINE init_XDMFTopologyList()
  ! Setup param list for cell type conversions
  ! id is XDMF topology id,
  ! n is number of vertices,
  ! multiple valid names exist for the same topology, ex: Tri_6 == Triangle_6
  CALL XDMFTopologyList%add('Topology->Triangle->id'            , 4_SNK)
  CALL XDMFTopologyList%add('Topology->Triangle->n'             , 3_SNK)
  CALL XDMFTopologyList%add('Topology->Triangle_6->id'          ,36_SNK)
  CALL XDMFTopologyList%add('Topology->Triangle_6->n'           , 6_SNK)
  CALL XDMFTopologyList%add('Topology->Tri_6->id'               ,36_SNK)
  CALL XDMFTopologyList%add('Topology->Tri_6->n'                , 6_SNK)
  CALL XDMFTopologyList%add('Topology->Quadrilateral->id'       , 5_SNK)
  CALL XDMFTopologyList%add('Topology->Quadrilateral->n'        , 4_SNK)
  CALL XDMFTopologyList%add('Topology->Quadrilateral_8->id'     ,37_SNK)
  CALL XDMFTopologyList%add('Topology->Quadrilateral_8->n'      , 8_SNK)
  CALL XDMFTopologyList%add('Topology->Quad_8->id'              ,37_SNK)
  CALL XDMFTopologyList%add('Topology->Quad_8->n'               , 8_SNK)
  CALL XDMFTopologyList%add('XDMFID->4' ,'Triangle'       )
  CALL XDMFTopologyList%add('XDMFID->36','Triangle_6'     )
  CALL XDMFTopologyList%add('XDMFID->5' ,'Quadrilateral'  )
  CALL XDMFTopologyList%add('XDMFID->37','Quadrilateral_8')
ENDSUBROUTINE init_XDMFTopologyList
!
!-------------------------------------------------------------------------------
!> @brief Returns the hdf5 group where heavy data is stored
!> @param xmle XML element storing the file path and group to the heavy data
!> @returns group The hdf5 group where the heavy data is stored
!>
FUNCTION getH5GroupFromXMLContent(xmle) RESULT(group)
  TYPE(XMLElementType),INTENT(IN) :: xmle
  TYPE(StringType) :: content, group
  TYPE(StringType), ALLOCATABLE :: segments(:)
  ! Content of xmle should be the h5 filename with the
  ! path to the data
  content=xmle%getContent()
  ! Split file from group data
  segments=content%split(':')
  ! Just grab the group data
  group=segments(2)%substr(2,LEN(segments(2)))
  ! Replace newline char
  group = group%replace(NEW_LINE("A"),"")
ENDFUNCTION getH5GroupFromXMLContent
!
!-------------------------------------------------------------------------------
!> @brief Create the XDMF mesh object
!> @param mesh the parent mesh
!> @param xmle the child XML element
!> @param h5 the HDF5 fil containing mesh data
!>
RECURSIVE SUBROUTINE create_XDMFMesh_from_file(mesh, xmle, h5)
  CHARACTER(LEN=*),PARAMETER :: myName='create_XDMFMesh_from_file'
  TYPE(XDMFMeshType),TARGET, INTENT(INOUT)  :: mesh
  TYPE(XMLElementType), INTENT(INOUT) :: xmle
  TYPE(HDF5FileType), INTENT(INOUT) :: h5
  TYPE(XMLElementType), POINTER :: xmle_children(:)
  TYPE(StringType) :: strIn, strOut
  INTEGER(SIK) :: i, nGrid, iMesh

  ! If this xml element has children
  IF(xmle%hasChildren()) THEN
    ! Determine the number or XML children that are grids
    CALL xmle%getChildren(xmle_children)
    nGrid = 0
    DO i=1,SIZE(xmle_children)
      IF(xmle_children(i)%name%upper() == 'GRID') nGrid = nGrid + 1
    ENDDO

    ! If some children are grids, it is not a leaf
    IF(nGrid > 0) THEN
      ! Allocate children for current mesh and create entities
      ALLOCATE(mesh%children(nGrid))
      iMesh=1
      ! Recursively create children of each grid
      DO i=1,SIZE(xmle_children)
        IF(xmle_children(i)%name%upper() == 'GRID') THEN
          strIn='Name'
          CALL xmle_children(i)%getAttributeValue(strIn,strOut)
          mesh%children(iMesh)%name = strOut
          mesh%children(iMesh)%parent => mesh
          CALL create_XDMFMesh_from_file(mesh%children(iMesh), xmle_children(i), h5)
          iMesh = iMesh + 1
        ENDIF
      ENDDO
    ! If this mesh does not have grid children it is a leaf on the tree.
    ! Add vertices, cells, etc.
    ELSE
      CALL setup_leaf_XDMFMesh_from_file(mesh, xmle, h5)
    ENDIF
  ELSE
    CALL eXDMF%raiseError(modName//'::'//myName// &
      ' - Expected the XML element '//mesh%name//' to have children.')
  ENDIF
ENDSUBROUTINE create_XDMFMesh_from_file
!
!-------------------------------------------------------------------------------
!> @brief Setup the leaf mesh objects which contain vertices, cells, etc.
!> @param mesh the parent mesh
!> @param xmle the child XML element
!> @param h5 the HDF5 file containing mesh data
!>
SUBROUTINE setup_leaf_XDMFMesh_from_file(mesh, xmle, h5)
  CHARACTER(LEN=*),PARAMETER :: myName='setup_leaf_XDMFMesh_from_file'
  TYPE(XDMFMeshType),INTENT(INOUT),TARGET  :: mesh
  TYPE(XMLElementType), INTENT(INOUT) :: xmle
  TYPE(HDF5FileType), INTENT(INOUT) :: h5
  TYPE(XMLElementType), POINTER :: xmle_children(:), ele_children(:)
  TYPE(StringType) :: elname, strIn, strOut, group, dtype, toponame, &
    xdmf_id_str
  TYPE(StringType),ALLOCATABLE :: strArray(:)
  INTEGER(SLK) :: nverts, ncells,ivert,i,j
  INTEGER(SNK) :: ncell_sets, xdmf_id, verts_per_cell
  INTEGER(SNK),ALLOCATABLE :: dshape(:)
  REAL(SSK),ALLOCATABLE :: vals4_2d(:,:)
  REAl(SDK),ALLOCATABLE :: vals8_2d(:,:)
  INTEGER(SNK),ALLOCATABLE :: ivals4_1d(:),ivals4_2d(:,:)
  INTEGER(SLK),ALLOCATABLE :: ivals8_1d(:),ivals8_2d(:,:)
  TYPE(XDMFCellSet), ALLOCATABLE :: cell_sets_temp(:)

  IF(.NOT.xmle%hasChildren())  CALL eXDMF%raiseError(modName//'::'//myName// &
    ' - XML element for '//CHAR(mesh%name)//' should have children.')
  CALL xmle%getChildren(xmle_children)
  ! Each XML element has a type of information.
  ! Handle each with a CASE
  DO i=1,SIZE(xmle_children)
    elname=xmle_children(i)%name%upper()
    SELECTCASE(ADJUSTL(elname))
    CASE("GEOMETRY")
      ! GeometryType
      strIn='GeometryType'
      CALL xmle_children(i)%getAttributeValue(strIn,strOut)
      IF(strOut /= 'XYZ') THEN
        CALL eXDMF%raiseError(modName//'::'//myName// &
          ' - GeometryType only supports XYZ right now.')
      ENDIF
      ! Format
      CALL xmle_children(i)%getChildren(ele_children)
      IF(SIZE(ele_children) /= 1) CALL eXDMF%raiseError(modName//'::'//myName//&
        ' - Expected Geometry to have only one child.')
      strIn='Format'
      CALL ele_children(1)%getAttributeValue(strIn,strOut)
      IF(strOut /= 'HDF') THEN
        CALL eXDMF%raiseError(modName//'::'//myName// &
          ' - only supports HDF5 geometry data right now.')
      ENDIF
      ! Vertex Data
      strIn='Dimensions'
      CALL ele_children(1)%getAttributeValue(strIn,strOut)
      strArray=strOut%split()
      IF(strArray(2) /= '3') CALL eXDMF%raiseError(modName//'::'//myName//&
        ' - Expected vertex data to be 3 dimensional.')
      nverts=strArray(1)%stoi()
      group=getH5GroupFromXMLContent(ele_children(1))
      ! Make sure the h5 path exists
      IF(.NOT.h5%pathExists(CHAR(group)))THEN
        CALL eXDMF%raiseError(modName//'::'//myName//&
          ' - HDF5 group containing vertex data does not exist in h5 file.')
      ENDIF
      group = group%replace("/", "->")
      ! Data shape
      dshape=h5%getDataShape(CHAR(group))
      IF(.NOT.(dshape(1) == 3 .AND. dshape(2) == nverts))THEN
        CALL eXDMF%raiseError(modName//'::'//myName//&
          ' - HDF5 vertex data shape does not match XDMF vertex data shape.')
      ENDIF
      ! Data type
      dtype=h5%getDataType(CHAR(group))
      IF(dtype == 'SSK') THEN
        CALL h5%fread(CHAR(group),vals4_2d)
      ELSE
        CALL h5%fread(CHAR(group),vals8_2d)
      ENDIF
      ALLOCATE(mesh%vertices(3,nverts))
      IF(dtype == 'SSK') THEN
        mesh%vertices=vals4_2d
        DEALLOCATE(vals4_2d)
      ELSE
        mesh%vertices=vals8_2d
        DEALLOCATE(vals8_2d)
      ENDIF
    CASE("TOPOLOGY")
      ! TopologyType
      strIn='TopologyType'
      CALL xmle_children(i)%getAttributeValue(strIn,toponame)
      IF(toponame%upper() == 'MIXED') THEN
        ! Mixed topology
        ! Format
        CALL xmle_children(i)%getChildren(ele_children)
        IF(SIZE(ele_children) /= 1) CALL eXDMF%raiseError(modName//'::'//myName//&
          ' - Expected Topology to have only one child.')
        strIn='Format'
        CALL ele_children(1)%getAttributeValue(strIn,strOut)
        IF(strOut /= 'HDF') THEN
          CALL eXDMF%raiseWarning(modName//'::'//myName// &
            ' - only supports HDF5 topology data right now.')
        ENDIF
        ! Topology Data
        strIn='NumberOfElements'
        CALL xmle_children(i)%getAttributeValue(strIn,strOut)
        ncells=strOut%stoi()
        group=getH5GroupFromXMLContent(ele_children(1))
        ! Make sure the h5 path exists
        IF(.NOT.h5%pathExists(CHAR(group)))THEN
          CALL eXDMF%raiseError(modName//'::'//myName//&
            ' - HDF5 group containing topology data does not exist in h5 file.')
        ENDIF
        group = group%replace("/", "->")
        ! Data shape
        dshape=h5%getDataShape(CHAR(group))
        IF(SIZE(dshape) /= 1)THEN
          CALL eXDMF%raiseError(modName//'::'//myName//&
            ' - HDF5 mixed topology data shape does not match XDMF data shape.')
        ENDIF
        ! Data type
        dtype=h5%getDataType(CHAR(group))
        IF(dtype == 'SNK') THEN
          CALL h5%fread(CHAR(group),ivals4_1d)
        ELSE
          CALL h5%fread(CHAR(group),ivals8_1d)
        ENDIF
        ALLOCATE(mesh%cells(ncells))
        ivert = 1
        IF(dtype == 'SNK') THEN
          DO j=1,ncells
            xdmf_id = ivals4_1d(ivert)
            xdmf_id_str = xdmf_id
            IF(.NOT.XDMFTopologyList%has('XDMFID->'//ADJUSTL(xdmf_id_str))) THEN
              CALL eXDMF%raiseError(modName//'::'//myName//&
                ' - Topology type '//TRIM(xdmf_id_str)//' not supported')
            ELSE
              CALL XDMFTopologyList%get('XDMFID->'//ADJUSTL(xdmf_id_str), toponame)
              CALL XDMFTopologyList%get(ADJUSTL(toponame)//'->n', nverts)
            ENDIF
            ALLOCATE(mesh%cells(j)%vertex_list(nverts+1))
            mesh%cells(j)%vertex_list(1) = xdmf_id
            mesh%cells(j)%vertex_list(2:nverts+1) = ivals4_1d(ivert:ivert+nverts) + 1
            ivert = ivert + nverts
          ENDDO
          DEALLOCATE(ivals4_1d)
        ELSE
          DO j=1,ncells
            ! This number should be well below HUGE(ivals8_1d), so narrowing
            ! is not a concern
            xdmf_id = ivals8_1d(ivert)
            xdmf_id_str = xdmf_id
            IF(.NOT.XDMFTopologyList%has('XDMFID->'//ADJUSTL(xdmf_id_str))) THEN
              CALL eXDMF%raiseError(modName//'::'//myName//&
                ' - Topology type '//TRIM(xdmf_id_str)//' not supported')
            ELSE
              CALL XDMFTopologyList%get('XDMFID->'//ADJUSTL(xdmf_id_str), toponame)
              CALL XDMFTopologyList%get(ADJUSTL(toponame)//'->n', &
                verts_per_cell)
            ENDIF
            ALLOCATE(mesh%cells(j)%vertex_list(verts_per_cell+1))
            mesh%cells(j)%vertex_list(1) = xdmf_id
            mesh%cells(j)%vertex_list(2:verts_per_cell+1) = &
              ivals8_1d(ivert+1:ivert+verts_per_cell) + 1
            ivert = ivert + verts_per_cell + 1
          ENDDO
          DEALLOCATE(ivals8_1d)
        ENDIF
      ELSE
        ! Single topology
        IF(.NOT.XDMFTopologyList%has(CHAR(toponame))) CALL eXDMF%raiseError(modName// &
          '::'//myName//' - Topology type '//TRIM(strOut)//' not supported')
        ! XDMF ID
        CALL XDMFTopologyList%get(CHAR(toponame)//'->id', xdmf_id)
        ! Format
        CALL xmle_children(i)%getChildren(ele_children)
        IF(SIZE(ele_children) /= 1) CALL eXDMF%raiseError(modName//'::'//myName//&
          ' - Expected Topology to have only one child.')
        strIn='Format'
        CALL ele_children(1)%getAttributeValue(strIn,strOut)
        IF(strOut /= 'HDF') THEN
          CALL eXDMF%raiseWarning(modName//'::'//myName// &
            ' - only supports HDF5 topology data right now.')
        ENDIF
        ! Topology Data
        strIn='NumberOfElements'
        CALL xmle_children(i)%getAttributeValue(strIn,strOut)
        ncells=strOut%stoi()
        strIn='NodesPerElement'
        CALL xmle_children(i)%getAttributeValue(strIn,strOut)
        nverts=strOut%stoi()
        group=getH5GroupFromXMLContent(ele_children(1))
        ! Make sure the h5 path exists
        IF(.NOT.h5%pathExists(CHAR(group)))THEN
          CALL eXDMF%raiseError(modName//'::'//myName//&
            ' - HDF5 group containing topology data does not exist in h5 file.')
        ENDIF
        group = group%replace("/", "->")
        ! Data shape
        dshape=h5%getDataShape(CHAR(group))
        IF(.NOT.(dshape(1) == nverts .AND. dshape(2) == ncells))THEN
          CALL eXDMF%raiseError(modName//'::'//myName//&
            ' - HDF5 mixed topology data shape does not match XDMF data shape.')
        ENDIF
        ! Data type
        dtype=h5%getDataType(CHAR(group))
        IF(dtype == 'SNK') THEN
          CALL h5%fread(CHAR(group),ivals4_2d)
        ELSE
          CALL h5%fread(CHAR(group),ivals8_2d)
        ENDIF
        ALLOCATE(mesh%cells(ncells))
        IF(dtype == 'SNK') THEN
          DO j=1,ncells
            ALLOCATE(mesh%cells(j)%vertex_list(nverts + 1))
            mesh%cells(j)%vertex_list(1) = xdmf_id
            ! Account for 0 based to 1 based index switch
            mesh%cells(j)%vertex_list(2:) = ivals4_2d(:, j) + 1
          ENDDO
          DEALLOCATE(ivals4_2d)
        ELSE
          DO j=1,ncells
            ALLOCATE(mesh%cells(j)%vertex_list(nverts + 1))
            mesh%cells(j)%vertex_list(1) = xdmf_id
            ! Account for 0 based to 1 based index switch
            mesh%cells(j)%vertex_list(2:) = ivals8_2d(:, j) + 1
          ENDDO
          DEALLOCATE(ivals8_2d)
        ENDIF
        mesh%singleTopology = .TRUE.
      ENDIF
    CASE("ATTRIBUTE")
      strIn='Name'
      CALL xmle_children(i)%getAttributeValue(strIn,strOut)
      IF(strOut%upper() == 'MATERIALID') THEN
        ! Format
        CALL xmle_children(i)%getChildren(ele_children)
        IF(SIZE(ele_children) /= 1) CALL eXDMF%raiseError(modName//'::'//myName//&
          ' - Expected Attribute to have only one child.')
        strIn='Format'
        CALL ele_children(1)%getAttributeValue(strIn,strOut)
        IF(strOut /= 'HDF') THEN
          CALL eXDMF%raiseWarning(modName//'::'//myName// &
            ' - only supports HDF5 material data right now.')
        ENDIF
        ! Material Data
        strIn='Dimensions'
        CALL ele_children(1)%getAttributeValue(strIn,strOut)
        ncells=strOut%stoi()
        IF(.NOT.(ALLOCATED(mesh%cells) .AND. ncells == SIZE(mesh%cells))) THEN
          CALL eXDMF%raiseError(modName//'::'//myName//&
            ' - material data is before topology data, or is the wrong size.')
        ENDIF
        group=getH5GroupFromXMLContent(ele_children(1))
        ! Make sure the h5 path exists
        IF(.NOT.h5%pathExists(CHAR(group)))THEN
          CALL eXDMF%raiseError(modName//'::'//myName//&
            ' - HDF5 group containing material data does not exist in h5 file.')
        ENDIF
        group = group%replace("/", "->")
        ! Data shape
        dshape=h5%getDataShape(CHAR(group))
        IF(.NOT.(SIZE(dshape) == 1 .AND. dshape(1) == ncells)) THEN
          CALL eXDMF%raiseError(modName//'::'//myName//&
            ' - material data in h5 file is the wrong size or shape.')
        ENDIF
        ! Data type
        dtype=h5%getDataType(CHAR(group))
        IF(dtype == 'SNK') THEN
          CALL h5%fread(CHAR(group),ivals4_1d)
        ELSE
          CALL h5%fread(CHAR(group),ivals8_1d)
        ENDIF
        ALLOCATE(mesh%material_ids(ncells))
        IF(dtype == 'SNK') THEN
          ! Account for 0 based to 1 based index switch
          mesh%material_ids = ivals4_1d + 1
          DEALLOCATE(ivals4_1d)
        ELSE
          ! Account for 0 based to 1 based index switch
          ! material ids will not exceed MAX(INTEGER(4)),
          ! so narrowing will not occur.
          mesh%material_ids = ivals8_1d + 1
          DEALLOCATE(ivals8_1d)
        ENDIF
      ELSE
        CALL eXDMF%raiseWarning(modName//'::'//myName//' - mesh attribute '//&
          TRIM(strOut)//' not supported')
      ENDIF
    CASE("SET")
      ! SetType
      strIn='SetType'
      CALL xmle_children(i)%getAttributeValue(strIn,strOut)
      IF(strOut /= 'Cell') THEN
        CALL eXDMF%raiseWarning(modName//'::'//myName// &
          ' - only supports SetType="Cell" right now.')
      ENDIF
      ! SetName
      strIn='Name'
      CALL xmle_children(i)%getAttributeValue(strIn,elname)
      ! Format
      CALL xmle_children(i)%getChildren(ele_children)
      IF(SIZE(ele_children) /= 1) CALL eXDMF%raiseError(modName//'::'//myName//&
        ' - Expected Set to have only one child.')
      strIn='Format'
      strIn='Format'
      CALL ele_children(1)%getAttributeValue(strIn,strOut)
      IF(strOut /= 'HDF') THEN
        CALL eXDMF%raiseWarning(modName//'::'//myName// &
          ' - only supports HDF5 cell set data right now.')
      ENDIF
      ! Cell Set Data
      strIn='Dimensions'
      CALL ele_children(1)%getAttributeValue(strIn,strOut)
      ncells=strOut%stoi()
      IF(.NOT.(ALLOCATED(mesh%cells) .AND. ncells <= SIZE(mesh%cells))) THEN
        CALL eXDMF%raiseError(modName//'::'//myName//&
          ' - material data is before topology data, or is too big.')
      ENDIF
      group=getH5GroupFromXMLContent(ele_children(1))
      ! Make sure the h5 path exists
      IF(.NOT.h5%pathExists(CHAR(group)))THEN
        CALL eXDMF%raiseError(modName//'::'//myName//&
          ' - HDF5 group containing set data does not exist in h5 file.')
      ENDIF
      group = group%replace("/", "->")
      ! Data shape
      dshape=h5%getDataShape(CHAR(group))
      IF(.NOT.(SIZE(dshape) == 1 .AND. dshape(1) == ncells)) THEN
        CALL eXDMF%raiseError(modName//'::'//myName//&
          ' - set data in h5 file is the wrong size or shape.')
      ENDIF
      ! Data type
      dtype=h5%getDataType(CHAR(group))
      IF(dtype == 'SNK') THEN
        CALL h5%fread(CHAR(group),ivals4_1d)
      ELSE
        CALL h5%fread(CHAR(group),ivals8_1d)
      ENDIF
      ! Resize cell sets if needed
      ! This is expected to happen infrequently
      IF(ALLOCATED(mesh%cell_sets)) THEN
        ! Copy current sets to temp, deallocate current sets
        ncell_sets = SIZE(mesh%cell_sets)
        ALLOCATE(cell_sets_temp(ncell_sets))
        DO j=1, ncell_sets
          ALLOCATE(cell_sets_temp(j)%cell_list(SIZE(mesh%cell_sets(j)%cell_list)))
          cell_sets_temp(j)%cell_list = mesh%cell_sets(j)%cell_list
          cell_sets_temp(j)%name = mesh%cell_sets(j)%name
          DEALLOCATE(mesh%cell_sets(j)%cell_list)
        ENDDO
        DEALLOCATE(mesh%cell_sets)
        ! Reallocate cell sets to be on bigger and copy all old sets over
        ALLOCATE(mesh%cell_sets(ncell_sets+1))
        DO j=1, ncell_sets
          ALLOCATE(mesh%cell_sets(j)%cell_list(SIZE(cell_sets_temp(j)%cell_list)))
          mesh%cell_sets(j)%cell_list = cell_sets_temp(j)%cell_list
          mesh%cell_sets(j)%name = cell_sets_temp(j)%name
          DEALLOCATE(cell_sets_temp(j)%cell_list)
          CALL cell_sets_temp(j)%name%clear()
        ENDDO
        DEALLOCATE(cell_sets_temp)
      ELSE
        ncell_sets = 0
        ALLOCATE(mesh%cell_sets(1))
      ENDIF
      ! Add the one new cell set
      mesh%cell_sets(ncell_sets + 1)%name = elname
      ALLOCATE(mesh%cell_sets(ncell_sets + 1)%cell_list(ncells))
      IF(dtype == 'SNK') THEN
        ! Account for 0 based to 1 based index switch
        mesh%cell_sets(ncell_sets + 1)%cell_list = ivals4_1d + 1
        DEALLOCATE(ivals4_1d)
      ELSE
        mesh%cell_sets(ncell_sets + 1)%cell_list = ivals8_1d + 1
        DEALLOCATE(ivals8_1d)
      ENDIF

    CASE DEFAULT
      CALL eXDMF%raiseWarning(modName//'::'//myName// &
        ' - Unsupported data in XDMF file '//CHAR(elname))
    ENDSELECT
  ENDDO
ENDSUBROUTINE setup_leaf_XDMFMesh_from_file
!
!-------------------------------------------------------------------------------
!> @brief Imports the mesh data in the file to a mesh object.
!> @param strpath the string holding the path to the XDMF file
!> @param mesh the XDMF mesh object
!>
SUBROUTINE importXDMFMesh(strpath, mesh)
  CHARACTER(LEN=*),PARAMETER :: myName='importXDMFMesh'
  TYPE(StringType),INTENT(INOUT) :: strpath
  TYPE(XDMFMeshType),INTENT(OUT),TARGET  :: mesh
  TYPE(XMLFileType) :: xml
  TYPE(HDF5FileType) :: h5
  TYPE(XMLElementType),POINTER :: xmle, children(:)
  TYPE(StringType) :: strIn, strOut
  INTEGER(SIK) :: i, gridIdx
  CHARACTER(LEN=200) :: charpath

  ! Initialize the XDMFTopologyList if it has not been
  IF(.NOT.XDMFTopologyList%has('Topology')) CALL init_XDMFTopologyList()

  !H5
  ! NOTE: it is assumed that the h5 and xml files have the same name.
  i = LEN_TRIM(strpath)
  charpath = CHAR(strpath)
  CALL h5%init(charpath(1:i-4)//"h5",'READ')
  CALL h5%fopen()

  !XML
  CALL xml%importFromDisk(ADJUSTL(strpath))
  xmle => xml%root
  IF(.NOT.ASSOCIATED(xmle)) CALL eXDMF%raiseError(modName//'::'//myName// &
    ' - XML data import encountered an error. Pointer to root not associated.')
  IF(.NOT.xmle%name%upper() == 'XDMF') CALL eXDMF%raiseError(modName//'::'//&
    myName//' - Expected XDMF XML element to be the root element.')

  ! Version
  strIn='Version'
  CALL xmle%getAttributeValue(strIn,strOut)
  IF(strOut /= '3.0') THEN
    CALL eXDMF%raiseError(modName//'::'//myName// &
      ' - Currently only supports XDMF version 3.0')
  ENDIF

  ! Domain
  CALL xmle%getChildren(children)
  IF(.NOT.(SIZE(children) == 1 .AND. children(1)%name%upper() == 'DOMAIN'))THEN
    CALL eXDMF%raiseError(modName//'::'//myName// &
    ' - Expected XDMF XML element to have one child (Domain).')
  ENDIF

  ! Information
  ! NOTE: It is assumed that material information is before any grids
  ! and that all grids are contained in one overall grid.
  CALL children(1)%getChildren(children)
  IF (SIZE(children) == 2) THEN
    IF(.NOT.children(1)%name%upper() == 'INFORMATION')THEN
      CALL eXDMF%raiseError(modName//'::'//myName// &
        ' - Expected Domain XML element to have one Information before Grid.')
    ENDIF
    IF(.NOT.children(2)%name%upper() == 'GRID')THEN
      CALL eXDMF%raiseError(modName//'::'//myName// &
        ' - Expected Domain XML element to have Grid after Information.')
    ENDIF
    gridIdx = 2
  ELSE IF(SIZE(children) == 1) THEN
    IF(.NOT.children(1)%name%upper() == 'GRID')THEN
      CALL eXDMF%raiseError(modName//'::'//myName// &
        ' - Expected Domain XML element to have Grid child.')
    ENDIF
    gridIdx = 1
  ELSE
    CALL eXDMF%raiseError(modName//'::'//myName// &
      ' - Expecting information and grid elements only.')
  ENDIF

  ! Init root mesh
  strIn="Name"
  CALL children(gridIdx)%getAttributeValue(strIn,strOut)
  mesh%name = strOut

  ! Create grids
  CALL create_XDMFMesh_from_file(mesh, children(gridIdx), h5)

  ! Setup bounding boxes
  CALL mesh%recomputeBoundingBox()

ENDSUBROUTINE importXDMFMesh
!
!-------------------------------------------------------------------------------
!> @brief Clears the XDMF mesh
!> @param thismesh the XDMF mesh object
!>
RECURSIVE SUBROUTINE clear_XDMFMeshType(thismesh)
  CLASS(XDMFMeshType), INTENT(INOUT) :: thismesh
  INTEGER(SNK) :: i

  CALL thismesh%name%clear()
  thismesh%singleTopology = .FALSE.
  thismesh%boundingBox = 0.0_SDK
  IF(ASSOCIATED(thismesh%parent)) thismesh%parent => NULL()
  IF(ASSOCIATED(thismesh%children)) THEN
    DO i=1,SIZE(thismesh%children)
      CALL thismesh%children(i)%clear()
    ENDDO
    thismesh%children => NULL()
  ENDIF
  IF( ALLOCATED(thismesh%vertices) ) DEALLOCATE(thismesh%vertices)
  IF( ALLOCATED(thismesh%cells) ) THEN
    DO i=1, SIZE(thismesh%cells)
      DEALLOCATE(thismesh%cells(i)%vertex_list)
    ENDDO
    DEALLOCATE(thismesh%cells)
  ENDIF
  IF( ALLOCATED(thismesh%material_ids) ) DEALLOCATE(thismesh%material_ids)
  IF( ALLOCATED(thismesh%cell_sets) ) THEN
    DO i=1, SIZE(thismesh%cell_sets)
      CALL thismesh%cell_sets(i)%name%clear()
      DEALLOCATE(thismesh%cell_sets(i)%cell_list)
    ENDDO
    DEALLOCATE(thismesh%cell_sets)
  ENDIF
ENDSUBROUTINE clear_XDMFMeshType
!
!-------------------------------------------------------------------------------
!> @brief Gets the number of grid levels to a leaf node
!> @param thismesh the XDMF mesh object
!> @returns n the number of levels to a leaf node
!>
RECURSIVE FUNCTION distanceToLeaf_XDMFMeshType(thismesh) RESULT(n)
  CLASS(XDMFMeshType), INTENT(INOUT) :: thismesh
  INTEGER(SNK) :: n
  ! It is assumed that all siblings have children if any sibling has children
  ! Hence it is sufficient to assess thismesh%children(1)'s depth recursively
  ! Ex:     Possible                Not Possible
  !           L1                        L1
  !          /  \                      /  \
  !      L2_1   L2_2               L2_1   L2_2 <-- Should have children.
  !     /  \     |                /  \
  ! L3_1  L3_2  L3_3           L3_1  L3_2
  IF(ASSOCIATED(thismesh%children))THEN
    n = thismesh%children(1)%distanceToLeaf() + 1
  ELSE
    n = 0
  ENDIF
ENDFUNCTION distanceToLeaf_XDMFMeshType
!
!-------------------------------------------------------------------------------
!> @brief Recompute the bounding box for this mesh and all children.
!> @param thismesh the XDMF mesh object
!>
RECURSIVE SUBROUTINE recomputeBoundingBox_XDMFMeshType(thismesh)
  CLASS(XDMFMeshType), INTENT(INOUT) :: thismesh
  REAL(SDK) :: xmin, xmax, ymin, ymax
  INTEGER(SLK) :: i
  xmin = HUGE(xmin)
  xmax = -HUGE(xmax)
  ymin = HUGE(ymin)
  ymax = -HUGE(ymax)
  IF(ASSOCIATED(thismesh%children))THEN
    DO i = 1, SIZE(thismesh%children)
      CALL thismesh%children(i)%recomputeBoundingBox()
    ENDDO
    DO i = 1, SIZE(thismesh%children)
      IF(thismesh%children(i)%boundingBox(1) < xmin) &
        xmin = thismesh%children(i)%boundingBox(1)
      IF(thismesh%children(i)%boundingBox(2) > xmax) &
        xmax = thismesh%children(i)%boundingBox(2)
      IF(thismesh%children(i)%boundingBox(3) < ymin) &
        ymin = thismesh%children(i)%boundingBox(3)
      IF(thismesh%children(i)%boundingBox(4) > ymax) &
        ymax = thismesh%children(i)%boundingBox(4)
    ENDDO
  ELSE
    DO i = 1, SIZE(thismesh%vertices, DIM=2)
      IF(thismesh%vertices(1,i) < xmin) xmin = thismesh%vertices(1,i)
      IF(thismesh%vertices(1,i) > xmax) xmax = thismesh%vertices(1,i)
      IF(thismesh%vertices(2,i) < ymin) ymin = thismesh%vertices(2,i)
      IF(thismesh%vertices(2,i) > ymax) ymax = thismesh%vertices(2,i)
    ENDDO
  ENDIF
  thismesh%boundingBox = (/xmin, xmax, ymin, ymax/)
ENDSUBROUTINE recomputeBoundingBox_XDMFMeshType
!
!-------------------------------------------------------------------------------
!> @brief Assigns an XDMF mesh type to another
!> @param thismesh the XDMF mesh object being assigned to
!> @param thatmesh the XDMF mesh object being assigned from
!>
RECURSIVE SUBROUTINE assign_XDMFMeshType(thismesh, thatmesh)
  TYPE(XDMFMeshType), INTENT(INOUT) :: thismesh
  TYPE(XDMFMeshType), INTENT(IN) :: thatmesh
  INTEGER(SNK) :: i

  thismesh%name = thatmesh%name
  thismesh%singleTopology = thatmesh%singleTopology
  thismesh%boundingBox = thatmesh%boundingBox
  IF(ASSOCIATED(thatmesh%parent)) thismesh%parent => thatmesh%parent
  ! NOTE: Children cannot be recursively cleared without risk of
  ! modify other mesh objects due to the pointer to other meshes.
  ! Therefore, it is assumed that one will manually clear a mesh
  ! if the children are to be deleted.
  IF(ASSOCIATED(thatmesh%children)) THEN
    ALLOCATE(thismesh%children(SIZE(thatmesh%children)))
    thismesh%children => thatmesh%children
  ENDIF
  IF( ALLOCATED(thatmesh%vertices) ) THEN
    IF(ALLOCATED(thismesh%vertices)) DEALLOCATE(thismesh%vertices)
    ALLOCATE(thismesh%vertices(3, SIZE(thatmesh%vertices, DIM=2)))
    thismesh%vertices = thatmesh%vertices
  ENDIF
  IF( ALLOCATED(thatmesh%cells) ) THEN
    IF(ALLOCATED(thismesh%cells))THEN
      DO i=1, SIZE(thismesh%cells)
        DEALLOCATE(thismesh%cells(i)%vertex_list)
      ENDDO
      DEALLOCATE(thismesh%cells)
    ENDIF
    ALLOCATE(thismesh%cells(SIZE(thatmesh%cells)))
    DO i = 1, SIZE(thatmesh%cells)
      ALLOCATE(thismesh%cells(i)%vertex_list(SIZE(thatmesh%cells(i)%vertex_list)))
      thismesh%cells(i)%vertex_list = thatmesh%cells(i)%vertex_list
    ENDDO
  ENDIF
  IF( ALLOCATED(thatmesh%material_ids) )THEN
    IF(ALLOCATED(thismesh%material_ids)) DEALLOCATE(thismesh%material_ids)
    ALLOCATE(thismesh%material_ids(SIZE(thatmesh%material_ids)))
    thismesh%material_ids = thatmesh%material_ids
  ENDIF
  IF( ALLOCATED(thatmesh%cell_sets) ) THEN
    IF(ALLOCATED(thismesh%cell_sets))THEN
      DO i=1, SIZE(thismesh%cell_sets)
        DEALLOCATE(thismesh%cell_sets(i)%cell_list)
      ENDDO
      DEALLOCATE(thismesh%cell_sets)
    ENDIF
    ALLOCATE(thismesh%cell_sets(SIZE(thatmesh%cell_sets)))
    DO i = 1, SIZE(thatmesh%cell_sets)
      ALLOCATE(thismesh%cell_sets(i)%cell_list(SIZE(thatmesh%cell_sets(i)%cell_list)))
      thismesh%cell_sets(i)%cell_list = thatmesh%cell_sets(i)%cell_list
      thismesh%cell_sets(i)%name = thatmesh%cell_sets(i)%name
    ENDDO
  ENDIF
ENDSUBROUTINE assign_XDMFMeshType
!
!-------------------------------------------------------------------------------
!> @brief Export the leaf nodes of the mesh hierarchy
!> @param mesh the mesh
!> @param xmle the XML element
!> @param strpath the string holding the path to the XDMF file
!> @param h5 the HDF5 file
!>
RECURSIVE SUBROUTINE export_leaf_XDMF(mesh, xmle, strpath, h5)
  TYPE(XDMFMeshType),INTENT(IN)  :: mesh
  TYPE(XMLElementType),TARGET,INTENT(INOUT) :: xmle
  TYPE(StringType),INTENT(IN) :: strpath
  TYPE(HDF5FileType), INTENT(INOUT) :: h5
  TYPE(XMLElementType),POINTER :: current_xml, child_xml, children(:), &
    children2(:)
  TYPE(StringType) :: str_name, str_value, str1, str2, toponame, xdmf_id_str
  INTEGER(SNK) :: nchildren, ichild, verts_in_cell
  INTEGER(SLK) :: xdmf_id, nverts, ncells, i, j, ivert
  CHARACTER(LEN=200) :: charpath
  INTEGER(SLK), ALLOCATABLE :: vertex_list_2d(:, :), vertex_list_1d(:), cell_list_1d(:)
  INTEGER(SNK),PARAMETER :: GEOMETRY_IDX=1
  INTEGER(SNK),PARAMETER :: TOPOLOGY_IDX=2

  REQUIRE(ALLOCATED(mesh%vertices) .AND. ALLOCATED(mesh%cells))

  ! Create HDF5 group
  CALL h5%mkdir(CHAR(mesh%name))

  ! Determine number of xmle children
  ! Geometry, Topology, Material, Cell sets
  ! Assumes mesh has geometry and topology
  nchildren = 2
  IF(ALLOCATED(mesh%material_ids)) nchildren = nchildren + 1
  IF(ALLOCATED(mesh%cell_sets)) nchildren = nchildren + SIZE(mesh%cell_sets)
  ALLOCATE(children(nchildren))
  CALL xmle%setChildren(children)

  ! GEOMETRY
  current_xml => children(GEOMETRY_IDX)
  str_name="Geometry"
  CALL current_xml%setName(str_name)
  CALL current_xml%setParent(xmle)

  str_name= "GeometryType"
  str_value = "XYZ"
  CALL current_xml%setAttribute(str_name, str_value)

  children => NULL()
  ALLOCATE(children(1))
  CALL current_xml%setChildren(children)
  child_xml => children(1)
  str_name="DataItem"
  CALL child_xml%setName(str_name)
  CALL child_xml%setParent(current_xml)

  str_name= "DataType"
  str_value = "Float"
  CALL child_xml%setAttribute(str_name, str_value)

  nverts=SIZE(mesh%vertices, DIM=2)
  str_name="Dimensions"
  str1 = nverts
  str2 = "3"
  str_value = str1//" "//str2
  CALL child_xml%setAttribute(str_name, str_value)

  str_name= "Format"
  str_value = "HDF"
  CALL child_xml%setAttribute(str_name, str_value)

  str_name= "Precision"
  str_value = "8"
  CALL child_xml%setAttribute(str_name, str_value)

  i = LEN_TRIM(strpath)
  charpath = CHAR(strpath)
  child_xml%content = charpath(1:i-4)//"h5:/"//mesh%name//"/vertices"

  CALL h5%fwrite(CHAR(mesh%name)//'->vertices',mesh%vertices)

  children => NULL()

  ! TOPOLOGY
  CALL xmle%getChildren(children)
  current_xml => children(TOPOLOGY_IDX)
  str_name="Topology"
  CALL current_xml%setName(str_name)
  CALL current_xml%setParent(xmle)

  ! Single topology
  IF(mesh%singleTopology)THEN
    str_name= "TopologyType"
    xdmf_id = mesh%cells(1)%vertex_list(1)
    xdmf_id_str = xdmf_id
    CALL XDMFTopologyList%get('XDMFID->'//ADJUSTL(xdmf_id_str), toponame)
    CALL XDMFTopologyList%get(ADJUSTL(toponame)//'->n', verts_in_cell)
    str_value = toponame
    CALL current_xml%setAttribute(str_name, str_value)

    str_name= "NumberOfElements"
    ncells = SIZE(mesh%cells)
    str_value = ncells
    CALL current_xml%setAttribute(str_name, str_value)

    str_name= "NodesPerElement"
    str_value = verts_in_cell
    CALL current_xml%setAttribute(str_name, str_value)

    children => NULL()
    ALLOCATE(children(1))

    CALL current_xml%setChildren(children)
    child_xml => children(1)
    str_name="DataItem"
    CALL child_xml%setName(str_name)
    CALL child_xml%setParent(current_xml)

    str_name= "DataType"
    str_value = "Int"
    CALL child_xml%setAttribute(str_name, str_value)

    str_name="Dimensions"
    str1 = ncells
    str2 = verts_in_cell
    str_value = str1//" "//str2
    CALL child_xml%setAttribute(str_name, str_value)

    str_name= "Format"
    str_value = "HDF"
    CALL child_xml%setAttribute(str_name, str_value)

    str_name= "Precision"
    str_value = "8"
    CALL child_xml%setAttribute(str_name, str_value)

    i = LEN_TRIM(strpath)
    charpath = CHAR(strpath)
    child_xml%content = charpath(1:i-4)//"h5:/"//mesh%name//"/cells"

    ALLOCATE(vertex_list_2d(verts_in_cell, ncells))
    DO i = 1, ncells
      ! Convert 1 based to 0 based index
      vertex_list_2d(:,i) = mesh%cells(i)%vertex_list(2:) - 1
    ENDDO
    CALL h5%fwrite(CHAR(mesh%name)//'->cells',vertex_list_2d)
    DEALLOCATE(vertex_list_2d)
  ! Mixed topology
  ELSE
    str_name= "TopologyType"
    str_value = "Mixed"
    CALL current_xml%setAttribute(str_name, str_value)

    str_name= "NumberOfElements"
    ncells = SIZE(mesh%cells)
    str_value = ncells
    CALL current_xml%setAttribute(str_name, str_value)

    children => NULL()
    ALLOCATE(children(1))
    CALL current_xml%setChildren(children)
    child_xml => children(1)
    str_name="DataItem"
    CALL child_xml%setName(str_name)
    CALL child_xml%setParent(current_xml)

    str_name= "DataType"
    str_value = "Int"
    CALL child_xml%setAttribute(str_name, str_value)

    nverts = 0
    DO i = 1, ncells
      nverts = nverts + SIZE(mesh%cells(i)%vertex_list)
    ENDDO
    str_name="Dimensions"
    str_value = nverts
    CALL child_xml%setAttribute(str_name, str_value)

    str_name= "Format"
    str_value = "HDF"
    CALL child_xml%setAttribute(str_name, str_value)

    str_name= "Precision"
    str_value = "8"
    CALL child_xml%setAttribute(str_name, str_value)

    i = LEN_TRIM(strpath)
    charpath = CHAR(strpath)
    child_xml%content = charpath(1:i-4)//"h5:/"//mesh%name//"/cells"

    ALLOCATE(vertex_list_1d(nverts))
    ivert = 1
    DO i = 1, ncells
      nverts = SIZE(mesh%cells(i)%vertex_list)
      ! Convert 1 based to 0 based index
      vertex_list_1d(ivert) = mesh%cells(i)%vertex_list(1)
      vertex_list_1d(ivert + 1 : ivert + nverts - 1) = mesh%cells(i)%vertex_list(2:) - 1
      ivert = ivert + nverts
    ENDDO
    CALL h5%fwrite(CHAR(mesh%name)//'->cells',vertex_list_1d)
    DEALLOCATE(vertex_list_1d)
  ENDIF

  ichild = 3

  ! MATERIAL ID
  IF(ALLOCATED(mesh%material_ids))THEN
    children => NULL()
    CALL xmle%getChildren(children)
    current_xml => children(ichild)
    str_name="Attribute"
    CALL current_xml%setName(str_name)
    CALL current_xml%setParent(xmle)

    str_name= "Center"
    str_value = "Cell"
    CALL current_xml%setAttribute(str_name, str_value)

    str_name= "Name"
    str_value = "MaterialID"
    CALL current_xml%setAttribute(str_name, str_value)

    children => NULL()
    ALLOCATE(children(1))
    CALL current_xml%setChildren(children)
    child_xml => children(1)

    str_name="DataItem"
    CALL child_xml%setName(str_name)
    CALL child_xml%setParent(current_xml)

    str_name= "DataType"
    str_value = "Int"
    CALL child_xml%setAttribute(str_name, str_value)

    str_name="Dimensions"
    ncells = SIZE(mesh%cells)
    str_value = ncells
    CALL child_xml%setAttribute(str_name, str_value)

    str_name= "Format"
    str_value = "HDF"
    CALL child_xml%setAttribute(str_name, str_value)

    str_name= "Precision"
    str_value = "4"
    CALL child_xml%setAttribute(str_name, str_value)

    i = LEN_TRIM(strpath)
    charpath = CHAR(strpath)
    child_xml%content = charpath(1:i-4)//"h5:/"//mesh%name//"/material_id"

    ALLOCATE(vertex_list_1d(ncells))
    ! Convert 1 based to 0 based index
    vertex_list_1d = mesh%material_ids - 1
    CALL h5%fwrite(CHAR(mesh%name)//'->material_id',vertex_list_1d)
    DEALLOCATE(vertex_list_1d)
    ichild = ichild + 1
  ENDIF

  ! CELL SETS
  IF(ALLOCATED(mesh%cell_sets))THEN
    children => NULL()
    CALL xmle%getChildren(children)
    DO i=ichild, nchildren
      current_xml => children(i)
      str_name="Set"
      CALL current_xml%setName(str_name)
      CALL current_xml%setParent(xmle)

      str_name= "Name"
      str_value = mesh%cell_sets(i - ichild + 1)%name
      CALL current_xml%setAttribute(str_name, str_value)

      str_name= "SetType"
      str_value = "Cell"
      CALL current_xml%setAttribute(str_name, str_value)

      ALLOCATE(children2(1))
      CALL current_xml%setChildren(children2)
      child_xml => children2(1)
      str_name="DataItem"
      CALL child_xml%setName(str_name)
      CALL child_xml%setParent(current_xml)

      str_name= "DataType"
      str_value = "Int"
      CALL child_xml%setAttribute(str_name, str_value)

      str_name="Dimensions"
      ncells = SIZE(mesh%cell_sets(i-ichild+1)%cell_list)
      str_value = ncells
      CALL child_xml%setAttribute(str_name, str_value)

      str_name= "Format"
      str_value = "HDF"
      CALL child_xml%setAttribute(str_name, str_value)

      str_name= "Precision"
      str_value = "8"
      CALL child_xml%setAttribute(str_name, str_value)

      j = LEN_TRIM(strpath)
      charpath = CHAR(strpath)
      child_xml%content = charpath(1:j-4)//"h5:/"//mesh%name//"/"//mesh%cell_sets(i-ichild+1)%name

      ALLOCATE(cell_list_1d(ncells))
      ! Convert 1 based to 0 based index
      cell_list_1d = mesh%cell_sets(i-ichild+1)%cell_list - 1
      CALL h5%fwrite(CHAR(mesh%name)//'->'//CHAR(mesh%cell_sets(i-ichild+1)%name),cell_list_1d)
      DEALLOCATE(cell_list_1d)
    ENDDO
  ENDIF
ENDSUBROUTINE export_leaf_XDMF
!
!-------------------------------------------------------------------------------
!> @brief Create the xml hierarchy for the mesh
!> @param mesh the mesh
!> @param xmle the XML element
!> @param strpath the string holding the path to the XDMF file
!> @param h5 the HDF5 file
!>
RECURSIVE SUBROUTINE create_xml_hierarchy_XDMF(mesh, xmle, strpath, h5)
  TYPE(XDMFMeshType),INTENT(INOUT)  :: mesh
  TYPE(XMLElementType),TARGET,INTENT(INOUT) :: xmle
  TYPE(StringType),INTENT(INOUT) :: strpath
  TYPE(HDF5FileType), INTENT(INOUT) :: h5
  TYPE(XMLElementType), POINTER :: children(:)
  INTEGER(SNK) :: i
  TYPE(StringType) :: str_name, str_value

  ! If this mesh has children
  IF(ASSOCIATED(mesh%children)) THEN
    ! Add XML element children
    ALLOCATE(children(SIZE(mesh%children)))
    CALL xmle%setChildren(children)
    DO i=1,SIZE(mesh%children)
      ! Set attributes then recurse
      str_name="Grid"
      CALL children(i)%setName(str_name)
      CALL children(i)%setParent(xmle)
      str_name='Name'
      str_value = mesh%children(i)%name
      CALL children(i)%setAttribute(str_name, str_value)
      str_name='GridType'
      IF(ASSOCIATED(mesh%children(i)%children))THEN
        str_value = 'Tree'
      ELSE
        str_value = 'Uniform'
      ENDIF
      CALL children(i)%setAttribute(str_name, str_value)

      CALL create_xml_hierarchy_XDMF(mesh%children(i), children(i), strpath, h5)
    ENDDO
  ELSE
    CALL export_leaf_XDMF(mesh, xmle, strpath, h5)
  ENDIF
ENDSUBROUTINE create_xml_hierarchy_XDMF
!
!-------------------------------------------------------------------------------
!> @brief Exports mesh data to an XDMF file.
!> @param strpath the string holding the path to the XDMF file
!> @param mesh the XDMF mesh object
!>
SUBROUTINE exportXDMFMesh(strpath, mesh)
  CHARACTER(LEN=*),PARAMETER :: myName='exportXDMFMesh'
  TYPE(StringType),INTENT(INOUT) :: strpath
  TYPE(XDMFMeshType),INTENT(INOUT)  :: mesh
  TYPE(XMLFileType) :: xml
  TYPE(HDF5FileType) :: h5
  TYPE(XMLElementType),POINTER :: xmle, children(:), children2(:)
  TYPE(StringType) :: str_name, str_value
  INTEGER(SNK) :: i
  CHARACTER(LEN=200) :: charpath

  ! Create HDF5 file
  i = LEN_TRIM(strpath)
  charpath = CHAR(strpath)
  CALL h5%init(charpath(1:i-4)//"h5",'NEW')
  CALL h5%fopen()

  ! Create XML file
  CALL xml%init(ADJUSTL(strpath),.FALSE.)
  xmle => xml%root
  IF(.NOT.ASSOCIATED(xmle)) CALL eXDMF%raiseError(modName//'::'//myName// &
    ' - XML data init encountered an error. Pointer to root not associated.')
  !   Set Xdmf
  str_name='Xdmf'
  CALL xmle%setName(str_name)
  str_name='Version'
  str_value = '3.0'
  CALL xmle%setAttribute(str_name, str_value)
  !   Set Domain
  ALLOCATE(children(1))
  CALL xmle%setChildren(children)
  str_name='Domain'
  CALL children(1)%setName(str_name)
  CALL children(1)%setParent(xml%root)
  ! Setup the grid that contains everything
  xmle => children(1)

  ALLOCATE(children2(1))
  CALL xmle%setChildren(children2)
  str_name="Grid"
  CALL children2(1)%setName(str_name)
  CALL children2(1)%setParent(xmle)
  str_name='Name'
  str_value = mesh%name
  CALL children2(1)%setAttribute(str_name, str_value)
  str_name='GridType'
  IF(.NOT.ASSOCIATED(mesh%children))THEN
    str_value = 'Uniform'
  ELSE
    str_value = 'Tree'
  ENDIF
  CALL children2(1)%setAttribute(str_name, str_value)

  ! Recursively add xml elements for each grid. Only the leaves have vertices,
  ! so only the leaves have HDF5 groups/data.
  xmle => children2(1)
  CALL create_xml_hierarchy_XDMF(mesh, xmle, strpath, h5)

  ! Finish up
  CALL xml%exportToDisk(CHAR(strpath))
  CALL h5%fclose()

ENDSUBROUTINE exportXDMFMesh
#endif
ENDMODULE XDMFMesh

max_stars_repo_name: picmc/Futility

max_issues_repo_name: picmc/Futility

max_forks_repo_name: picmc/Futility


text:
Describe Users/JaneBF here.

Longtime Davis resident, kids in town, work local. My main reason for a page is my Top Tens

Topten favorite foods in Davis
   1. Crab cakes at Mustard Seed for lunch.
   2. Fish shabushabu in spicy broth at Red Orchid
   3. Greek crepe at Crepeville
   4. Curried fish at Kathmandu Kitchen Kathmandu
   5. Sadly, the Burger King BK fish
   6. Feta prawns at Symposium
   7. Red snapper at Sophias
   8. Salmon burrito at Dos Coyotes
   9. Virtually everthing at Tucos
   10. ?

Local hiking and dog walking
   1. Putah Creek on UC property west of campus
   2. BLM land about (left and right) the Stebbins Cold Canyon Reserve Cold Canyon reserve
   3. Point Isabel in Berkeley

Places that fuss about dogs for absolutely no reason except they can (you realize an energetic dog gets virtually no exercise if hes on a leash).
   1. Davis cemetery
   2. Davis greenbelts
   3. Slide Hill Park
   4. parts of UC davis property where dogs are supposed to be on leashes
   5. Community park
   6. Vic Fazio wildlife area
   7. UC Davis arboretum
   8.

Great lists!  Theres no reason to be ashamed about liking the BK Big Fish.  If it appeals to your tastes, dont worry about perception.  Enjoy what you enjoy!  Users/JabberWokky Evan JabberWokky Edwards

Welcome to the wiki!  I gave you some links for your top tens.  It seems that you and I share interests in eating out and in offleash activities with our dogs!   Users/CovertProfessor

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      subroutine addmesh_amr(cmoc,cmob,ierror)
C
C ######################################################################
C
C     PURPOSE -
C
C        This routine adds two meshes together to preserve the AMR
C           data structures.
C
C     INPUT ARGUMENTS -
C
C        cmoc    - The course mesh_object (source1).
C        cmob    - The fine mesh_object (source2).
C
C     OUTPUT ARGUMENTS -
C
C        cmoc   - The resulting mesh_object.
C        ierror - Error flag.
C
C     CHANGE HISTORY -
C
C    $Log: addmesh_amr.f,v $
C    Revision 2.00  2007/11/05 19:45:45  spchu
C    Import to CVS
C
CPVCS
CPVCS       Rev 1.0   21 Mar 2002 10:07:38   dcg
CPVCS    Initial revision.
CPVCS
CPVCS       Rev 1.89   Wed Apr 05 13:35:20 2000   nnc
CPVCS    Minor source modifications required by the Absoft compiler.
CPVCS
CPVCS       Rev 1.88   Wed Sep 01 13:36:34 1999   dcg
CPVCS    define asmallnumber
CPVCS
CPVCS       Rev 1.87   Tue Sep 01 08:37:06 1998   dcg
CPVCS    remove unused subroutines
C
C ######################################################################
C
      implicit none
C
      include 'local_element.h'
C
      character*132 logmess
C
C ######################################################################
C
      character*(*) cmoc, cmob
      integer ierror
C
C ######################################################################
C
      pointer (ipitetpar, itetpar)
      pointer (ipitetkid, itetkid)
      pointer (ipitetlev, itetlev)
      integer itetpar(1000000),
     *        itetkid(1000000),
     *        itetlev(1000000)
C
C ######################################################################
C
      integer npointsa, numteta
      integer npointsb, numtetb
      integer npointsc, numtetc
C
      pointer (ipitettyp, itettyp)
      pointer (ipitetoff, itetoff)
      pointer (ipjtetoff, jtetoff)
      integer itettyp(1000000), itetoff(1000000), jtetoff(1000000)
      pointer (ipitet, itet1)
      pointer (ipjtet, jtet1)
      integer itet1(1000000), jtet1(1000000)
C
      pointer (ipxic, xic)
      pointer (ipyic, yic)
      pointer (ipzic, zic)
      real*8 xic(1000000), yic(1000000), zic(1000000)
C
C ######################################################################
C
      pointer (ipielmb, ielmb)
      integer ielmb(1000000)
      pointer (ipxelmb, xelmb)
      pointer (ipyelmb, yelmb)
      pointer (ipzelmb, zelmb)
      real*8 xelmb(1000000), yelmb(1000000), zelmb(1000000)
C
      real*8 xavg, yavg, zavg
C
C ######################################################################
C
      integer icscode, lenout, length, icmotp, i, i1, it, itb, itc
      integer icharlnf
C
      character*32 isubname, cmoattnam
      character*8192 cbuff
C
C ######################################################################
C
      isubname='addmesh_amr'
C
      call cmo_select(cmoc,icscode)
      if(icscode.ne.0) call x3d_error(isubname,'cmo_select')
C
      call cmo_get_info('nnodes',cmoc,npointsa,length,icmotp,icscode)
      call cmo_get_info('nelements',cmoc,numteta,length,icmotp,icscode)
C
      cmoattnam='itetpar'
      call mmfindbk(cmoattnam,cmoc,ipitetpar,lenout,icscode)
      if(icscode.ne.0) then
         cbuff='cmo/addatt/' //
     *         cmoc(1:icharlnf(cmoc)) //
     *         '/' //
     *         cmoattnam(1:icharlnf(cmoattnam)) //
     *         '/VINT' //
     *         '/scalar/nelements/linear/permanent/gxa/0.0' //
     *         ' ; finish '
         call dotaskx3d(cbuff,icscode)
         call mmfindbk(cmoattnam,cmoc,ipitetpar,lenout,icscode)
         do i1=1,numteta
            itetpar(i1)=0
         enddo
      endif
      cmoattnam='itetkid'
      call mmfindbk(cmoattnam,cmoc,ipitetkid,lenout,icscode)
      if(icscode.ne.0) then
         cbuff='cmo/addatt/' //
     *         cmoc(1:icharlnf(cmoc)) //
     *         '/' //
     *         cmoattnam(1:icharlnf(cmoattnam)) //
     *         '/VINT' //
     *         '/scalar/nelements/linear/permanent/gxa/0.0' //
     *         ' ; finish '
         call dotaskx3d(cbuff,icscode)
         call mmfindbk(cmoattnam,cmoc,ipitetkid,lenout,icscode)
         do i1=1,numteta
            itetkid(i1)=0
         enddo
      endif
      cmoattnam='itetlev'
      call mmfindbk(cmoattnam,cmoc,ipitetlev,lenout,icscode)
      if(icscode.ne.0) then
         cbuff='cmo/addatt/' //
     *         cmoc(1:icharlnf(cmoc)) //
     *         '/' //
     *         cmoattnam(1:icharlnf(cmoattnam)) //
     *         '/VINT' //
     *         '/scalar/nelements/linear/permanent/gxa/0.0' //
     *         ' ; finish '
         call dotaskx3d(cbuff,icscode)
         call mmfindbk(cmoattnam,cmoc,ipitetlev,lenout,icscode)
         do i1=1,numteta
            itetlev(i1)=0
         enddo
      endif
C
      call cmo_get_info('nnodes',cmob,npointsb,length,icmotp,icscode)
      call cmo_get_info('nelements',cmob,numtetb,length,icmotp,icscode)
      call cmo_get_info('xic',cmob,ipxic,length,icmotp,icscode)
      call cmo_get_info('yic',cmob,ipyic,length,icmotp,icscode)
      call cmo_get_info('zic',cmob,ipzic,length,icmotp,icscode)
      call cmo_get_info('itettyp',cmob,ipitettyp,length,icmotp,icscode)
      call cmo_get_info('itetoff',cmob,ipitetoff,length,icmotp,icscode)
      call cmo_get_info('jtetoff',cmob,ipjtetoff,length,icmotp,icscode)
      call cmo_get_info('itet',cmob,ipitet,length,icmotp,icscode)
      call cmo_get_info('jtet',cmob,ipjtet,length,icmotp,icscode)
C
      length=numtetb
      call mmgetblk('ielmb',isubname,ipielmb,length,2,icscode)
      call mmgetblk('xelmb',isubname,ipxelmb,length,2,icscode)
      call mmgetblk('yelmb',isubname,ipyelmb,length,2,icscode)
      call mmgetblk('zelmb',isubname,ipzelmb,length,2,icscode)
C
      do it=1,numtetb
         xavg=0.0
         yavg=0.0
         zavg=0.0
         do i=1,nelmnen(itettyp(it))
            i1=itet1(itetoff(it)+i)
            xavg=xavg+xic(i1)
            yavg=yavg+yic(i1)
            zavg=zavg+zic(i1)
         enddo
         xelmb(it)=xavg/nelmnen(itettyp(it))
         yelmb(it)=yavg/nelmnen(itettyp(it))
         zelmb(it)=zavg/nelmnen(itettyp(it))
      enddo
C
      call table_element(cmoc,
     &                   ipxelmb,ipyelmb,ipzelmb,numtetb,
     &                   ipielmb,
     &                   icscode)
C
      call addmesh_merge(cmoc,cmob,icscode)
C
      call cmo_get_info('nnodes',cmoc,npointsc,length,icmotp,icscode)
      call cmo_get_info('nelements',cmoc,numtetc,length,icmotp,icscode)
C
      call cmo_get_info('itetpar',cmoc,ipitetpar,length,icmotp,icscode)
      call cmo_get_info('itetkid',cmoc,ipitetkid,length,icmotp,icscode)
      call cmo_get_info('itetlev',cmoc,ipitetlev,length,icmotp,icscode)
C
      do it=1,numtetb
         itb=numteta+it
         itc=ielmb(it)
         if(itc.gt.0) then
            itetpar(itb)=itc
            itetkid(itc)=itb
            itetlev(itb)=itetlev(itc)+1
         else
            write(logmess,'(a,i10)') 'Not in any element: ',itb
            call writloga('default',0,logmess,0,icscode)
         endif
      enddo
C
      cbuff='filter/1,0,0 ; geniee ; finish'
      call dotaskx3d(cbuff,icscode)
C
      goto 9999
 9999 continue
      call mmrelprt(isubname,icscode)
      return
      end

max_forks_repo_name: daniellivingston/LaGriT


text:
c     Introduced by Ben Allanach to provide a wrapper to C class
      subroutine dilogwrap(a, b, c, d)

      implicit none

      double precision a, b, c, d
      complex*16 x, l, pCSPEN

      x = dcmplx(a, b)
      l = pCSPEN(x)

      c = dble(l)
      d = dimag(l)

      return
      end

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     SPENCE-FUNKTION KOMPLEX, FREI NACH HOLLIK                     C
C---------------------------------------------------------------------C
C     20.07.83    LAST CHANGED 10.05.89        ANSGAR DENNER        C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      FUNCTION pCSPEN(ZZ)

      COMPLEX*16 pCSPEN,W,SUM,ZZ,Z,U
      DOUBLE PRECISION RZ,AZ,A1
      DOUBLE PRECISION B(9)
C     BEACHTE:                 B(N)=B2N
C     B(1)=1./6.
C     B(2)=-1./30.
C     B(3)=1./42.
C     B(4)=-1./30.
C     B(5)=5./66.
C     B(6)=-691./2730.
C     B(7)=7./6.
C     B(8)=-3617./510.
C     B(9)=43867./798.
C     PI=3.1415926535897932384
C     PI*PI/6.=1.6449..., PI*PI/3=3.28986...

      B(1)=0.1666666666666666666666666667d0
      B(2)=-0.0333333333333333333333333333d0
      B(3)=0.0238095238095238095238095238d0
      B(4)=-0.0333333333333333333333333333d0
      B(5)=0.0757575757575757575757575758d0
      B(6)=-0.2531135531135531135531135531d0
      B(7)=1.1666666666666666666666666667d0
      B(8)=-7.09215686274509804d0
      B(9)=54.97117794486215539d0

c     write(*,*) 'z:',z
      Z =ZZ*DCMPLX(1D0)
      RZ=DREAL(Z)
      AZ=CDABS(Z)
      A1=CDABS(1D0-Z)
c     write(*,*)'z, rz, az, a1:',z,rz,az,a1
C     IF((SNGL(RZ) .EQ. 0.0) .AND. (SNGL(DIMAG(Z)) .EQ. 0.0)) THEN
C     ---> CHANGED  10.5.89
      IF(AZ .LT. 1D-20) THEN
         pCSPEN=-CDLOG(1D0-Z)
c     write(*,*) 'cspen:', cspen
         RETURN
      END IF
      IF((SNGL(RZ) .EQ. 1.0) .AND. (SNGL(DIMAG(Z)) .EQ. 0.0)) THEN
         pCSPEN=1.64493406684822643D0
c     write(*,*) 'cspen:', cspen
         RETURN
      END IF
      IF(RZ.GT.5D-1) GOTO 20
      IF(AZ.GT.1D0) GOTO 10
      W=-CDLOG(1D0-Z)
      SUM=W-0.25D0*W*W
      U=W
      IF(CDABS(U).LT.1D-10) GOTO 2
c     write(*,*) 'u:',u
c     write(*,*) 'sum:',sum
      DO 1 K=1,9
         U=U*W*W/DFLOAT(2*K*(2*K+1))
         IF(CDABS(U*B(K)/SUM).LT.1D-20) GOTO 2
         SUM=SUM+U*B(K)
 1    CONTINUE
 2    pCSPEN=SUM
c     write(*,*) 'cspen:', cspen
      RETURN
 10   W=-CDLOG(1D0-1D0/Z)
      SUM=W-0.25D0*W*W
      U=W
      IF(CDABS(U).LT.1D-10) GOTO 12

      DO 11 K=1,9
         U=U*W*W/DFLOAT(2*K*(2*K+1))
         IF(CDABS(B(K)*U/SUM).LT.1D-20) GOTO 12
         SUM=SUM+U*B(K)
 11   CONTINUE
 12   pCSPEN=-SUM-1.64493406684822643D0-.5D0*CDLOG(-Z)**2
c     write(*,*) 'cspen:', cspen
      RETURN
 20   IF(A1.GT.1D0) GOTO 30
      W=-CDLOG(Z)
      SUM=W-0.25D0*W*W
      U=W
      IF(CDABS(U).LT.1D-10) GOTO 22
      DO 21 K=1,9
         U=U*W*W/DFLOAT(2*K*(2*K+1))
         IF(CDABS(U*B(K)/SUM).LT.1D-20) GOTO 22
         SUM=SUM+U*B(K)
 21   CONTINUE
 22   pCSPEN=-SUM+1.64493406684822643D0-CDLOG(Z)*CDLOG(1D0-Z)
c     write(*,*) 'cspen:', cspen
      RETURN
 30   W=CDLOG(1D0-1D0/Z)
      SUM=W-0.25D0*W*W
      U=W
      IF(CDABS(U).LT.1D-10) GOTO 32
      DO 31 K=1,9
         U=U*W*W/DFLOAT(2*K*(2*K+1))
         IF(CDABS(U*B(K)/SUM).LT.1D-20) GOTO 32
         SUM=SUM+U*B(K)
 31   CONTINUE
 32   pCSPEN=SUM+3.28986813369645287D0
     *     +.5D0*CDLOG(Z-1D0)**2-CDLOG(Z)*CDLOG(1D0-Z)
c     write(*,*) 'cspen:', cspen
      END


text:
      double complex Lla(4,4,4,4),Lal(4,4,4,4)
      double complex Rla(4,4,4,4),Ral(4,4,4,4)
      common/basic/Lla,Lal,Rla,Ral

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text:
      implicit none
      integer ns,nt,nd
      double precision, allocatable :: source(:,:),targ(:,:)
      double complex, allocatable :: charge(:,:),dipstr(:,:)
      real *8, allocatable :: dipvec(:,:,:)
      double complex, allocatable :: pot(:,:)
      double complex, allocatable :: pottarg(:,:)

      double precision eps
      double complex eye,zk
      integer i,j,k,idim,ier
      double precision hkrand,pi,thet,phi,t1,t2,omp_get_wtime

      data eye/(0.0d0,1.0d0)/

      pi = 4.0d0*atan(1.0d0)

c
cc      initialize printing routine
c
      call prini(6,13)
      write(*,*)
      write(*,*)
      write(*,*) "================================="
      call prin2("This code is an example fortran driver*",i,0)
      call prin2("On output, the code prints sample pot,pottarg*",i,0)
      write(*,*)
      write(*,*)

      zk = 2.2d0

      ns = 1600
      nt = 62500

      nd = 1

      allocate(source(2,ns))
      allocate(targ(2,nt))
      allocate(charge(nd,ns))
      allocate(dipstr(nd,ns))
      allocate(dipvec(nd,2,ns))
      allocate(pot(nd,ns))
      allocate(pottarg(nd,nt))

      eps = 0.51d-6

      write(*,*) "=========================================="

c
c
c       example demonstrating use of
c        source to source+targ, charges+dipoles, pot
c

c
cc      generate sources uniformly on the sphere
c
c
      do i=1,ns
        thet = hkrand(0)*pi
        phi = hkrand(0)*2*pi
        source(1,i) = sin(thet)*cos(phi)
        source(2,i) = sin(thet)*sin(phi)

        do idim=1,nd
          charge(idim,i) = hkrand(0) + eye*hkrand(0)
          dipstr(idim,i) = hkrand(0) + eye*hkrand(0)
          dipvec(idim,1,i) = hkrand(0)
          dipvec(idim,2,i) = hkrand(0)
          pot(idim,i) = 0
        enddo
      enddo

      do i=1,nt
        targ(1,i) = -3.0d0 + hkrand(0)*6.0d0
        targ(2,i) = -3.0d0 + hkrand(0)*6.0d0

        do idim=1,nd
          pottarg(idim,i) = 0
        enddo
      enddo

       call cpu_time(t1)
C$       t1 = omp_get_wtime()

       call hfmm2d_t_cd_p(eps,zk,ns,source,charge,
     1      dipstr,dipvec,nt,targ,pottarg,ier)

       call cpu_time(t2)
C$       t2 = omp_get_wtime()

       call prin2('time taken=*',t2-t1,1)
       call prin2("potential at targets=*",pottarg,12)

      stop
      end
c----------------------------------------------------------
c
cc
c
c


text:
      SUBROUTINE DINTM (NDIMI,ANSWER,WORK,NWORK,IOPT)
c Copyright (c) 1996 California Institute of Technology, Pasadena, CA.
c ALL RIGHTS RESERVED.
c Based on Government Sponsored Research NAS7-03001.
C>> 2009-11-03 DINTM  Krogh  Initialized lots of variables.
C>> 1996-03-31 DINTM  Krogh  Removed unused variable in common.
c>> 1995-11-20 DINTM  Krogh  Converted from SFTRAN to Fortran 77.
C>> 1994-11-14 DINTM  Krogh  Declared all vars.
c>> 1994-10-19 DINTM  Krogh  Changes to use M77CON
c>> 1994-07-07 DINTM  Snyder set up for CHGTYP.
C>> 1993-05-18 DINTM  Krogh -- Changed "END" to "END PROGRAM"
C>> 1994-05-02 DINTM  Snyder corrected some comments
C>> 1991-09-20 DINTM  Krogh converted '(1)' dimensioning to '(*)'.
C>> 1987-11-19 DINTM  Snyder  Initial code.
c
c--D replaces "?": ?INTA,?intc,?intec,?INTF,?INTM,?INTMA,?INTNC,?INTOP
C
C     ******************************************************************
C
C     THIS SUBROUTINE ATTEMPTS TO CALCULATE THE INTEGRAL OF A FUNCTION
C     PROVIDED BY THE USER.  THE FUNCTION IS PROVIDED BY THE USER VIA A
C     SUBROUTINE REFERENCED BY CALL DINTF(F,X,IFLAG) OR BY REVERSE
C     COMMUNICATION.  ALL ABSCISSAE ARE WITHIN THE REGION OF
C     INTEGRATION.
C
C     THE RESULT IS OBTAINED USING QUADRATURE FORMULAE DUE TO
C     T. N. L. PATTERSON, MATHEMATICS OF COMPUTATION, VOLUME 22,
C     PAGES 847-856, 1968.
C
C     *****    WARNING   ***********************************************
C
C     THE RELIABILITY AND EFFICIENCY OF THIS PROGRAM ARE STRONGLY
C     INFLUENCED BY DISCONTINUITIES IN THE FUNCTION OR IT'S
C     DERIVATIVES, INTEGRABLE SINGULARITIES IN THE REGION OF
C     INTEGRATION, AND NON-INTEGRABLE SINGULARITIES NEAR THIS REGION.
C     (INCLUDING COMPLEX VALUES).  THE EFFICIENCY AND RELIABILITY
C     OF INTEGRATING SUCH FUNCTIONS MAY BE GREATLY IMPROVED BY MANUALLY
C     SUBDIVIDING THE REGION OF INTEGRATION AT THE DISCONTINUITY,
C     SINGULARITY OR CLOSEST POINT TO A POLE, AND SUMMING THE ANSWERS.
C     A CHANGE OF VARIABLE TO ELIMINATE OR REDUCE THE STRENGTH OF THE
C     SINGULARITY WILL SIGNIFICANTLY IMPROVE PERFORMANCE.
C
C     *****    FORMAL ARGUMENTS    *************************************
C
C NDIMI   IS THE NUMBER OF DIMENSIONS OF INTEGRATION.
C ANSWER  IS THE ESTIMATE OF THE INTEGRAL.  WHEN REVERSE COMMUNICATION
C         IS SPECIFIED IT IS USED TO PASS FUNCTION VALUES FROM
C         THE USER PROGRAM TO DINTA OR DINTMA.
      DOUBLE PRECISION ANSWER
C WORK    CONTAINS THE LIMITS, WORKING SPACE AND MAY BE REFERENCED BY
C         THE OPTION VECTOR (SEE IOPT BELOW).  WHEN THE INTEGRATION
C         IS COMPLETE, WORK(1) CONTAINS THE ESTIMATED ABSOLUTE ERROR.
C         WHEN REVERSE COMMUNICATION IS SPECIFIED, WORK(1) IS USED TO
C         PASS ABSCISSAE FROM DINTA OR DINTMA TO THE USER PROGRAM.
C         WORK(1), ..., WORK(NDIMI) ARE COMPONENTS OF THE VECTOR OF
C         ABSCISSAE.  WORK(NDIMI+1), ..., WORK(2*NDIMI) ARE LOWER
C         LIMITS.  WORK(2*NDIMI+1), ..., WORK(3*NDIMI) ARE UPPER LIMITS.
C         THE ABSCISSAE AND LIMITS ARE STORED INNERMOST FIRST.  DINTF
C         IS CALLED TO REQUEST THE LIMITS OF EVERY DIMENSION, BUT
C         CONSTANT LIMITS (THE LIMITS OF THE OUTER DIMENSION ARE ALWAYS
C         CONSTANT) MAY INSTEAD BE STORED BEFORE INTEGRATION BEGINS.
C         WORK(3*NDIMI+1), ..., WORK(3*NDIMI+KWORK*(NDIMI-1)) ARE
C         WORKING STORAGE, WHERE KWORK DEPENDS ON THE MACHINE AND THE
C         PRECISION OF THE PROGRAM VERSION.  KWORK MUST SPECIFY ENOUGH
C         STORAGE FOR 4 VARIABLES WHICH ARE ALWAYS DOUBLE PRECISION,
C         139 VARIABLES WHICH ARE DOUBLE PRECISION IN THE DOUBLE
C         PRECISION VERSION, 30 VARIABLES WHICH ARE DOUBLE PRECISION
C         IN THE DOUBLE PRECISION PROGRAM IF THE LARGEST REPRESENTABLE
C         NUMBERS IN SINGLE AND DOUBLE PRECISION ARE DIFFERENT BY MORE
C         THAN A FACTOR OF 4, 29 INTEGER VARIABLES AND 11 LOGICAL
C         VARIABLES.  REMEMBER TO CONSIDER THE DATA TYPE OF WORK WHEN
C         SPECIFYING THE VALUE OF KWORK.  EXAMPLES ARE
C-- Begin mask code changes
C           IF PRECISION IS REAL, AND INTEGER AND LOGICAL VARIABLES
C           OCCUPY THE SAME SPACE AS REAL VARIABLES, KWORK=217.
C           IF PRECISION IS REAL, AND INTEGER AND LOGICAL VARIABLES
C           OCCUPY HALF AS MUCH SPACE AS REAL VARIABLES, KWORK=?
C           IF PRECISION IS PARTLY DOUBLE, AND INTEGER AND LOGICAL
C           VARIABLES OCCUPY THE SAME SPACE AS REAL VARIABLES,
C           KWORK=?
C           IF PRECISION IS PARTLY DOUBLE, AND INTEGER AND LOGICAL
C           VARIABLES OCCUPY HALF AS MUCH SPACE AS REAL VARIABLES,
C           KWORK=?
C           IF PRECISION IS DOUBLE, AND INTEGER AND LOGICAL VARIABLES
C           OCCUPY THE SAME SPACE AS REAL VARIABLES, KWORK=193.
C           IF PRECISION IS DOUBLE, AND INTEGER AND LOGICAL VARIABLES
C           OCCUPY HALF AS MUCH SPACE AS REAL VARIABLES, KWORK=?
C           FOR PORTABILITY, ASSUME THE WORST CASE.  THAT IS, IF THE
C           PRECISION IS SINGLE, KWORK=217.  IF THE PRECISION IS DOUBLE,
C           KWORK=193.
C         WORK(3*NDIMI+KWORK*(NDIMI-1)+1), ... ARE AVAILABLE FOR
C         REFERENCE THROUGH THE OPTION VECTOR.
C-- End mask code changes
      DOUBLE PRECISION WORK(*)
C NWORK   IS THE NUMBER OF ELEMENTS ALLOCATED FOR WORK.
C         NWORK MUST BE GREATER THAN 3*NDIMI+KWORK*(NDIMI-1).
C IOPT    IS A VECTOR OF INTEGERS USED TO RETURN THE STATUS AND TO
C         SPECIFY OPTIONS.
      INTEGER IOPT(*)
C         IOPT(1) IS USED TO RETURN A STATUS INDICATOR TO THE USER.
C         VALUES OF THIS INDICATOR ARE
C         -NDIM   - NORMAL TERMINATION, EITHER THE ABSOLUTE OR THE
C                   RELATIVE ERROR CRITERIA ARE SATISFIED.
C         -NDIM-1 - NORMAL TERMINATION, NEITHER THE ABSOLUTE NOR THE
C                   RELATIVE ERROR CRITERIA ARE SATISFIED, BUT THE
C                   ERROR RELATIVE TO THE OBTAINABLE PRECISION
C                   PRECISION CRITERION IS SATISFIED.
C         -NDIM-2 - NORMAL TERMINATION, BUT NONE OF THE ERROR CRITERIA
C                   ARE SATISFIED.
C         -NDIM-3 - NOT ENOUGH SPACE IN WORK, NWORK IS TOO SMALL.
C         -NDIM-4 - BAD IOPT VALUE.
C         -NDIM-5 - TOO MANY FUNCTION VALUES NEEDED.
C         -NDIM-KDIM-5 - APPARENT NON-INTEGRABLE SINGULARITY IN
C                   DIMENSION KDIM.  ANSWER CONTAINS THE APPROXIMATE
C                   ABSCISSA OF THE SINGULARITY IN THE KDIM-TH
C                   DIMENSION, WORK(KDIM+1) THROUGH WORK(NDIM)
C                   CONTAIN THE ABSCISSAE OF EXTERIOR INTEGRALS.
C         -2*NDIM-KDIM-5 - INCORRECT INNER INTEGRAL DIMENSIONALITY
C                   SPECIFIED DURING EVALUATION OF INTEGRAL AT
C                   DIMENSION KDIM, THAT IS, IABS(IFLAG(1)) = KDIM.
C         ENTRIES IN IOPT STARTING WITH IOPT(2) ARE DESCRIBED BELOW.
C         IOPT(I)   ENTRY IN IOPT(I) MEANS
C           0       NO MORE OPTIONS, BEGIN INTEGRATION.
C           1       IOPT(I+1) CONTAINS THE UNIT NUMBER FOR OUTPUT.
C           2       IOPT(I+1) IS AN NDIMI DIGIT INTEGER, WHERE
C                   EACH DIGIT IS THE DIAGNOSTIC PRINT LEVEL FOR
C                   ONE DIMENSION.  THE LOW ORDER DIGIT IS THE PRINT
C                   LEVEL FOR THE INNER DIMENSION.
C                   0 - NO PRINTING
C                   1 - MINIMUM PRINTING - ERROR MESSAGES (DEFAULT)
C                   2 - PANEL BOUNDARIES AND ANSWERS
C                   3 - ERROR ESTIMATES FOR EACH QUADRATURE FORMULA
C                   4 - DETAILED OUTPUT (DIFFERENCE LINES, ETC).
C           3       WORK(IOPT(I+1)) CONTAINS THE ABSOLUTE TOLERANCE,
C                   WORK(IOPT(I+1)+1) CONTAINS TOLERANCE RELATIVE TO
C                   THE LOCALLY OBTAINABLE PRECISION, AND
C                   WORK(IOPT(I+1)+2) CONTAINS TOLERANCE RELATIVE TO
C                   THE VALUE OF THE INTEGRAL.  THE TOLERANCE RELATIVE
C                   TO THE LOCALLY OBTAINABLE PRECISION IS SPECIFIED AS
C                   THE FRACTION OF LOCALLY OBTAINABLE DIGITS THAT ARE
C                   PERMITTED TO BE WRONG, AND IS INTERNALLY BOUNDED
C                   BETWEEN 0.0 AND 1.0.  IF ANY OF THE ERROR CONTROL
C                   CRITERIA ARE SATISFIED, THE ANSWER IS ACCEPTED.  IF
C                   THIS OPTION IS NOT USED, THE ABSOLUTE AND RELATIVE
C                   TOLERANCES ARE SET TO ZERO, AND THE TOLERANCE
C                   RELATIVE TO THE LOCALLY OBTAINABLE PRECISION IS SET
C                   TO 0.25.
C           4       WORK(IOPT(I+1)) CONTAINS ABSOLUTE ERROR COMMITTED
C                   COMPUTING F.  WORK(IOPT(I+1)) MAY BE CHANGED DURING
C                   THE INTEGRATION.
C           5       WORK(IOPT(I+1)) CONTAINS RELATIVE ERROR EXPECTED TO
C                   BE COMMITTED COMPUTING F.  CHANGES TO THIS VALUE
C                   DURING INTEGRATION WILL NOT BE DETECTED.
C           6       USE REVERSE COMMUNICATION -
C                   CALL DINTM (NDIMI,ANSWER,WORK,NWORK,IOPT)
C                 1 CALL DINTMA (ANSWER,WORK,IOPT(1))
C                   IF (IOPT(1)) 3,2,4
C                 2 COMPUTE THE INNERMOST INTEGRAND -
C                        ANSWER=F(WORK(K)), K=1,NDIMI
C                   THEN EITHER
C                        GO TO 1
C                   OR (FOR BETTER EFFICIENCY)
C                        CALL DINTA (ANSWER,WORK,IOPT(1))
C                        IF (IOPT(1)) 1,2,5
C                   VALUES OF IOPT(1) PRODUCED BY DINTA AND DINTMA
C                   HAVE DIFFERENT MEANINGS.
C                 3 IF (IOPT(1)+NDIMI.LE.0) GO TO 6
C                   THE INTEGRAL OVER THE (-IOPT(1))TH DIMENSION IS
C                   CONTAINED IN ANSWER.  IF A TRANSFORMATION OF
C                   THE INTEGRAL IS NOT NECESSARY BEFORE IT IS USED
C                   AS AN INTEGRAND FOR THE (1-IOPT(1))TH DIMENSION
C                        GO TO 1
C                   ELSE COMPUTE ANSWER=G(ANSWER,WORK(K)),
C                   K=1-IOPT(1),NDIMI AND G IS THE TRANSFORMATION.
C                   COMPUTE ALSO WORK(1)=PARTIAL DERIVATIVE OF G
C                   WITH RESPECT TO ANSWER THEN
C                        GO TO 1
C                 4 COMPUTE WORK(NDIMI+IOPT(1))=A(WORK(K)),
C                   K=IOPT(1)+1,NDIMI, WORK(2*NDIMI+IOPT(1))=B(WORK(K))
C                   K=IOPT(1)+1,NDIM, WHERE A IS THE LOWER LIMIT OF
C                   INTEGRATION FOR THE (IOPT(1))TH DIMENSION, AND B IS
C                   THE UPPER LIMIT.  IF A TRANSFORMATION AS DISCUSSED
C                   ABOVE IS TO BE APPLIED WHEN THE INTEGRATION OF THIS
C                   DIMENSION IS COMPLETE, COMPUTE WORK(1)=AN ESTIMATE
C                   OF THE UPPER BOUND OF THE PARTIAL DERIVATIVE OF G
C                   WITH RESPECT TO THE INTEGRAL.  IF G IS LINEAR IN
C                   THE INTEGRAL, THIS DERIVATIVE CAN BE COMPUTED
C                   WITHOUT KNOWING THE INTEGRAL, RATHER THAN
C                   ESTIMATED.  IF THERE IS A SINGULARITY IN THE
C                   INTEGRAND FOR THIS DIMENSION (PERHAPS INTRODUCED
C                   BY THE BOUNDARY OR THE TRANSFORMATION G), SET
C                   IOPT(1) TO THE LOCATION IN WORK OF THE ABSCISSA IN
C                   THIS DIMENSION OF THE SINGULARITY, AND WITH THE
C                   SIGN SET AS DESCRIBED FOR OPTION 11 BELOW.  THE
C                   MAGNITUDE OF IOPT(1) MUST BE GREATER THAN NDIM FOR
C                   THIS REQUEST TO BE DETECTED.  IN PARTICULAR, IF
C                   THERE IS A SINGULARITY AT ONE OF THE LIMITS, SET
C                   IOPT(1) TO IOPT(1)+NDIM OR -(IOPT(1)+2*NDIM).  THEN
C                        GO TO 1.
C                 5 IOPT(1)=-(IOPT(1)+NDIM)
C                 6 CONTINUE
C                   AT THIS POINT, IERR CONTAINS THE STATUS DESCRIBED
C                   BELOW, FOR OPTION 8.
C           7       SET MINIMUM INDEX OF QUADRATURE FORMULA TO USE
C                   BEFORE SUBDIVISION TO IOPT(I+1).
C           8       NO EFFECT.  IOPT(I+1) MAY BE USED TO PASS
C                   INFORMATION TO DINTF.  INCREMENT I BY 2.
C           9       IOPT(I+1) IS THE MAXIMUM NUMBER OF FUNCTION VALUES
C                   TO USE TO INTEGRATE THE ENTIRE PROBLEM.  IF
C                   IOPT(I+1) .LE. 0, THE NUMBER OF FUNCTION VALUES
C                   IS NOT CONTROLLED.
C          10       IOPT(I+1) IS USED TO RETURN THE NUMBER OF FUNCTION
C                   VALUES USED TO INTEGRATE THE ENTIRE PROBLEM.
C          11       WORK(IABS(IOPT(I+1))) IS THE LOCATION OF A
C                   SINGULARITY OR DISCONTINUITY IN THE OUTERMOST
C                   INTEGRAL.  IF THE LOCATION IS INSIDE THE INTERVAL,
C                   THE INTERVAL IS IMMEDIATELY SUBDIVIDED.  IF
C                   IOPT(I+1) .GT. 0, A T**2 SUBSTITUTION BASED AT
C                   WORK(...) WILL BE USED.  IF IOPT(I+1) .LT. 0,
C                   A T**4 SUBSTITUTION WILL BE USED.  IF IOPT(I+1)
C                   .EQ. 0, NO SUBSTITUTION WILL BE USED.  TO NOTIFY
C                   THE PROGRAM OF SINGULARITIES OR DISCONTINUITIES IN
C                   INTERIOR INTEGRANDS, SET IFLAG NEGATIVE WHEN
C                   ASKED TO COMPUTE LIMITS.  SEE THE DESCRIPTION
C                   OF REVERSE COMMUNICATION ABOVE, OR THE DESCRIPTION
C                   OF DINTF BELOW.
C          12       WORK(IOPT(I+1)) CONTAINS THE ABSOLUTE ERROR IN THE
C                   LOWER LIMIT OF THE CURRENT DIMENSION.  THE ERROR IN
C                   THE UPPER LIMIT IS IN WORK(IOPT(I+1)+1).
C          13       IOPT(I+1) IS THE LOCATION IN IOPT IN WHICH THE USER
C                   IS TO STORE NON-STANDARD CHANGES TO THE DIMENSION OF
C                   INNER INTEGRALS.  THE DEFAULT VALUE IS 1, ALLOWING
C                   SUCH NONSTANDARD CHANGES TO BE STORED IN THE CELL OF
C                   IFLAG (=IOPT) USUALLY USED FOR COMMUNICATION WITH
C                   DINTF, IN THE CASE WHEN IT IS NOT SIMULTANEOUSLY
C                   NECESSARY TO INDICATE A SINGULARITY.  SEE THE
C                   DESCRIPTION OF THE INTERFACE TO DINTF BELOW.
C
C     ALL OPTIONS ARE SET TO THEIR NOMINAL VALUES BEFORE THE OPTION
C     VECTOR IS PROCESSED.
C
C    *****     DINTF     ***********************************************
C
C     DINTF IS REFERENCED VIA CALL DINTF (ANSWER,WORK,IFLAG).
C     WORK AND ANSWER ARE AS DESCRIBED ABOVE.  VALUES OF IFLAG ARE
C     IFLAG .LT. 0
C                    THE INTEGRAL OVER THE (-IFLAG)TH DIMENSION IS
C                    CONTAINED IN ANSWER.  IF A TRANSFORMATION OF
C                    THE INTEGRAL IS NOT NECESSARY BEFORE IT IS USED
C                    AS AN INTEGRAND FOR THE (1-IFLAG)TH DIMENSION
C                         RETURN
C                    ELSE IF THE INTEGRAND DEPENDS ON ONLY ONE INNER
C                    INTEGRAL COMPUTE ANSWER=G(ANSWER,WORK(K)),
C                    ALSO MULTIPLY WORK(1) BY THE PARTIAL DERIVATIVE
C                    OF G WITH RESPECT TO ANSWER THEN
C                         RETURN
C                    ELSE IF THE INTEGRAND DEPENDS ON MORE INTEGRALS
C                    SAVE ANSWER AND WORK(1), AND SET IFLAG(IXKDIM) TO
C                    THE DIMENSION OF INTEGRAL NEXT TO BE COMPUTED
C                    (IXKDIM IS THE VALUE SPECIFIED BY OPTION 12, OR 1
C                    IF OPTION 12 IS NOT SELECTED), THEN
C                         RETURN
C                    ELSE (ALL INNER INTEGRALS UPON WHICH THE INNER
C                    INTEGRAND DEPENDS HAVE BEEN EVALUATED) CALCULATE
C                    ANSWER=G(ANSWER,WORK(K),SAVED INTEGRALS)
C                    K=1-IFLAG, NDIMI (G IS THE TRANSFORMATION),
C                    MULTIPLY WORK(1) BY THE PARTIAL DERIVATIVE OF G
C                    WITH RESPECT TO THE LAST INTEGRAL, AND ADD THE
C                    PRODUCT OF EACH OF THE SAVED VALUES OF WORK(1) AND
C                    PARTIAL DERIVATIVE OF G WITH RESPECT TO THE
C                    CORRESPONDING SAVED INTEGRAL ONTO WORK(1).  THAT
C                    IS, WORK(1) IS THE TOTAL ERROR IN G, CALCULATED AS
C                    THE INNER PRODUCT OF THE ERRORS IN THE INDIVIDUAL
C                    ARGUMENTS, AND THE PARTIAL DERIVATIVES OF G WITH
C                    RESPECT TO THOSE ARGUMENTS.
C     IFLAG .EQ. 0
C                    COMPUTE THE INNERMOST INTEGRAND -
C                         ANSWER=F(WORK(K)), K=1,NDIMI
C                    AND
C                         RETURN
C     IFLAG .GT. 0
C                    COMPUTE WORK(NDIMI+IFLAG)=A(WORK(K)),
C                    K=IFLAG+1,NDIMI, WORK(2*NDIMI+IFLAG)=B(WORK(K)),
C                    K=IFLAG+1,NDIM, WHERE A IS THE LOWER LIMIT OF
C                    INTEGRATION FOR THE (IFLAG)TH DIMENSION, AND B IS
C                    THE UPPER LIMIT.  IF A TRANSFORMATION AS DISCUSSED
C                    ABOVE IS TO BE APPLIED WHEN THE INTEGRATION OF THIS
C                    DIMENSION IS COMPLETE, COMPUTE WORK(1)=AN ESTIMATE
C                    OF THE UPPER BOUND OF THE PARTIAL DERIVATIVE
C                    OF G WITH RESPECT TO THE INTEGRAL.  IF G IS
C                    LINEAR IN THE INTEGRAL, THIS DERIVATIVE MAY BE
C                    COMPUTED WITHOUT KNOWING THE INTEGRAL, RATHER
C                    THAN ESTIMATED.  IF THERE IS A SINGULARITY IN
C                    THE INTEGRAND FOR THIS DIMENSION (PERHAPS
C                    INTRODUCED BY THE BOUNDARY OR THE TRANSFORMATION),
C                    SET IFLAG TO THE LOCATION IN WORK OF THE ABSCISSA
C                    IN THIS DIMENSION OF THE SINGULARITY, AND SET THE
C                    SIGN AS FOR OPTION 11 ABOVE.  NOTE THAT THE
C                    MAGNITUDE OF IFLAG MUST BE GREATER THAN NDIM FOR
C                    THIS REQUEST TO BE DETECTED.  IN PARTICULAR, IF
C                    THERE IS A SINGULARITY AT ONE OF THE LIMITS, SET
C                    THE MAGNITUDE OF IFLAG TO IFLAG+NDIM OR
C                    IFLAG+2*NDIM, WITH THE SIGN AS DESCRIBED FOR
C                    OPTION 11.  THEN
C                         RETURN.
C                    IF THE FIRST INTEGRAL TO BE COMPUTED DOES NOT HAVE
C                    DIMENSION IFLAG(1)-1, SET IFLAG(IXKDIM) TO THE
C                    DIMENSION OF INTEGRAL TO BE NEXT COMPUTED.  SEE
C                    OPTION 12 AND THE DISCUSSION OF IFLAG .LT. 0 ABOVE
C                    FOR INFORMATION ON IXKDIM.
C
C     *****    EXTERNAL REFERENCES     *********************************
C
C DINTMA  TO DO THE INTEGRATION.
C DINTOP  TO SET OPTIONS.
C
C     *****    LOCAL VARIABLES     *************************************
C
      INTEGER NDIMI, NWORK
C
C     *****    COMMON STORAGE     **************************************
C
C     COMMON /DINTNC/ CONTAINS VARIABLES NOT SEPARATELY SAVED FOR
C     EACH DIMENSION OF A MULTIPLE QUADRATURE.  COMMON /DINTC/
C     CONTAINS VARIABLES THAT MUST BE SAVED FOR EACH DIMENSION OF THE
C     QUADRATURE.  THE VARIABLES IN EACH COMMON BLOCK ARE STORED IN THE
C     ORDER - ALWAYS DOUBLE, DOUBLE IF DOUBLE PRECISION PROGRAM, DOUBLE
C     IF DOUBLE PRECISION PROGRAM AND EXPONENT RANGE OF DOUBLE AND
C     SINGLE VERY DIFFERENT, SINGLE, INTEGER, LOGICAL.  A PAD OF LOGICAL
C     VARIABLES IS INCLUDED AT THE END OF /DINTC/.  THE DIMENSION OF
C     THE PAD MAY NEED TO BE VARIED SO THAT NO VARIABLES BEYOND THE END
C     OF THE COMMON BLOCK ARE ALTERED.
C
C     DECLARATIONS OF COMMON /DINTNC/ VARIABLES.
C
      DOUBLE PRECISION AINIT, BINIT, FNCVAL, S, TP
      DOUBLE PRECISION FER, FER1, RELOBT, TPS, XJ, XJP
      INTEGER     FEA,       FEA1,      INC,       INC2,      IPRINT,
     1 ISTOP(2,2),JPRINT,    KDIM,      KK,        KMAXF,     NDIM,
     2 NFINDX,    NFMAX,     NFMAXM,    RELTOL,    REVERM,    REVERS,
     3 WHEREM
      LOGICAL NEEDH
C
C     DECLARATIONS OF COMMON /DINTC/ VARIABLES.
C
c--D Next line special: S => D, X => Q, D => D, P => D
      DOUBLE PRECISION ACUM, PACUM, RESULT(2)
C     139 $.TYPE.$ VARIABLES
      DOUBLE PRECISION
     1 AACUM,     ABSCIS,    DELMIN,    DELTA,     DIFF,      DISCX(2),
     2 END(2),    ERRINA,    ERRINB,    FAT(2),    FSAVE,
     3 FUNCT(24), F1,        F2,        LOCAL(4),  PAACUM,    PF1,
     4 PF2,       PHISUM,    PHTSUM,    PX,        SPACE(6),
     5 STEP(2),   START(2),  SUM,       T,         TA,        TASAVE,
     6 TB,        TEND,      WORRY(2),  X,         X1,
     7 X2,        XT(17),    FT(17),    PHI(34)
c Note XT, FT, and PHI above are last, because they must be in adjacent
c locations in DINTC.
C     30 $DSTYP$ VARIABLES
      DOUBLE PRECISION
     1 ABSDIF,    COUNT,     EDUE2A,    EDUE2B,    EP,        EPNOIZ,
     2 EPS,       EPSMAX,    EPSMIN,    EPSO,      EPSR,      EPSS,
     3 ERR,       ERRAT(2),  ERRC,      ERRF,      ERRI,      ERRT(2),
     4 ESOLD,     EXTRA,     PEPSMN,    RE,        RELEPS,    REP,
     5 REPROD,    RNDC,      TLEN,      XJUMP
C     29 INTEGER VARIABLES
      INTEGER     DISCF,     DISCHK,    ENDPTS,    I,         INEW,
     1 IOLD,      IP,        IXKDIM,    J,         J1,        J1OLD,
     2 J2,        J2OLD,     K,         KAIMT,     KMAX,      KMIN,
     3 L,         LENDT,     NFEVAL,    NFJUMP,    NSUB,      NSUBSV,
     4 NXKDIM,    PART,      SEARCH,    TALOC,     WHERE,     WHERE2
C     11 TO 18 LOGICALS (7 ARE PADDING).
      LOGICAL     DID1,      FAIL,      FATS(2),   FSAVED,    HAVDIF,
     1 IEND,      INIT,      ROUNDF,    XCDOBT(2), PAD(7)
C
C     THE COMMON BLOCKS.
C
      COMMON /DINTNC/
c        1       2       3     4        5       6       7        8
     W AINIT,  BINIT,  FNCVAL, S,      TP,     FER,    FER1,   RELOBT,
c       9      10       11      12      13       1       2        3
     X TPS,    XJ,     XJP,    FEA,    FEA1,   KDIM,    INC,    INC2,
c     4 (2,2)    8       9     10       11      12       13      14
     Y ISTOP,  JPRINT, IPRINT, KK,     KMAXF,  NDIM,   NFINDX, NFMAX,
c        15     16       17      18      19      20
     Z NFMAXM, RELTOL, REVERM, REVERS, WHEREM, NEEDH
      COMMON /DINTC/
     1 ACUM,   PACUM,  RESULT
      COMMON /DINTC/
c        1     2 (4)     6      7        8       9      10     11 (2)
     1 AACUM,  LOCAL,  ABSCIS, TA,     DELTA,  DELMIN, DIFF,   DISCX,
c     13 (2)     15      16    17 (2)   19     20 (24) 44
     2 END,    ERRINA, ERRINB, FAT,    FSAVE,  FUNCT,  F2,
c       45      46     47       48      49     50      51 (6)
     3 PAACUM, PF1,    PF2,    PHISUM, PHTSUM, PX,     SPACE,
c      57 (2)  59 (2)   61     62        63    64       65
     4 STEP,   START,  SUM,    T,      TASAVE, TB,     TEND,
c      66 (2)  68      69      70      71       72
     5 WORRY,  X1,     X2,     X,      F1,     COUNT,
c      73 (17) 90 (17) 107 (34)
     6 XT,     FT,     PHI
      COMMON /DINTC/
c       141     142    143     144      145     146
     1 ABSDIF, EDUE2A, EDUE2B, EP,     EPNOIZ, EPSMAX,
c       147     148     149    150 (2)  152     153
     2 EPSO,   EPSR,   EPSS,   ERRAT,  ERRC,   ERRF,
c     154 (2)   156     157     158     159    160
     3 ERRT,   ESOLD,  EXTRA,  PEPSMN, RELEPS, REP,
c       161     162     163
     4 RNDC,   TLEN,   XJUMP,
c       164    165      166    167    168       169
     5 ERRI,   ERR,    EPSMIN, EPS,    RE,     REPROD
      COMMON /DINTC/
c       170     171     172
     1 DISCF,  DISCHK, ENDPTS, INEW,   IOLD,   IP,     IXKDIM,
     2 J,      J1,     J1OLD,  J2,     J2OLD,  KMAX,   KMIN,
     3 L,      LENDT,  NFJUMP, NSUBSV, NXKDIM, TALOC,  WHERE2,
c      1       2          3      4       5         6      7       8
     4 I,      K,      KAIMT,  NSUB,   PART,   SEARCH, WHERE, NFEVAL
      COMMON /DINTC/
     1 DID1,   FAIL,   FATS,   FSAVED, HAVDIF, IEND,   INIT,   ROUNDF,
     2 XCDOBT, PAD
      SAVE /DINTNC/, /DINTC/
C
C     THE VARIABLES HERE DEFINE THE MACHINE ENVIRONMENT.  ALL ARE SET
C     IN DINTOP.  THE MEANING ATTACHED TO THESE VARIABLES CAN BE
C     FOUND BY LOOKING AT THE DEFINITIONS IN DINTOP.
      DOUBLE PRECISION
     1  EMEPS,  EEPSM8, EDELM2, EDELM3, ESQEPS, ERSQEP, ERSQE6, EMINF,
     2  ESMALL, ENZER,  EDELM1, ENINF
      COMMON /DINTEC/
     1  EMEPS,  EEPSM8, EDELM2, EDELM3, ESQEPS, ERSQEP, ERSQE6, EMINF,
     2  ESMALL, ENZER,  EDELM1, ENINF
      SAVE /DINTEC/

c For initialization of common blocks
      integer KREAL, KINT, KLOG, ITEMP
      parameter (KINT=29, KREAL=169, KLOG=11)
      logical LMOVE(KLOG)
      integer IMOVE(KINT)
      double precision RMOVE(KREAL)
      equivalence (AACUM, RMOVE)
      equivalence (DISCF, IMOVE)
      equivalence (DID1, LMOVE)
C
C
C     *****    PROCEDURES     ******************************************
C

c Initialize common blocks to avoid references to undefined variables.

      do 5 ITEMP = 1, KLOG
        LMOVE(ITEMP) = .false.
 5    continue
      do 10 ITEMP = 1, KREAL
        RMOVE(ITEMP) = 0.D0
 10   continue
      do 15 ITEMP = 1, KINT
        IMOVE(ITEMP) = 0
 15   continue

      WHEREM=0
      NFEVAL=0
      NDIM=NDIMI
      KDIM=NWORK
      X = 0.D0

C     KDIM IS TEMPORARILY USED IN DINTMA TO CHECK THE DIMENSION OF
C     WORK.
      CALL DINTOP (IOPT,WORK)
      IF (IOPT(1).EQ.0) THEN
C
C     CALL DINTMA TO DO THE INTEGRATION.
C
         IF (REVERM.EQ.0) CALL DINTMA (ANSWER,WORK,IOPT)
      ELSE
         IOPT(1)=-NDIM-4
      END IF
      RETURN
      END