Tests.f90

Tests.f90 provides Fortran code examples for many XMDF APIs.

MODULE TestDatasets

USE Xmdf

CHARACTER(LEN=*), PARAMETER :: DATASETS_LOCATION = 'Datasets'
CHARACTER(LEN=*), PARAMETER :: SCALAR_A_LOCATION = 'Scalars/ScalarA'
CHARACTER(LEN=*), PARAMETER :: SCALAR_A_LOCATION_FULL = 'Datasets/Scalars/ScalarA'
CHARACTER(LEN=*), PARAMETER :: SCALAR_B_LOCATION = 'Scalars/ScalarB'
CHARACTER(LEN=*), PARAMETER :: VECTOR2D_A_LOCATION = 'Vectors/Vector2D_A'
CHARACTER(LEN=*), PARAMETER :: VECTOR2D_B_LOCATION = 'Vectors/Vector2D_B'

CONTAINS
! ---------------------------------------------------------------------------
! FUNCTION  tdEditScalarAValues
! PURPOSE   
! NOTES     
! ---------------------------------------------------------------------------
SUBROUTINE TD_EDIT_SCALAR_A_VALUES(a_Filename, a_Compression, error)
  CHARACTER(LEN=*), INTENT(IN) :: a_Filename
  INTEGER, INTENT(IN) :: a_Compression
  INTEGER, INTENT(OUT) :: error
  INTEGER xFileId, xScalarId
  INTEGER, PARAMETER :: editNumValues = 3
  INTEGER  editTimestep
  INTEGER, DIMENSION(editNumValues) :: indices
  REAL, DIMENSION(editNumValues) :: new_values

  CALL TD_WRITE_SCALAR_A(a_Filename, a_Compression, error)
  if (error < 0) then
    return
  endif

    ! open the file and edit the values  
  CALL XF_OPEN_FILE(a_Filename, .FALSE., xFileId, error)
  if (error < 0) then
    return
  endif

  CALL XF_OPEN_GROUP(xFileId, SCALAR_A_LOCATION_FULL, xScalarId, error);
  if (error < 0) then
    CALL XF_CLOSE_FILE(xFileId, error)
    return
  endif
  
    ! Edit values in timestep 1, make index 1 = 4, index 5 = 40,
    ! and index 10 = 400
  editTimestep = 1
  indices(1) = 1;
  indices(2) = 5;
  indices(3) = 10;
  new_values(1) = 4.0;
  new_values(2) = 40.0;
  new_values(3) = 400.0;

  CALL XF_CHANGE_SCALAR_VALUES_TIMESTEP_FLOAT(xScalarId, editTimestep, editNumValues, &
                               indices, new_values, error)
  if (error < 0) then
    CALL XF_CLOSE_GROUP(xScalarId, error)
    CALL XF_CLOSE_FILE(xFileId, error)
    return
  endif

!  RDJ - Technically we shouldn't have to close and reopen a file that we are 
!        editing but it seems to fix a crash that the tests are having.
!  CALL XF_CLOSE_GROUP(xScalarId, error)
!  CALL XF_CLOSE_FILE(xFileId, error)
!    ! open the file and edit the values  
!  CALL XF_OPEN_FILE(a_Filename, .FALSE., xFileId, error)
!  if (error < 0) then
!    return
!  endif
!
!  CALL XF_OPEN_GROUP(xFileId, SCALAR_A_LOCATION_FULL, xScalarId, error);
!  if (error < 0) then
!    CALL XF_CLOSE_FILE(xFileId, error)
!    return
!  endif

    ! Edit values in timestep 2, make index 2 = 6, index 3 = 60, and 
    ! index 9 = 600
  editTimestep = 2
  indices(1) = 2
  indices(2) = 3
  indices(3) = 9
  new_values(1) = -6.0
  new_values(2) = 60.0
  new_values(3) = 6000.0

  CALL XF_CHANGE_SCALAR_VALUES_TIMESTEP_FLOAT(xScalarId, editTimestep, editNumValues, &
                               indices, new_values, error)
  if (error < 0) then
    CALL XF_CLOSE_GROUP(xScalarId, error)
    CALL XF_CLOSE_FILE(xFileId, error)
    return
  endif

  CALL XF_CLOSE_GROUP(xScalarId, error)
  CALL XF_CLOSE_FILE(xFileId, error)

  return
END SUBROUTINE
! tdEditScalarAValues

! ---------------------------------------------------------------------------
! FUNCTION  tdReadScalarAIndices
! PURPOSE   
! NOTES     
! ---------------------------------------------------------------------------
SUBROUTINE TD_READ_SCALAR_A_INDICES (a_Filename, a_nIndices, a_indices, error)
  CHARACTER(LEN=*), INTENT(IN) :: a_Filename
  INTEGER, INTENT(IN) :: a_nIndices
  INTEGER, DIMENSION(*), INTENT(IN) :: a_indices
  INTEGER, INTENT(OUT) :: error
  INTEGER     xFileId, xDsetsId, xScalarAId
  INTEGER     nTimesteps
  REAL, ALLOCATABLE, DIMENSION(:) ::  fValues
  INTEGER     nValues
  INTEGER     id, i, j

  ! open the file
  CALL XF_OPEN_FILE(a_Filename, .TRUE., xFileId, error)
  if (error < 0) then
    return
  endif

  ! open the dataset group
  CALL XF_OPEN_GROUP(xFileId, DATASETS_LOCATION, xDsetsId, error)
  if (error >= 0) then
    CALL XF_OPEN_GROUP(xDsetsId, SCALAR_A_LOCATION, xScalarAId, error)
  endif
  if (error < 0) then
    return
  endif

  ! Find out the number of timesteps in the file
  CALL XF_GET_DATASET_NUM_TIMES(xScalarAId, nTimesteps, error)
  if (error < 0) then
    return
  endif
  if (nTimesteps < 1) then
    error = -1
    return
  endif

  ! Read the values for the index
  nValues = nTimesteps*a_nIndices
  allocate(fValues(nValues))
  CALL XF_READ_SCALAR_VALUES_AT_INDICES_FLOAT(xScalarAId, a_nIndices, a_indices, 1, &
                                            nTimesteps, fValues, error)
  if (error < 0) then
    return
  endif

  ! output the data
  WRITE(*,*) ''
  WRITE(*,*) 'Reading scalar A indices'
  id = 1;
  do i = 1, nTimesteps
    WRITE(*,*) 'Timestep: ', i
    do j = 1, a_nIndices
      WRITE(*,*) 'index:', a_indices(j), ' value: ', fValues(id)
      id = id + 1
    enddo
  enddo
  WRITE(*,*) ''

  deallocate(fValues)

  return

END SUBROUTINE
! TD_READ_SCALARA_INDICES

! --------------------------------------------------------------------------
! FUNCTION tdReadDatasets
! PURPOSE  Read a dataset group from an XMDF file and output information to
!          to a text file
! NOTES    
! --------------------------------------------------------------------------
RECURSIVE SUBROUTINE TD_READ_DATASETS (a_xGroupId, a_FileUnit, error)
INTEGER, INTENT(IN) :: a_xGroupId
INTEGER, INTENT(IN)        :: a_FileUnit
INTEGER, INTENT(OUT)       :: error
INTEGER                   nPaths, nMaxPathLength, j
CHARACTER, ALLOCATABLE, DIMENSION(:) :: Paths
CHARACTER(LEN=500)       IndividualPath
INTEGER                   nStatus, i
INTEGER            xScalarId, xVectorId, xMultiId
INTEGER                   nMultiDatasets

xScalarId = NONE
xVectorId = NONE

nMultiDatasets = 0
nPaths = 0
nMaxPathLength = 0

  ! Look for scalar datasets
call XF_GET_SCALAR_DATASETS_INFO(a_xGroupId, nPaths, nMaxPathLength, nStatus)
if (nStatus >= 0 .AND. nPaths > 0) then
  allocate(Paths(nPaths*nMaxPathLength))
  call XF_GET_SCALAR_DATASET_PATHS(a_xGroupId, nPaths, nMaxPathLength, Paths, &
                                                                         error)
endif
if (nStatus < 0) then
  error = -1
  return
endif

  ! Output number and paths to scalar datasets
WRITE(a_FileUnit,*) 'Number of Scalars ', nPaths
do i=2, nPaths
  IndividualPath = ''
  do j=1, nMaxPathLength-1
    IndividualPath(j:j) = Paths((i-1)*nMaxPathLength+j)
  enddo
  WRITE(a_FileUnit,*) 'Reading scalar: ', IndividualPath(1:nMaxPathLength-1)
  call XF_OPEN_GROUP(a_xGroupId, IndividualPath(1:nMaxPathLength-1), &
                                                           xScalarId, nStatus)
  if (nStatus < 0) then
    error = -1
  return
  endif

  call TDI_READ_SCALAR(xScalarId, a_FileUnit, nStatus)
  call XF_CLOSE_GROUP(xScalarId, error)
  if (nStatus < 0) then
    WRITE(*,*) 'Error reading scalar dataset.'
    error = -1
  return
  endif
enddo

if (allocated(Paths)) deallocate(Paths)
  ! Look for vector datasets
call XF_GET_VECTOR_DATASETS_INFO(a_xGroupId, nPaths, nMaxPathLength, nStatus)
if (nStatus >= 0 .AND. nPaths > 0) then
  allocate(Paths(nPaths*nMaxPathLength))
  call XF_GET_VECTOR_DATASET_PATHS(a_xGroupId, nPaths, nMaxPathLength, Paths, error)
endif
if (nStatus < 0) then
  error = -1
  return
endif

  ! Output number and paths to scalar datasets
WRITE(a_FileUnit,*) 'Number of Vectors ', nPaths
do i=2, nPaths
  do j=1, nMaxPathLength-1
    IndividualPath(j:j) = Paths((i-1)*nMaxPathLength+j)
  enddo
  WRITE(a_FileUnit,*) 'Reading Vector: ', &
                      IndividualPath(1:nMaxPathLength-1)
  call XF_OPEN_GROUP(a_xGroupId, IndividualPath(1:nMaxPathLength-1), &
                                                          xVectorId, nStatus)
  if (nStatus < 0) then
    error = -1
  return
  endif
  call TDI_READ_VECTOR(xVectorId, a_FileUnit, nStatus)
  call XF_CLOSE_GROUP(xVectorId, error)
  if (nStatus < 0) then
    WRITE(*,*) 'Error reading vector dataset.'
    error = -1
  return
  endif
enddo

if (allocated(Paths)) deallocate(Paths)

! find multidataset folders
call XF_GET_GRP_PTHS_SZ_MLT_DSETS(a_xGroupId, nMultiDatasets, &
                                                      nMaxPathLength, nStatus)
if (nStatus >= 0 .AND. nMultiDatasets > 0) then
  allocate(Paths(nMultiDatasets*nMaxPathLength))
  call XF_GET_ALL_GRP_PATHS_MLT_DSETS(a_xGroupId, nMultiDatasets, &
                                                 nMaxPathLength, Paths, error)
  if (nStatus < 0) then
    error = -1
    return
  endif

  ! Output number and paths to multidatasets
  WRITE(a_FileUnit,*) 'Number of Multidatasets ', nMultiDatasets
  do i=2, nMultiDatasets
    IndividualPath = ''
    do j=1, nMaxPathLength-1
      IndividualPath(j:j) = Paths((i-1)*nMaxPathLength+j)
    enddo
    WRITE(a_FileUnit,*) 'Reading multidataset: ', &
                                             IndividualPath(1:nMaxPathLength-1)
    call XF_OPEN_GROUP(a_xGroupId, IndividualPath(1:nMaxPathLength-1), &
                                                           xMultiId, nStatus)
    if (nStatus < 0) then
      error = -1
    return
    endif

    call TD_READ_DATASETS(xMultiId, a_FileUnit, nStatus)
    call XF_CLOSE_GROUP(xMultiId, error)
    if (nStatus < 0) then
      WRITE(*,*) 'Error reading multidatasets.'
      error = -1
    return
    endif
  enddo
endif
if (allocated(Paths)) deallocate(Paths)

error = 1
return

END SUBROUTINE
!tdReadDatasets
! --------------------------------------------------------------------------
! FUNCTION tdReadActivityScalarAIndex
! PURPOSE  Read all timestep values for a particular index
! NOTES    
! --------------------------------------------------------------------------
SUBROUTINE TD_READ_ACTIVITY_SCALAR_A_INDEX(a_Filename, a_Index, error)
CHARACTER(LEN=*), INTENT(IN) :: a_Filename
INTEGER, INTENT(IN)   :: a_Index
INTEGER, INTENT(OUT)  :: error
INTEGER                  status
INTEGER           xFileId, xDsetsId, xScalarAId
INTEGER                  nTimesteps, i
INTEGER, ALLOCATABLE  :: bActive(:)

xFileId = NONE
xDsetsId = NONE
xScalarAId = NONE

  ! open the file
call XF_OPEN_FILE(a_Filename, .TRUE., xFileId, status)
if (status < 0) then
  error = -1
  return
endif

  ! open the dataset group
call XF_OPEN_GROUP(xFileId, DATASETS_LOCATION, xDsetsId, status)
if (status >= 0) then
  call XF_OPEN_GROUP(xDsetsId, SCALAR_A_LOCATION, xScalarAId, status)
endif
if (status < 0) then
  error = status
  return
endif

  ! Find out the number of timesteps in the file
CALL XF_GET_DATASET_NUM_TIMES(xScalarAId, nTimesteps, status)
if (status < 0) then
  error = status
  return
endif

if (nTimesteps < 1) then
  error = -1
  return
endif

  ! Read the values for the index
allocate(bActive(nTimesteps))
call XF_READ_ACTIVE_VALS_AT_INDEX(xScalarAId, a_Index, 1, nTimesteps, &
                                       bActive, status)
  ! output the data
WRITE(*,*) ''
WRITE(*,*) 'Reading activity for scalar A slice at index: ', a_Index
do i=1, nTimesteps
  WRITE(*,*) bActive(i), ' '
enddo

deallocate(bActive)

error = status
return

END SUBROUTINE
! tdReadActivityScalarAIndex

! --------------------------------------------------------------------------
! FUNCTION tdReadScalarAIndex
! PURPOSE  Read all timestep values for a particular index
! NOTES    
! --------------------------------------------------------------------------
SUBROUTINE TD_READ_SCALAR_A_INDEX (a_Filename, a_Index, error)
CHARACTER(LEN=*), INTENT(IN) :: a_Filename
INTEGER, INTENT(IN)   :: a_Index
INTEGER, INTENT(OUT)  :: error
INTEGER              status
INTEGER       xFileId, xDsetsId, xScalarAId
INTEGER              nTimesteps, i
REAL, ALLOCATABLE :: fValues(:)

xFileId = NONE
xDsetsId = NONE
xScalarAId = NONE

  ! open the file
call XF_OPEN_FILE(a_Filename, .TRUE., xFileId, status)
if (status < 0) then
  error = -1
  return
endif

  ! open the dataset group
call XF_OPEN_GROUP(xFileId, DATASETS_LOCATION, xDsetsId, status)
if (status >= 0) then
  call XF_OPEN_GROUP(xDsetsId, SCALAR_A_LOCATION, xScalarAId, status)
endif
if (status < 0) then
  error = status
  return
endif

  ! Find out the number of timesteps in the file
call XF_GET_DATASET_NUM_TIMES(xScalarAId, nTimesteps, status)
if (status < 0) then
  error = status
  return
endif

if (nTimesteps < 1) then
  error = -1
  return
endif

  ! Read the values for the index
allocate (fValues(nTimesteps))
call XF_READ_SCALAR_VALUES_AT_INDEX(xScalarAId, a_Index, 1, nTimesteps, &
                                     fValues, status)

  ! output the data
WRITE(*,*) ''
WRITE(*,*) 'Reading scalar A slice at index: ', a_Index
do i=1, nTimesteps
  WRITE(*,*) fValues(i), ' '
enddo

deallocate(fValues)

error = status
return

END SUBROUTINE
! tdReadScalarAtIndex

! --------------------------------------------------------------------------
! FUNCTION tdWriteScalarA
! PURPOSE  Write scalar Dataset to an HDF5 File
! NOTES    This tests dynamic data sets, and activity
!          This dataset is dynamic concentrations (mg/L) with output times
!          in minutes.
!          Dataset is for a mesh and so nActive is the number of elements
!          which is not the same as the nValues which would be number of nodes
!          reads/writes a reference time in julian days
! --------------------------------------------------------------------------
SUBROUTINE TD_WRITE_SCALAR_A (a_Filename, a_Compression, error)
  CHARACTER(LEN=*), INTENT(IN) :: a_Filename
  INTEGER, INTENT(IN)          :: a_Compression
  INTEGER, INTENT(OUT)         :: error
  INTEGER      xFileId, xDsetsId, xScalarAId, xCoordId
  INTEGER      nValues, nTimes, nActive
  REAL(DOUBLE) dTime, dJulianReftime
  INTEGER      iTimestep, iActive, iHpgnZone
  REAL         fValues(10) ! nValues
  INTEGER*1    bActivity(10) ! activity
  INTEGER      i, status

  ! initialize the data
  nValues = 10
  nTimes = 3
  nActive = 8
  dTime = 0.0

  ! 5th item in data set is always inactive, others active
  do iActive = 1, nActive
    bActivity(iActive) = 1
  enddo 
  bActivity(6) = 0


  ! create the file
  call XF_CREATE_FILE(a_Filename, .TRUE., xFileId, error)
  if (error .LT. 0) then
      ! close the file
    call XF_CLOSE_FILE(xFileId, error)
    return
  endif

  ! create the group where we will put all the datasets 
  call XF_CREATE_GENERIC_GROUP(xFileId, DATASETS_LOCATION, xDsetsId, status)
  if (status < 0) then
    call XF_CLOSE_FILE(xFileId, error)
    error = -1
  return
  endif

  ! Create the scalar A dataset group
  call XF_CREATE_SCALAR_DATASET(xDsetsId, SCALAR_A_LOCATION, 'mg/L', &
              TS_HOURS, a_Compression, xScalarAId, status)
  if (status .LT. 0) then
      ! close the dataset
    call XF_CLOSE_GROUP(xScalarAId, error)
    call XF_CLOSE_GROUP(xDsetsId, error)
    call XF_CLOSE_FILE(xFileId, error)
    error = status
    return 
  endif

  ! Add in a reftime.  This is a julian day for:
  ! noon July 1, 2003
  dJulianReftime = 2452822.0;
  call XF_WRITE_REFTIME(xScalarAId, dJulianReftime, status)
  if (status < 0) then
    call XF_CLOSE_GROUP(xScalarAId, error)
    call XF_CLOSE_GROUP(xDsetsId, error)
    call XF_CLOSE_FILE(xFileId, error)
  endif

  ! Loop through timesteps adding them to the file
  do iTimestep = 1, nTimes
    ! We will have an 0.5 hour timestep
    dTime = iTimestep * 0.5

    fValues(1) = dTime
    do i = 2, nValues
      fValues(i) = fValues(i-1)*2.5
    end do

    ! write the dataset array values
    call XF_WRITE_SCALAR_TIMESTEP(xScalarAId, dTime, nValues, fValues, error)
    if (error .GE. 0) then
      ! write activity array
      call XF_WRITE_ACTIVITY_TIMESTEP(xScalarAId, nActive, bActivity, error)
    end if 
  enddo

  ! Write Coordinate file - for ScalarA, we will set the coordinate system
  !   to be Geographic HPGN, with HPGN settings written to the file.
  call XF_CREATE_COORDINATE_GROUP(xFileId, xCoordId, status)
  if (status < 0) then
    call XF_CLOSE_GROUP(xScalarAId, error)
    call XF_CLOSE_GROUP(xDsetsId, error)
    call XF_CLOSE_FILE(xFileId, error)
    error = -1
  return
  endif

    ! set HPGN Zone for test
  iHpgnZone = 29   ! Utah
    ! Write Coordinate Information to file
  call XF_SET_HORIZ_DATUM(xCoordId, HORIZ_DATUM_GEOGRAPHIC_HPGN, error)
  call XF_SET_HORIZ_UNITS(xCoordId, COORD_UNITS_METERS, error)
  call XF_SET_VERT_DATUM(xCoordId, VERT_DATUM_LOCAL, error)
  call XF_SET_VERT_UNITS(xCoordId, COORD_UNITS_METERS, error)

    ! write additional information
  call XF_SET_HPGN_AREA(xCoordId, iHpgnZone, error)

  call XF_CLOSE_GROUP(xCoordId, error)
  xCoordId = 0;

  ! close the dataset
  call XF_CLOSE_GROUP(xScalarAId, error)
  call XF_CLOSE_GROUP(xDsetsId, error)
  call XF_CLOSE_FILE(xFileId, error)

  return
END SUBROUTINE
! tdWriteScalarA

! --------------------------------------------------------------------------
! FUNCTION tdWriteScalarAPIECES
! PURPOSE  Write scalar Dataset to an HDF5 File
! NOTES    This tests dynamic data sets, and activity
!          This dataset is dynamic concentrations (mg/L) with output times
!          in minutes.
!          Dataset is for a mesh and so nActive is the number of elements
!          which is not the same as the nValues which would be number of nodes
!          reads/writes a reference time in julian days
! --------------------------------------------------------------------------
SUBROUTINE TD_WRITE_SCALAR_A_PIECES (a_Filename, a_Compression, error)
  CHARACTER(LEN=*), INTENT(IN) :: a_Filename
  INTEGER, INTENT(IN)          :: a_Compression
  INTEGER, INTENT(OUT)         :: error
  INTEGER      xFileId, xDsetsId, xScalarAId, xCoordId
  INTEGER      nValues, nTimes, nActive
  REAL(DOUBLE) dTime, dJulianReftime
  INTEGER      iTimestep, iActive, iHpgnZone
  REAL         fValues(10) ! nValues
  INTEGER*1    bActivity(10) ! activity
  INTEGER      i, status
  REAL         minvalue, maxvalue
  INTEGER      timestepId
  INTEGER      activeTs

  ! initialize the data
  nValues = 10
  nTimes = 3
  nActive = 8
  dTime = 0.0

  ! 5th item in data set is always inactive, others active
  do iActive = 1, nActive
    bActivity(iActive) = 1
  enddo 
  bActivity(6) = 0


  ! create the file
  call XF_CREATE_FILE(a_Filename, .TRUE., xFileId, error)
  if (error .LT. 0) then
      ! close the file
    call XF_CLOSE_FILE(xFileId, error)
    return
  endif

  ! create the group where we will put all the datasets 
  call XF_CREATE_GENERIC_GROUP(xFileId, DATASETS_LOCATION, xDsetsId, status)
  if (status < 0) then
    call XF_CLOSE_FILE(xFileId, error)
    error = -1
  return
  endif

  ! Create the scalar A dataset group
  call XF_CREATE_SCALAR_DATASET(xDsetsId, SCALAR_A_LOCATION, 'mg/L', &
              TS_HOURS, a_Compression, xScalarAId, status)
  if (status .LT. 0) then
      ! close the dataset
    call XF_CLOSE_GROUP(xScalarAId, error)
    call XF_CLOSE_GROUP(xDsetsId, error)
    call XF_CLOSE_FILE(xFileId, error)
    error = status
    return 
  endif

  ! Add in a reftime.  This is a julian day for:
  ! noon July 1, 2003
  dJulianReftime = 2452822.0;
  call XF_WRITE_REFTIME(xScalarAId, dJulianReftime, status)
  if (status < 0) then
    call XF_CLOSE_GROUP(xScalarAId, error)
    call XF_CLOSE_GROUP(xDsetsId, error)
    call XF_CLOSE_FILE(xFileId, error)
  endif

  ! Loop through timesteps adding them to the file
  do iTimestep = 1, nTimes
    ! We will have an 0.5 hour timestep
    dTime = iTimestep * 0.5

    fValues(1) = dTime
    minvalue = fValues(1)
    maxvalue = fValues(1)
    do i = 2, nValues
      fValues(i) = fValues(i-1)*2.5
      minvalue = min(minvalue, fValues(i))
      maxvalue = max(maxvalue, fValues(i))
    end do

    ! write the dataset array values
    call XF_INITIALIZE_SCALAR_TIMESTEP(xScalarAId, dTime, nValues, minvalue, &
            maxvalue, timestepId, error)
    
    ! write data in pairs
    do i = 1, nValues, +2
      call XF_WRITE_SCALAR_TIMESTEP_PORTION(xScalarAId, timestepId, 2, i, &
                                        fValues(i), error)
    enddo

    if (error .GE. 0) then
      ! write activity array
      call XF_INITIALIZE_ACTIVITY_TIMESTEP(xScalarAId, nActive, activeTs, error)
      
      do i = 1, nActive, +2
        call XF_WRITE_ACTIVITY_TIMESTEP_PORTION(xScalarAId, activeTs, 2, &
                i, bActivity(i), error)
      enddo
      
    end if 
  enddo

  ! Write Coordinate file - for ScalarA, we will set the coordinate system
  !   to be Geographic HPGN, with HPGN settings written to the file.
  call XF_CREATE_COORDINATE_GROUP(xFileId, xCoordId, status)
  if (status < 0) then
    call XF_CLOSE_GROUP(xScalarAId, error)
    call XF_CLOSE_GROUP(xDsetsId, error)
    call XF_CLOSE_FILE(xFileId, error)
    error = -1
  return
  endif

    ! set HPGN Zone for test
  iHpgnZone = 29   ! Utah
    ! Write Coordinate Information to file
  call XF_SET_HORIZ_DATUM(xCoordId, HORIZ_DATUM_GEOGRAPHIC_HPGN, error)
  call XF_SET_HORIZ_UNITS(xCoordId, COORD_UNITS_METERS, error)
  call XF_SET_VERT_DATUM(xCoordId, VERT_DATUM_LOCAL, error)
  call XF_SET_VERT_UNITS(xCoordId, COORD_UNITS_METERS, error)

    ! write additional information
  call XF_SET_HPGN_AREA(xCoordId, iHpgnZone, error)

  call XF_CLOSE_GROUP(xCoordId, error)
  xCoordId = 0;

  ! close the dataset
  call XF_CLOSE_GROUP(xScalarAId, error)
  call XF_CLOSE_GROUP(xDsetsId, error)
  call XF_CLOSE_FILE(xFileId, error)

  return
END SUBROUTINE
! tdWriteScalarAPieces

! --------------------------------------------------------------------------
! FUNCTION TD_WRITE_SCALAR_B
! PURPOSE  Write scalar Dataset to an HDF5 File
! NOTES    This tests dynamic data sets, and activity
!          This dataset is dynamic concentrations (mg/L) with output times
!          in minutes.
!          Dataset is for a mesh and so nActive is the number of elements
!          which is not the same as the nValues which would be number of nodes
!          reads/writes a reference time in julian days
! --------------------------------------------------------------------------
SUBROUTINE TD_WRITE_SCALAR_B (a_Filename, a_Compression, a_Overwrite, error)
  CHARACTER(LEN=*), INTENT(IN) :: a_Filename
  INTEGER, INTENT(IN)          :: a_Compression
  LOGICAL, INTENT(IN)          :: a_Overwrite
  INTEGER, INTENT(OUT)         :: error
  INTEGER      xFileId, xDsetsId, xScalarBId, xCoordId
  INTEGER      nValues, nTimes, nActive
  REAL(DOUBLE) dTime, dJulianReftime
  INTEGER      iTimestep, iActive
  REAL         fValues(10) ! nValues
  INTEGER*1    bActivity(10) ! activity
  INTEGER      i, status

  ! initialize the data
  nValues = 10
  nTimes = 3
  nActive = 8
  dTime = 0.0
  i = 0

  ! 5th item in data set is always inactive, others active
  do iActive = 1, nActive
    bActivity(iActive) = 1
  enddo 
  bActivity(6) = 0

  if (a_Overwrite) then
      ! open the already-existing file
    call XF_OPEN_FILE(a_Filename, .FALSE., xFileId, status)
    if (status < 0) then
      error = -1
      return
    endif
      ! open the group where we have all the datasets
    call XF_OPEN_GROUP(xFileId, DATASETS_LOCATION, xDsetsId, status)
    if (status < 0) then
      call XF_CLOSE_FILE(xFileId, error)
      error = -1
      return
    endif
  else
      ! create the file
    call XF_CREATE_FILE(a_Filename, .TRUE., xFileId, error)
    if (error .LT. 0) then
        ! close the file
      call XF_CLOSE_FILE(xFileId, error)
      return
    endif

      ! create the group where we will put all the datasets 
    call XF_CREATE_GENERIC_GROUP(xFileId, DATASETS_LOCATION, xDsetsId, status)
    if (status < 0) then
      call XF_CLOSE_FILE(xFileId, error)
      error = -1
    return
    endif
  endif

  ! Create/Overwrite the scalar B dataset group
  call XF_CREATE_SCALAR_DATASET(xDsetsId, SCALAR_B_LOCATION, 'mg/L', &
              TS_HOURS, a_Compression, xScalarBId, status)
  if (status < 0) then
      ! close the dataset
    call XF_CLOSE_GROUP(xScalarBId, error)
    call XF_CLOSE_GROUP(xDsetsId, error)
    call XF_CLOSE_FILE(xFileId, error)
    error = status
    return 
  endif

  ! Add in a reftime.  This is a julian day for:
  ! noon July 1, 2003
  dJulianReftime = 2452822.0;
  call XF_WRITE_REFTIME(xScalarBId, dJulianReftime, status)
  if (status < 0) then
    call XF_CLOSE_GROUP(xScalarBId, error)
    call XF_CLOSE_GROUP(xDsetsId, error)
    call XF_CLOSE_FILE(xFileId, error)
  endif

  if (.NOT. a_Overwrite) then
      ! Loop through timesteps adding them to the file
    do iTimestep = 1, nTimes
        ! We will have an 0.5 hour timestep
      dTime = iTimestep * 0.5

      fValues(1) = dTime
      do i = 2, nValues
        fValues(i) = fValues(i-1)*2.5
      end do

        ! write the dataset array values
      call XF_WRITE_SCALAR_TIMESTEP(xScalarBId, dTime, nValues, fValues, error)
      if (error .GE. 0) then
          ! write activity array
        call XF_WRITE_ACTIVITY_TIMESTEP(xScalarBId, nActive, bActivity, error)
      end if
      if (error < 0) then
          call XF_CLOSE_GROUP(xScalarBId, error)
          call XF_CLOSE_GROUP(xDsetsId, error)
          call XF_CLOSE_FILE(xFileId, error)
      endif
    enddo
  else
      ! Loop through timesteps adding them to the file
    do iTimestep = 1, nTimes
        ! We will have an 1.5 hour timestep
      dTime = iTimestep * 1.5

      fValues(1) = dTime
      do i = 2, nValues
        fValues(i) = fValues(i-1)*1.5
      end do

        ! write the dataset array values
      call XF_WRITE_SCALAR_TIMESTEP(xScalarBId, dTime, nValues, fValues, error)
      if (error .GE. 0) then
          ! write activity array
        call XF_WRITE_ACTIVITY_TIMESTEP(xScalarBId, nActive, bActivity, error)
      end if
      if (error < 0) then
          call XF_CLOSE_GROUP(xScalarBId, error)
          call XF_CLOSE_GROUP(xDsetsId, error)
          call XF_CLOSE_FILE(xFileId, error)
      endif
    enddo
  endif

  if (.NOT. a_Overwrite) then
    ! Write Coordinate file - for ScalarB, we will set the coordinate system
    !   to be UTM, with UTM Zone settings written to the file.
    call XF_CREATE_COORDINATE_GROUP(xFileId, xCoordId, status)
    if (status < 0) then
      call XF_CLOSE_GROUP(xScalarBId, error)
      call XF_CLOSE_GROUP(xDsetsId, error)
      call XF_CLOSE_FILE(xFileId, error)
    error = -1
    return
    endif

     ! Write Coord Info to file
    call XF_SET_HORIZ_DATUM(xCoordId, HORIZ_DATUM_UTM, error)
    call XF_SET_HORIZ_UNITS(xCoordId, COORD_UNITS_METERS, error)

    call XF_SET_VERT_DATUM(xCoordId, VERT_DATUM_LOCAL, error)
    call XF_SET_VERT_UNITS(xCoordId, COORD_UNITS_METERS, error)

      ! write additional information - we'll use the max value for this test
    call XF_SET_UTM_ZONE(xCoordId, UTM_ZONE_MAX, error)

    call XF_CLOSE_GROUP(xCoordId, error)
    xCoordId = 0
  endif

  ! close the dataset
  call XF_CLOSE_GROUP(xScalarBId, error)
  call XF_CLOSE_GROUP(xDsetsId, error)
  call XF_CLOSE_FILE(xFileId, error)

  error = 1
  return
END SUBROUTINE
! tdWriteScalarB
!------------------------------------------------------------------------------
!  FUNCTION TD_WRITE_COORDS_TO_MULTI
!  PURPOSE  Write coordinate system to a multidataset file
!  NOTES
!------------------------------------------------------------------------------
SUBROUTINE TD_WRITE_COORDS_TO_MULTI (a_xFileId, error)
INTEGER, INTENT(IN) :: a_xFileId
INTEGER, INTENT(OUT)       :: error
INTEGER    xCoordId
INTEGER           status

  ! Write Coordinate file - for Multidatasets, we will set the coordinate system
  !   to be UTM, with UTM Zone settings written to the file.
  call XF_CREATE_COORDINATE_GROUP(a_xFileId, xCoordId, status)
  if (status < 0) then
    error = status
  return
  endif

    ! Write Coord Info to file
  call XF_SET_HORIZ_DATUM(xCoordId, HORIZ_DATUM_UTM, error)
  call XF_SET_HORIZ_UNITS(xCoordId, COORD_UNITS_METERS, error)

  call XF_SET_VERT_DATUM(xCoordId, VERT_DATUM_LOCAL, error)
  call XF_SET_VERT_UNITS(xCoordId, COORD_UNITS_METERS, error)

    ! write additional information - we'll use the max value for this test
  call XF_SET_UTM_ZONE(xCoordId, UTM_ZONE_MAX, error)

  call XF_CLOSE_GROUP(xCoordId, error)
  xCoordId = 0

  return
END SUBROUTINE

! --------------------------------------------------------------------------
! FUNCTION tdWriteScalarAToMulti
! PURPOSE  Write scalar Dataset to an HDF5 File
! NOTES    This tests dynamic data sets, and activity
!          This dataset is dynamic concentrations (mg/L) with output times
!          in minutes.
!          Dataset is for a mesh and so nActive is the number of elements
!          which is not the same as the nValues which would be number of nodes
!          reads/writes a reference time in julian days
! --------------------------------------------------------------------------
SUBROUTINE TD_WRITE_SCALAR_A_TO_MULTI (a_GroupID, status)
 ! CHARACTER(LEN=*), INTENT(IN) :: a_Filename
 ! INTEGER, INTENT(IN)          :: a_Compression
 ! INTEGER, INTENT(OUT)         :: error
  INTEGER      xFileId, xDsetsId, xScalarAId
  INTEGER      a_GroupID
  INTEGER      nValues, nTimes, nActive
  REAL(DOUBLE) dTime, dJulianReftime
  INTEGER      iTimestep, iActive
  REAL         fValues(10) ! nValues
  INTEGER*1    bActivity(10) ! activity
  INTEGER      i, status

  ! initialize the data
  nValues = 10
  nTimes  = 3
  nActive = 8
  dTime   = 0.0

  ! 5th item in data set is always inactive, others active
  do iActive = 1, nActive
    bActivity(iActive) = 1
  enddo 
  bActivity(6) = 0

  ! Create the scalar A dataset group
  call XF_CREATE_SCALAR_DATASET(a_GroupID, SCALAR_A_LOCATION, 'mg/L', &
              TS_HOURS, NONE, xScalarAId, status)
  if (status .LT. 0) then
      ! close the dataset
    call XF_CLOSE_GROUP(xScalarAId, status)
    call XF_CLOSE_GROUP(xDsetsId, status)
    call XF_CLOSE_FILE(xFileId, status)
    return 
  endif

  ! Add in a reftime.  This is a julian day for:
  ! noon July 1, 2003
  dJulianReftime = 2452822.0;
  call XF_WRITE_REFTIME(xScalarAId, dJulianReftime, status)
  if (status < 0) then
    call XF_CLOSE_GROUP(xScalarAId, status)
    call XF_CLOSE_GROUP(xDsetsId, status)
    call XF_CLOSE_FILE(xFileId, status)
  endif

  ! Loop through timesteps adding them to the file
  do iTimestep = 1, nTimes
    ! We will have an 0.5 hour timestep
    dTime = iTimestep * 0.5

    fValues(1) = dTime
    do i = 2, nValues
      fValues(i) = fValues(i-1)*2.5
    end do

    ! write the dataset array values
    call XF_WRITE_SCALAR_TIMESTEP(xScalarAId, dTime, nValues, fValues, status)
    if (status .GE. 0) then
      ! write activity array
      call XF_WRITE_ACTIVITY_TIMESTEP(xScalarAId, nActive, bActivity, status)
    end if 
  enddo

  ! close the dataset
  call XF_CLOSE_GROUP(xScalarAId, status)
  !call XF_CLOSE_GROUP(a_GroupID, status)
  !call XF_CLOSE_FILE(a_FileID, status)

  return
END SUBROUTINE
! tdWriteScalarAToMulti
! --------------------------------------------------------------------------
! FUNCTION tdReadVector2DAIndex
! PURPOSE  Read all timestep values for a particular index
! NOTES    
! --------------------------------------------------------------------------
SUBROUTINE TD_READ_VECTOR2D_A_INDEX (a_Filename, a_Index, error)
CHARACTER(LEN=*), INTENT(IN) :: a_Filename
INTEGER, INTENT(IN)   :: a_Index
INTEGER, INTENT(OUT)  :: error
INTEGER           status
INTEGER    xFileId, xDsetsId, xVector2DA
INTEGER           nTimesteps, i
REAL, ALLOCATABLE     :: fValues(:)

xFileId = NONE
xDsetsId = NONE
xVector2DA = NONE

  ! open the file
call XF_OPEN_FILE(a_Filename, .TRUE., xFileId, status)
if (status < 0) then
  error = -1
  return
endif

  ! open the dataset group
call XF_OPEN_GROUP(xFileId, DATASETS_LOCATION, xDsetsId, status)
if (status >= 0) then
  call XF_OPEN_GROUP(xDsetsId, VECTOR2D_A_LOCATION, xVector2DA, status)
endif
if (status < 0) then
  error = status
  return
endif

  ! Find out the number of timesteps in the file
call XF_GET_DATASET_NUM_TIMES(xVector2DA, nTimesteps, status)
if (status < 0) then
  error = status
  return
endif

if (nTimesteps < 1) then
  error = -1
  return
endif

  ! Read the values for the index
allocate(fValues(nTimesteps*2))
call XF_READ_VECTOR_VALUES_AT_INDEX(xVector2DA, a_Index, 1, nTimesteps, 2, &
                                     fValues, status)

  ! output the data
WRITE(*,*) ''
WRITE(*,*) 'Reading vector 2D A slice at index: ', a_Index
do i=1, nTimesteps
  WRITE(*,*) fValues(i*2-1), ' ', fValues(i*2)
enddo
WRITE(*,*) ''

deallocate(fValues)

error = status
return

END SUBROUTINE
!tdReadVector2DAIndex

! --------------------------------------------------------------------------
! FUNCTION tdWriteVector2D_A
! PURPOSE  Write scalar Dataset to an HDF5 File
! NOTES    This tests dynamic data sets, and activity
!          This dataset is dynamic concentrations (mg/L) with output times
!          in minutes.
!          Dataset is for a mesh and so nActive is the number of elements
!          which is not the same as the nValues which would be number of nodes
!          reads/writes a reference time in julian days
! --------------------------------------------------------------------------
SUBROUTINE TD_WRITE_VECTOR2D_A (a_Filename, a_Compression, error)
  CHARACTER(LEN=*), INTENT(IN) :: a_Filename
  INTEGER, INTENT(IN)          :: a_Compression
  INTEGER, INTENT(OUT)         :: error
  INTEGER      xFileId, xDsetsId, xVector2D_A, xCoordId
  INTEGER      nValues, nTimes, nComponents, nActive
  REAL(DOUBLE) dTime
  INTEGER      iTimestep, iActive
  REAL, DIMENSION(2, 100) :: fValues ! nComponents, nValues
  INTEGER*1    bActivity(100) ! activity
  INTEGER      i, j, status
  INTEGER      iHpgnZone

  ! initialize the data
  nComponents = 2
  nValues = 100
  nTimes = 6
  nActive = 75
  dTime = 0.0

  ! 5th item in data set is always inactive, others active
  bActivity(1) = 0
  do iActive = 2, nActive
    if (mod(iActive-1, 3) == 0) then
      bActivity(iActive) = 0
    else
    bActivity(iActive) = 1
  endif
  enddo

  ! create the file
  call XF_CREATE_FILE(a_Filename, .TRUE., xFileId, error)
  if (error .LT. 0) then
    ! close the dataset
    call XF_CLOSE_FILE(xFileId, error)
    return
  endif

  ! create the group where we will put all the datasets 
  call XF_CREATE_GENERIC_GROUP(xFileId, DATASETS_LOCATION, xDsetsId, status)
  if (status < 0) then
    call XF_CLOSE_FILE(xFileId, error)
    error = -1
  return
  endif

  ! Create the vector dataset group
  call XF_CREATE_VECTOR_DATASET(xDsetsId, VECTOR2D_A_LOCATION, 'ft/s', &
              TS_SECONDS, a_Compression, xVector2D_A, status)
  if (status .LT. 0) then
      ! close the dataset
    call XF_CLOSE_GROUP(xVector2D_A, error)
    call XF_CLOSE_GROUP(xDsetsId, error)
    call XF_CLOSE_FILE(xFileId, error)
    error = status
    return 
  endif

  ! Loop through timesteps adding them to the file
  do iTimestep = 1, nTimes
    ! We will have an 0.5 hour timestep
    dTime = iTimestep * 0.5

    do i = 1, nValues
      do j = 1, nComponents
        fValues(j,i) = ((i-1)*nComponents + j)*dTime
      end do
    end do

    ! write the dataset array values
    call XF_WRITE_VECTOR_TIMESTEP(xVector2D_A, dTime, nValues, nComponents, &
                                  fValues, error)
    if (error .GE. 0) then
      ! write activity array
      call XF_WRITE_ACTIVITY_TIMESTEP(xVector2D_A, nActive, bActivity, error)
    end if 
  enddo

  ! Write Coordinate file - for Vector2D_A, we will set the coordinate system
  !   to be Geographic HPGN, with HPGN settings written to the file.
  call XF_CREATE_COORDINATE_GROUP(xFileId, xCoordId, status)
  if (status < 0) then
    call XF_CLOSE_GROUP(xVector2D_A, error)
    call XF_CLOSE_GROUP(xDsetsId, error)
    call XF_CLOSE_FILE(xFileId, error)
  error = -1
  return
  endif

    ! set HPGN info for test
  iHpgnZone = 29   ! Utah

  call XF_SET_HORIZ_DATUM(xCoordId, HORIZ_DATUM_GEOGRAPHIC_HPGN, error)
  call XF_SET_HORIZ_UNITS(xCoordId, COORD_UNITS_METERS, error)
  call XF_SET_VERT_DATUM(xCoordId, VERT_DATUM_LOCAL, error)
  call XF_SET_VERT_UNITS(xCoordId, COORD_UNITS_METERS, error)

    ! write additional information
  call XF_SET_HPGN_AREA(xCoordId, iHpgnZone, error)

  call XF_CLOSE_GROUP(xCoordId, error)
  xCoordId = 0

  ! close the dataset
  call XF_CLOSE_GROUP(xVector2D_A, error)
  call XF_CLOSE_GROUP(xDsetsId, error)
  call XF_CLOSE_FILE(xFileId, error)

  return
END SUBROUTINE
! tdWriteVector2D_A

! --------------------------------------------------------------------------
! FUNCTION tdWriteVector2D_A_PIECES
! PURPOSE  Write scalar Dataset to an HDF5 File
! NOTES    This tests dynamic data sets, and activity
!          This dataset is dynamic concentrations (mg/L) with output times
!          in minutes.
!          Dataset is for a mesh and so nActive is the number of elements
!          which is not the same as the nValues which would be number of nodes
!          reads/writes a reference time in julian days
! --------------------------------------------------------------------------
SUBROUTINE TD_WRITE_VECTOR2D_A_PIECES (a_Filename, a_Compression, error)
  CHARACTER(LEN=*), INTENT(IN) :: a_Filename
  INTEGER, INTENT(IN)          :: a_Compression
  INTEGER, INTENT(OUT)         :: error
  INTEGER      xFileId, xDsetsId, xVector2D_A, xCoordId
  INTEGER      nValues, nTimes, nComponents, nActive
  REAL(DOUBLE) dTime
  INTEGER      iTimestep, iActive
  REAL*4, DIMENSION(2, 100) :: fValues ! nComponents, nValues
  INTEGER*1    bActivity(100) ! activity
  INTEGER      i, j, status
  INTEGER      iHpgnZone
  INTEGER      nValuesToWrite, nComponentsToWrite, startComponent
  REAL*4       minvalue, maxvalue
  INTEGER      timeId
  REAL*8       mag

  ! initialize the data
  nComponents = 2
  nValues = 100
  nTimes = 6
  nActive = 75
  dTime = 0.0

  ! 5th item in data set is always inactive, others active
  bActivity(1) = 0
  do iActive = 2, nActive
    if (mod(iActive-1, 3) == 0) then
      bActivity(iActive) = 0
    else
    bActivity(iActive) = 1
  endif
  enddo

  ! create the file
  call XF_CREATE_FILE(a_Filename, .TRUE., xFileId, error)
  if (error .LT. 0) then
    ! close the dataset
    call XF_CLOSE_FILE(xFileId, error)
    return
  endif

  ! create the group where we will put all the datasets 
  call XF_CREATE_GENERIC_GROUP(xFileId, DATASETS_LOCATION, xDsetsId, status)
  if (status < 0) then
    call XF_CLOSE_FILE(xFileId, error)
    error = -1
  return
  endif

  ! Create the vector dataset group
  call XF_CREATE_VECTOR_DATASET(xDsetsId, VECTOR2D_A_LOCATION, 'ft/s', &
              TS_SECONDS, a_Compression, xVector2D_A, status)
  if (status .LT. 0) then
      ! close the dataset
    call XF_CLOSE_GROUP(xVector2D_A, error)
    call XF_CLOSE_GROUP(xDsetsId, error)
    call XF_CLOSE_FILE(xFileId, error)
    error = status
    return 
  endif

  ! Loop through timesteps adding them to the file
  do iTimestep = 1, nTimes
    ! We will have an 0.5 hour timestep
    dTime = iTimestep * 0.5

    do i = 1, nValues
      do j = 1, nComponents
        fValues(j,i) = ((i-1)*nComponents + j)*dTime
      end do
    end do

    minvalue = 99999.0
    maxvalue = 0.0
    do i = 1, nValues
      mag = 0.0
      do j = 1, nComponents
        mag = mag + fValues(j, i)**2
      end do
      mag = mag**0.5
     
      minvalue = min(minvalue, mag)
      maxvalue = max(maxvalue, mag)
    end do

    ! write the dataset array values
    call XF_INITIALIZE_VECTOR_TIMESTEP(xVector2D_A, dTime, nValues, nComponents, &
            minvalue, maxvalue, timeId, error)
    
    nValuesToWrite = 2
    nComponentsToWrite = 2
    startComponent = 1
    
    do i = 1, nValues, +2
      call XF_WRITE_VECTOR_TIMESTEP_PORTION(xVector2D_A, timeId, nValuesToWrite, &
              nComponentsToWrite, i, startComponent, fValues(1, i), error)
    enddo

    if (error .GE. 0) then
      ! write activity array
      call XF_WRITE_ACTIVITY_TIMESTEP(xVector2D_A, nActive, bActivity, error)
    end if 
  enddo

  ! Write Coordinate file - for Vector2D_A, we will set the coordinate system
  !   to be Geographic HPGN, with HPGN settings written to the file.
  call XF_CREATE_COORDINATE_GROUP(xFileId, xCoordId, status)
  if (status < 0) then
    call XF_CLOSE_GROUP(xVector2D_A, error)
    call XF_CLOSE_GROUP(xDsetsId, error)
    call XF_CLOSE_FILE(xFileId, error)
  error = -1
  return
  endif

    ! set HPGN info for test
  iHpgnZone = 29   ! Utah

  call XF_SET_HORIZ_DATUM(xCoordId, HORIZ_DATUM_GEOGRAPHIC_HPGN, error)
  call XF_SET_HORIZ_UNITS(xCoordId, COORD_UNITS_METERS, error)
  call XF_SET_VERT_DATUM(xCoordId, VERT_DATUM_LOCAL, error)
  call XF_SET_VERT_UNITS(xCoordId, COORD_UNITS_METERS, error)

    ! write additional information
  call XF_SET_HPGN_AREA(xCoordId, iHpgnZone, error)

  call XF_CLOSE_GROUP(xCoordId, error)
  xCoordId = 0

  ! close the dataset
  call XF_CLOSE_GROUP(xVector2D_A, error)
  call XF_CLOSE_GROUP(xDsetsId, error)
  call XF_CLOSE_FILE(xFileId, error)

  return
END SUBROUTINE
! tdWriteVector2D_A_PIECES

! --------------------------------------------------------------------------
! FUNCTION TD_WRITE_VECTOR2D_B
! PURPOSE  Write scalar Dataset to an HDF5 File
! NOTES    This tests dynamic data sets, and activity
!          This dataset is dynamic concentrations (mg/L) with output times
!          in minutes.
!          Dataset is for a mesh and so nActive is the number of elements
!          which is not the same as the nValues which would be number of nodes
!          reads/writes a reference time in julian days
! --------------------------------------------------------------------------
SUBROUTINE TD_WRITE_VECTOR2D_B (a_Filename, a_Compression, a_Overwrite, error)
  CHARACTER(LEN=*), INTENT(IN) :: a_Filename
  INTEGER, INTENT(IN)          :: a_Compression
  LOGICAL, INTENT(IN)          :: a_Overwrite
  INTEGER, INTENT(OUT)         :: error
  INTEGER      xFileId, xDsetsId, xVector2D_B, xCoordId
  INTEGER      nValues, nTimes, nComponents, nActive
  REAL(DOUBLE) dTime
  INTEGER      iTimestep, iActive
  REAL, DIMENSION(2, 100) :: fValues
  INTEGER*1    bActivity(100)
  INTEGER      i, j, status

    ! initialize the data
  nComponents = 2
  nValues = 100
  nTimes = 6
  nActive = 75
  dTime = 0.0

    ! 5th item in data set is always inactive, others active
  bActivity(1) = 0
  do iActive = 2, nActive
    if (mod(iActive-1, 3) == 0) then
      bActivity(iActive) = 0
    else
      bActivity(iActive) = 1
    endif
  enddo

  if (a_Overwrite) then
      ! open the already-existing file
    call XF_OPEN_FILE(a_Filename, .FALSE., xFileId, status)
    if (status < 0) then
      error = -1
      return
    endif
      ! open the group where we have all the datasets
    call XF_OPEN_GROUP(xFileId, DATASETS_LOCATION, xDsetsId, status)
    if (status < 0) then
      call XF_CLOSE_FILE(xFileId, error)
      error = -1
      return
    endif
  else
      ! create the file
    call XF_CREATE_FILE(a_Filename, .TRUE., xFileId, error)
    if (error .LT. 0) then
        ! close the dataset
      call XF_CLOSE_FILE(xFileId, error)
      return
    endif

      ! create the group where we will put all the datasets 
    call XF_CREATE_GENERIC_GROUP(xFileId, DATASETS_LOCATION, xDsetsId, status)
    if (status < 0) then
      call XF_CLOSE_FILE(xFileId, error)
      error = -1
      return
    endif
  endif

    ! Create/Overwrite the vector dataset group
  call XF_CREATE_VECTOR_DATASET(xDsetsId, VECTOR2D_B_LOCATION, 'ft/s', &
                                TS_SECONDS, a_Compression, xVector2D_B, status)
  if (status .LT. 0) then
      ! close the dataset
    call XF_CLOSE_GROUP(xVector2D_B, error)
    call XF_CLOSE_GROUP(xDsetsId, error)
    call XF_CLOSE_FILE(xFileId, error)
    error = status
    return 
  endif

  if (.NOT. a_Overwrite) then
      ! Loop through timesteps adding them to the file
    do iTimestep = 1, nTimes
        ! We will have an 0.5 hour timestep
      dTime = iTimestep * 0.5
      do i = 1, nValues
        do j = 1, nComponents
          fValues(j,i) = ((i-1)*nComponents + j)*dTime
        end do
      end do
        ! write the dataset array values
      call XF_WRITE_VECTOR_TIMESTEP(xVector2D_B, dTime, nValues, nComponents, &
                                    fValues, error)
      if (error .GE. 0) then
          ! write activity array
        call XF_WRITE_ACTIVITY_TIMESTEP(xVector2D_B, nActive, bActivity, error)
      end if
      if (error < 0) then
        call XF_CLOSE_GROUP(xVector2D_B, error)
        call XF_CLOSE_GROUP(xDsetsId, error)
        call XF_CLOSE_FILE(xFileId, error)
      endif
    enddo
  else
      ! Loop through timesteps adding them to the file
    do iTimestep = 1, nTimes
        ! We will have an 1.5 hour timestep
      dTime = iTimestep * 1.5
      do i = 1, nValues
        do j = 1, nComponents
          fValues(j,i) = ((i-1)*nComponents + j)*dTime
        end do
      end do
        ! write the dataset array values
      call XF_WRITE_VECTOR_TIMESTEP(xVector2D_B, dTime, nValues, nComponents, &
                                    fValues, error)
      if (error .GE. 0) then
          ! write activity array
        call XF_WRITE_ACTIVITY_TIMESTEP(xVector2D_B, nActive, bActivity, error)
      end if
      if (error < 0) then
        call XF_CLOSE_GROUP(xVector2D_B, error)
        call XF_CLOSE_GROUP(xDsetsId, error)
        call XF_CLOSE_FILE(xFileId, error)
      endif
    enddo
  endif

  if (.NOT. a_Overwrite) then
    ! Write Coordinate file - for ScalarB, we will set the coordinate system
    !   to be UTM, with UTM Zone settings written to the file.
    call XF_CREATE_COORDINATE_GROUP(xFileId, xCoordId, status)
    if (status < 0) then
      call XF_CLOSE_GROUP(xVector2D_B, error)
      call XF_CLOSE_GROUP(xDsetsId, error)
      call XF_CLOSE_FILE(xFileId, error)
    error = -1
    return
    endif

      ! write the coordinate data to the file
    call XF_SET_HORIZ_DATUM(xCoordId, HORIZ_DATUM_UTM, error)
    call XF_SET_HORIZ_UNITS(xCoordId, COORD_UNITS_METERS, error)
    call XF_SET_VERT_DATUM(xCoordId, VERT_DATUM_LOCAL, error)
    call XF_SET_VERT_UNITS(xCoordId, COORD_UNITS_METERS, error)

      ! write additional information - we'll use the max UTM zone for the test
    call XF_SET_UTM_ZONE(xCoordId, UTM_ZONE_MAX, error)

    call XF_CLOSE_GROUP(xCoordId, error)
    xCoordId = 0
  endif

  ! close the dataset
  call XF_CLOSE_GROUP(xVector2D_B, error)
  call XF_CLOSE_GROUP(xDsetsId, error)
  call XF_CLOSE_FILE(xFileId, error)

  return
END SUBROUTINE
! tdWriteVector2D_B

! --------------------------------------------------------------------------
! FUNCTION tdWriteVector2D_AToMulti
! PURPOSE  Write scalar Dataset to an HDF5 File
! NOTES    This tests dynamic data sets, and activity
!          This dataset is dynamic concentrations (mg/L) with output times
!          in minutes.
!          Dataset is for a mesh and so nActive is the number of elements
!          which is not the same as the nValues which would be number of nodes
!          reads/writes a reference time in julian days
! --------------------------------------------------------------------------
SUBROUTINE TD_WRITE_VECTOR2D_A_TO_MULTI (a_FileID, a_GroupID, status)
  INTEGER      xVector2D_A
  INTEGER      a_FileID, a_GroupID
  INTEGER      nValues, nTimes, nComponents, nActive
  REAL(DOUBLE) dTime
  INTEGER      iTimestep, iActive
  REAL, DIMENSION(2, 100) :: fValues ! nComponents, nValues
  INTEGER*1    bActivity(100) ! activity
  INTEGER      i, j, status

  ! initialize the data
  nComponents = 2
  nValues = 100
  nTimes = 6
  nActive = 75
  dTime = 0.0

  ! 5th item in data set is always inactive, others active
  bActivity(1) = 0
  do iActive = 2, nActive
    if (mod(iActive-1, 3) == 0) then
      bActivity(iActive) = 0
    else
    bActivity(iActive) = 1
  endif
  enddo

  ! Create the vector dataset group
  call XF_CREATE_VECTOR_DATASET(a_GroupID, VECTOR2D_A_LOCATION, 'ft/s', &
              TS_SECONDS, NONE, xVector2D_A, status)
  if (status .LT. 0) then
      ! close the dataset
    call XF_CLOSE_GROUP(xVector2D_A, status)
    call XF_CLOSE_GROUP(a_GroupID, status)
    call XF_CLOSE_FILE(a_FileID, status)
    return 
  endif

  ! Loop through timesteps adding them to the file
  do iTimestep = 1, nTimes
    ! We will have an 0.5 hour timestep
    dTime = iTimestep * 0.5

    do i = 1, nValues
      do j = 1, nComponents
        fValues(j,i) = ((i-1)*nComponents + j)*dTime
      end do
    end do

    ! write the dataset array values
    call XF_WRITE_VECTOR_TIMESTEP(xVector2D_A, dTime, nValues, nComponents, &
                                  fValues, status)
    if (status .GE. 0) then
      ! write activity array
      call XF_WRITE_ACTIVITY_TIMESTEP(xVector2D_A, nActive, bActivity, status)
    end if 
  enddo

  ! close the dataset
  call XF_CLOSE_GROUP(xVector2D_A, status)
  return
END SUBROUTINE
! tdWriteVector2D_AToMulti
! --------------------------------------------------------------------------
! FUNCTION tdiReadScalar
! PURPOSE  Read a scalar from an XMDF file and output information to
!          to a text file
! NOTES    
! --------------------------------------------------------------------------
SUBROUTINE TDI_READ_SCALAR (a_xScalarId, FileUnit, error)
  INTEGER, INTENT(IN) ::  a_xScalarId
  INTEGER, INTENT(IN) ::         FileUnit
  INTEGER, INTENT(OUT) :: error
  INTEGER             nTimes, nValues, nActive
  LOGICAL*2             bUseReftime
  INTEGER             iTime
  CHARACTER(LEN=100)   TimeUnits
  REAL(DOUBLE), ALLOCATABLE :: Times(:)
  REAL, ALLOCATABLE         :: Values(:), Minimums(:), Maximums(:)
  INTEGER, ALLOCATABLE      :: Active(:)
  REAL(DOUBLE)                 Reftime
nTimes = NONE
nValues = NONE
nActive = None

  ! read the time units
  call XF_GET_DATASET_TIME_UNITS(a_xScalarId, TimeUnits, error)
  if (error < 0) return

  WRITE(FileUnit,*) 'Time units: ', TimeUnits(1:LEN_TRIM(TimeUnits))

  ! see if we are using a reftime
  call XF_USE_REFTIME (a_xScalarId, bUseReftime, error)
  if (error < 0) then
    return
  endif
  if (bUseReftime) then
    call XF_READ_REFTIME (a_xScalarId, Reftime, error)
    if (error < 0) then
      return
  endif
    WRITE(FileUnit,*) 'Reftime: ', Reftime
  endif

  ! read in the number of values and number of active values
  call XF_GET_DATASET_NUMVALS(a_xScalarId, nValues, error)
  if (error .GE. 0) then
    call XF_GET_DATASET_NUMACTIVE(a_xScalarId, nActive, error)
  endif
  if (error .LT. 0) return 

  if (nValues <= 0) then
    WRITE(FileUnit, *) 'No data to read in.'
    error = -1
    return 
  endif

  ! read in the number of times
  call XF_GET_DATASET_NUM_TIMES(a_xScalarId, nTimes, error)
  if (error < 0) then
    return 
  endif

  ! Read in the individual time values
  allocate(Times(nTimes))

  call XF_GET_DATASET_TIMES(a_xScalarId, nTimes, Times, error)
  if (error < 0) return 

  ! Read in the minimum and maximum values
  allocate(Minimums(nTimes))
  allocate(Maximums(nTimes))

  call XF_GET_DATASET_MINS(a_xScalarId, nTimes, Minimums, error)
  if (error >= 0) then
    call XF_GET_DATASET_MAXS(a_xScalarId, nTimes, Maximums, error)
  endif
  if (error < 0) then
    deallocate(Times)
    deallocate(Minimums)
    deallocate(Maximums)
    return
  endif

  allocate(Values(nValues))
  if (nActive .GT. 0) then
    allocate(Active(nActive))
  endif

  WRITE(FileUnit,*) 'Number Timesteps: ', nTimes
  WRITE(FileUnit,*) 'Number Values: ', nValues
  WRITE(FileUnit,*) 'Number Active: ', nActive
  WRITE(FileUnit,*) ''

  ! loop through the timesteps, read the values and active values and write
  ! them to the text file
  do iTime = 1, nTimes
    call XF_READ_SCALAR_VALUES_TIMESTEP(a_xScalarId, iTime, nValues, Values, error)
    if (error >= 0 .AND. nActive > 0) then
      call XF_READ_ACTIVITY_TIMESTEP(a_xScalarId, iTime, nActive, Active, error)
    endif

    ! Write the time, min, max, values and active values to the text output
    ! file.
    WRITE(FileUnit,*) 'Timestep at  ', Times(iTime)
    WRITE(FileUnit,*) 'Min: ', Minimums(iTime)
    WRITE(FileUnit,*) 'Max: ', Maximums(iTime)

    WRITE(FileUnit,*) 'Values:'
    WRITE(FileUnit,*) Values(1:nValues)
    WRITE(FileUnit,*) ''

    WRITE(FileUnit,*) 'Activity:'
    WRITE(FileUnit,*) Active(1:nActive)
    WRITE(FileUnit,*) ''
  end do

  if (allocated(Times)) then
    deallocate(Times)
  endif
  
  if (allocated(Minimums)) then
    deallocate(Minimums)
  endif

  if (allocated(Maximums)) then
    deallocate(Maximums)
  endif

  if (allocated(Values)) then
    deallocate(Values)
  endif

  if (allocated(Active)) then
    deallocate(Active)
  endif

  return
END SUBROUTINE
! tdiReadScalar

! --------------------------------------------------------------------------
! FUNCTION TDI_READ_VECTOR
! PURPOSE  Read a vector from an XMDF file and output information to
!          to a text file
! NOTES    
! --------------------------------------------------------------------------
SUBROUTINE TDI_READ_VECTOR (a_xVectorId, FileUnit, error)
  INTEGER, INTENT(IN) ::  a_xVectorId
  INTEGER, INTENT(IN) ::         FileUnit
  INTEGER, INTENT(OUT) :: error
  INTEGER             nTimes, nValues, nActive, nComponents
  INTEGER             iTime, i
  LOGICAL*2            bUseReftime
  CHARACTER(LEN=100)   TimeUnits
  REAL(DOUBLE), ALLOCATABLE :: Times(:)
  REAL, ALLOCATABLE, DIMENSION (:, :) :: Values
  REAL, ALLOCATABLE         :: Minimums(:), Maximums(:)
  INTEGER, ALLOCATABLE      :: Active(:)
  REAL(DOUBLE)                 Reftime

nTimes = NONE
nValues = NONE
nActive = NONE
nComponents = NONE

  ! read the time units
  call XF_GET_DATASET_TIME_UNITS(a_xVectorId, TimeUnits, error)
  if (error < 0) return

  WRITE(FileUnit,*) 'Time units: ', TimeUnits(1:LEN_TRIM(TimeUnits))

  ! see if we are using a reftime
  call XF_USE_REFTIME (a_xVectorId, bUseReftime, error)
  if (error < 0) then
    return
  endif
  if (bUseReftime) then
    call XF_READ_REFTIME (a_xVectorId, Reftime, error)
    if (error < 0) then
      return
  endif
    WRITE(FileUnit,*) 'Reftime: ', Reftime
  endif

  ! read in the number of values and number of active values
  call XF_GET_DATASET_NUMVALS(a_xVectorId, nValues, error)
  if (error .GE. 0) then
    call XF_GET_DATASET_NUMCOMPONENTS(a_xVectorId, nComponents, error)
    if (error .GE. 0) then
      call XF_GET_DATASET_NUMACTIVE(a_xVectorId, nActive, error)
    endif
  endif
  if (error .LT. 0) return 

  if (nValues <= 0) then
    WRITE(FileUnit, *) 'No data to read in.'
    error = -1
    return 
  endif

  ! read in the number of times
  call XF_GET_DATASET_NUM_TIMES(a_xVectorId, nTimes, error)
  if (error < 0) then
    return 
  endif

  ! Read in the individual time values
  allocate(Times(nTimes))

  call XF_GET_DATASET_TIMES(a_xVectorId, nTimes, Times, error)
  if (error < 0) return 

  ! Read in the minimum and maximum values
  allocate(Minimums(nTimes))
  allocate(Maximums(nTimes))

  call XF_GET_DATASET_MINS(a_xVectorId, nTimes, Minimums, error)
  if (error >= 0) then
    call XF_GET_DATASET_MAXS(a_xVectorId, nTimes, Maximums, error)
  endif
  if (error < 0) then
    deallocate(Times)
    deallocate(Minimums)
    deallocate(Maximums)
    return
  endif

  allocate(Values(nComponents, nValues))
  if (nActive .GT. 0) then
    allocate(Active(nActive))
  endif

  WRITE(FileUnit,*) 'Number Timesteps: ', nTimes
  WRITE(FileUnit,*) 'Number Values: ', nValues
  WRITE(FileUnit,*) 'Number Components: ', nComponents
  WRITE(FileUnit,*) 'Number Active: ', nActive

  ! loop through the timesteps, read the values and active values and write
  ! them to the text file
  do iTime = 1, nTimes
    call XF_READ_VECTOR_VALUES_TIMESTEP(a_xVectorId, iTime, nValues, &
                                        nComponents, Values, error)
    if (error >= 0 .AND. nActive > 0) then
      call XF_READ_ACTIVITY_TIMESTEP(a_xVectorId, iTime, nActive, Active, error)
    endif

    ! Write the time, min, max, values and active values to the text output
    ! file.
  WRITE(FileUnit,*) ''
    WRITE(FileUnit,*) 'Timestep at  ', Times(iTime)
    WRITE(FileUnit,*) 'Min: ', Minimums(iTime)
    WRITE(FileUnit,*) 'Max: ', Maximums(iTime)

    WRITE(FileUnit,*) 'Values:'
    do i=1, nValues
      WRITE(FileUnit,*) Values(1:nComponents,i:i)
    enddo
    WRITE(FileUnit,*) ''

    WRITE(FileUnit,*) 'Activity:'
    WRITE(FileUnit,*) Active(1:nActive)
    WRITE(FileUnit,*) ''    
  WRITE(FileUnit,*) ''    

  end do

  if (allocated(Times)) then
    deallocate(Times)
  endif
  
  if (allocated(Minimums)) then
    deallocate(Minimums)
  endif

  if (allocated(Maximums)) then
    deallocate(Maximums)
  endif

  if (allocated(Values)) then
    deallocate(Values)
  endif

  if (allocated(Active)) then
    deallocate(Active)
  endif

  return
END SUBROUTINE
! tdiReadVector

END MODULE TestDatasets

TestDatasets.f90 tests datasets

MODULE TestGeomPaths

USE TestDatasets
USE Xmdf
USE ErrorDefinitions
USE XmdfDefs

CONTAINS

SUBROUTINE TM_READ_TEST_PATHS(Filename, OutFilename, error)
  CHARACTER(LEN=*), INTENT(IN) :: Filename, OutFilename
  INTEGER, INTENT(OUT)         :: error
  INTEGER                         OutUnit  
  INTEGER               :: xFileId, xGroupId
  INTEGER                         nGroups, nMaxPathLength, nDims, nPaths, nTimes
  INTEGER                         i, j, nStatus
  CHARACTER,ALLOCATABLE        :: Paths(:)
  REAL(DOUBLE), ALLOCATABLE    :: times(:), locs(:)
  REAL(DOUBLE)                    Nullval
  CHARACTER(LEN=BIG_STRING_SIZE) :: IndividualPath
  INTEGER                           StartLoc
  INTEGER, DIMENSION(2)        :: PathIndices
  INTEGER                         iTime, iPath

  ! open the XMDF file 
  CALL XF_OPEN_FILE(Filename, .FALSE., xFileId, nStatus)
  if (nStatus < 0) then
   WRITE(*) 'Unable to open XMDF file TM_READ_TEST_PATHS.'
    error = nStatus;
    return
  endif

  ! open the Output file 
  OutUnit = 79
  OPEN(UNIT=OutUnit, FILE=OutFilename, STATUS='REPLACE', ACTION='WRITE', &
     IOSTAT = error)
  if (OutUnit == 0) then
    call XF_CLOSE_FILE(xFileId, error)
    error = -1
    return
  endif

  ! find the geomotric path groups
  ! Get the number and paths of datasets in the file.
  CALL XF_GRP_PTHS_SZ_FOR_GEOM_PTHS(xFileId, nGroups,             &
                                              nMaxPathLength, nStatus)
  if (nStatus >= 0 .AND. nGroups > 0) then
    allocate (Paths(nMaxPathLength*nGroups))
    CALL XF_GRP_PTHS_FOR_GEOM_PTHS(xFileId, nGroups,                &
                                         nMaxPathLength, Paths, nStatus)
  endif
  if (nStatus < 0) then
    CALL XF_CLOSE_FILE(xFileId, nStatus)
  error = -1
    return
  endif
  ! Report the number and paths to individual geom paths groups in the file.
  WRITE(OutUnit,*) 'Number of geometric paths in file: ', nGroups

  do i = 1, nGroups 
    StartLoc = (i-1)*nMaxPathLength + 1
  IndividualPath = ''
  do j = 1, nMaxPathLength - 1
    IndividualPath(j:j) = Paths(StartLoc+j-1)
  enddo
    WRITE(OutUnit,*) 'Reading particles in group: ', &
                                     IndividualPath(1:LEN_TRIM(IndividualPath))

    CALL XF_OPEN_GROUP(xFileId, IndividualPath, xGroupId, nStatus)
    if (nStatus >= 0) then
      ! read the dimensionality of the paths
      CALL XF_GET_PATH_DIMENSIONALITY(xGroupId, nDims, nStatus)
      if (nStatus >= 0) then
        WRITE(OutUnit,*) 'Group dimensionality:', nDims
        ! read the number of paths
        CALL XF_GET_NUMBER_OF_PATHS(xGroupId, nPaths, nStatus)
        if (nStatus >= 0) then
          WRITE(OutUnit,*) 'Number of paths in group:', nPaths
          ! read the number of timesteps
          CALL XF_GET_NUMBER_OF_TIMES(xGroupId, nTimes, nStatus)
          if (nStatus >= 0) then
            WRITE(OutUnit,*) 'Number of timesteps in group:', nTimes
            ! allocate times array
      allocate(times(nTimes))
            CALL XF_GET_PATH_TIMES_ARRAY(xGroupId, nTimes, times, nStatus)
            if (nStatus >= 0) then
                ! get space for the data location values
              allocate(locs(nPaths*nDims))
                ! read null value 
              CALL XF_GET_PATH_NULL_VAL(xGroupId, NullVal, nStatus)
              if (nStatus >= 0) then
                WRITE(OutUnit,*) 'Null value:', NullVal
                do iTime=1, nTimes
                  WRITE(OutUnit,*) 'Timestep: ', times(iTime)
                    ! read the data for this timestep
                  CALL XF_READ_PATH_LOCATIONS_AT_TIME(xGroupId, iTime,    &
                                                     1, nPaths, locs, nStatus)
                  if (nStatus >= 0) then
                      ! write  the data for this timestep
                    WRITE(OutUnit,*) '  X        Y'
                    if (nDims == 3) then
                      WRITE(OutUnit,*) '       Z'
                    endif
                    WRITE(OutUnit,*) ''
                    do j=1, nPaths 
                      if (locs(j*nDims) == NullVal) then
                        WRITE(OutUnit,*) 'Particle not active yet'
                      else 
                        WRITE(OutUnit,*) locs((j-1)*nDims+1), ' ', locs((j-1)*nDims+2)
                        if (nDims == 3) then
                          WRITE(OutUnit,*) ' ', locs((j-1)*nDims+3)
                        endif
                        WRITE(OutUnit,*) ''
                      endif
                    enddo ! Loop through the paths
                  endif ! if timestep read
                enddo ! Loop through timesteps
              endif ! if null value read
              if (allocated(locs)) deallocate (locs)
            endif ! if we get get the times array
                    ! get space for the data location values - 1 particle 
                    ! all times
        allocate(locs(nTimes*nDims))
            do iPath=1, nPaths
                ! read the values for this particle for all timesteps
              CALL XF_READ_PATH_LOCS_FOR_PART(xGroupId, iPath, &
                                   1, nTimes, locs, nStatus)
              if (nStatus >= 0) then
                ! write  the data for this path
                WRITE(OutUnit,*) 'Time       X        Y'
                if (nDims == 3) then
                  WRITE(OutUnit,*) '        Z'
                endif
                WRITE(OutUnit, *) ''
                do j=1, nTimes
                  if (locs((j-1)*nDims+1) .NE. NullVal) then
                    WRITE(OutUnit,*) times(j), ' ', locs((j-1)*nDims + 1), ' ', &
                            locs((j-1)*nDims+2)
                    if (nDims == 3) then
                      WRITE(OutUnit,*) ' ', locs((j-1)*nDims+3)
                    endif
                    WRITE(OutUnit,*) ''
                  endif
                enddo ! Loop through the times
              endif ! if read path locations for particle
            enddo ! Loop through the paths
            if (allocated(locs)) deallocate(locs)
          endif ! if allocation of locs suceeded
                ! get space for the data location values - 2 particle 
                ! all times
      allocate(locs(nTimes*nDims*2))
          PathIndices(1) = 1
          PathIndices(2) = nPaths                
          CALL XF_READ_PATH_LOCS_FOR_PARTS(xGroupId, 2,    &
                              PathIndices, 1, nTimes, locs, nStatus)
          if (nStatus >= 0) then
            ! write  the data for these 2 paths
            if (nDims == 3) then
              WRITE(OutUnit,*) 'Timestep       X1        Y1        Z1        Xn        Yn        Zn'
            else
              WRITE(OutUnit,*) 'Timestep       X1        Y1        Xn        Yn'
            endif
            do j=1, nTimes
              if (nDims == 3) then
                WRITE(OutUnit,*) times(j), ' ', locs((j-1)*2*nDims+1), ' ',  &
             locs((j-1)*2*nDims+2),' ', locs((j-1)*2*nDims+3),' ',       &
           locs((j-1)*2*nDims+4), locs((j-1)*2*nDims+5), ' ',locs((j-1)*2*nDims+6)
              else
                WRITE(OutUnit,*) times(j),' ', locs((j-1)*2*nDims + 1), ' ',         &
                           locs((j-1)*2*nDims+2), ' ',locs((j-1)*2*nDims+3)
              endif
      enddo
        else
          WRITE(*,*) 'Error reading path locations for multiple particles'
        endif
      endif
          if (allocated(locs)) deallocate(locs)
          if (allocated(times)) deallocate(times)
          CALL XF_CLOSE_GROUP(xGroupId, nStatus)
        if (nStatus < 0) then
          WRITE(*)'Error reading geometric paths..'
        endif
      endif
  endif
  enddo ! loop through groups
    ! free the paths
  if (allocated(Paths)) deallocate(Paths)
    ! close the files
  close(OutUnit)
  CALL XF_CLOSE_FILE(xFileId, nStatus)
  error = nStatus;
  return

END SUBROUTINE
! tmReadTestPaths 

SUBROUTINE TM_WRITE_TEST_PATHS(Filename, Compression, error)
  CHARACTER(LEN=*), INTENT(IN) :: Filename
  INTEGER, INTENT(IN)         :: Compression
  INTEGER, INTENT(OUT)         :: error 
  INTEGER                      nPaths
  REAL(DOUBLE), DIMENSION(6)   :: pLocs
  REAL, DIMENSION(2)           :: pSpeeds
  REAL(DOUBLE)                 NullVal

  REAL                         NullValSinglePrec
  INTEGER               xFileId, xPathGroupId, xSpeedId, xPropGroupId
  INTEGER                      status
  REAL(DOUBLE)                 Timeval

  nPaths = 0
  NullVal = -99999.9d0
  NullValSinglePrec = -99999.9

  ! create the file
  call XF_CREATE_FILE(Filename, .TRUE., xFileId, status)
  if (status < 0) then
    error = -1
    return
  endif

  ! create the group to store the particle paths 
  call XF_CREATE_GEOMETRIC_PATH_GROUP(xFileId, "particles", &
           'abcdefglasdfjoaieur', Compression, xPathGroupId, NullVal, status)
  if (status < 0) then
    error = -1
    return
  endif

  ! create the data set to store the speed 
  call XF_CREATE_SCALAR_DSET_EXTNDBL(xPathGroupId, 'Vmag', 'm/s', &
             TS_SECONDS, NullValSinglePrec, Compression, xSpeedId, status)
  if (status < 0) then
    error = -1
    return
  endif

  call XF_CREATE_PROPERTY_GROUP(xSpeedId, xPropGroupId, status)
  if (status < 0) then
    call XF_CLOSE_GROUP(xSpeedId, status)
    error = -1
    return
  endif
  call XF_WRITE_PROPERTY_FLOAT(xPropGroupId, PROP_NULL_VALUE, 1, &
                               NullValSinglePrec, -1, status)
  call XF_CLOSE_GROUP(xPropGroupId, status)

  ! Setup the arrays for the path group at timestep 0 
  ! particle location at the first timestep
  nPaths = 1
  pLocs(1) = 1.0
  pLocs(2) = 2.0
  pLocs(3) = 3.0
  
    ! store the particles for the first timestep
  Timeval = 0.0
  call XF_WRITE_PARTICLE_TIMESTEP(xPathGroupId, 3, Timeval, nPaths, pLocs, status)
  if (status < 0) then
    error = -1
    return 
  endif
   ! set up and store the speed at timestep 0
  pSpeeds(1) = 1.1
  Timeval = 0.0
  call XF_WRITE_SCALAR_TIMESTEP(xSpeedId, Timeval, 1, pSpeeds, status) 
  if (status < 0) then
    error = -1
    return
  endif

  ! Setup the arrays for the path group at timestep 1 
  ! particle location at the first timestep
  pLocs(1) = 4.0
  pLocs(2) = 5.0
  pLocs(3) = 6.0
  
  ! store the particles for the second timestep
  Timeval = 1.0
  call XF_WRITE_PARTICLE_TIMESTEP(xPathGroupId, 3, Timeval, nPaths, pLocs, status)
  if (status < 0) then
    error = -1
    return
  endif
  
  ! set up and store the speed at timestep 2
  pSpeeds(1) = 1.2
  call XF_WRITE_SCALAR_TIMESTEP(xSpeedId, Timeval, 1, pSpeeds, status) 
  if (status < 0) then
    error = -1
    return
  endif

  ! Setup the arrays for the path group at timestep 3-add a particle
  ! particle location at the first timestep
  nPaths = 2
  pLocs(1) = 7.0
  pLocs(2) = 8.0
  pLocs(3) = 9.0
  pLocs(4) = -1.0
  pLocs(5) = -2.0
  pLocs(6) = -3.0
  
  ! store the particles for the timestep 3
  Timeval = 2.0
  call XF_WRITE_PARTICLE_TIMESTEP(xPathGroupId, 3, Timeval, nPaths, pLocs, status)
  if (status < 0) then
    error = -1
    return
  endif

  ! extend the data set for speed
  call XF_EXTEND_SCALAR_DATASET(xSpeedId, 2, status)
  if (status < 0) then
    error = -1
    return
  endif

  ! set up and store the speed at timestep 4
  pSpeeds(1) = 1.3
  pSpeeds(2) = 2.1
  call XF_WRITE_SCALAR_TIMESTEP(xSpeedId, Timeval, 2, pSpeeds, status) 
  if (status < 0) then
    error = -1
    return
  endif

  ! Setup the arrays for the path group at timestep 3-inactive particle(static)
  ! particle location at the first timestep
  pLocs(1) = 7.0
  pLocs(2) = 8.0
  pLocs(3) = 9.0
  pLocs(4) = -4.0
  pLocs(5) = -5.0
  pLocs(6) = -6.0
  
  ! store the particles for timestep 4
  Timeval = 3.0
  call XF_WRITE_PARTICLE_TIMESTEP(xPathGroupId, 3, Timeval, nPaths, pLocs, status)
  if (status < 0) then
    error = -1
    return
  endif

  ! set up and store the speed at timestep 4
  pSpeeds(1) = NullVal
  pSpeeds(2) = 2.2
  call XF_WRITE_SCALAR_TIMESTEP(xSpeedId, Timeval, 2, pSpeeds, status) 
  if (status < 0) then
    error = -1
    return
  endif

  ! close the resources
  call XF_CLOSE_GROUP(xSpeedId, status)
  call XF_CLOSE_GROUP(xPathGroupId, status)
  call XF_CLOSE_FILE(xFileId, status)

  error = 1
  return
END SUBROUTINE

END MODULE TestGeomPaths

TestGeomPaths.f90 tests geometry paths

MODULE TestGrid

USE ErrorDefinitions
USE XmdfDefs
USE Xmdf

CHARACTER(LEN=*),PARAMETER :: GRID_CART2D_GROUP_NAME = 'Grid Cart2D Group'
CHARACTER(LEN=*),PARAMETER :: GRID_CURV2D_GROUP_NAME = 'Grid Curv2D Group'
CHARACTER(LEN=*),PARAMETER :: GRID_CART3D_GROUP_NAME = 'Grid Cart3D Group'

CONTAINS

!****************************
! ---------------------------------------------------------------------------
! FUNCTION  TG_READ_GRID
! PURPOSE   Read a grid and write data to a text file
! NOTES     
! ---------------------------------------------------------------------------
SUBROUTINE  TG_READ_GRID(a_Id, a_Outfile, error)
INTEGER, INTENT(IN) :: a_Id
INTEGER, INTENT(IN)        :: a_Outfile
INTEGER, INTENT(OUT)       :: error
INTEGER nGridType, nExtrudeType, nDims, nCellsI, nCellsJ
INTEGER nCellsK, nLayers, nOrientation, nValsI, nValsJ
INTEGER nValsK, nExtrudeVals, nCompOrigin, nUDir
CHARACTER(265) strGridType, strExtrudeType
LOGICAL*2 bDefined
REAL(DOUBLE) dOrigin(3), dBearing, dDip, dRoll
REAL(DOUBLE), ALLOCATABLE :: dExtrudeVals(:), dCoordI(:), dCoordJ(:)
REAL(DOUBLE), ALLOCATABLE :: dCoordK(:)
INTEGER i

nGridType = 0
nExtrudeType = 0
nDims = 0
nCellsI = 0
nCellsJ = 0
nCellsK = 0
nLayers = 0
nOrientation = 0
nValsI = 0
nValsJ = 0
nValsK = 0
nExtrudeVals = 0
nCompOrigin = 1
nUDir = 1
bDefined = .FALSE.
dBearing = 0.0
dDip = 0.0
dRoll =0.0
i = 0
error = 1

  ! Grid type
call XF_GET_GRID_TYPE(a_Id, nGridType, error)
if (error < 0) then
  return
endif
SELECT CASE (nGridType) 
  CASE (GRID_TYPE_CARTESIAN)
    strGridType = 'Cartesian'
  CASE (GRID_TYPE_CURVILINEAR)
    strGridType = 'Curvilinear'
  CASE (GRID_TYPE_CARTESIAN_EXTRUDED)
    strGridType = 'Cartesian extruded'
  CASE (GRID_TYPE_CURVILINEAR_EXTRUDED)
    strGridType = 'Curvilinear extruded'
  CASE DEFAULT
    WRITE(*,*) 'Invalid grid type'
    error = -1
END SELECT
WRITE(a_Outfile,*) 'The grid type is: ', strGridType(1:LEN_TRIM(strGridType))

  ! Number of dimensions
call XF_GET_NUMBER_OF_DIMENSIONS(a_Id, nDims, error)
if (error .LT. 0) then
  return
endif
if (nDims == 2) then
  WRITE(a_Outfile,*) 'The grid is two-dimensional'
elseif (nDims == 3) then
  WRITE(a_Outfile,*) 'The grid is three-dimensional'
else
  WRITE(*,*) 'The grid dimensions are invalid'
  error = -1
endif

  ! Extrusion type if applicable
if (nGridType .EQ. GRID_TYPE_CARTESIAN_EXTRUDED .OR. &
    nGridType .EQ. GRID_TYPE_CURVILINEAR_EXTRUDED) then
  call XF_GET_EXTRUSION_TYPE(a_Id, nExtrudeType, error)
  if (error < 0) then
    return
  endif
  SELECT CASE (nExtrudeType)
    case (EXTRUDE_SIGMA)
      strExtrudeType = 'Sigma stretch'
    case (EXTRUDE_CARTESIAN)
      strExtrudeType = 'Cartesian'
    case (EXTRUDE_CURV_AT_CORNERS)
      strExtrudeType = 'Curvilinear at Corners'
    case (EXTRUDE_CURV_AT_CELLS)
      strExtrudeType = 'Curvilinear at Cells'
  END SELECT
  WRITE(a_Outfile,*) 'The grid is extruding using: ', &
                     strExtrudeType(1:LEN_TRIM(strExtrudeType))
endif

  ! Origin
call XF_ORIGIN_DEFINED(a_Id, bDefined, error)
if (error < 0) then
  return
endif
if (bDefined) then
  call XF_GET_ORIGIN(a_Id, dOrigin(1), dOrigin(2), dOrigin(3), error)
  if (error < 0) then
    return
  endif
  WRITE(a_Outfile,*) 'The grid origin is ', dOrigin(1), ' ',&
                      dOrigin(2), ' ', dOrigin(3)
endif
  
  ! Orientation
call XF_GET_ORIENTATION(a_Id, nOrientation, error)
if (error < 0) then
  return
endif
if (nOrientation == ORIENTATION_RIGHT_HAND) then
  WRITE(a_Outfile,*) 'The grid has a right hand orientation'
elseif (nOrientation == ORIENTATION_LEFT_HAND) then
  WRITE(a_Outfile,*) 'The grid has a left hand orientation'
else 
  WRITE(*,*) 'Invalid grid orientation';
  error = -1
  return
endif

  ! Bearing
call XF_BEARING_DEFINED(a_Id, bDefined, error)
if (error < 0) then
  return
endif
if (bDefined) then
  call XF_GET_BEARING(a_Id, dBearing, error)
  if (error < 0) then
    return
  endif
  WRITE(a_Outfile,*) 'The grid bearing is ', dBearing
endif

  ! Dip
call XF_DIP_DEFINED(a_Id, bDefined, error)
if (error < 0) then
  return
endif
if (bDefined) then
  call XF_GET_DIP(a_Id, dDip, error)
  if (error < 0) then
    return
  endif
  WRITE(a_Outfile,*) 'The grid Dip is ', dDip
endif

if (nDims == 3) then
  ! Roll
  call XF_ROLL_DEFINED(a_Id, bDefined, error)
  if (error < 0) then
    return
  endif
  if (bDefined) then
    call XF_GET_ROLL(a_Id, dRoll, error)
    if (error < 0) then
      return
    endif
    WRITE(a_Outfile,*) 'The grid Roll is ', dRoll
  endif
endif

  ! Computational origin
call XF_COMPUTATIONAL_ORIGIN_DEFINED(a_Id, bDefined, error)
if (error < 0) then
  return
endif
if (bDefined) then
  call XF_GET_COMPUTATIONAL_ORIGIN(a_Id, nCompOrigin, error)
  if (error < 0) then
    return
  endif
  WRITE(a_Outfile,*) 'The grid Computational Origin is ', nCompOrigin
else 
  WRITE(a_Outfile,*) 'The grid Computational Origin is not defined';
endif


  ! U Direction
call XF_GET_U_DIRECTION_DEFINED(a_Id, bDefined, error)
if (error < 0) then
  return 
endif
if (bDefined) then
  call XF_GET_U_DIRECTION(a_Id, nUDir, error)
  if (error < 0) then
    return
  endif
    WRITE(a_Outfile,*) 'The grid U Direction is ', nUDir
else 
  WRITE(a_Outfile,*) 'The grid U Direction is not defined'
endif

  ! number of cells in each direction
call XF_GET_NUMBER_CELLS_IN_I(a_Id, nCellsI, error)
if (error >= 0) then
  call XF_GET_NUMBER_CELLS_IN_J(a_Id, nCellsJ, error)
  if ((error >= 0) .AND. (nDims == 3)) then 
    call XF_GET_NUMBER_CELLS_IN_K(a_Id, nCellsK, error)
  endif
endif
if (error < 0) then
  return
endif
WRITE(a_Outfile,*) 'Number of cells in I ', nCellsI
WRITE(a_Outfile,*) 'Number of cells in J ', nCellsJ
if (nDims == 3) then
  WRITE(a_Outfile,*) 'Number of cells in K ', nCellsK
endif

  ! Grid coordinates 
if (nGridType == GRID_TYPE_CARTESIAN .OR. &
    nGridType == GRID_TYPE_CARTESIAN_EXTRUDED) then
  nValsI = nCellsI
  nValsJ = nCellsJ
  if (nDims == 3) then
    nValsK = nCellsK
  endif
elseif (nGridType == GRID_TYPE_CURVILINEAR .OR. &
        nGridType == GRID_TYPE_CURVILINEAR_EXTRUDED) then
  if (nDims == 3) then
      ! three dimensions
    nValsK = (nCellsI + 1) * (nCellsJ + 1) * (nCellsK + 1)
    nValsJ = nValsK
    nValsI = nValsJ
  else
      ! two dimensions
    nValsJ = (nCellsI + 1) * (nCellsJ + 1)
    nValsI = nValsJ
  endif
else
  WRITE(*,*) 'Invalid grid type'
  error = -1
  return
endif

ALLOCATE(dCoordI(nValsI))
ALLOCATE(dCoordJ(nValsJ))
if (nDims == 3) then
  ALLOCATE(dCoordK(nValsK))
endif

call XF_GET_GRID_COORDS_I(a_Id, nValsI, dCoordI, error)
if (error >= 0) then
  call XF_GET_GRID_COORDS_J(a_Id, nValsJ, dCoordJ, error)
  if ((error >= 0) .AND. (nDims == 3)) then
    call XF_GET_GRID_COORDS_K(a_Id, nValsK, dCoordK, error)
  endif
endif
if (error < 0) then
  WRITE(*,*) 'Error reading coordinates'
  error = -1
  return
endif

WRITE(a_Outfile,*) 'The Coordinates in direction I:'
do i = 1, nValsI 
  if (mod(i,5) == 0) then
    WRITE(a_Outfile,*) ''
  endif
  WRITE(a_Outfile,*) dCoordI(i)
enddo
WRITE(a_Outfile,*) ''

WRITE(a_Outfile,*) 'The Coordinates in direction J:'
do i = 1, nValsJ
  if (mod(i,5) == 0) then
    WRITE(a_Outfile,*) ''
  endif
  WRITE(a_Outfile,*) dCoordJ(i)
enddo
WRITE(a_Outfile,*) ''

if (nDims == 3) then
  WRITE(a_Outfile,*) 'The Coordinates in direction K:'
  do i = 1, nValsK
    if (mod(i,5) == 0) then
      WRITE(a_Outfile,*) ''
    endif
    WRITE(a_Outfile,*) dCoordK(i)
  enddo
endif
WRITE(a_Outfile,*) ''

if (ALLOCATED(dCoordI)) DEALLOCATE(dCoordI)
if (ALLOCATED(dCoordJ)) DEALLOCATE(dCoordJ)
if (ALLOCATED(dCoordK)) DEALLOCATE(dCoordK)

!  // Extrude data
if (nGridType .EQ. GRID_TYPE_CARTESIAN_EXTRUDED .OR. &
    nGridType .EQ. GRID_TYPE_CURVILINEAR_EXTRUDED) then
  call XF_GET_EXTRUDE_NUM_LAYERS(a_Id, nLayers, error)
  if (error < 0) then
    return
  endif

  SELECT CASE(nExtrudeType)
    case (EXTRUDE_SIGMA)
      nExtrudeVals = nLayers
    case (EXTRUDE_CURV_AT_CORNERS)
      nExtrudeVals = (nCellsI + 1) * (nCellsJ + 1) * nLayers
    case (EXTRUDE_CURV_AT_CELLS)
      nExtrudeVals = nCellsI * nCellsJ * nLayers
  END SELECT

  ALLOCATE(dExtrudeVals(nExtrudeVals))

  call XF_GET_EXTRUDE_VALUES(a_Id, nExtrudeVals, dExtrudeVals, error)
  if (error < 0) then
    return
  endif

  WRITE(*,*) 'The extrude values are:'
  do i = 1, nExtrudeVals
    if (mod(i,5) == 0) then
      WRITE(a_Outfile,*) ''
    endif
    WRITE(a_Outfile,*) dExtrudeVals(i)
  enddo
  if (ALLOCATED(dExtrudeVals)) DEALLOCATE(dExtrudeVals)
endif

return

END SUBROUTINE TG_READ_GRID

!****************************
!----------------------------------------------------------------------------
! FUNCTION   TG_WRITE_TEST_GRID_CART_2D
! PURPOSE    Write a file that contains data for a 2D Cartesian Grid
! NOTES      A picture of the grid is in the file (TestGridCart2D.gif)
!            returns TRUE on success and FALSE on failure
!----------------------------------------------------------------------------
SUBROUTINE TG_WRITE_TEST_GRID_CART_2D(Filename, error)
CHARACTER(LEN=*), INTENT(IN) :: Filename
INTEGER, INTENT(OUT) :: error
INTEGER        nDimensions
INTEGER        nCellsI, nCellsJ
INTEGER        nGridType
INTEGER        nCompOrigin, nUDir
REAL(DOUBLE)   dOriginX, dOriginY, dOriginZ
INTEGER        nOrientation
REAL(DOUBLE)   dBearing
REAL(DOUBLE)   PlanesI(5), PlanesJ(5)
INTEGER        i, j, iSpcZone
INTEGER xFileId, xGridId, xCoordId
INTEGER        status
INTEGER        tmpOut1, tmpOut2

  nDimensions = 2
  nCellsI = 5
  nCellsJ = 5
  nGridType = GRID_TYPE_CARTESIAN
  nCompOrigin = 4
  nUDir = -2
  dOriginX = 10.0
  dOriginY = 10.0
  dOriginZ = 0.0
  nOrientation = ORIENTATION_RIGHT_HAND
  dBearing = 45.0
  xFileId = NONE
  xGridId = NONE

    ! Fill in the grid plane data with a constant size of 30
  do i = 1, nCellsI
    PlanesI(i) = i*30.0
  enddo
  do j = 1, nCellsJ
    PlanesJ(j) = j*30.0
  enddo

    ! create the file
  call XF_CREATE_FILE(Filename, .TRUE., xFileId, status)
  if (status < 0) then
    error = -1  
    return
  endif

    ! create the group to store the grid
  call XF_CREATE_GROUP_FOR_GRID(xFileId, GRID_CART2D_GROUP_NAME, xGridId, status)
  if (status < 0) then
    call XF_CLOSE_FILE(xFileId, error)
    error = -1
    return
  endif

    ! Write the grid information to the file
  call XF_SET_GRID_TYPE(xGridId, nGridType, tmpOut1)
  call XF_SET_NUMBER_OF_DIMENSIONS(xGridId, nDimensions, tmpOut2)

  if ((tmpOut1 < 0) .OR. (tmpOut2 < 0)) then
    call XF_CLOSE_GROUP(xGridId, error)
    call XF_CLOSE_FILE(xFileId, error)
    error = -1
    return
  endif

    ! set origin and orientation
  call XF_SET_ORIGIN(xGridId, dOriginX, dOriginY, dOriginZ, tmpOut1) 
  call XF_SET_ORIENTATION(xGridId, nOrientation, tmpOut2)

  if ((tmpOut1 < 0) .OR. (tmpOut2 .LT. 0)) then
    call XF_CLOSE_GROUP(xGridId, error)
    call XF_CLOSE_FILE(xFileId, error)
    error = -1
    return
  endif

    ! Set bearing
  call XF_SET_BEARING(xGridId, dBearing, tmpOut1)
  if (tmpOut1 < 0) then
    call XF_CLOSE_GROUP(xGridId, error)
    call XF_CLOSE_FILE(xFileId, error)
    error = -1
    return
  endif

   ! Set computational origin
  call XF_SET_COMPUTATIONAL_ORIGIN(xGridId, nCompOrigin, tmpOut1)
  if (tmpOut1 < 0) then
    call XF_CLOSE_GROUP(xGridId, error)
    call XF_CLOSE_FILE(xFileId, error)
    error = -1
    return
  endif

    ! Set u direction
  call XF_SET_U_DIRECTION(xGridId, nUDir, tmpOut1)
  if (tmpOut1 < 0) then
    call XF_CLOSE_GROUP(xGridId, error)
    call XF_CLOSE_FILE(xFileId, error)
    error = -1
    return
  endif

    ! Write the grid geometry to the file
    ! Set the number of cells in each direction
  call XF_SET_NUMBER_CELLS_IN_I(xGridId, nCellsI, tmpOut1)
  call XF_SET_NUMBER_CELLS_IN_J(xGridId, nCellsJ, tmpOut2)
  if ((tmpOut1 < 0) .OR. (tmpOut2 < 0)) then
    call XF_CLOSE_GROUP(xGridId, error)
    call XF_CLOSE_FILE(xFileId, error)
    error = -1
    return
  endif

    ! Set the grid plane locations
  call XF_SET_GRID_COORDS_I(xGridId, nCellsI, PlanesI, tmpOut1)
  call XF_SET_GRID_COORDS_J(xGridId, nCellsJ, PlanesJ, tmpOut2)  
  if ((tmpOut1 < 0) .OR. (tmpOut2 < 0)) then
    call XF_CLOSE_GROUP(xGridId, error)
    call XF_CLOSE_FILE(xFileId, error)
    error = -1
    return
  endif

  ! Write Coordinate file - for GridCart2D, we will set the coordinate system
  !   to be State Plane NAD27.
  call XF_CREATE_COORDINATE_GROUP(xFileId, xCoordId, status)
  if (status < 0) then
    call XF_CLOSE_GROUP(xGridId, error)
    call XF_CLOSE_FILE(xFileId, error)
  error = status
  return
  endif

  iSpcZone = 3601  ! Oregon North

  call XF_SET_HORIZ_DATUM(xCoordId, HORIZ_DATUM_STATE_PLANE_NAD27, error)
  call XF_SET_HORIZ_UNITS(xCoordId, COORD_UNITS_US_FEET, error)

  call XF_SET_VERT_DATUM(xCoordId, VERT_DATUM_LOCAL, error)
  call XF_SET_VERT_UNITS(xCoordId, COORD_UNITS_US_FEET, error)

    ! write additional information
  call XF_SET_SPC_ZONE(xCoordId, iSpcZone, error)

  call XF_CLOSE_GROUP(xCoordId, error)
  xCoordId = 0

  ! release memory
call XF_CLOSE_GROUP(xGridId, error)
call XF_CLOSE_FILE(xFileId, error)
return

END SUBROUTINE TG_WRITE_TEST_GRID_CART_2D

!****************************
!------------------------------------------------------------------------------
! FUNCTION   tgWriteTestGridCurv2D
! PURPOSE    Write a file that contains data for a 2D Curvilinear Grid
! NOTES      A picture of the grid is TestGridCurv2D.gif
!            returns TRUE on success and FALSE on failure
!------------------------------------------------------------------------------
SUBROUTINE TG_WRITE_TEST_GRID_CURV_2D(Filename, Compression, error)
CHARACTER(LEN=*), INTENT(IN) :: Filename
INTEGER, INTENT(IN) :: Compression
INTEGER, INTENT(OUT) :: error
INTEGER        nDimensions
INTEGER        nCompOrigin, nUDir
INTEGER        nCellsI, nCellsJ
INTEGER        nCells
INTEGER        nCorners
INTEGER        nGridType
REAL(DOUBLE)   xVals(16), yVals(16) !There are 16 corners
INTEGER        i
INTEGER xFileId, xGridId, xPropId, xDatasetsId, xScalarId
INTEGER xCoordId
REAL(DOUBLE)   dNullValue(1)
INTEGER        nOrientation
REAL           fDsetCellVals(6) !For cell-centered dataset
REAL           fDsetCornerVals(12) !For corner-centered dataset
REAL(DOUBLE)   tempdouble, dCppLat, dCppLon
INTEGER*1      bDsetCellActive(6)
INTEGER*1      bDsetCornerActive(12)
INTEGER        status
INTEGER        tmpOut1, tmpOut2

  nDimensions = 2
  nCompOrigin = 1
  nUDir = 1;
  nCellsI = 2
  nCellsJ = 3
  nCells = nCellsI*nCellsJ
  nCorners = (nCellsI + 1)*(nCellsJ + 1)
  nGridType = GRID_TYPE_CURVILINEAR
  xFileId = NONE
  xGridId = NONE
  xPropId = NONE
  xDatasetsId = NONE
  xScalarId = NONE
  dNullValue(1) = -999.0
  nOrientation = ORIENTATION_RIGHT_HAND

    ! There is no cell in the top right corner so we have a NullValue for
    ! the top right corner

    ! xValues row by row
  xVals(1) = 0.0
  xVals(2) = 7.5
  xVals(3) = 15.0
  xVals(4) = 2.5
  xVals(5) = 10.0
  xVals(6) = 17.5
  xVals(7) = 3.5
  xVals(8) = 11.0
  xVals(9) = 18.5
  xVals(10) = 0.0
  xVals(11) = 7.5
  xVals(12) = dNullValue(1)

    ! yValues row by row
  yVals(1) = 0.0
  yVals(2) = 0.0
  yVals(3) = 0.0
  yVals(4) = 10.0
  yVals(5) = 10.0
  yVals(6) = 10.0
  yVals(7) = 20.0
  yVals(8) = 20.0
  yVals(9) = 20.0
  yVals(10) = 30.0
  yVals(11) = 30.0
  yVals(12) = dNullValue(1)

    ! cell centered velocity dataset values
  fDsetCellVals(1) = 2.1
  fDsetCellVals(2) = 2.0
  fDsetCellVals(3) = 1.9
  fDsetCellVals(4) = 2.3
  fDsetCellVals(5) = 2.5
  fDsetCellVals(6) = dNullValue(1)

    ! all are active except the last value
  do i = 1, nCells
    bDsetCellActive(i) = 1
  enddo
  bDsetCellActive(nCells) = 0

    ! corner centered elevation dataset values
  fDsetCornerVals(1) = 1.0
  fDsetCornerVals(2) = 0.8 
  fDsetCornerVals(3) = 1.2
  fDsetCornerVals(4) = 1.4
  fDsetCornerVals(5) = 1.8
  fDsetCornerVals(6) = 2.2
  fDsetCornerVals(7) = 1.8
  fDsetCornerVals(8) = 1.4
  fDsetCornerVals(9) = 2.0
  fDsetCornerVals(10) = 1.0
  fDsetCornerVals(11) = 1.8
  fDsetCornerVals(12) = 2.2

    ! all are active except the last value
  do i = 1, nCorners
    bDsetCornerActive(i) = 1
  enddo
  bDsetCornerActive(nCorners) = 0

    ! create the file
  call XF_CREATE_FILE(Filename, .TRUE., xFileId, status)
  if (status < 0) then
    error = -1
    return
  endif

    ! create the group to store the grid
  call XF_CREATE_GROUP_FOR_GRID(xFileId, GRID_CURV2D_GROUP_NAME, xGridId, status)
  if (status < 0) then
    call XF_CLOSE_FILE(xFileId, error)
    error = -1
    return
  endif

    ! Write the grid information to the file
  call XF_SET_GRID_TYPE(xGridId, nGridType, tmpOut1)
  call XF_SET_NUMBER_OF_DIMENSIONS(xGridId, nDimensions, tmpOut2)
  if ((tmpOut1 < 0) .OR. (tmpOut2 < 0)) then
    call XF_CLOSE_GROUP(xGridId, error)
    call XF_CLOSE_FILE(xFileId, error)
    error = -1
    return
  endif

    ! set orientation
  call XF_SET_ORIENTATION(xGridId, nOrientation, tmpOut1)  
  if (tmpOut1 < 0 ) then
    call XF_CLOSE_GROUP(xGridId, error)
    call XF_CLOSE_FILE(xFileId, error)
    error = -1
    return
  endif

    ! Set computational origin
  call XF_SET_COMPUTATIONAL_ORIGIN(xGridId, nCompOrigin, tmpOut1)  
  if (tmpOut1 < 0 ) then
    call XF_CLOSE_GROUP(xGridId, error)
    call XF_CLOSE_FILE(xFileId, error)
    error = -1
    return
  endif

    ! Set u direction
  call XF_SET_U_DIRECTION(xGridId, nUDir, tmpOut1)  
  if (tmpOut1 < 0 ) then
    call XF_CLOSE_GROUP(xGridId, error)
    call XF_CLOSE_FILE(xFileId, error)
    error = -1
    return
  endif

    ! Write the grid geometry to the file
    ! Set the number of cells in each direction
  call XF_SET_NUMBER_CELLS_IN_I(xGridId, nCellsI, tmpOut1)
  call XF_SET_NUMBER_CELLS_IN_J(xGridId, nCellsJ, tmpOut2)
  if ((tmpOut1 < 0) .OR. (tmpOut2 < 0)) then
    call XF_CLOSE_GROUP(xGridId, error)
    call XF_CLOSE_FILE(xFileId, error)
    error = -1
    return
  endif
    ! Set a NullValue.  This is used to identify locations in the grid that are
    ! not being used.  In our case no geometry is defined for the top right
    ! corner.
  call XF_CREATE_GRID_PROPERTY_GROUP(xGridId, xPropId, tmpOut1)
  if (xPropId < 0) then
    call XF_CLOSE_GROUP(xGridId, error)
    call XF_CLOSE_FILE(xFileId, error)
    error = -1
    return
  endif

  call XF_WRITE_PROPERTY_DOUBLE(xPropId, PROP_NULL_VALUE, 1, dNullValue, NONE, &
                                tmpOut1)  
  if (tmpOut1 < 0) then
    call XF_CLOSE_GROUP(xPropId, error)
    call XF_CLOSE_GROUP(xGridId, error)
    call XF_CLOSE_FILE(xFileId, error)
    error = -1
    return
  endif
  call XF_CLOSE_GROUP(xPropId, error)

    ! Set the grid plane locations
  call XF_SET_GRID_COORDS_I(xGridId, nCorners, xVals, tmpOut1)
  call XF_SET_GRID_COORDS_J(xGridId, nCorners, yVals, tmpOut2)
  if ((tmpOut1 < 0) .OR. (tmpOut2 < 0)) then
    call XF_CLOSE_GROUP(xGridId, error)
    call XF_CLOSE_FILE(xFileId, error)
    error = -1
    return
  endif

    ! Create the datasets group
  call XF_CREATE_GENERIC_GROUP(xGridId, 'Datasets', xDatasetsId, tmpOut1)  
  if (tmpOut1 < 0) then
    call XF_CLOSE_GROUP(xGridId, error)
    call XF_CLOSE_FILE(xFileId, error)
    error = -1
    return
  endif

    ! Create the cell-centered dataset
  call XF_CREATE_SCALAR_DATASET(xDatasetsId, 'Velocity Mag', 'ft/s', TS_MINUTES, &
                                Compression, xScalarId, tmpOut1)  
  if (tmpOut1 < 0) then
    call XF_CLOSE_GROUP(xDatasetsId, error)
    call XF_CLOSE_GROUP(xGridId, error)
    call XF_CLOSE_FILE(xFileId, error)
    error = -1
    return
  endif

    ! specify that the dataset is cell-centered
  call XF_SCALAR_DATA_LOCATION(xScalarId, GRID_LOC_CENTER, tmpOut1)
  if (tmpOut1 < 0) then
    call XF_CLOSE_GROUP(xScalarId, error)
    call XF_CLOSE_GROUP(xDatasetsId, error)
    call XF_CLOSE_GROUP(xGridId, error)
    call XF_CLOSE_FILE(xFileId, error)
    error = -1
    return
  endif

    ! Write the data
    ! set a temporary variable to pass into the function
  tempdouble = 0.0
  call XF_WRITE_SCALAR_TIMESTEP(xScalarId, tempdouble, nCells, fDsetCellVals, tmpOut1)
  call XF_WRITE_ACTIVITY_TIMESTEP(xScalarId, nCells, bDsetCellActive, tmpOut2)
  if ((tmpOut1 < 0) .OR. (tmpOut1 < 0)) then
    call XF_CLOSE_GROUP(xScalarId, error)
    call XF_CLOSE_GROUP(xDatasetsId, error)
    call XF_CLOSE_GROUP(xGridId, error)
    call XF_CLOSE_FILE(xFileId, error)
    error = -1
    return
  endif

    ! close the cell-centered dataset
  call XF_CLOSE_GROUP(xScalarId, error)

    ! Create the corner-centered dataset
  call XF_CREATE_SCALAR_DATASET(xDatasetsId, 'elevation', 'ft', TS_MINUTES, &
                                Compression, xScalarId, tmpOut1)  
  if (tmpOut1 < 0) then
    call XF_CLOSE_GROUP(xDatasetsId, error)
    call XF_CLOSE_GROUP(xGridId, error)
    call XF_CLOSE_FILE(xFileId, error)
    error = -1
    return
  endif

    ! specify that the dataset is corner-centered
  call XF_SCALAR_DATA_LOCATION(xScalarId, GRID_LOC_CORNER, tmpOut1)  
  if (tmpOut1 < 0) then
    call XF_CLOSE_GROUP(xScalarId, error)
    call XF_CLOSE_GROUP(xDatasetsId, error)
    call XF_CLOSE_GROUP(xGridId, error)
    call XF_CLOSE_FILE(xFileId, error)
    error = -1
    return
  endif

    ! Write the data
    ! set a temporary variable to pass into the function
  tempdouble = 0.0
  call XF_WRITE_SCALAR_TIMESTEP(xScalarId, tempdouble, nCorners, fDsetCornerVals, tmpOut1)
  call XF_WRITE_ACTIVITY_TIMESTEP(xScalarId, nCorners, bDsetCornerActive, tmpOut2)
  if ((tmpOut1 < 0) .OR. (tmpOut2 < 0)) then
    call XF_CLOSE_GROUP(xScalarId, error)
    call XF_CLOSE_GROUP(xDatasetsId, error)
    call XF_CLOSE_GROUP(xGridId, error)
    call XF_CLOSE_FILE(xFileId, error)
    error = -1
    return
  endif

  ! Write Coordinate file - for GridCurv2D, we will set the coordinate system
  !   to be CPP, with CPP Latitude and CPP Longitude settings written 
  !   to the file.
  call XF_CREATE_COORDINATE_GROUP(xFileId, xCoordId, status)
  if (status < 0) then
    call XF_CLOSE_GROUP(xScalarId, error)
    call XF_CLOSE_GROUP(xDatasetsId, error)
    call XF_CLOSE_GROUP(xGridId, error)
    call XF_CLOSE_FILE(xFileId, error)
    error = status
    return
  endif

  dCppLat = 56.0;  ! Made-up value
  dCppLon = 23.0;  ! Made-up value

  call XF_SET_HORIZ_DATUM(xCoordId, HORIZ_DATUM_CPP, error)
  call XF_SET_HORIZ_UNITS(xCoordId, COORD_UNITS_METERS, error)

  call XF_SET_VERT_DATUM(xCoordId, VERT_DATUM_LOCAL, error)
  call XF_SET_VERT_UNITS(xCoordId, COORD_UNITS_METERS, error)

    ! write additional information
  call XF_SET_CPP_LAT(xCoordId, dCppLat, error)
  call XF_SET_CPP_LON(xCoordId, dCppLon, error)

  call XF_CLOSE_GROUP(xCoordId, error)
  xCoordId = 0

  ! release memory
call XF_CLOSE_GROUP(xScalarId, error)
call XF_CLOSE_GROUP(xDatasetsId, error)
call XF_CLOSE_GROUP(xGridId, error)
call XF_CLOSE_FILE(xFileId, error)
return

END SUBROUTINE TG_WRITE_TEST_GRID_CURV_2D

!****************************
!------------------------------------------------------------------------------
! FUNCTION   tgWriteTestGridCart3D
! PURPOSE    Write a file that contains data for a 2D Cartesian Grid
! NOTES      A picture of the grid is in the file (TestGridCart2D.gif)
!            returns TRUE on success and FALSE on failure
!------------------------------------------------------------------------------
SUBROUTINE TG_WRITE_TEST_GRID_CART_3D(Filename, Compression, error)
CHARACTER(LEN=*), INTENT(IN) :: Filename
INTEGER, INTENT(IN) :: Compression
INTEGER, INTENT(OUT) :: error
INTEGER        nDimensions
INTEGER        nCompOrigin, nUDir
INTEGER        nCellsI, nCellsJ, nCellsK
INTEGER        nGridType
REAL(DOUBLE)   dOriginX, dOriginY, dOriginZ
INTEGER        nOrientation
REAL(DOUBLE)   dBearing, dDip, dRoll
REAL(DOUBLE)   PlanesI(5), PlanesJ(5), PlanesK(3)
INTEGER        i, j, status, iSpcZone
INTEGER xFileId, xGridId, xPropId, xCoordId
INTEGER        nCells
INTEGER        Active(75)
INTEGER        tmpOut1, tmpOut2, tmpOut3

nDimensions = 3
nCompOrigin = 8
nUDir = -2
nCellsI = 5
nCellsJ = 5
nCellsK = 3
nGridType = GRID_TYPE_CARTESIAN
dOriginX = 10.0
dOriginY = 10.0
dOriginZ = 0.0
nOrientation = ORIENTATION_RIGHT_HAND
dBearing = 45.0
dDip = 0.0
dRoll = 0.0
xFileId = NONE
xGridId = NONE
xPropId = NONE
nCells = nCellsI * nCellsJ * nCellsK

  ! Fill in the grid plane data with a constant size of 30
do i = 1, nCellsI
  PlanesI(i) = i*30.0
enddo
do j = 1, nCellsJ
  PlanesJ(j) = j*30.0
enddo
do j = 1,  nCellsK
  PlanesK(j) = j*30.0
enddo

  ! fill in the activity array
  ! default array to active
do i = 1, nCells
  Active(i) = 1
enddo

  ! two cells are inactive (identified by array index)
  ! i = 0, j = 0, k = 0  and i = 4, j = 4, k = 0
Active(1) = 0
Active(4*nCellsJ*nCellsK+4*nCellsK+1) = 0

  ! create the file
call XF_CREATE_FILE(Filename, .TRUE., xFileId, tmpOut1)  
if (tmpOut1 < 0) then
  error = -1
  return
endif

  ! create the group to store the grid
call XF_CREATE_GROUP_FOR_GRID(xFileId, GRID_CART3D_GROUP_NAME, xGridId, tmpOut1)
if (tmpOut1 < 0) then
  call XF_CLOSE_FILE(xFileId, error)
  error = -1
  return
endif

  ! Write the grid information to the file
call XF_SET_GRID_TYPE(xGridId, nGridType, tmpOut1)
call XF_SET_NUMBER_OF_DIMENSIONS(xGridId, nDimensions, tmpOut2)
if ((tmpOut1 < 0) .OR. (tmpOut2 < 0)) then
  call XF_CLOSE_GROUP(xGridId, error)
  call XF_CLOSE_FILE(xFileId, error)
  error = -1
  return
endif

  ! set origin and orientation
call XF_SET_ORIGIN(xGridId, dOriginX, dOriginY, dOriginZ, tmpOut1)
call XF_SET_ORIENTATION(xGridId, nOrientation, tmpOut2)
if ((tmpOut1 < 0) .OR. (tmpOut2 < 0)) then
  call XF_CLOSE_GROUP(xGridId, error)
  call XF_CLOSE_FILE(xFileId, error)
  error = -1
  return
endif
 
  ! Set bearing, dip and roll
call XF_SET_BEARING(xGridId, dBearing, tmpOut1)
call XF_SET_DIP(xGridId, dDip, tmpOut2)
call XF_SET_ROLL(xGridId, dRoll, tmpOut3)
if ((tmpOut1 < 0) .OR. (tmpOut2 < 0) .OR. (tmpOut3 < 0)) then
  call XF_CLOSE_GROUP(xGridId, error)
  call XF_CLOSE_FILE(xFileId, error)
  error = -1
  return
endif

 ! Set computational origin
call XF_SET_COMPUTATIONAL_ORIGIN(xGridId, nCompOrigin, tmpOut1)
if (tmpOut1 < 0) then
  call XF_CLOSE_GROUP(xGridId, error)
  call XF_CLOSE_FILE(xFileId, error)
  error = -1
  return
endif

  ! Set u direction
call XF_SET_U_DIRECTION(xGridId, nUDir, tmpOut1)
if (tmpOut1 < 0) then
  call XF_CLOSE_GROUP(xGridId, error)
  call XF_CLOSE_FILE(xFileId, error)
  error = -1
  return
endif

  ! Write the grid geometry to the file
  ! Set the number of cells in each direction
call XF_SET_NUMBER_CELLS_IN_I(xGridId, nCellsI, tmpOut1)
call XF_SET_NUMBER_CELLS_IN_J(xGridId, nCellsJ, tmpOut2)
call XF_SET_NUMBER_CELLS_IN_K(xGridId, nCellsK, tmpOut3)  
if ((tmpOut1 < 0) .OR. (tmpOut2 < 0) .OR. (tmpOut3 < 0)) then
  call XF_CLOSE_GROUP(xGridId, error)
  call XF_CLOSE_FILE(xFileId, error)
  error = -1
  return
endif

  ! Set the grid plane locations
call XF_SET_GRID_COORDS_I(xGridId, nCellsI, PlanesI, tmpOut1)
call XF_SET_GRID_COORDS_J(xGridId, nCellsJ, PlanesJ, tmpOut2)
call XF_SET_GRID_COORDS_K(xGridId, nCellsK, PlanesK, tmpOut3)  

if ((tmpOut1 < 0) .OR. (tmpOut2 < 0) .OR. (tmpOut3 < 0)) then
  call XF_CLOSE_GROUP(xGridId, error)
  call XF_CLOSE_FILE(xFileId, error)
  error = -1
  return
endif

  ! Write the activity array
call XF_CREATE_GRID_CELL_PROP_GRP(xGridId, xPropId, tmpOut1)  
if (xPropId < 0) then
  call XF_CLOSE_GROUP(xGridId, error)
  call XF_CLOSE_FILE(xFileId, error)
  error = -1
  return
endif

call XF_WRITE_PROPERTY_INT(xPropId, PROP_ACTIVITY, nCells, Active, &
                           Compression, tmpOut1)

if (tmpOut1 < 0) then
  call XF_CLOSE_GROUP(xPropId, error)
  call XF_CLOSE_GROUP(xGridId, error)
  call XF_CLOSE_FILE(xFileId, error)
  error = -1
  return
endif
  
call XF_CLOSE_GROUP(xPropId, error)

  ! Write Coordinate file - for GridCart3D, we will set the coordinate system
  !   to be State Plane NAD27.
  call XF_CREATE_COORDINATE_GROUP(xFileId, xCoordId, status)
  if (status < 0) then
    call XF_CLOSE_GROUP(xGridId, error)
    call XF_CLOSE_FILE(xFileId, error)
    error = status
  endif

  iSpcZone = 3601; ! Oregon North

  call XF_SET_HORIZ_DATUM(xCoordId, HORIZ_DATUM_STATE_PLANE_NAD27, error)
  call XF_SET_HORIZ_UNITS(xCoordId, COORD_UNITS_US_FEET, error)

  call XF_SET_VERT_DATUM(xCoordId, VERT_DATUM_LOCAL, error)
  call XF_SET_VERT_UNITS(xCoordId, COORD_UNITS_US_FEET, error)

    ! write additional information
  call XF_SET_SPC_ZONE(xCoordId, iSpcZone, error)

  call XF_CLOSE_GROUP(xCoordId, error)
  xCoordId = 0

  ! release memory
call XF_CLOSE_GROUP(xGridId, error)
call XF_CLOSE_FILE(xFileId, error)
return

END SUBROUTINE TG_WRITE_TEST_GRID_CART_3D

END MODULE TestGrid

TestGrid.f90 tests grids

MODULE TestMesh

USE TestDatasets
USE Xmdf
USE ErrorDefinitions
USE XmdfDefs

CHARACTER(LEN=*), PARAMETER :: MESH_A_GROUP_NAME = 'MeshA Group'
CHARACTER(LEN=*), PARAMETER :: MESH_B_GROUP_NAME = 'MeshB Group'


CONTAINS

!****************************
! ---------------------------------------------------------------------------
! FUNCTION  TM_READ_MESH
! PURPOSE   Read a mesh file and write a text output file
! NOTES     
! ---------------------------------------------------------------------------
SUBROUTINE TM_READ_MESH (xGroupId, a_OutFile, error)
INTEGER, INTENT(IN) :: xGroupId
INTEGER, INTENT(IN)        :: a_OutFile
INTEGER, INTENT(OUT)       :: error
INTEGER           nElems, nNodes, nNodesPerElem, nElemType, nNodeId
LOGICAL*2           ElementsOneType
INTEGER           status, i, j
INTEGER           StrType, UIntType, IntType, DblType, FloatType
INTEGER    xPropGroupId
INTEGER, ALLOCATABLE      :: ElemTypes(:), NodesInElem(:)
REAL(DOUBLE), ALLOCATABLE :: XNodeLocs(:), YNodeLocs(:), ZNodeLocs(:)

    ! Get the number of elements, nodes, and Maximum number of nodes per element
  call XF_GET_NUMBER_OF_ELEMENTS (xGroupId, nElems, status)
  if (status >= 0) then
    call XF_GET_NUMBER_OF_NODES (xGroupId, nNodes, status)
    if (status >= 0) then
      call XF_GET_MAX_NODES_IN_ELEM (xGroupId, nNodesPerElem, status)
    endif
  endif

  if (status < 0) then
    error = -1
    return
  endif

    ! Do Element information first
  WRITE(a_OutFile,*) 'Number of Elements: ', nElems

    ! Element types
  call XF_ARE_ALL_ELEMS_SAME_TYPE (xGroupId, ElementsOneType, status)
  if (status < 0) then
    error = -1
    return
  endif

  if (ElementsOneType) then
    call XF_READ_ELEM_TYPES_SINGLE_VALUE (xGroupId, nElemType, status)
    WRITE(a_OutFile,*) 'All elements are type ', nElemType
  else
    allocate (ElemTypes(nElems))
    call XF_READ_ELEM_TYPES (xGroupId, nElems, ElemTypes, status)
    if (status < 0) then
      error = -1
    return
    endif
    WRITE(a_OutFile,*) 'Element Types:'
    do i=1, nElems
      WRITE(a_OutFile,*) 'Elem ', i, ', ', 'Type ', ElemTypes(i)
    enddo
    deallocate (ElemTypes)
  endif

    ! Nodes in each element
  allocate (NodesInElem(nElems*nNodesPerElem))
  call XF_READ_ELEM_NODE_IDS (xGroupId, nElems, nNodesPerElem, NodesInElem, error)
  if (error < 0) then
    WRITE (*,*) 'Error reading mesh'
    return
  endif

  do i=1, nElems
    WRITE(a_OutFile,*) 'Elem: ', i, ' - '
    do j=1, nNodesPerElem
      nNodeId = NodesInElem((i-1)*nNodesPerElem + j)
    if (nNodeId > 0) then  ! -1 is for unused array locations
        WRITE(a_OutFile,*) nNodeId, ' '
      endif
    enddo
    WRITE (a_OutFile,*) ''
  enddo

  if (allocated(NodesInElem)) deallocate (NodesInElem)

    ! NodeLocations
  allocate (XNodeLocs(nNodes))
  allocate (YNodeLocs(nNodes))
  allocate (ZNodeLocs(nNodes))

  call XF_READ_X_NODE_LOCATIONS (xGroupId, nNodes, XNodeLocs, status)
  if (status >= 0) then
    call XF_READ_Y_NODE_LOCATIONS (xGroupId, nNodes, YNodeLocs, status)
    if (status >= 0) then
      call XF_READ_Z_NODE_LOCATIONS (xGroupId, nNodes, ZNodeLocs, status)
    endif
  endif

  WRITE(a_OutFile,*) 'Node Locations:'
  do i=1, nNodes
    WRITE(a_OutFile,*) 'Node: ', i, ' Location: ', XNodeLocs(i), ' ', &
                        YNodeLocs(i), ' ', ZNodeLocs(i)
  enddo

    deallocate (XNodeLocs)
    deallocate (YNodeLocs)
    deallocate (ZNodeLocs)

    ! Open Properties Group
  call XF_OPEN_GROUP(xGroupId, 'PROPERTIES', xPropGroupId, status)
  if (status < 0) then
    WRITE(a_OutFile,*) ''
    WRITE(a_OutFile,*) 'Properties Group not found'
    WRITE(a_OutFile,*) ''
    error = -1
    return
  endif

  call XF_GET_PROPERTY_TYPE(xPropGroupId, 'String', StrType, status)
  call XF_GET_PROPERTY_TYPE(xPropGroupId, 'UInt', UIntType, status)
  call XF_GET_PROPERTY_TYPE(xPropGroupId, 'Int', IntType, status)
  call XF_GET_PROPERTY_TYPE(xPropGroupId, 'Double', DblType, status)
  call XF_GET_PROPERTY_TYPE(xPropGroupId, 'Float', FloatType, status)

    ! Property Types:
  WRITE(a_OutFile,*) ''
  if (StrType == XF_TYPE_STRING) then
    WRITE(a_OutFile,*) 'String Property Type Read Correctly'
  else
    WRITE(a_OutFile,*) 'Error in Getting String Property Type'
  endif
  if (UIntType == XF_TYPE_UINT) then
    WRITE(a_OutFile,*) 'Unsigned Integer Property Type Read Correctly'
  else
    WRITE(a_OutFile,*) 'Error in Getting Unsigned Integer Property Type'
  endif
  if (IntType == XF_TYPE_INT) then
    WRITE(a_OutFile,*) 'Integer Property Type Read Correctly'
  else
    WRITE(a_OutFile,*) 'Error in Getting Integer Property Type'
  endif
  if (DblType == XF_TYPE_DOUBLE) then
    WRITE(a_OutFile,*) 'Double Property Type Read Correctly'
  else
    WRITE(a_OutFile,*) 'Error in Getting Double Property Type'
  endif
  if (FloatType == XF_TYPE_FLOAT) then
    WRITE(a_OutFile,*) 'Float Property Type Read Correctly'
  else
    WRITE(a_OutFile,*) 'Error in Getting Float Property Type'
  endif
  WRITE(a_OutFile,*) ''

  error = 0
  return

END SUBROUTINE

!****************************
!------------------------------------------------------------------------------
! FUNCTION   TM_WRITE_TEST_MESH_A
! PURPOSE    Write a file that contains data for an all tri mesh
! NOTES      A picture of the mesh is in the file (TestMeshA.gif)
!            error equals TRUE on success and FALSE on failure
!------------------------------------------------------------------------------
SUBROUTINE TM_WRITE_TEST_MESH_A (Filename, Compression, error)
CHARACTER(LEN=*), INTENT(IN) :: Filename
INTEGER, INTENT(IN)          :: Compression
INTEGER, INTENT(OUT)         :: error
INTEGER           nElements, nNodes
INTEGER    xFileId, xMeshId, xPropGroupId, xCoordId
REAL(DOUBLE), DIMENSION(5) :: dNodeLocsX
REAL(DOUBLE), DIMENSION(5) :: dNodeLocsY
REAL(DOUBLE), DIMENSION(5) :: dNodeLocsZ
INTEGER, DIMENSION(3,3)    :: iElementNodes
INTEGER           status, ElemType, propint(1), iEllipse
CHARACTER(LEN=BIG_STRING_SIZE)  propstring
INTEGER           propuint(1)
REAL(DOUBLE)      propdouble(1), dMajorR, dMinorR
REAL              propfloat

    ! set Element type to EL_TYPE_TRI_LINEAR
  ElemType = EL_TYPE_TRI_LINEAR

  nElements = 3
  nNodes = 5
  xFileId = NONE
  xMeshId = NONE

    ! Setup the arrays for the mesh data
    ! nodes
  dNodeLocsX(1) = 0.0
  dNodeLocsX(2) = 5.0
  dNodeLocsX(3) = 0.0
  dNodeLocsX(4) = 5.0
  dNodeLocsX(5) = 7.5

  dNodeLocsY(1) = 5.0
  dNodeLocsY(2) = 5.0
  dNodeLocsY(3) = 0.0
  dNodeLocsY(4) = 0.0
  dNodeLocsY(5) = 2.5

  dNodeLocsZ(1) = 0.0
  dNodeLocsZ(2) = 0.0
  dNodeLocsZ(3) = 0.0
  dNodeLocsZ(4) = 0.0
  dNodeLocsZ(5) = 0.0

    ! nodes for each element must be counter-clockwize
  iElementNodes(1,1) = 1
  iElementNodes(2,1) = 3
  iElementNodes(3,1) = 2

  iElementNodes(1,2) = 2
  iElementNodes(2,2) = 3
  iElementNodes(3,2) = 4

  iElementNodes(1,3) = 5
  iElementNodes(2,3) = 2
  iElementNodes(3,3) = 4

    ! create the file
  call XF_CREATE_FILE (Filename, .TRUE., xFileId, status)
  if (status < 0) then
      ! close the resources
    call XF_CLOSE_FILE (xFileId, error)
    error = -1
    return
  endif

    ! create the group to store the mesh
  call XF_CREATE_GROUP_FOR_MESH (xFileId, MESH_A_GROUP_NAME, xMeshId, status)
  if (status < 0) then
      ! close the resources
    call XF_CLOSE_GROUP (xMeshId, error)
    call XF_CLOSE_FILE (xFileId, error)
    error = -1
    return
  endif

    ! Element types - all are linear triangles
  call XF_SET_ALL_ELEMS_SAME_TYPE (xMeshId, ElemType, status)
  if (status < 0) then
      ! close the resources
    call XF_CLOSE_GROUP (xMeshId, error)
    call XF_CLOSE_FILE (xFileId, error)
    error = -1
    return
  endif

    ! node information
  call XF_SET_NUMBER_OF_NODES (xMeshId, nNodes, status)
  if (status < 0) then
      ! close the resources
    call XF_CLOSE_GROUP (xMeshId, error)
    call XF_CLOSE_FILE (xFileId, error)
    error = -1
    return
  endif

  call XF_WRITE_X_NODE_LOCATIONS (xMeshId, nNodes, dNodeLocsX, Compression, status)
  if (status < 0) then
      ! close the resources
    call XF_CLOSE_GROUP (xMeshId, error)
    call XF_CLOSE_FILE (xFileId, error)
    error = -1
    return
  endif

  call  XF_WRITE_Y_NODE_LOCATIONS (xMeshId, nNodes, dNodeLocsY, status)
  if (status < 0) then
      ! close the resources
    call XF_CLOSE_GROUP (xMeshId, error)
    call XF_CLOSE_FILE (xFileId, error)
    error = -1
    return
  endif

  call XF_WRITE_Z_NODE_LOCATIONS (xMeshId, nNodes, dNodeLocsZ, status)
  if (status < 0) then
      ! close the resources
    call XF_CLOSE_GROUP (xMeshId, error)
    call XF_CLOSE_FILE (xFileId, error)
    error = -1
    return
  endif

    ! element information
  call XF_SET_NUMBER_OF_ELEMENTS (xMeshId, nElements, status)
  if (status < 0) then
      ! close the resources
    call XF_CLOSE_GROUP (xMeshId, error)
    call XF_CLOSE_FILE (xFileId, error)
    error = -1
    return
  endif
 
    ! Write the node array ids for the nodes in each element
  call XF_WRITE_ELEM_NODE_IDS (xMeshId, nElements, 3, iElementNodes, &
                               Compression, status)
  if (status < 0) then
      ! close the resources
    call XF_CLOSE_GROUP (xMeshId, error)
    call XF_CLOSE_FILE (xFileId, error)
    error = -1
    return
  endif

    ! Write the Property File
  call XF_CREATE_MESH_PROPERTY_GROUP(xMeshId, xPropGroupId, status)
  if (status < 0) then
    call XF_CLOSE_GROUP(xMeshId, error)
    call XF_CLOSE_FILE(xFileId, error)
    error = -1
    return
  endif

  propstring = 'Property String'
  propuint = 5
  propint = -5
  propdouble = 5.6789012345d0
  propfloat = 5.6789

  call XF_WRITE_PROPERTY_STRING(xPropGroupId, 'String', propstring, status)
  call XF_WRITE_PROPERTY_UINT(xPropGroupId, 'UInt', 1, propuint, NONE, status)
  call XF_WRITE_PROPERTY_INT(xPropGroupId, 'Int', 1, propint, NONE, status)
  call XF_WRITE_PROPERTY_DOUBLE(xPropGroupId, 'Double', 1, &
                                propdouble, NONE, status)
  call XF_WRITE_PROPERTY_FLOAT(xPropGroupId, 'Float', 1, &
                                propfloat, NONE, status)

  ! Write Coordinate file - for MeshA, we will set the coordinate system to be
  !   Geogrpahic, with Latitude, Longitude, and user-defined ellipsoid settings
  !   written to the file.
  call XF_CREATE_COORDINATE_GROUP(xFileId, xCoordId, status)
  if (status < 0) then
    call XF_CLOSE_GROUP (xPropGroupId, error)
    call XF_CLOSE_GROUP (xMeshId, error)
    call XF_CLOSE_FILE (xFileId, error)
    error = status
    return
  endif

    ! set coordinate values
  iEllipse = 32    ! User defined
  dMajorR = 45.0   ! Made up
  dMinorR = 32.0   ! Made up

  call XF_SET_HORIZ_DATUM(xCoordId, HORIZ_DATUM_GEOGRAPHIC, error)
  call XF_SET_HORIZ_UNITS(xCoordId, COORD_UNITS_US_FEET, error)

  call XF_SET_VERT_DATUM(xCoordId, VERT_DATUM_LOCAL, error)
  call XF_SET_VERT_UNITS(xCoordId, COORD_UNITS_US_FEET, error)

    ! write additional information
  call XF_SET_ELLIPSE(xCoordId, iEllipse, error)
  call XF_SET_LAT(xCoordId, LATITUDE_NORTH, error)
  call XF_SET_LON(xCoordId, LONGITUDE_EAST, error)
  call XF_SET_MAJOR_R(xCoordId, dMajorR, error)
  call XF_SET_MINOR_R(xCoordId, dMinorR, error)

  call XF_CLOSE_GROUP(xCoordId, error)
  xCoordId = 0

    ! close the resources
  call XF_CLOSE_GROUP (xPropGroupId, error)
  call XF_CLOSE_GROUP (xMeshId, error)
  call XF_CLOSE_FILE (xFileId, error)

  return

END SUBROUTINE

!****************************
!------------------------------------------------------------------------------
! FUNCTION   TM_WRITE_TEST_MESH_B
! PURPOSE    Write a file that contains data for an mixed quad/tri linear mesh
! NOTES      A picture of the mesh is in the file (TestMeshB.gif)
!            error equals TRUE on success and FALSE on failure
!------------------------------------------------------------------------------
SUBROUTINE TM_WRITE_TEST_MESH_B (Filename, Compression, error)
CHARACTER(LEN=*), INTENT(IN) :: Filename
INTEGER, INTENT(IN)          :: Compression
INTEGER, INTENT(OUT)         :: error
INTEGER           nElements, nNodes, nMaxNodePerElem
INTEGER    xFileId, xMeshId, xPropGroupId, xCoordId
REAL(DOUBLE), DIMENSION(5)   :: dNodeLocsX
REAL(DOUBLE), DIMENSION(5)   :: dNodeLocsY
REAL(DOUBLE), DIMENSION(5)   :: dNodeLocsZ
INTEGER, DIMENSION(4,2)      :: iElementNodes
INTEGER, DIMENSION(2)        :: iElementTypes
INTEGER           status, propint(1), iEllipse
CHARACTER(LEN=BIG_STRING_SIZE)  propstring
INTEGER           propuint(1)
REAL(DOUBLE)      propdouble(1)
REAL              propfloat

  nElements = 2
  nNodes = 5
  nMaxNodePerElem = 4
  xFileId = NONE
  xMeshId = NONE

    ! Setup the arrays for the mesh data
    ! nodes
  dNodeLocsX(1) = 0.0
  dNodeLocsX(2) = 5.0
  dNodeLocsX(3) = 0.0
  dNodeLocsX(4) = 5.0
  dNodeLocsX(5) = 7.5

  dNodeLocsY(1) = 5.0
  dNodeLocsY(2) = 5.0
  dNodeLocsY(3) = 0.0
  dNodeLocsY(4) = 0.0
  dNodeLocsY(5) = 2.5

  dNodeLocsZ(1) = 0.0
  dNodeLocsZ(2) = 0.0
  dNodeLocsZ(3) = 0.0
  dNodeLocsZ(4) = 0.0
  dNodeLocsZ(5) = 0.0
  
    ! nodes for each element must be counter-clockwize
  iElementNodes(1,1) = 1
  iElementNodes(2,1) = 3
  iElementNodes(3,1) = 4
  iElementNodes(4,1) = 2

  iElementNodes(1,2) = 2
  iElementNodes(2,2) = 4
  iElementNodes(3,2) = 5
  iElementNodes(4,2) = NONE;

  iElementTypes(1) = EL_TYPE_QUADRILATERAL_LINEAR;
  iElementTypes(2) = EL_TYPE_TRI_LINEAR;

    ! create the file
  call XF_CREATE_FILE (Filename, .TRUE., xFileId, status)
  if (status < 0) then
      ! close the resources
    call XF_CLOSE_FILE  (xFileId, error)
    error = -1
    return
  endif

    ! create the group to store the mesh
  call XF_CREATE_GROUP_FOR_MESH (xFileId, MESH_B_GROUP_NAME, xMeshId, status)
  if (status < 0) then
      ! close the resources
    call XF_CLOSE_GROUP (xMeshId, error)
    call XF_CLOSE_FILE  (xFileId, error)
    error = -1
    return
  endif

    ! node information
  call XF_SET_NUMBER_OF_NODES (xMeshId, nNodes, status)
  if (status < 0) then
      ! close the resources
    call XF_CLOSE_GROUP (xMeshId, error)
    call XF_CLOSE_FILE  (xFileId, error)
    error = -1
    return
  endif
  
  call XF_WRITE_X_NODE_LOCATIONS (xMeshId, nNodes, dNodeLocsX, &
                                  Compression, status)
  if (status < 0) then
      ! close the resources
    call XF_CLOSE_GROUP (xMeshId, error)
    call XF_CLOSE_FILE  (xFileId, error)
    error = -1
    return
  endif

  call XF_WRITE_Y_NODE_LOCATIONS (xMeshId, nNodes, dNodeLocsY, status)
  if (status < 0) then
      ! close the resources
    call XF_CLOSE_GROUP (xMeshId, error)
    call XF_CLOSE_FILE  (xFileId, error)
    error = -1
    return
  endif

  call XF_WRITE_Z_NODE_LOCATIONS (xMeshId, nNodes, dNodeLocsZ, status)
  if (status < 0) then
      ! close the resources
    call XF_CLOSE_GROUP (xMeshId, error)
    call XF_CLOSE_FILE  (xFileId, error)
    error = -1
    return
  endif
 
    ! element information
  call XF_SET_NUMBER_OF_ELEMENTS (xMeshId, nElements, status)
  if (status < 0) then
      ! close the resources
    call XF_CLOSE_GROUP (xMeshId, error)
    call XF_CLOSE_FILE  (xFileId, error)
    error = -1
    return
  endif

    ! Element types
  call XF_WRITE_ELEM_TYPES (xMeshId, nElements, iElementTypes, &
                            Compression, status)
  if (status < 0) then
      ! close the resources
    call XF_CLOSE_GROUP (xMeshId, error)
    call XF_CLOSE_FILE  (xFileId, error)
    error = -1
    return
  endif
 
    ! Write the node array ids for the nodes in each element
  call XF_WRITE_ELEM_NODE_IDS (xMeshId, nElements, nMaxNodePerElem, &
                               iElementNodes, Compression, status)
  if (status < 0) then
      ! close the resources
    call XF_CLOSE_GROUP (xMeshId, error)
    call XF_CLOSE_FILE  (xFileId, error)
    error = -1
    return
  endif

    ! Write the Property File
  call XF_CREATE_MESH_PROPERTY_GROUP(xMeshId, xPropGroupId, status)
  if (status < 0) then
    call XF_CLOSE_GROUP(xMeshId, error)
    call XF_CLOSE_GROUP(xFileId, error)
    error = -1
    return
  endif

  propstring = 'String Property'
  propuint = 2
  propint = -2
  propdouble = 2.3456789012d0
  propfloat = 2.3456

  call XF_WRITE_PROPERTY_STRING(xPropGroupId, 'String', propstring, status)
  call XF_WRITE_PROPERTY_UINT(xPropGroupId, 'UInt', 1, propuint, NONE, status)
  call XF_WRITE_PROPERTY_INT(xPropGroupId, 'Int', 1, propint, NONE, status)
  call XF_WRITE_PROPERTY_DOUBLE(xPropGroupId, 'Double', 1, &
                                propdouble, NONE, status)
  call XF_WRITE_PROPERTY_FLOAT(xPropGroupId, 'Float', 1, &
                               propfloat, NONE, status)

  ! Write Coordinate file - for MeshB, we will set the coordinate system to be
  !   Geogrpahic, with Latitude, Longitude, and standard ellipsoid settings
  !   written to the file.
  call XF_CREATE_COORDINATE_GROUP(xFileId, xCoordId, status)
  if (status < 0) then
    call XF_CLOSE_GROUP (xPropGroupId, error)
    call XF_CLOSE_GROUP (xMeshId, error)
    call XF_CLOSE_FILE  (xFileId, error)
    error = status
    return
  endif

    ! set coordinate values
  iEllipse = 21   ! International 1924

  call XF_SET_HORIZ_DATUM(xCoordId, HORIZ_DATUM_GEOGRAPHIC, error)
  call XF_SET_HORIZ_UNITS(xCoordId, COORD_UNITS_METERS, error)

  call XF_SET_VERT_DATUM(xCoordId, VERT_DATUM_NGVD_88, error)
  call XF_SET_VERT_UNITS(xCoordId, COORD_UNITS_METERS, error)

    ! write additional information
  call XF_SET_ELLIPSE(xCoordId, iEllipse, error)
  call XF_SET_LAT(xCoordId, LATITUDE_SOUTH, error)
  call XF_SET_LON(xCoordId, LONGITUDE_WEST, error)

  call XF_CLOSE_GROUP(xCoordId, error)
  xCoordId = 0

    ! close the resources
  call XF_CLOSE_GROUP (xPropGroupId, error)
  call XF_CLOSE_GROUP (xMeshId, error)
  call XF_CLOSE_FILE  (xFileId, error)

  return

END SUBROUTINE

END MODULE TestMesh

TestMesh.f90 tests meshes

MODULE TestTimestep

USE Xmdf

CHARACTER(LEN=*), PARAMETER :: DATASETS_LOCATION = 'Datasets'
CHARACTER(LEN=*), PARAMETER :: SCALAR_A_LOCATION = 'Scalars/ScalarA'
CHARACTER(LEN=*), PARAMETER :: SCALAR_B_LOCATION = 'Scalars/ScalarB'
CHARACTER(LEN=*), PARAMETER :: VECTOR2D_A_LOCATION = 'Vectors/Vector2D_A'
CHARACTER(LEN=*), PARAMETER :: VECTOR2D_B_LOCATION = 'Vectors/Vector2D_B'

CONTAINS

! --------------------------------------------------------------------------
! FUNCTION ttiTestNumTimes
! PURPOSE  Change the NumTimes to truncate timesteps
! NOTES    
! --------------------------------------------------------------------------
RECURSIVE SUBROUTINE TTI_Test_Num_Times (a_DatasetId, a_Itimestep, error)
INTEGER, INTENT(IN) :: a_DatasetId
INTEGER, INTENT(IN)        :: a_Itimestep
INTEGER, INTENT(OUT)       :: error
INTEGER                    NumTimes
INTEGER                    Itimestep

error = 1

! Fortran loops are 1 to 3 but C is 0 to 2, etc
Itimestep = a_Itimestep - 1

  ! truncate just written timestep and test error conditions
if (1 == Itimestep .OR. 3 == Itimestep .OR. 5 == Itimestep) then
      ! Test setting NumTimes after end of dataset
      call XF_SET_DATASET_NUM_TIMES( a_DatasetId, Itimestep + 2, error )
      if (error >= 0) then
         WRITE(*,*) 'ERROR1: XF_SET_DATASET_NUM_TIMES must return ERROR.'
      endif

      if (1 == Itimestep) then
         Itimestep = 1;
      endif
      if (3 == Itimestep) then
         Itimestep = 2;
      endif
      if (5 == Itimestep) then
         Itimestep = 3;
      endif

      ! Write actual NumTimes
      call XF_SET_DATASET_NUM_TIMES( a_DatasetId, Itimestep, error )
      if (error < 0) then
         WRITE(*,*) 'ERROR2: xfSetDatasetNumTimes must NOT return error.'
      endif

      ! Test setting NumTimes after end step.
      call XF_SET_DATASET_NUM_TIMES( a_DatasetId, Itimestep + 1, error )
      if (error >= 0) then
         WRITE(*,*) 'ERROR3: xfSetDatasetNumTimes must return ERROR.'
      endif

      ! Test reading NumTimes
      call XF_GET_DATASET_NUM_TIMES( a_DatasetId, NumTimes, error )
      if (error < 0) then
         WRITE(*,*) 'ERROR4: xfSetDatasetNumTimes must NOT return error.'
      endif
      if (NumTimes .NE. Itimestep) then
         WRITE(*,*) 'ERROR5: xfGetDatasetNumTimes must return CORRECT NumTimes.'
      endif
endif

return

END SUBROUTINE
!ttiTestNumTimes
! --------------------------------------------------------------------------
! FUNCTION ttReadDatasets
! PURPOSE  Read a dataset group from an XMDF file and output information to
!          to a text file
! NOTES    
! --------------------------------------------------------------------------
RECURSIVE SUBROUTINE TT_READ_DATASETS (a_xGroupId, a_FileUnit, error)
INTEGER, INTENT(IN) :: a_xGroupId
INTEGER, INTENT(IN)        :: a_FileUnit
INTEGER, INTENT(OUT)       :: error
INTEGER                   nPaths, nMaxPathLength, j
CHARACTER, ALLOCATABLE, DIMENSION(:) :: Paths
CHARACTER(LEN=500)       IndividualPath
INTEGER                   nStatus, i
INTEGER            xScalarId, xVectorId, xMultiId
INTEGER                   nMultiDatasets

xScalarId = NONE
xVectorId = NONE

nMultiDatasets = 0
nPaths = 0
nMaxPathLength = 0

  ! Look for scalar datasets
call XF_GET_SCALAR_DATASETS_INFO(a_xGroupId, nPaths, nMaxPathLength, nStatus)
if (nStatus >= 0 .AND. nPaths > 0) then
  allocate(Paths(nPaths*nMaxPathLength))
  call XF_GET_SCALAR_DATASET_PATHS(a_xGroupId, nPaths, nMaxPathLength, Paths, &
                                                                         error)
endif
if (nStatus < 0) then
  error = -1
  return
endif

  ! Output number and paths to scalar datasets
WRITE(a_FileUnit,*) 'Number of Scalars ', nPaths
do i=2, nPaths
  IndividualPath = ''
  do j=1, nMaxPathLength-1
    IndividualPath(j:j) = Paths((i-1)*nMaxPathLength+j)
  enddo
  WRITE(a_FileUnit,*) 'Reading scalar: ', IndividualPath(1:nMaxPathLength-1)
  call XF_OPEN_GROUP(a_xGroupId, IndividualPath(1:nMaxPathLength-1), &
                                                           xScalarId, nStatus)
  if (nStatus < 0) then
    error = -1
  return
  endif

  call TTI_READ_SCALAR(xScalarId, a_FileUnit, nStatus)
  call XF_CLOSE_GROUP(xScalarId, error)
  if (nStatus < 0) then
    WRITE(*,*) 'Error reading scalar dataset.'
    error = -1
  return
  endif
enddo

if (allocated(Paths)) deallocate(Paths)
  ! Look for vector datasets
call XF_GET_VECTOR_DATASETS_INFO(a_xGroupId, nPaths, nMaxPathLength, nStatus)
if (nStatus >= 0 .AND. nPaths > 0) then
  allocate(Paths(nPaths*nMaxPathLength))
  call XF_GET_VECTOR_DATASET_PATHS(a_xGroupId, nPaths, nMaxPathLength, Paths, error)
endif
if (nStatus < 0) then
  error = -1
  return
endif

  ! Output number and paths to scalar datasets
WRITE(a_FileUnit,*) 'Number of Vectors ', nPaths
do i=2, nPaths
  do j=1, nMaxPathLength-1
    IndividualPath(j:j) = Paths((i-1)*nMaxPathLength+j)
  enddo
  WRITE(a_FileUnit,*) 'Reading Vector: ', &
                      IndividualPath(1:nMaxPathLength-1)
  call XF_OPEN_GROUP(a_xGroupId, IndividualPath(1:nMaxPathLength-1), &
                                                          xVectorId, nStatus)
  if (nStatus < 0) then
    error = -1
  return
  endif
  call TTI_READ_VECTOR(xVectorId, a_FileUnit, nStatus)
  call XF_CLOSE_GROUP(xVectorId, error)
  if (nStatus < 0) then
    WRITE(*,*) 'Error reading vector dataset.'
    error = -1
  return
  endif
enddo

if (allocated(Paths)) deallocate(Paths)

! find multidataset folders
call XF_GET_GRP_PTHS_SZ_MLT_DSETS(a_xGroupId, nMultiDatasets, &
                                                      nMaxPathLength, nStatus)
if (nStatus >= 0 .AND. nMultiDatasets > 0) then
  allocate(Paths(nMultiDatasets*nMaxPathLength))
  call XF_GET_ALL_GRP_PATHS_MLT_DSETS(a_xGroupId, nMultiDatasets, &
                                                 nMaxPathLength, Paths, error)
  if (nStatus < 0) then
    error = -1
    return
  endif

  ! Output number and paths to multidatasets
  WRITE(a_FileUnit,*) 'Number of Multidatasets ', nMultiDatasets
  do i=2, nMultiDatasets
    IndividualPath = ''
    do j=1, nMaxPathLength-1
      IndividualPath(j:j) = Paths((i-1)*nMaxPathLength+j)
    enddo
    WRITE(a_FileUnit,*) 'Reading multidataset: ', &
                                             IndividualPath(1:nMaxPathLength-1)
    call XF_OPEN_GROUP(a_xGroupId, IndividualPath(1:nMaxPathLength-1), &
                                                           xMultiId, nStatus)
    if (nStatus < 0) then
      error = -1
    return
    endif

    call TT_READ_DATASETS(xMultiId, a_FileUnit, nStatus)
    call XF_CLOSE_GROUP(xMultiId, error)
    if (nStatus < 0) then
      WRITE(*,*) 'Error reading multidatasets.'
      error = -1
    return
    endif
  enddo
endif
if (allocated(Paths)) deallocate(Paths)

error = 1
return

END SUBROUTINE
!ttReadDatasets
! --------------------------------------------------------------------------
! FUNCTION ttReadActivityScalarAIndex
! PURPOSE  Read all timestep values for a particular index
! NOTES    
! --------------------------------------------------------------------------
SUBROUTINE TT_READ_ACTIVITY_SCALAR_A_INDEX(a_Filename, a_Index, error)
CHARACTER(LEN=*), INTENT(IN) :: a_Filename
INTEGER, INTENT(IN)   :: a_Index
INTEGER, INTENT(OUT)  :: error
INTEGER                  status
INTEGER           xFileId, xDsetsId, xScalarAId
INTEGER                  nTimesteps, i
INTEGER, ALLOCATABLE  :: bActive(:)

xFileId = NONE
xDsetsId = NONE
xScalarAId = NONE

  ! open the file
call XF_OPEN_FILE(a_Filename, .TRUE., xFileId, status)
if (status < 0) then
  error = -1
  return
endif

  ! open the dataset group
call XF_OPEN_GROUP(xFileId, DATASETS_LOCATION, xDsetsId, status)
if (status >= 0) then
  call XF_OPEN_GROUP(xDsetsId, SCALAR_A_LOCATION, xScalarAId, status)
endif
if (status < 0) then
  error = status
  return
endif

  ! Find out the number of timesteps in the file
CALL XF_GET_DATASET_NUM_TIMES(xScalarAId, nTimesteps, status)
if (status < 0) then
  error = status
  return
endif

if (nTimesteps < 1) then
  error = -1
  return
endif

  ! Read the values for the index
allocate(bActive(nTimesteps))
call XF_READ_ACTIVE_VALS_AT_INDEX(xScalarAId, a_Index, 1, nTimesteps, &
                                       bActive, status)
  ! output the data
WRITE(*,*) ''
WRITE(*,*) 'Reading activity for scalar A slice at index: ', a_Index
do i=1, nTimesteps
  WRITE(*,*) bActive(i), ' '
enddo

deallocate(bActive)

error = status
return

END SUBROUTINE
! ttReadActivityScalarAIndex

! --------------------------------------------------------------------------
! FUNCTION ttReadScalarAIndex
! PURPOSE  Read all timestep values for a particular index
! NOTES    
! --------------------------------------------------------------------------
SUBROUTINE TT_READ_SCALAR_A_INDEX (a_Filename, a_Index, error)
CHARACTER(LEN=*), INTENT(IN) :: a_Filename
INTEGER, INTENT(IN)   :: a_Index
INTEGER, INTENT(OUT)  :: error
INTEGER              status
INTEGER       xFileId, xDsetsId, xScalarAId
INTEGER              nTimesteps, i
REAL, ALLOCATABLE :: fValues(:)

xFileId = NONE
xDsetsId = NONE
xScalarAId = NONE

  ! open the file
call XF_OPEN_FILE(a_Filename, .TRUE., xFileId, status)
if (status < 0) then
  error = -1
  return
endif

  ! open the dataset group
call XF_OPEN_GROUP(xFileId, DATASETS_LOCATION, xDsetsId, status)
if (status >= 0) then
  call XF_OPEN_GROUP(xDsetsId, SCALAR_A_LOCATION, xScalarAId, status)
endif
if (status < 0) then
  error = status
  return
endif

  ! Find out the number of timesteps in the file
call XF_GET_DATASET_NUM_TIMES(xScalarAId, nTimesteps, status)
if (status < 0) then
  error = status
  return
endif

if (nTimesteps < 1) then
  error = -1
  return
endif

  ! Read the values for the index
allocate (fValues(nTimesteps))
call XF_READ_SCALAR_VALUES_AT_INDEX(xScalarAId, a_Index, 1, nTimesteps, &
                                     fValues, status)

  ! output the data
WRITE(*,*) ''
WRITE(*,*) 'Reading scalar A slice at index: ', a_Index
do i=1, nTimesteps
  WRITE(*,*) fValues(i), ' '
enddo

deallocate(fValues)

error = status
return

END SUBROUTINE
! ttReadScalarAtIndex

! --------------------------------------------------------------------------
! FUNCTION ttWriteScalarA
! PURPOSE  Write scalar Dataset to an HDF5 File
! NOTES    This tests dynamic data sets, and activity
!          This dataset is dynamic concentrations (mg/L) with output times
!          in minutes.
!          Dataset is for a mesh and so nActive is the number of elements
!          which is not the same as the nValues which would be number of nodes
!          reads/writes a reference time in julian days
! --------------------------------------------------------------------------
SUBROUTINE TT_WRITE_SCALAR_A (a_Filename, a_Compression, error)
  CHARACTER(LEN=*), INTENT(IN) :: a_Filename
  INTEGER, INTENT(IN)          :: a_Compression
  INTEGER, INTENT(OUT)         :: error
  INTEGER      xFileId, xDsetsId, xScalarAId, xCoordId
  INTEGER      nValues, nTimes, nActive
  REAL(DOUBLE) dTime, dJulianReftime
  INTEGER      iTimestep, iActive, iHpgnZone
  REAL         fValues(10) ! nValues
  INTEGER*1    bActivity(10) ! activity
  INTEGER      i, status

  ! initialize the data
  nValues = 10
  nTimes = 3
  nActive = 8
  dTime = 0.0

  ! 5th item in data set is always inactive, others active
  do iActive = 1, nActive
    bActivity(iActive) = 1
  enddo 
  bActivity(6) = 0


  ! create the file
  call XF_CREATE_FILE(a_Filename, .TRUE., xFileId, error)
  if (error .LT. 0) then
      ! close the file
    call XF_CLOSE_FILE(xFileId, error)
    return
  endif

  ! create the group where we will put all the datasets 
  call XF_CREATE_GENERIC_GROUP(xFileId, DATASETS_LOCATION, xDsetsId, status)
  if (status < 0) then
    call XF_CLOSE_FILE(xFileId, error)
    error = -1
  return
  endif

  ! Create the scalar A dataset group
  call XF_CREATE_SCALAR_DATASET(xDsetsId, SCALAR_A_LOCATION, 'mg/L', &
              TS_HOURS, a_Compression, xScalarAId, status)
  if (status .LT. 0) then
      ! close the dataset
    call XF_CLOSE_GROUP(xScalarAId, error)
    call XF_CLOSE_GROUP(xDsetsId, error)
    call XF_CLOSE_FILE(xFileId, error)
    error = status
    return 
  endif

  ! Add in a reftime.  This is a julian day for:
  ! noon July 1, 2003
  dJulianReftime = 2452822.0;
  call XF_WRITE_REFTIME(xScalarAId, dJulianReftime, status)
  if (status < 0) then
    call XF_CLOSE_GROUP(xScalarAId, error)
    call XF_CLOSE_GROUP(xDsetsId, error)
    call XF_CLOSE_FILE(xFileId, error)
  endif

  ! Loop through timesteps adding them to the file
  do iTimestep = 1, nTimes
    ! We will have an 0.5 hour timestep
    dTime = iTimestep * 0.5

    fValues(1) = dTime
    do i = 2, nValues
      fValues(i) = fValues(i-1)*2.5
    end do

    ! write the dataset array values
    call XF_WRITE_SCALAR_TIMESTEP(xScalarAId, dTime, nValues, fValues, error)
    if (error .GE. 0) then
      ! write activity array
      call XF_WRITE_ACTIVITY_TIMESTEP(xScalarAId, nActive, bActivity, error)
    end if 

    call TTI_Test_Num_Times(xScalarAid, iTimestep, error)
  enddo

  ! Write Coordinate file - for ScalarA, we will set the coordinate system
  !   to be Geographic HPGN, with HPGN settings written to the file.
  call XF_CREATE_COORDINATE_GROUP(xFileId, xCoordId, status)
  if (status < 0) then
    call XF_CLOSE_GROUP(xScalarAId, error)
    call XF_CLOSE_GROUP(xDsetsId, error)
    call XF_CLOSE_FILE(xFileId, error)
    error = -1
  return
  endif

    ! set HPGN Zone for test
  iHpgnZone = 29   ! Utah
    ! Write Coordinate Information to file
  call XF_SET_HORIZ_DATUM(xCoordId, HORIZ_DATUM_GEOGRAPHIC_HPGN, error)
  call XF_SET_HORIZ_UNITS(xCoordId, COORD_UNITS_METERS, error)
  call XF_SET_VERT_DATUM(xCoordId, VERT_DATUM_LOCAL, error)
  call XF_SET_VERT_UNITS(xCoordId, COORD_UNITS_METERS, error)

    ! write additional information
  call XF_SET_HPGN_AREA(xCoordId, iHpgnZone, error)

  call XF_CLOSE_GROUP(xCoordId, error)
  xCoordId = 0;

  ! close the dataset
  call XF_CLOSE_GROUP(xScalarAId, error)
  call XF_CLOSE_GROUP(xDsetsId, error)
  call XF_CLOSE_FILE(xFileId, error)

  return
END SUBROUTINE
! ttWriteScalarA

! --------------------------------------------------------------------------
! FUNCTION TT_WRITE_SCALAR_B
! PURPOSE  Write scalar Dataset to an HDF5 File
! NOTES    This tests dynamic data sets, and activity
!          This dataset is dynamic concentrations (mg/L) with output times
!          in minutes.
!          Dataset is for a mesh and so nActive is the number of elements
!          which is not the same as the nValues which would be number of nodes
!          reads/writes a reference time in julian days
! --------------------------------------------------------------------------
SUBROUTINE TT_WRITE_SCALAR_B (a_Filename, a_Compression, a_Overwrite, error)
  CHARACTER(LEN=*), INTENT(IN) :: a_Filename
  INTEGER, INTENT(IN)          :: a_Compression
  LOGICAL, INTENT(IN)          :: a_Overwrite
  INTEGER, INTENT(OUT)         :: error
  INTEGER      xFileId, xDsetsId, xScalarBId, xCoordId
  INTEGER      nValues, nTimes, nActive
  REAL(DOUBLE) dTime, dJulianReftime
  INTEGER      iTimestep, iActive
  REAL         fValues(10) ! nValues
  INTEGER*1    bActivity(10) ! activity
  INTEGER      i, status

  ! initialize the data
  nValues = 10
  nTimes = 3
  nActive = 8
  dTime = 0.0
  i = 0

  ! 5th item in data set is always inactive, others active
  do iActive = 1, nActive
    bActivity(iActive) = 1
  enddo 
  bActivity(6) = 0

  if (a_Overwrite) then
      ! open the already-existing file
    call XF_OPEN_FILE(a_Filename, .FALSE., xFileId, status)
    if (status < 0) then
      error = -1
      return
    endif
      ! open the group where we have all the datasets
    call XF_OPEN_GROUP(xFileId, DATASETS_LOCATION, xDsetsId, status)
    if (status < 0) then
      call XF_CLOSE_FILE(xFileId, error)
      error = -1
      return
    endif
  else
      ! create the file
    call XF_CREATE_FILE(a_Filename, .TRUE., xFileId, error)
    if (error .LT. 0) then
        ! close the file
      call XF_CLOSE_FILE(xFileId, error)
      return
    endif

      ! create the group where we will put all the datasets 
    call XF_CREATE_GENERIC_GROUP(xFileId, DATASETS_LOCATION, xDsetsId, status)
    if (status < 0) then
      call XF_CLOSE_FILE(xFileId, error)
      error = -1
    return
    endif
  endif

  ! Create/Overwrite the scalar B dataset group
  call XF_CREATE_SCALAR_DATASET(xDsetsId, SCALAR_B_LOCATION, 'mg/L', &
              TS_HOURS, a_Compression, xScalarBId, status)
  if (status < 0) then
      ! close the dataset
    call XF_CLOSE_GROUP(xScalarBId, error)
    call XF_CLOSE_GROUP(xDsetsId, error)
    call XF_CLOSE_FILE(xFileId, error)
    error = status
    return 
  endif

  ! Add in a reftime.  This is a julian day for:
  ! noon July 1, 2003
  dJulianReftime = 2452822.0;
  call XF_WRITE_REFTIME(xScalarBId, dJulianReftime, status)
  if (status < 0) then
    call XF_CLOSE_GROUP(xScalarBId, error)
    call XF_CLOSE_GROUP(xDsetsId, error)
    call XF_CLOSE_FILE(xFileId, error)
  endif

  if (.NOT. a_Overwrite) then
      ! Loop through timesteps adding them to the file
    do iTimestep = 1, nTimes
        ! We will have an 0.5 hour timestep
      dTime = iTimestep * 0.5

      fValues(1) = dTime
      do i = 2, nValues
        fValues(i) = fValues(i-1)*2.5
      end do

        ! write the dataset array values
      call XF_WRITE_SCALAR_TIMESTEP(xScalarBId, dTime, nValues, fValues, error)
      if (error .GE. 0) then
          ! write activity array
        call XF_WRITE_ACTIVITY_TIMESTEP(xScalarBId, nActive, bActivity, error)
      end if
      if (error < 0) then
          call XF_CLOSE_GROUP(xScalarBId, error)
          call XF_CLOSE_GROUP(xDsetsId, error)
          call XF_CLOSE_FILE(xFileId, error)
      endif
    enddo
  else
      ! Loop through timesteps adding them to the file
    do iTimestep = 1, nTimes
        ! We will have an 1.5 hour timestep
      dTime = iTimestep * 1.5

      fValues(1) = dTime
      do i = 2, nValues
        fValues(i) = fValues(i-1)*1.5
      end do

        ! write the dataset array values
      call XF_WRITE_SCALAR_TIMESTEP(xScalarBId, dTime, nValues, fValues, error)
      if (error .GE. 0) then
          ! write activity array
        call XF_WRITE_ACTIVITY_TIMESTEP(xScalarBId, nActive, bActivity, error)
      end if
      if (error < 0) then
          call XF_CLOSE_GROUP(xScalarBId, error)
          call XF_CLOSE_GROUP(xDsetsId, error)
          call XF_CLOSE_FILE(xFileId, error)
      endif
    enddo
  endif

  if (.NOT. a_Overwrite) then
    ! Write Coordinate file - for ScalarB, we will set the coordinate system
    !   to be UTM, with UTM Zone settings written to the file.
    call XF_CREATE_COORDINATE_GROUP(xFileId, xCoordId, status)
    if (status < 0) then
      call XF_CLOSE_GROUP(xScalarBId, error)
      call XF_CLOSE_GROUP(xDsetsId, error)
      call XF_CLOSE_FILE(xFileId, error)
    error = -1
    return
    endif

     ! Write Coord Info to file
    call XF_SET_HORIZ_DATUM(xCoordId, HORIZ_DATUM_UTM, error)
    call XF_SET_HORIZ_UNITS(xCoordId, COORD_UNITS_METERS, error)

    call XF_SET_VERT_DATUM(xCoordId, VERT_DATUM_LOCAL, error)
    call XF_SET_VERT_UNITS(xCoordId, COORD_UNITS_METERS, error)

      ! write additional information - we'll use the max value for this test
    call XF_SET_UTM_ZONE(xCoordId, UTM_ZONE_MAX, error)

    call XF_CLOSE_GROUP(xCoordId, error)
    xCoordId = 0
  endif

  ! close the dataset
  call XF_CLOSE_GROUP(xScalarBId, error)
  call XF_CLOSE_GROUP(xDsetsId, error)
  call XF_CLOSE_FILE(xFileId, error)

  error = 1
  return
END SUBROUTINE
! ttWriteScalarB
!------------------------------------------------------------------------------
!  FUNCTION TT_WRITE_COORDS_TO_MULTI
!  PURPOSE  Write coordinate system to a multidataset file
!  NOTES
!------------------------------------------------------------------------------
SUBROUTINE TT_WRITE_COORDS_TO_MULTI (a_xFileId, error)
INTEGER, INTENT(IN) :: a_xFileId
INTEGER, INTENT(OUT)       :: error
INTEGER    xCoordId
INTEGER           status

  ! Write Coordinate file - for Multidatasets, we will set the coordinate system
  !   to be UTM, with UTM Zone settings written to the file.
  call XF_CREATE_COORDINATE_GROUP(a_xFileId, xCoordId, status)
  if (status < 0) then
    error = status
  return
  endif

    ! Write Coord Info to file
  call XF_SET_HORIZ_DATUM(xCoordId, HORIZ_DATUM_UTM, error)
  call XF_SET_HORIZ_UNITS(xCoordId, COORD_UNITS_METERS, error)

  call XF_SET_VERT_DATUM(xCoordId, VERT_DATUM_LOCAL, error)
  call XF_SET_VERT_UNITS(xCoordId, COORD_UNITS_METERS, error)

    ! write additional information - we'll use the max value for this test
  call XF_SET_UTM_ZONE(xCoordId, UTM_ZONE_MAX, error)

  call XF_CLOSE_GROUP(xCoordId, error)
  xCoordId = 0

  return
END SUBROUTINE

! --------------------------------------------------------------------------
! FUNCTION ttWriteScalarAToMulti
! PURPOSE  Write scalar Dataset to an HDF5 File
! NOTES    This tests dynamic data sets, and activity
!          This dataset is dynamic concentrations (mg/L) with output times
!          in minutes.
!          Dataset is for a mesh and so nActive is the number of elements
!          which is not the same as the nValues which would be number of nodes
!          reads/writes a reference time in julian days
! --------------------------------------------------------------------------
SUBROUTINE TT_WRITE_SCALAR_A_TO_MULTI (a_GroupID, status)
 ! CHARACTER(LEN=*), INTENT(IN) :: a_Filename
 ! INTEGER, INTENT(IN)          :: a_Compression
 ! INTEGER, INTENT(OUT)         :: error
  INTEGER      xFileId, xDsetsId, xScalarAId
  INTEGER      a_GroupID
  INTEGER      nValues, nTimes, nActive
  REAL(DOUBLE) dTime, dJulianReftime
  INTEGER      iTimestep, iActive
  REAL         fValues(10) ! nValues
  INTEGER*1    bActivity(10) ! activity
  INTEGER      i, status

  ! initialize the data
  nValues = 10
  nTimes  = 3
  nActive = 8
  dTime   = 0.0

  ! 5th item in data set is always inactive, others active
  do iActive = 1, nActive
    bActivity(iActive) = 1
  enddo 
  bActivity(6) = 0

  ! Create the scalar A dataset group
  call XF_CREATE_SCALAR_DATASET(a_GroupID, SCALAR_A_LOCATION, 'mg/L', &
              TS_HOURS, NONE, xScalarAId, status)
  if (status .LT. 0) then
      ! close the dataset
    call XF_CLOSE_GROUP(xScalarAId, status)
    call XF_CLOSE_GROUP(xDsetsId, status)
    call XF_CLOSE_FILE(xFileId, status)
    return 
  endif

  ! Add in a reftime.  This is a julian day for:
  ! noon July 1, 2003
  dJulianReftime = 2452822.0;
  call XF_WRITE_REFTIME(xScalarAId, dJulianReftime, status)
  if (status < 0) then
    call XF_CLOSE_GROUP(xScalarAId, status)
    call XF_CLOSE_GROUP(xDsetsId, status)
    call XF_CLOSE_FILE(xFileId, status)
  endif

  ! Loop through timesteps adding them to the file
  do iTimestep = 1, nTimes
    ! We will have an 0.5 hour timestep
    dTime = iTimestep * 0.5

    fValues(1) = dTime
    do i = 2, nValues
      fValues(i) = fValues(i-1)*2.5
    end do

    ! write the dataset array values
    call XF_WRITE_SCALAR_TIMESTEP(xScalarAId, dTime, nValues, fValues, status)
    if (status .GE. 0) then
      ! write activity array
      call XF_WRITE_ACTIVITY_TIMESTEP(xScalarAId, nActive, bActivity, status)
    end if 

    call TTI_Test_Num_Times(xScalarAId, iTimestep, status)
  enddo

  ! close the dataset
  call XF_CLOSE_GROUP(xScalarAId, status)
  !call XF_CLOSE_GROUP(a_GroupID, status)
  !call XF_CLOSE_FILE(a_FileID, status)

  return
END SUBROUTINE
! ttWriteScalarAToMulti
! --------------------------------------------------------------------------
! FUNCTION ttReadVector2DAIndex
! PURPOSE  Read all timestep values for a particular index
! NOTES    
! --------------------------------------------------------------------------
SUBROUTINE TT_READ_VECTOR2D_A_INDEX (a_Filename, a_Index, error)
CHARACTER(LEN=*), INTENT(IN) :: a_Filename
INTEGER, INTENT(IN)   :: a_Index
INTEGER, INTENT(OUT)  :: error
INTEGER           status
INTEGER    xFileId, xDsetsId, xVector2DA
INTEGER           nTimesteps, i
REAL, ALLOCATABLE     :: fValues(:)

xFileId = NONE
xDsetsId = NONE
xVector2DA = NONE

  ! open the file
call XF_OPEN_FILE(a_Filename, .TRUE., xFileId, status)
if (status < 0) then
  error = -1
  return
endif

  ! open the dataset group
call XF_OPEN_GROUP(xFileId, DATASETS_LOCATION, xDsetsId, status)
if (status >= 0) then
  call XF_OPEN_GROUP(xDsetsId, VECTOR2D_A_LOCATION, xVector2DA, status)
endif
if (status < 0) then
  error = status
  return
endif

  ! Find out the number of timesteps in the file
call XF_GET_DATASET_NUM_TIMES(xVector2DA, nTimesteps, status)
if (status < 0) then
  error = status
  return
endif

if (nTimesteps < 1) then
  error = -1
  return
endif

  ! Read the values for the index
allocate(fValues(nTimesteps*2))
call XF_READ_VECTOR_VALUES_AT_INDEX(xVector2DA, a_Index, 1, nTimesteps, 2, &
                                     fValues, status)

  ! output the data
WRITE(*,*) ''
WRITE(*,*) 'Reading vector 2D A slice at index: ', a_Index
do i=1, nTimesteps
  WRITE(*,*) fValues(i*2-1), ' ', fValues(i*2)
enddo
WRITE(*,*) ''

deallocate(fValues)

error = status
return

END SUBROUTINE
!ttReadVector2DAIndex

! --------------------------------------------------------------------------
! FUNCTION ttWriteVector2D_A
! PURPOSE  Write scalar Dataset to an HDF5 File
! NOTES    This tests dynamic data sets, and activity
!          This dataset is dynamic concentrations (mg/L) with output times
!          in minutes.
!          Dataset is for a mesh and so nActive is the number of elements
!          which is not the same as the nValues which would be number of nodes
!          reads/writes a reference time in julian days
! --------------------------------------------------------------------------
SUBROUTINE TT_WRITE_VECTOR2D_A (a_Filename, a_Compression, error)
  CHARACTER(LEN=*), INTENT(IN) :: a_Filename
  INTEGER, INTENT(IN)          :: a_Compression
  INTEGER, INTENT(OUT)         :: error
  INTEGER      xFileId, xDsetsId, xVector2D_A, xCoordId
  INTEGER      nValues, nTimes, nComponents, nActive
  REAL(DOUBLE) dTime
  INTEGER      iTimestep, iActive
  REAL, DIMENSION(2, 100) :: fValues ! nComponents, nValues
  INTEGER*1    bActivity(100) ! activity
  INTEGER      i, j, status
  INTEGER      iHpgnZone

  ! initialize the data
  nComponents = 2
  nValues = 100
  nTimes = 6
  nActive = 75
  dTime = 0.0

  ! 5th item in data set is always inactive, others active
  bActivity(1) = 0
  do iActive = 2, nActive
    if (mod(iActive-1, 3) == 0) then
      bActivity(iActive) = 0
    else
    bActivity(iActive) = 1
  endif
  enddo

  ! create the file
  call XF_CREATE_FILE(a_Filename, .TRUE., xFileId, error)
  if (error .LT. 0) then
    ! close the dataset
    call XF_CLOSE_FILE(xFileId, error)
    return
  endif

  ! create the group where we will put all the datasets 
  call XF_CREATE_GENERIC_GROUP(xFileId, DATASETS_LOCATION, xDsetsId, status)
  if (status < 0) then
    call XF_CLOSE_FILE(xFileId, error)
    error = -1
  return
  endif

  ! Create the vector dataset group
  call XF_CREATE_VECTOR_DATASET(xDsetsId, VECTOR2D_A_LOCATION, 'ft/s', &
              TS_SECONDS, a_Compression, xVector2D_A, status)
  if (status .LT. 0) then
      ! close the dataset
    call XF_CLOSE_GROUP(xVector2D_A, error)
    call XF_CLOSE_GROUP(xDsetsId, error)
    call XF_CLOSE_FILE(xFileId, error)
    error = status
    return 
  endif

  ! Loop through timesteps adding them to the file
  do iTimestep = 1, nTimes
    ! We will have an 0.5 hour timestep
    dTime = iTimestep * 0.5

    do i = 1, nValues
      do j = 1, nComponents
        fValues(j,i) = ((i-1)*nComponents + j)*dTime
      end do
    end do

    ! write the dataset array values
    call XF_WRITE_VECTOR_TIMESTEP(xVector2D_A, dTime, nValues, nComponents, &
                                  fValues, error)
    if (error .GE. 0) then
      ! write activity array
      call XF_WRITE_ACTIVITY_TIMESTEP(xVector2D_A, nActive, bActivity, error)
    end if 

    call TTI_Test_Num_Times(xVector2d_A, iTimestep, error)
  enddo

  ! Write Coordinate file - for Vector2D_A, we will set the coordinate system
  !   to be Geographic HPGN, with HPGN settings written to the file.
  call XF_CREATE_COORDINATE_GROUP(xFileId, xCoordId, status)
  if (status < 0) then
    call XF_CLOSE_GROUP(xVector2D_A, error)
    call XF_CLOSE_GROUP(xDsetsId, error)
    call XF_CLOSE_FILE(xFileId, error)
  error = -1
  return
  endif

    ! set HPGN info for test
  iHpgnZone = 29   ! Utah

  call XF_SET_HORIZ_DATUM(xCoordId, HORIZ_DATUM_GEOGRAPHIC_HPGN, error)
  call XF_SET_HORIZ_UNITS(xCoordId, COORD_UNITS_METERS, error)
  call XF_SET_VERT_DATUM(xCoordId, VERT_DATUM_LOCAL, error)
  call XF_SET_VERT_UNITS(xCoordId, COORD_UNITS_METERS, error)

    ! write additional information
  call XF_SET_HPGN_AREA(xCoordId, iHpgnZone, error)

  call XF_CLOSE_GROUP(xCoordId, error)
  xCoordId = 0

  ! close the dataset
  call XF_CLOSE_GROUP(xVector2D_A, error)
  call XF_CLOSE_GROUP(xDsetsId, error)
  call XF_CLOSE_FILE(xFileId, error)

  return
END SUBROUTINE
! ttWriteVector2D_A

! --------------------------------------------------------------------------
! FUNCTION TT_WRITE_VECTOR2D_B
! PURPOSE  Write scalar Dataset to an HDF5 File
! NOTES    This tests dynamic data sets, and activity
!          This dataset is dynamic concentrations (mg/L) with output times
!          in minutes.
!          Dataset is for a mesh and so nActive is the number of elements
!          which is not the same as the nValues which would be number of nodes
!          reads/writes a reference time in julian days
! --------------------------------------------------------------------------
SUBROUTINE TT_WRITE_VECTOR2D_B (a_Filename, a_Compression, a_Overwrite, error)
  CHARACTER(LEN=*), INTENT(IN) :: a_Filename
  INTEGER, INTENT(IN)          :: a_Compression
  LOGICAL, INTENT(IN)          :: a_Overwrite
  INTEGER, INTENT(OUT)         :: error
  INTEGER      xFileId, xDsetsId, xVector2D_B, xCoordId
  INTEGER      nValues, nTimes, nComponents, nActive
  REAL(DOUBLE) dTime
  INTEGER      iTimestep, iActive
  REAL, DIMENSION(2, 100) :: fValues
  INTEGER*1    bActivity(100)
  INTEGER      i, j, status

    ! initialize the data
  nComponents = 2
  nValues = 100
  nTimes = 6
  nActive = 75
  dTime = 0.0

    ! 5th item in data set is always inactive, others active
  bActivity(1) = 0
  do iActive = 2, nActive
    if (mod(iActive-1, 3) == 0) then
      bActivity(iActive) = 0
    else
      bActivity(iActive) = 1
    endif
  enddo

  if (a_Overwrite) then
      ! open the already-existing file
    call XF_OPEN_FILE(a_Filename, .FALSE., xFileId, status)
    if (status < 0) then
      error = -1
      return
    endif
      ! open the group where we have all the datasets
    call XF_OPEN_GROUP(xFileId, DATASETS_LOCATION, xDsetsId, status)
    if (status < 0) then
      call XF_CLOSE_FILE(xFileId, error)
      error = -1
      return
    endif
  else
      ! create the file
    call XF_CREATE_FILE(a_Filename, .TRUE., xFileId, error)
    if (error .LT. 0) then
        ! close the dataset
      call XF_CLOSE_FILE(xFileId, error)
      return
    endif

      ! create the group where we will put all the datasets 
    call XF_CREATE_GENERIC_GROUP(xFileId, DATASETS_LOCATION, xDsetsId, status)
    if (status < 0) then
      call XF_CLOSE_FILE(xFileId, error)
      error = -1
      return
    endif
  endif

    ! Create/Overwrite the vector dataset group
  call XF_CREATE_VECTOR_DATASET(xDsetsId, VECTOR2D_B_LOCATION, 'ft/s', &
                                TS_SECONDS, a_Compression, xVector2D_B, status)
  if (status .LT. 0) then
      ! close the dataset
    call XF_CLOSE_GROUP(xVector2D_B, error)
    call XF_CLOSE_GROUP(xDsetsId, error)
    call XF_CLOSE_FILE(xFileId, error)
    error = status
    return 
  endif

  if (.NOT. a_Overwrite) then
      ! Loop through timesteps adding them to the file
    do iTimestep = 1, nTimes
        ! We will have an 0.5 hour timestep
      dTime = iTimestep * 0.5
      do i = 1, nValues
        do j = 1, nComponents
          fValues(j,i) = ((i-1)*nComponents + j)*dTime
        end do
      end do
        ! write the dataset array values
      call XF_WRITE_VECTOR_TIMESTEP(xVector2D_B, dTime, nValues, nComponents, &
                                    fValues, error)
      if (error .GE. 0) then
          ! write activity array
        call XF_WRITE_ACTIVITY_TIMESTEP(xVector2D_B, nActive, bActivity, error)
      end if
      if (error < 0) then
        call XF_CLOSE_GROUP(xVector2D_B, error)
        call XF_CLOSE_GROUP(xDsetsId, error)
        call XF_CLOSE_FILE(xFileId, error)
      endif
    enddo
  else
      ! Loop through timesteps adding them to the file
    do iTimestep = 1, nTimes
        ! We will have an 1.5 hour timestep
      dTime = iTimestep * 1.5
      do i = 1, nValues
        do j = 1, nComponents
          fValues(j,i) = ((i-1)*nComponents + j)*dTime
        end do
      end do
        ! write the dataset array values
      call XF_WRITE_VECTOR_TIMESTEP(xVector2D_B, dTime, nValues, nComponents, &
                                    fValues, error)
      if (error .GE. 0) then
          ! write activity array
        call XF_WRITE_ACTIVITY_TIMESTEP(xVector2D_B, nActive, bActivity, error)
      end if
      if (error < 0) then
        call XF_CLOSE_GROUP(xVector2D_B, error)
        call XF_CLOSE_GROUP(xDsetsId, error)
        call XF_CLOSE_FILE(xFileId, error)
      endif
    enddo
  endif

  if (.NOT. a_Overwrite) then
    ! Write Coordinate file - for ScalarB, we will set the coordinate system
    !   to be UTM, with UTM Zone settings written to the file.
    call XF_CREATE_COORDINATE_GROUP(xFileId, xCoordId, status)
    if (status < 0) then
      call XF_CLOSE_GROUP(xVector2D_B, error)
      call XF_CLOSE_GROUP(xDsetsId, error)
      call XF_CLOSE_FILE(xFileId, error)
    error = -1
    return
    endif

      ! write the coordinate data to the file
    call XF_SET_HORIZ_DATUM(xCoordId, HORIZ_DATUM_UTM, error)
    call XF_SET_HORIZ_UNITS(xCoordId, COORD_UNITS_METERS, error)
    call XF_SET_VERT_DATUM(xCoordId, VERT_DATUM_LOCAL, error)
    call XF_SET_VERT_UNITS(xCoordId, COORD_UNITS_METERS, error)

      ! write additional information - we'll use the max UTM zone for the test
    call XF_SET_UTM_ZONE(xCoordId, UTM_ZONE_MAX, error)

    call XF_CLOSE_GROUP(xCoordId, error)
    xCoordId = 0
  endif

  ! close the dataset
  call XF_CLOSE_GROUP(xVector2D_B, error)
  call XF_CLOSE_GROUP(xDsetsId, error)
  call XF_CLOSE_FILE(xFileId, error)

  return
END SUBROUTINE
! ttWriteVector2D_B

! --------------------------------------------------------------------------
! FUNCTION ttWriteVector2D_AToMulti
! PURPOSE  Write scalar Dataset to an HDF5 File
! NOTES    This tests dynamic data sets, and activity
!          This dataset is dynamic concentrations (mg/L) with output times
!          in minutes.
!          Dataset is for a mesh and so nActive is the number of elements
!          which is not the same as the nValues which would be number of nodes
!          reads/writes a reference time in julian days
! --------------------------------------------------------------------------
SUBROUTINE TT_WRITE_VECTOR2D_A_TO_MULTI (a_FileID, a_GroupID, status)
  INTEGER      xVector2D_A
  INTEGER      a_FileID, a_GroupID
  INTEGER      nValues, nTimes, nComponents, nActive
  REAL(DOUBLE) dTime
  INTEGER      iTimestep, iActive
  REAL, DIMENSION(2, 100) :: fValues ! nComponents, nValues
  INTEGER*1    bActivity(100) ! activity
  INTEGER      i, j, status

  ! initialize the data
  nComponents = 2
  nValues = 100
  nTimes = 6
  nActive = 75
  dTime = 0.0

  ! 5th item in data set is always inactive, others active
  bActivity(1) = 0
  do iActive = 2, nActive
    if (mod(iActive-1, 3) == 0) then
      bActivity(iActive) = 0
    else
    bActivity(iActive) = 1
  endif
  enddo

  ! Create the vector dataset group
  call XF_CREATE_VECTOR_DATASET(a_GroupID, VECTOR2D_A_LOCATION, 'ft/s', &
              TS_SECONDS, NONE, xVector2D_A, status)
  if (status .LT. 0) then
      ! close the dataset
    call XF_CLOSE_GROUP(xVector2D_A, status)
    call XF_CLOSE_GROUP(a_GroupID, status)
    call XF_CLOSE_FILE(a_FileID, status)
    return 
  endif

  ! Loop through timesteps adding them to the file
  do iTimestep = 1, nTimes
    ! We will have an 0.5 hour timestep
    dTime = iTimestep * 0.5

    do i = 1, nValues
      do j = 1, nComponents
        fValues(j,i) = ((i-1)*nComponents + j)*dTime
      end do
    end do

    ! write the dataset array values
    call XF_WRITE_VECTOR_TIMESTEP(xVector2D_A, dTime, nValues, nComponents, &
                                  fValues, status)
    if (status .GE. 0) then
      ! write activity array
      call XF_WRITE_ACTIVITY_TIMESTEP(xVector2D_A, nActive, bActivity, status)
    end if 

    call TTI_Test_Num_Times(xVector2D_A, iTimestep, status)
  enddo

  ! close the dataset
  call XF_CLOSE_GROUP(xVector2D_A, status)
  return
END SUBROUTINE
! ttWriteVector2D_AToMulti
! --------------------------------------------------------------------------
! FUNCTION ttiReadScalar
! PURPOSE  Read a scalar from an XMDF file and output information to
!          to a text file
! NOTES    
! --------------------------------------------------------------------------
SUBROUTINE TTI_READ_SCALAR (a_xScalarId, FileUnit, error)
  INTEGER, INTENT(IN) ::  a_xScalarId
  INTEGER, INTENT(IN) ::         FileUnit
  INTEGER, INTENT(OUT) :: error
  INTEGER             nTimes, nValues, nActive
  LOGICAL*2             bUseReftime
  INTEGER             iTime
  CHARACTER(LEN=100)   TimeUnits
  REAL(DOUBLE), ALLOCATABLE :: Times(:)
  REAL, ALLOCATABLE         :: Values(:), Minimums(:), Maximums(:)
  INTEGER, ALLOCATABLE      :: Active(:)
  REAL(DOUBLE)                 Reftime
nTimes = NONE
nValues = NONE
nActive = None

  ! read the time units
  call XF_GET_DATASET_TIME_UNITS(a_xScalarId, TimeUnits, error)
  if (error < 0) return

  WRITE(FileUnit,*) 'Time units: ', TimeUnits(1:LEN_TRIM(TimeUnits))

  ! see if we are using a reftime
  call XF_USE_REFTIME (a_xScalarId, bUseReftime, error)
  if (error < 0) then
    return
  endif
  if (bUseReftime) then
    call XF_READ_REFTIME (a_xScalarId, Reftime, error)
    if (error < 0) then
      return
  endif
    WRITE(FileUnit,*) 'Reftime: ', Reftime
  endif

  ! read in the number of values and number of active values
  call XF_GET_DATASET_NUMVALS(a_xScalarId, nValues, error)
  if (error .GE. 0) then
    call XF_GET_DATASET_NUMACTIVE(a_xScalarId, nActive, error)
  endif
  if (error .LT. 0) return 

  if (nValues <= 0) then
    WRITE(FileUnit, *) 'No data to read in.'
    error = -1
    return 
  endif

  ! read in the number of times
  call XF_GET_DATASET_NUM_TIMES(a_xScalarId, nTimes, error)
  if (error < 0) then
    return 
  endif

  ! Read in the individual time values
  allocate(Times(nTimes))

  call XF_GET_DATASET_TIMES(a_xScalarId, nTimes, Times, error)
  if (error < 0) return 

  ! Read in the minimum and maximum values
  allocate(Minimums(nTimes))
  allocate(Maximums(nTimes))

  call XF_GET_DATASET_MINS(a_xScalarId, nTimes, Minimums, error)
  if (error >= 0) then
    call XF_GET_DATASET_MAXS(a_xScalarId, nTimes, Maximums, error)
  endif
  if (error < 0) then
    deallocate(Times)
    deallocate(Minimums)
    deallocate(Maximums)
    return
  endif

  allocate(Values(nValues))
  if (nActive .GT. 0) then
    allocate(Active(nActive))
  endif

  WRITE(FileUnit,*) 'Number Timesteps: ', nTimes
  WRITE(FileUnit,*) 'Number Values: ', nValues
  WRITE(FileUnit,*) 'Number Active: ', nActive
  WRITE(FileUnit,*) ''

  ! loop through the timesteps, read the values and active values and write
  ! them to the text file
  do iTime = 1, nTimes
    call XF_READ_SCALAR_VALUES_TIMESTEP(a_xScalarId, iTime, nValues, Values, error)
    if (error >= 0 .AND. nActive > 0) then
      call XF_READ_ACTIVITY_TIMESTEP(a_xScalarId, iTime, nActive, Active, error)
    endif

    ! Write the time, min, max, values and active values to the text output
    ! file.
    WRITE(FileUnit,*) 'Timestep at  ', Times(iTime)
    WRITE(FileUnit,*) 'Min: ', Minimums(iTime)
    WRITE(FileUnit,*) 'Max: ', Maximums(iTime)

    WRITE(FileUnit,*) 'Values:'
    WRITE(FileUnit,*) Values(1:nValues)
    WRITE(FileUnit,*) ''

    WRITE(FileUnit,*) 'Activity:'
    WRITE(FileUnit,*) Active(1:nActive)
    WRITE(FileUnit,*) ''
  end do

  if (allocated(Times)) then
    deallocate(Times)
  endif
  
  if (allocated(Minimums)) then
    deallocate(Minimums)
  endif

  if (allocated(Maximums)) then
    deallocate(Maximums)
  endif

  if (allocated(Values)) then
    deallocate(Values)
  endif

  if (allocated(Active)) then
    deallocate(Active)
  endif

  return
END SUBROUTINE
! ttiReadScalar

! --------------------------------------------------------------------------
! FUNCTION TTI_READ_VECTOR
! PURPOSE  Read a vector from an XMDF file and output information to
!          to a text file
! NOTES    
! --------------------------------------------------------------------------
SUBROUTINE TTI_READ_VECTOR (a_xVectorId, FileUnit, error)
  INTEGER, INTENT(IN) ::  a_xVectorId
  INTEGER, INTENT(IN) ::         FileUnit
  INTEGER, INTENT(OUT) :: error
  INTEGER             nTimes, nValues, nActive, nComponents
  INTEGER             iTime, i
  LOGICAL*2            bUseReftime
  CHARACTER(LEN=100)   TimeUnits
  REAL(DOUBLE), ALLOCATABLE :: Times(:)
  REAL, ALLOCATABLE, DIMENSION (:, :) :: Values
  REAL, ALLOCATABLE         :: Minimums(:), Maximums(:)
  INTEGER, ALLOCATABLE      :: Active(:)
  REAL(DOUBLE)                 Reftime

nTimes = NONE
nValues = NONE
nActive = NONE
nComponents = NONE

  ! read the time units
  call XF_GET_DATASET_TIME_UNITS(a_xVectorId, TimeUnits, error)
  if (error < 0) return

  WRITE(FileUnit,*) 'Time units: ', TimeUnits(1:LEN_TRIM(TimeUnits))

  ! see if we are using a reftime
  call XF_USE_REFTIME (a_xVectorId, bUseReftime, error)
  if (error < 0) then
    return
  endif
  if (bUseReftime) then
    call XF_READ_REFTIME (a_xVectorId, Reftime, error)
    if (error < 0) then
      return
  endif
    WRITE(FileUnit,*) 'Reftime: ', Reftime
  endif

  ! read in the number of values and number of active values
  call XF_GET_DATASET_NUMVALS(a_xVectorId, nValues, error)
  if (error .GE. 0) then
    call XF_GET_DATASET_NUMCOMPONENTS(a_xVectorId, nComponents, error)
    if (error .GE. 0) then
      call XF_GET_DATASET_NUMACTIVE(a_xVectorId, nActive, error)
    endif
  endif
  if (error .LT. 0) return 

  if (nValues <= 0) then
    WRITE(FileUnit, *) 'No data to read in.'
    error = -1
    return 
  endif

  ! read in the number of times
  call XF_GET_DATASET_NUM_TIMES(a_xVectorId, nTimes, error)
  if (error < 0) then
    return 
  endif

  ! Read in the individual time values
  allocate(Times(nTimes))

  call XF_GET_DATASET_TIMES(a_xVectorId, nTimes, Times, error)
  if (error < 0) return 

  ! Read in the minimum and maximum values
  allocate(Minimums(nTimes))
  allocate(Maximums(nTimes))

  call XF_GET_DATASET_MINS(a_xVectorId, nTimes, Minimums, error)
  if (error >= 0) then
    call XF_GET_DATASET_MAXS(a_xVectorId, nTimes, Maximums, error)
  endif
  if (error < 0) then
    deallocate(Times)
    deallocate(Minimums)
    deallocate(Maximums)
    return
  endif

  allocate(Values(nComponents, nValues))
  if (nActive .GT. 0) then
    allocate(Active(nActive))
  endif

  WRITE(FileUnit,*) 'Number Timesteps: ', nTimes
  WRITE(FileUnit,*) 'Number Values: ', nValues
  WRITE(FileUnit,*) 'Number Components: ', nComponents
  WRITE(FileUnit,*) 'Number Active: ', nActive

  ! loop through the timesteps, read the values and active values and write
  ! them to the text file
  do iTime = 1, nTimes
    call XF_READ_VECTOR_VALUES_TIMESTEP(a_xVectorId, iTime, nValues, &
                                        nComponents, Values, error)
    if (error >= 0 .AND. nActive > 0) then
      call XF_READ_ACTIVITY_TIMESTEP(a_xVectorId, iTime, nActive, Active, error)
    endif

    ! Write the time, min, max, values and active values to the text output
    ! file.
  WRITE(FileUnit,*) ''
    WRITE(FileUnit,*) 'Timestep at  ', Times(iTime)
    WRITE(FileUnit,*) 'Min: ', Minimums(iTime)
    WRITE(FileUnit,*) 'Max: ', Maximums(iTime)

    WRITE(FileUnit,*) 'Values:'
    do i=1, nValues
      WRITE(FileUnit,*) Values(1:nComponents,i:i)
    enddo
    WRITE(FileUnit,*) ''

    WRITE(FileUnit,*) 'Activity:'
    WRITE(FileUnit,*) Active(1:nActive)
    WRITE(FileUnit,*) ''    
  WRITE(FileUnit,*) ''    

  end do

  if (allocated(Times)) then
    deallocate(Times)
  endif
  
  if (allocated(Minimums)) then
    deallocate(Minimums)
  endif

  if (allocated(Maximums)) then
    deallocate(Maximums)
  endif

  if (allocated(Values)) then
    deallocate(Values)
  endif

  if (allocated(Active)) then
    deallocate(Active)
  endif

  return
END SUBROUTINE
! ttiReadVector

END MODULE TestTimestep

TestTimestep.f90 tests timesteps

MODULE TestDefs

INTEGER, PARAMETER :: NUMTIMES1   = 5
INTEGER, PARAMETER :: NUMVALUES1  = 5
INTEGER, PARAMETER :: NUMACTIVE1  = 3
INTEGER, PARAMETER :: NUMTIMESADD = 1

CHARACTER(LEN=*), PARAMETER :: XMDF_VERSION_OUT_F = 'XMDF_Version_f.txt'
CHARACTER(LEN=*), PARAMETER :: MESH_A_FILE_F = 'mesh_a_file_f.h5'
CHARACTER(LEN=*), PARAMETER :: MESH_B_FILE_F = 'mesh_b_file_f.h5'
CHARACTER(LEN=*), PARAMETER :: MESH_A_OUT_F  = 'mesh_a_file_f.txt'
CHARACTER(LEN=*), PARAMETER :: MESH_B_OUT_F  = 'mesh_b_file_f.txt'

CHARACTER(LEN=*), PARAMETER :: GRID_CART2D_A_FILE_F = 'grid_cart2d_a_file_f.h5'
CHARACTER(LEN=*), PARAMETER :: GRID_CURV2D_A_FILE_F = 'grid_curv2d_a_file_f.h5'
CHARACTER(LEN=*), PARAMETER :: GRID_CART3D_A_FILE_F = 'grid_cart3d_a_file_f.h5'

CHARACTER(LEN=*), PARAMETER :: GRID_CART2D_A_OUT_F = 'grid_cart2d_a_out_f.txt'
CHARACTER(LEN=*), PARAMETER :: GRID_CURV2D_A_OUT_F = 'grid_curv2d_a_out_f.txt'
CHARACTER(LEN=*), PARAMETER :: GRID_CART3D_A_OUT_F = 'grid_cart3d_a_out_f.txt'

CHARACTER(LEN=*), PARAMETER :: MULTIDATASET_FILE_F = 'MultiDataSet_f.h5'
CHARACTER(LEN=*), PARAMETER :: MULTIDATASET_TEXT_F = 'MultiDataSet_f.txt'

CHARACTER(LEN=*), PARAMETER :: SCALAR_A_FILE_F   = 'ScalarA_f.h5'
CHARACTER(LEN=*), PARAMETER :: SCALAR_A_TEXT_F   = 'ScalarA_f.txt'
CHARACTER(LEN=*), PARAMETER :: SCALAR_A_PIECES_FILE_F   = 'ScalarA_Pieces_f.h5'
CHARACTER(LEN=*), PARAMETER :: SCALAR_A_EDITED_FILE_F   = 'ScalarA_edited_f.h5'
CHARACTER(LEN=*), PARAMETER :: SCALAR_A_EDITED_TEXT_F   = 'ScalarA_edited_f.txt'
CHARACTER(LEN=*), PARAMETER :: SCALAR_B_FILE_F   = 'ScalarB_f.h5'
CHARACTER(LEN=*), PARAMETER :: SCALAR_B_TEXT_F   = 'ScalarB_f.txt'
CHARACTER(LEN=*), PARAMETER :: VECTOR2D_A_FILE_F = 'Vector2D_A_f.h5'
CHARACTER(LEN=*), PARAMETER :: VECTOR2D_A_TEXT_F = 'Vector2D_A_f.txt'
CHARACTER(LEN=*), PARAMETER :: VECTOR2D_A_PIECES_FILE_F = 'Vector2D_A_Pieces_f.h5'
CHARACTER(LEN=*), PARAMETER :: VECTOR2D_B_FILE_F = 'Vector2D_B_f.h5'
CHARACTER(LEN=*), PARAMETER :: VECTOR2D_B_TEXT_F = 'Vector2D_B_f.txt'

CHARACTER(LEN=*), PARAMETER :: MESH_A_FILE_C = 'mesh_a_file_c.h5'
CHARACTER(LEN=*), PARAMETER :: MESH_B_FILE_C = 'mesh_b_file_c.h5'
CHARACTER(LEN=*), PARAMETER :: MESH_A_OUT_CF  = 'mesh_a_file_cf.txt'
CHARACTER(LEN=*), PARAMETER :: MESH_B_OUT_CF  = 'mesh_b_file_cf.txt'

CHARACTER(LEN=*), PARAMETER :: GRID_CART2D_A_FILE_C = 'grid_cart2d_a_file_c.h5'
CHARACTER(LEN=*), PARAMETER :: GRID_CURV2D_A_FILE_C = 'grid_curv2d_a_file_c.h5'
CHARACTER(LEN=*), PARAMETER :: GRID_CART3D_A_FILE_C = 'grid_cart3d_a_file_c.h5'

CHARACTER(LEN=*), PARAMETER :: GRID_CART2D_A_OUT_CF = 'grid_cart2d_a_out_cf.txt'
CHARACTER(LEN=*), PARAMETER :: GRID_CURV2D_A_OUT_CF = 'grid_curv2d_a_out_cf.txt'
CHARACTER(LEN=*), PARAMETER :: GRID_CART3D_A_OUT_CF = 'grid_cart3d_a_out_cf.txt'

CHARACTER(LEN=*), PARAMETER :: SCALAR_A_FILE_C   = 'ScalarA_c.h5'
CHARACTER(LEN=*), PARAMETER :: SCALAR_A_TEXT_CF   = 'ScalarA_cf.txt'
CHARACTER(LEN=*), PARAMETER :: SCALAR_B_FILE_C   = 'ScalarB_c.h5'
CHARACTER(LEN=*), PARAMETER :: SCALAR_B_TEXT_CF   = 'ScalarB_cf.txt'
CHARACTER(LEN=*), PARAMETER :: VECTOR2D_A_FILE_C = 'Vector2D_A_c.h5'
CHARACTER(LEN=*), PARAMETER :: VECTOR2D_A_TEXT_CF = 'Vector2D_A_cf.txt'
CHARACTER(LEN=*), PARAMETER :: VECTOR2D_B_FILE_C = 'Vector2D_B_c.h5'
CHARACTER(LEN=*), PARAMETER :: VECTOR2D_B_TEXT_CF = 'Vector2D_B_cf.txt'

CHARACTER(LEN=*), PARAMETER :: CALENDAR_OUT_F = 'Calendar_f.txt'

CHARACTER(LEN=*), PARAMETER :: GEOMPATH_A_FILE_F = 'Geompath_a_file_f.h5'
CHARACTER(LEN=*), PARAMETER :: GEOMPATH_A_FILE_F_OUT = 'Geompath_a_file_f_out.txt'

CHARACTER(LEN=*), PARAMETER :: TT_MULTIDATASET_FILE_F = 'TT_MultiDataSet_f.h5'
CHARACTER(LEN=*), PARAMETER :: TT_SCALAR_A_FILE_F   = 'TT_ScalarA_f.h5'
CHARACTER(LEN=*), PARAMETER :: TT_SCALAR_A_TEXT_F   = 'TT_ScalarA_f.txt'
CHARACTER(LEN=*), PARAMETER :: TT_VECTOR2D_A_FILE_F = 'TT_Vector2D_A_f.h5'

  ! Overwrite options in the function xfSetupToWriteDatasets
!INTEGER,PARAMETER :: XF_OVERWRITE_CLEAR_FILE          = 1
!INTEGER,PARAMETER :: XF_OVERWRITE_CLEAR_DATASET_GROUP = 2
!INTEGER,PARAMETER :: XF_OVERWRITE_NONE                = 3

END MODULE TestDefs

MODULE TestsModule

USE XMDF
USE XMDFDEFS
USE TestTimestep
USE TestDatasets
USE TestMesh
USE TestGrid
USE TestDefs
USE TestGeomPaths

CONTAINS

!**************************
!-----------------------------------------------------------------------------
! SUBROUTINE  TXI_TEST_TIMESTEPS
! PURPOSE     test to see if the code can read timestepfiles 
! NOTES
!------------------------------------------------------------------------------
SUBROUTINE TXI_TEST_TIMESTEPS(error)
INTEGER, INTENT(OUT) :: error
INTEGER    status, compression
INTEGER  MultiFileId, MultiGroupId
INTEGER         NumOpen

CHARACTER(LEN=37) SdoGuid

SdoGuid = '73289C80-6235-4fdc-9649-49E4F5AEB676'

status = 1
compression = NONE
! 'Path' should be able to be blank, but now it creates errors if it's blank
 call XF_SETUP_TO_WRITE_DATASETS(TT_MULTIDATASET_FILE_F, 'Multidatasets','', &
           SdoGuid, XF_OVERWRITE_CLEAR_FILE, MultiFileId, MultiGroupId, error)

  ! Write coordinates to multidatasets
call TT_WRITE_COORDS_TO_MULTI(MultiFileId, error)

  ! Write scalar A and Vector A to multidatasets.
call TT_WRITE_SCALAR_A_TO_MULTI(MultiGroupId, error)

call TT_WRITE_VECTOR2D_A_TO_MULTI(MultiFileId, MultiGroupId, error)

call XF_CLOSE_GROUP(MultiGroupId, status)
call XF_CLOSE_FILE(MultiFileId, status)

WRITE(*,*) 'Done writing multiple datasets...'

  ! scalar datasets
call TT_WRITE_SCALAR_A(TT_SCALAR_A_FILE_F, compression, status)
if (status < 0) then
  error = status
  return
endif
  
WRITE(*,*) 'Done writing scalar datasets...'

  ! vector datasets
call TT_WRITE_VECTOR2D_A(TT_VECTOR2D_A_FILE_F, compression, status)
if (status < 0) then
  WRITE(*,*) 'Error writing dataset vector2D_A.'
  error = status
  return
endif

WRITE(*,*) 'Done writing vector datasets...'

WRITE(*,*) 'Write edited scalar datasets...'
call TD_EDIT_SCALAR_A_VALUES(SCALAR_A_EDITED_FILE_F, compression, status);
if (status < 0) then
  error = status
  return
endif

WRITE(*,*) 'Done writing datasets...'

  ! Read the files back in
call TXI_READ_X_FORMAT_FILE(TT_SCALAR_A_FILE_F, SCALAR_A_TEXT_F, status)
if (status < 0) then
  WRITE(*,*) 'Error reading SCALAR A File (see TXI_READ_X_FORMAT_FILE)'
  error = status
  return
endif

call TXI_READ_X_FORMAT_FILE(SCALAR_A_EDITED_FILE_F, SCALAR_A_EDITED_TEXT_F, status)
if (status < 0) then
  WRITE(*,*) 'Error reading SCALAR A Edited File (see TXI_READ_X_FORMAT_FILE)'
  error = status
  return
endif

call  TXI_READ_X_FORMAT_FILE(TT_VECTOR2D_A_FILE_F, VECTOR2D_A_TEXT_F, status)
if (status < 0) then
  WRITE(*,*) 'Error reading VECTOR A Format File'
  error = status
  return
endif

call TXI_READ_X_FORMAT_FILE(TT_MULTIDATASET_FILE_F, MULTIDATASET_TEXT_F, status)
if (status < 0) then
  WRITE(*,*) 'Error reading Multidataset File (see TXI_READ_X_FORMAT_FILE)'
  error = status
  return
endif

WRITE(*,*) 'Done reading datasets...'

call XF_GET_NUM_OPEN_IDENTIFIERS(H5F_OBJ_ALL_F, NumOpen, error)

call XFI_CLOSE_OPEN_IDENTIFIERS(H5F_OBJ_ALL_F, error)

call XF_SETUP_TO_WRITE_DATASETS(TT_MULTIDATASET_FILE_F, 'Multidatasets','', &
     SdoGuid, XF_OVERWRITE_CLEAR_DATASET_GRP, MultiFileId, MultiGroupId, &
                                                                  error)

call TT_WRITE_SCALAR_A_TO_MULTI(MultiGroupId, error)

call XF_SETUP_TO_WRITE_DATASETS(TT_MULTIDATASET_FILE_F, 'Multidatasets','', &
     SdoGuid, XF_OVERWRITE_NONE, MultiFileId, MultiGroupId, error)

call TT_WRITE_VECTOR2D_A_TO_MULTI(MultiFileId, MultiGroupId, error)

  ! Test reading information at index for multiple timesteps
call TT_READ_SCALAR_A_INDEX(TT_SCALAR_A_FILE_F, 4, status)
if (status < 0) then
  error = status
  return
endif

WRITE(*,*) 'Done reading scalar data at index.'
  
call TT_READ_VECTOR2D_A_INDEX(TT_VECTOR2D_A_FILE_F, 6, status)
if (status < 0) then
  error = status
  return
endif

WRITE(*,*) 'Done reading vector data at index.'

call TT_READ_ACTIVITY_SCALAR_A_INDEX(TT_SCALAR_A_FILE_F, 6, status)
if (status < 0) then
  error = status
  return
endif

error = status
return

END SUBROUTINE


!**************************
!-----------------------------------------------------------------------------
! SUBROUTINE  TXI_TEST_DATASETS
! PURPOSE     test to see if the code can read datasetfiles 
! NOTES
!------------------------------------------------------------------------------
SUBROUTINE TXI_TEST_DATASETS(error)
INTEGER, INTENT(OUT) :: error
INTEGER    status, compression
INTEGER  MultiFileId, MultiGroupId
INTEGER         NumOpen
INTEGER, PARAMETER :: nIndices = 3
INTEGER, DIMENSION(nIndices) :: indices

CHARACTER(LEN=37) SdoGuid

SdoGuid = '73289C80-6235-4fdc-9649-49E4F5AEB676'

status = 1
compression = NONE
! 'Path' should be able to be blank, but now it creates errors if it's blank
call XF_SETUP_TO_WRITE_DATASETS(MULTIDATASET_FILE_F, 'Multidatasets','', &
           SdoGuid, XF_OVERWRITE_CLEAR_FILE, MultiFileId, MultiGroupId, error)

  ! Write coordinates to multidatasets
call TD_WRITE_COORDS_TO_MULTI(MultiFileId, error)

  ! Write scalar A and Vector A to multidatasets.
call TD_WRITE_SCALAR_A_TO_MULTI(MultiGroupId, error)

call TD_WRITE_VECTOR2D_A_TO_MULTI(MultiFileId, MultiGroupId, error)

call XF_CLOSE_GROUP(MultiGroupId, status)
call XF_CLOSE_FILE(MultiFileId, status)

WRITE(*,*) 'Done writing multiple datasets...'

  ! scalar datasets
call TD_WRITE_SCALAR_A(SCALAR_A_FILE_F, compression, status)
if (status < 0) then
  error = status
  return
endif

call TD_WRITE_SCALAR_A_PIECES(SCALAR_A_PIECES_FILE_F, compression, status)
if (status < 0) then
  error = status
  return
endif
  
WRITE(*,*) 'Done writing scalar datasets...'

  ! vector datasets
call TD_WRITE_VECTOR2D_A(VECTOR2D_A_FILE_F, compression, status)
if (status < 0) then
  WRITE(*,*) 'Error writing dataset vector2D_A.'
  error = status
  return
endif

call TD_WRITE_VECTOR2D_A_PIECES(VECTOR2D_A_PIECES_FILE_F, compression, status)
if (status < 0) then
  WRITE(*,*) 'Error writing dataset vector2D_A.'
  error = status
  return
endif

WRITE(*,*) 'Done writing vector datasets...'

WRITE(*,*) 'Done writing datasets...'

  ! Read the files back in
call TXI_READ_X_FORMAT_FILE(SCALAR_A_FILE_F, SCALAR_A_TEXT_F, status)
if (status < 0) then
  WRITE(*,*) 'Error reading SCALAR A File (see TXI_READ_X_FORMAT_FILE)'
  error = status
  return
endif

call  TXI_READ_X_FORMAT_FILE(VECTOR2D_A_FILE_F, VECTOR2D_A_TEXT_F, status)
if (status < 0) then
  WRITE(*,*) 'Error reading VECTOR A Format File'
  error = status
  return
endif

call TXI_READ_X_FORMAT_FILE(MULTIDATASET_FILE_F, MULTIDATASET_TEXT_F, status)
if (status < 0) then
  WRITE(*,*) 'Error reading Multidataset File (see TXI_READ_X_FORMAT_FILE)'
  error = status
  return
endif

WRITE(*,*) 'Done reading datasets...'

call XF_GET_NUM_OPEN_IDENTIFIERS(H5F_OBJ_ALL_F, NumOpen, error)

call XFI_CLOSE_OPEN_IDENTIFIERS(H5F_OBJ_ALL_F, error)

call XF_SETUP_TO_WRITE_DATASETS(MULTIDATASET_FILE_F, 'Multidatasets','', &
     SdoGuid, XF_OVERWRITE_CLEAR_DATASET_GRP, MultiFileId, MultiGroupId, &
                                                                  error)

call TD_WRITE_SCALAR_A_TO_MULTI(MultiGroupId, error)

call XF_SETUP_TO_WRITE_DATASETS(MULTIDATASET_FILE_F, 'Multidatasets','', &
     SdoGuid, XF_OVERWRITE_NONE, MultiFileId, MultiGroupId, error)

call TD_WRITE_VECTOR2D_A_TO_MULTI(MultiFileId, MultiGroupId, error)

  ! Test reading information at index for multiple timesteps
call TD_READ_SCALAR_A_INDEX(SCALAR_A_FILE_F, 4, status)
if (status < 0) then
  error = status
  return
endif

WRITE(*,*) 'Done reading scalar data at index.'
  
call TD_READ_VECTOR2D_A_INDEX(VECTOR2D_A_FILE_F, 6, status)
if (status < 0) then
  error = status
  return
endif

WRITE(*,*) 'Done reading vector data at index.'

call TD_READ_ACTIVITY_SCALAR_A_INDEX(SCALAR_A_FILE_F, 6, status)
if (status < 0) then
  error = status
  return
endif

    ! Test reading information at multiple indices
  indices(1) = 2;
  indices(2) = 3;
  indices(3) = 5;
  CALL TD_READ_SCALAR_A_INDICES(SCALAR_A_FILE_F, nIndices, indices, error)

error = status
return

END SUBROUTINE

!------------------------------------------------------------------------------
!  FUNCTION  TXI_TEST_OVERWRITE_DSETS
!  PURPOSE   Check to see if already-written datasets can be overwritten
!  NOTES
!------------------------------------------------------------------------------
SUBROUTINE TXI_TEST_OVERWRITE_DSETS(error)
INTEGER, INTENT(OUT) :: error
INTEGER    status, compression

  status = 1
  compression = NONE

    ! scalar datasets
  call TD_WRITE_SCALAR_B(SCALAR_B_FILE_F, compression, .FALSE., status)
  if (status < 0) then
    error = status
    return
  endif
    !overwrite scalar datasets
  call TD_WRITE_SCALAR_B(SCALAR_B_FILE_F, compression, .TRUE., status)
  if (status < 0) then
    error = status
    return
  endif
  
    ! vector datasets
  call TD_WRITE_VECTOR2D_B(VECTOR2D_B_FILE_F, compression, .FALSE., status)
  if (status < 0) then
    WRITE(*,*) 'Error writing dataset vector2D_B.'
    error = status
    return
  endif
    ! overwrite vector datasets
  call TD_WRITE_VECTOR2D_B(VECTOR2D_B_FILE_F, compression, .TRUE., status)
  if (status < 0) then
    WRITE(*,*) 'Error writing dataset vector2D_B.'
    error = status
    return
  endif

    ! Read the files back in
  call TXI_READ_X_FORMAT_FILE(SCALAR_B_FILE_F, SCALAR_B_TEXT_F, status)
  if (status < 0) then
    WRITE(*,*) 'Error reading SCALAR B File'
    error = status
    return
  endif

  call  TXI_READ_X_FORMAT_FILE(VECTOR2D_B_FILE_F, VECTOR2D_B_TEXT_F, status)
  if (status < 0) then
    WRITE(*,*) 'Error reading VECTOR B Format File'
    error = status
    return
  endif

  error = status
  return

END SUBROUTINE

!------------------------------------------------------------------------------
! FUNCTION TXI_TEST_GRIDS
! PURPOSE  Test to see if we can read and write grids
! NOTES
!------------------------------------------------------------------------------
SUBROUTINE TXI_TEST_GRIDS (error)
INTEGER, INTENT(OUT) :: error
INTEGER    compression

compression = NONE

WRITE(*,*) ''
WRITE(*,*) ''
WRITE(*,*) 'Writing grid data.'
WRITE(*,*) ''
 
call TG_WRITE_TEST_GRID_CART_2D(GRID_CART2D_A_FILE_F, error)
if (error < 0) then
  WRITE(*,*) 'Error writing grid Cartesian 2D A'
endif
WRITE(*,*) 'Finished writing grid Cartesian 2D A'

call TG_WRITE_TEST_GRID_CURV_2D(GRID_CURV2D_A_FILE_F, compression, error)
if (error < 0) then
  WRITE(*,*) 'Error writing grid Curvilinear 2D A'
endif
WRITE(*,*) 'Finished writing grid Curvilinear 2D A'

call TG_WRITE_TEST_GRID_CART_3D(GRID_CART3D_A_FILE_F, compression, error)
if (error < 0) then
  WRITE(*,*) 'Error writing grid Cartesian 3D A'
endif
WRITE(*,*) 'Finished writing grid Cartesian 3D A'
  ! read the files back in
call TXI_READ_X_FORMAT_FILE(GRID_CART2D_A_FILE_F, GRID_CART2D_A_OUT_F, error)
if (error < 0) then
  WRITE(*,*) 'Error reading grid Cartesian 2D A'
endif
WRITE(*,*) 'Finished reading grid Cartesian 2D A'

call TXI_READ_X_FORMAT_FILE(GRID_CURV2D_A_FILE_F, GRID_CURV2D_A_OUT_F, error)
if (error < 0) then
  WRITE(*,*) 'Error reading grid Curvilinear 2D A'
endif
WRITE(*,*) 'Finished reading grid Curvilinear 2D A'

call TXI_READ_X_FORMAT_FILE(GRID_CART3D_A_FILE_F, GRID_CART3D_A_OUT_F, error)
if (error < 0) then
  WRITE(*,*) 'Error reading grid Cartesian 3D A'
endif
WRITE(*,*) 'Finished reading grid Cartesian 3D A'
 
END SUBROUTINE

!**************************
! ---------------------------------------------------------------------------
! FUNCTION  TXI_TEST_MESHS
! PURPOSE   test to see if we can read and write meshes
! NOTES     
! ---------------------------------------------------------------------------
SUBROUTINE TXI_TEST_MESHS (error)
INTEGER, INTENT(OUT) :: error
INTEGER    status
INTEGER    compression

status = 1
compression = NONE

call TM_WRITE_TEST_MESH_A(MESH_A_FILE_F, compression, status)
if (status < 0) then
  WRITE(*,*) 'Error writing TestMeshA'
  error = status
  return
endif

call TM_WRITE_TEST_MESH_B(MESH_B_FILE_F, compression, status)
if (status < 0) then
  WRITE(*,*) 'Error writing TestMeshB'
  error = status
  return
endif

WRITE(*,*) 'Finished writing meshes.'

  ! read the files back in
call TXI_READ_X_FORMAT_FILE(MESH_A_FILE_F, MESH_A_OUT_F, status)
if (status < 0) then
  WRITE(*,*) 'Error reading TestMeshA'
  error = status
  return
endif

  ! read the files back in
call TXI_READ_X_FORMAT_FILE(MESH_B_FILE_F, MESH_B_OUT_F, status)
if (status < 0) then
  WRITE(*,*) 'Error reading TestMeshB'
  error = status
  return
endif

WRITE(*,*) 'Finished reading meshes.'

error = status
return

END SUBROUTINE

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

!---------------------------------------------------------------------------
! FUNCTION  txiTestC
! PURPOSE   test to see if fortran code can read file written with C.
! NOTES     
!---------------------------------------------------------------------------
SUBROUTINE TXI_TEST_C (error)
INTEGER, INTENT(OUT) :: error
INTEGER            nStatus
INTEGER     xFileId

error = 1

  !Check to see if files written with C exist
  ! Turn off error handling
!call H5Eset_auto_f(0, error)
  ! Try opening a file written with C to see if one exists.
call XF_OPEN_FILE(SCALAR_A_FILE_C, .TRUE., xFileId, nStatus)
  ! If the file written with C doesn't exist, return.
if (nStatus < 0) then
  call XF_CLOSE_FILE(xFileId, error)
    ! Restore previous error handler
  !call H5Eset_Auto_f(1, error)
  error = -1
  return
  ! If the file written with C does exist, assume all C files exist.
else
  call XF_CLOSE_FILE(xFileId, error)
    ! Restore previous error handler
  !call H5Eset_Auto_f(1, error)
endif

  ! Read the files back in
call TXI_READ_X_FORMAT_FILE(SCALAR_A_FILE_C, SCALAR_A_TEXT_CF, error)
if (error < 0) then
  return
endif
call TXI_READ_X_FORMAT_FILE(SCALAR_B_FILE_C, SCALAR_B_TEXT_CF, error)
if (error < 0) then
  return
endif

call TXI_READ_X_FORMAT_FILE(VECTOR2D_A_FILE_C, VECTOR2D_A_TEXT_CF, error)
if (error < 0) then
  return
endif
call TXI_READ_X_FORMAT_FILE(VECTOR2D_B_FILE_C, VECTOR2D_B_TEXT_CF, error)
if (error < 0) then
  return
endif
  
WRITE(*,*) 'Done reading C datasets...'

call TXI_READ_X_FORMAT_FILE(GRID_CART2D_A_FILE_C, GRID_CART2D_A_OUT_CF, error)
if (error < 0) then
  WRITE(*,*) 'Error reading C grid Cartesian 2D A'
endif
WRITE(*,*) 'Finished reading C grid Cartesian 2D A'

call TXI_READ_X_FORMAT_FILE(GRID_CURV2D_A_FILE_C, GRID_CURV2D_A_OUT_CF, error)
if (error < 0) then
  WRITE(*,*) 'Error reading C grid Curvilinear 2D A'
endif
WRITE(*,*) 'Finished reading C grid Curvilinear 2D A'

call TXI_READ_X_FORMAT_FILE(GRID_CART3D_A_FILE_C, GRID_CART3D_A_OUT_CF, error)
if (error < 0) then
  WRITE(*,*) 'Error reading C grid Cartesian 3D A'
endif
WRITE(*,*) 'Finished reading C grid Cartesian 3D A'

  ! read the files back in
call TXI_READ_X_FORMAT_FILE(MESH_A_FILE_C, MESH_A_OUT_CF, error)
if (error < 0) then
  WRITE(*,*) 'Error reading C TestMeshA'
  return
endif

  ! read the files back in
call TXI_READ_X_FORMAT_FILE(MESH_B_FILE_C, MESH_B_OUT_CF, error)
if (error < 0) then
  WRITE(*,*) 'Error reading C TestMeshB'
  return
endif

WRITE(*,*) 'Finished reading C meshes.'
  
return

END SUBROUTINE

!**************************
! --------------------------------------------------------------------------
! FUNCTION txiReadXFormatFile
! PURPOSE  Read a file using XMDF and write information about the data
!          contained in the file to a output file
! --------------------------------------------------------------------------
SUBROUTINE TXI_READ_X_FORMAT_FILE (a_XmdfFile, a_OutFile, error)
CHARACTER(LEN=*), INTENT(IN) :: a_XmdfFile
CHARACTER(LEN=*), INTENT(IN) :: a_OutFile
INTEGER, INTENT(OUT)  :: error
CHARACTER(LEN=BIG_STRING_SIZE) :: IndividualPath
CHARACTER,ALLOCATABLE :: Paths(:)
INTEGER     xFileId, xGroupId
INTEGER            nMeshGroups, nMaxPathLength, nGridGroups
INTEGER            FileUnit, StartLoc, nStatus, i, j
REAL               Version

xFileId  = NONE
xGroupId = NONE

  ! Open the XMDF file
call XF_OPEN_FILE(a_XmdfFile, .TRUE., xFileId, nStatus)
if (nStatus < 0) then
  call XF_CLOSE_FILE(xFileId, error)
  error = -1
  return
endif

  ! open the status file
FileUnit = 53
OPEN(UNIT=FileUnit, FILE=a_OutFile, STATUS='REPLACE', ACTION='WRITE', &
     IOSTAT = error)
if (FileUnit == 0) then
  call XF_CLOSE_FILE(xFileId, error)
  error = -1
  return
endif

WRITE(FileUnit,*) 'File ', a_XmdfFile, ' opened.'

  ! write the version number to the file
call XF_GET_LIBRARY_VERSION_FILE(xFileId, Version, error)
WRITE(FileUnit,*) 'XMDF Version: ', Version
WRITE(FileUnit,*) ''

  ! Read Coordinate System Informatioin to the .txt file if contained in
  ! file, if not skip
call TXI_TEST_COORD_SYSTEM(xFileId, FileUnit, nStatus)
WRITE(FileUnit,*) ''

  ! read all datasets not beneath a mesh, grid, or cross-sections
call TD_READ_DATASETS(xFileId,FileUnit, nStatus)
if (nStatus < 0) then
  call XF_CLOSE_FILE(xFileId, error)
  error = -1
  return
endif

  ! Get the number and paths of datasets in the file.
call XF_GRP_PTHS_SZ_FOR_MESHES(xFileId, nMeshGroups, &
                                       nMaxPathLength, error)
if (error >= 0 .AND. nMeshGroups > 0) then
  allocate (Paths(nMaxPathLength*nMeshGroups))
  call XF_GET_GROUP_PATHS_FOR_MESHES(xFileId, nMeshGroups, nMaxPathLength, &
                                     Paths, error)
endif

if (error < 0) then
  call XF_CLOSE_FILE(xFileId, error)
  error = -1
  return
endif

  ! Report the number and paths to individual meshes in the file.
WRITE(FileUnit,*) 'Number of meshes in file: ', nMeshGroups
WRITE(FileUnit,*) 'Paths:'
do i=1, nMeshGroups
  do j=1, nMaxPathLength-1
    IndividualPath(j:j) = Paths((i-1)*nMaxPathLength+j)
  enddo
  WRITE(FileUnit,*) IndividualPath(1:nMaxPathLength-1)
enddo

WRITE(FileUnit,*) ''

  ! Open each mesh group
!if (nMeshGroups > 0) allocate(IndividualPath(nMaxPathLength + 1))

do i=1, nMeshGroups
  ! copy the portion of the array where a single path is stored
  StartLoc = (i-1)*nMaxPathLength + 1
  IndividualPath = ''
  do j = 1, nMaxPathLength - 1
    IndividualPath(j:j) = Paths(StartLoc+j-1)
  enddo
     
  WRITE(FileUnit,*) 'Reading mesh in group: ', &
                     IndividualPath(1:nMaxPathLength-1)
  call XF_OPEN_GROUP(xFileId, IndividualPath(1:LEN_TRIM(IndividualPath)), &
            xGroupId, nStatus)
  if (nStatus >= 0) then
    call TM_READ_MESH(xGroupId, FileUnit, nStatus)
  endif
  if (nStatus < 0) then
    WRITE(*,*) 'Error reading mesh..'
  endif
enddo

if (allocated(Paths)) deallocate(Paths)
!if (allocated(IndividualPath)) deallocate(IndividualPath)

  ! Grid stuff
call XF_GRP_PTHS_SZ_FOR_GRIDS(xFileId, nGridGroups, &
                                       nMaxPathLength, nStatus)
if (nStatus >= 0 .AND. nGridGroups > 0) then
  allocate (Paths(nMaxPathLength*nGridGroups))
  call XF_GET_GROUP_PATHS_FOR_GRIDS(xFileId, nGridGroups, &
                                    nMaxPathLength, Paths, nStatus)
endif
if (nStatus < 0) then
  call XF_CLOSE_FILE(xFileId, error)
  error = -1
  return
endif

  ! Report the number and paths to individual meshes in the file.
WRITE(FileUnit,*) 'Number of grids in file: ', nGridGroups
WRITE(FileUnit,*) 'Paths:'
do i=1, nGridGroups
  do j=1, nMaxPathLength-1
    IndividualPath(j:j) = Paths((i-1)*nMaxPathLength+j)
  enddo
  WRITE(FileUnit,*) IndividualPath(1:LEN_TRIM(IndividualPath))
enddo

WRITE(FileUnit,*) ''

!if (nGridGroups > 0) allocate(IndividualPath(nMaxPathLength + 1))

  ! Open each grid group
do i=1, nGridGroups
  do j = 1, nMaxPathLength - 1
    IndividualPath(j:j) = Paths((i-1)*nMaxPathLength+j)
  enddo
  WRITE(FileUnit,*) 'Reading grid in group: ', &
                     IndividualPath(1:LEN_TRIM(IndividualPath))
  call XF_OPEN_GROUP(xFileId, IndividualPath(1:LEN_TRIM(IndividualPath)), &
                     xGroupId, nStatus)
  if (nStatus >= 0) then
    call TG_READ_GRID(xGroupId, FileUnit, nStatus)
  endif
  if (nStatus < 0) then
     WRITE(FileUnit,*) 'Error reading grid..'
  endif
enddo

if (allocated(Paths)) deallocate(Paths)
!if (allocated(IndividualPath)) deallocate(IndividualPath)
  
  ! TODO do grid, and cross-section stuff.
   
  ! close the files
call XF_CLOSE_FILE(xFileId, error)
CLOSE(FileUnit)
  
return

END SUBROUTINE

!-----------------------------------------------------------------------------
! SUBROUTINE TXI_TestCalendar
! PURPOSE    Check the Calculations of Julian date from calendar date or  
!            vice-versa.
! NOTES      era is #defined (use #defines): ERA_IS_BCE (BC), ERA_IS_CE (AD)
!-----------------------------------------------------------------------------
SUBROUTINE TXI_TEST_CALENDAR (error)
  INTEGER, INTENT(OUT) :: error
  INTEGER  era1, yr1, mo1, day1, hr1, min1, sec1
  INTEGER  era2, yr2, mo2, day2, hr2, min2, sec2
  INTEGER  era3, yr3, mo3, day3, hr3, min3, sec3, FileUnit
  INTEGER  era4, yr4, mo4, day4, hr4, min4, sec4, calendarworks
  DOUBLE PRECISION  julian1, julian2, julian3, julian4

  calendarworks = 0
  FileUnit = 53
  OPEN(UNIT=FileUnit, FILE=CALENDAR_OUT_F, STATUS='REPLACE', ACTION='WRITE', &
     IOSTAT = error)

  WRITE(FileUnit,*) 'Calendar conversion:'

  era1 = ERA_IS_BCE
  yr1  = 0
  mo1  = 0
  day1 = 0
  hr1  = 0
  min1 = 0
  sec1 = 0
  julian1 = 2655.5
  call XF_JULIAN_TO_CALENDAR(era1, yr1, mo1, day1, hr1, min1, sec1, julian1, error)

  era2 = ERA_IS_BCE
  yr2  = 4706
  mo2  = 4
  day2 = 10
  hr2  = 0
  min2 = 0
  sec2 = 0
  julian2 = 0.0
  call XF_CALENDAR_TO_JULIAN(era2, yr2, mo2, day2, hr2, min2, sec2, julian2, error)

  era3 = ERA_IS_CE
  yr3  = 2004
  mo3  = 6
  day3 = 3
  hr3  = 2
  min3 = 8
  sec3 = 32
  julian3 = 0.0
  call XF_CALENDAR_TO_JULIAN(era3, yr3, mo3, day3, hr3, min3, sec3, julian3, error)

  era4 = ERA_IS_BCE
  yr4  = 0
  mo4  = 0
  day4 = 0
  hr4  = 0
  min4 = 0
  sec4 = 0
  julian4 = 2453159.58926_double
  call XF_JULIAN_TO_CALENDAR(era4, yr4, mo4, day4, hr4, min4, sec4, julian4, error)

WRITE(FileUnit,*) ''
WRITE(FileUnit,*) 'Dates #1 & #2  were calculated with the same date:'
WRITE(FileUnit,*) ''
WRITE(FileUnit,*) era1, '/', yr1, '/', mo1, '/', day1
WRITE(FileUnit,*) '', hr1, ':', min1, ':',sec1, '--- julian =',julian1
WRITE(FileUnit,*) ''
WRITE(FileUnit,*) era2, '/', yr2, '/', mo2, '/', day2
WRITE(FileUnit,*) '', hr2, ':', min2, ':',sec2, '--- julian =',julian2
WRITE(FileUnit,*) ''
WRITE(FileUnit,*) 'Dates #3 & #4  were calculated with the same date:'
WRITE(FileUnit,*) ''
WRITE(FileUnit,*) era3, '/', yr3, '/', mo3, '/', day3
WRITE(FileUnit,*) '', hr3, ':', min3, ':',sec3, '--- julian =',julian3
WRITE(FileUnit,*) ''
WRITE(FileUnit,*) era4, '/', yr4, '/', mo4, '/', day4
WRITE(FileUnit,*) '', hr4, ':', min4, ':',sec4, '--- julian =',julian4

  if (era1==era2 .AND. era3==era4) then
    if (yr1==yr2 .AND. yr3==yr4) then
          if (mo1==mo2 .AND. mo3==mo4) then
            if (day1==day2 .AND. day3==day4) then
                  if (hr1==hr2 .AND. hr3==hr4) then
                    if (min1==min2 .AND. min3==min4) then
              if (julian1==julian2 .AND. (julian3-.0000001)<julian4 .AND. (julian3+.0000001)>julian4) then
          WRITE(*,*) ''
          WRITE(*,*) 'Calendar conversion works correctly.'
        calendarworks = 1
                      endif
            endif
          endif
        endif
      endif
    endif
  endif
  if (calendarworks .NE. 1) then
    WRITE(*,*) 'Calendar Stuff DOES NOT Work Correctly'
  error = -1
  else
    error = 1
  endif
  return

END SUBROUTINE
!------------------------------------------------------------------------------
! FUNCTION TXI_TEST_MULTI_DATASETS
! PURPOSE  To test the newly translated functions
! NOTES
!------------------------------------------------------------------------------
SUBROUTINE TXI_TEST_MULTI_DATASETS

!call XF_SETUP_TO_WRITE_DATASETS(a_Filename, a_MultiDatasetsGroupPath, &
!                        a_PathInMultiDatasetsGroup, a_SpatialDataObjectGuid, &
 !                       a_OverwriteOptions, a_FileId, a_GroupId, error)
!XF_OPEN_MULTI_DATASETS_GROUP
!XF_DELETE_GROUP

END SUBROUTINE

!------------------------------------------------------------------------------
! SUBROUTINE TXI_WRITE_XMDF_VERSION
! PURPOSE    Write the XMDF version number to the screen
! NOTES
!------------------------------------------------------------------------------
SUBROUTINE TXI_TEST_VERSION ()
INTEGER  error
REAL     Version

  WRITE(*,*) ''
  call XF_GET_LIBRARY_VERSION(Version, error)
  WRITE(*,*) 'The current version of XMDF is: ', Version
  WRITE(*,*) ''

  return

END SUBROUTINE

!------------------------------------------------------------------------------
! SUBROUTINE TXI_TEST_COORD_SYSTEM
! PURPOSE    Reads a file's Coordinate Group and prints coordinate data out
!            to each text file.
! NOTES
!------------------------------------------------------------------------------
SUBROUTINE TXI_TEST_COORD_SYSTEM (xFileId, a_OutFile, error)
INTEGER, INTENT(IN) :: xFileId
INTEGER, INTENT(IN)        :: a_OutFile
INTEGER, INTENT(OUT)       :: error
INTEGER   iHorizDatum, iHorizUnits, iVertDatum, iVertUnits
INTEGER   iLat, iLon, iUtmZone, iSpcZone, iHpgnArea, iEllipse
INTEGER   bHorizDatum, nStatus
REAL(DOUBLE)        dCppLat, dCppLon, dMajorR, dMinorR
INTEGER      xCoordId
CHARACTER(LEN=BIG_STRING_SIZE)    strHorizUnits, strVertDatum
CHARACTER(LEN=BIG_STRING_SIZE)    strVertUnits

  ! Coordinate stuff
  ! Read Coordinate Info
    ! Open the Coordinate group
  !call H5Eset_auto_f(0, error)
  call XF_OPEN_COORDINATE_GROUP(xFileId, xCoordId, nStatus)
  if (nStatus < 0) then
    WRITE(a_OutFile,*) ''
    WRITE(a_OutFile,*) 'Coordinate Group not found'
    WRITE(a_OutFile,*) ''
    error = nStatus
    return
  endif
  !call H5Eset_auto_f(1, error)

  WRITE(a_OutFile,*) ''
  WRITE(a_OutFile,*) 'Coordinate System:'

  call XF_GET_HORIZ_DATUM(xCoordId, iHorizDatum, bHorizDatum)
  call XF_GET_HORIZ_UNITS(xCoordId, iHorizUnits, error)
  call XF_GET_VERT_DATUM(xCoordId, iVertDatum, error)
  call XF_GET_VERT_UNITS(xCoordId, iVertUnits, error)

    ! set horizontal units
  if (iHorizUnits == 0) then
    strHorizUnits = 'Horizontal units = US Survey Feet (=0)'
  else if (iHorizUnits == 1) then
    strHorizUnits = 'Horizontal units = International Feet (=1)'
  else if (iHorizUnits == 2) then
    strHorizUnits = 'Horizontal units = Meters (=2)'
  else
    strHorizUnits = 'ERROR in reading Horizontal units'
  endif

    ! set vertical datum
  if (iVertDatum == 0) then
    strVertDatum = 'Vertical datum = Local (=0)'
  else if (iVertDatum == 1) then
    strVertDatum = 'Vertical datum = NGVD 29 (=1)'
  else if (iVertDatum == 2) then
    strVertDatum = 'Vertical datum = NGVD 88 (=2)'
  else
    strVertDatum = 'ERROR in reading the Vertical datum'
  endif

    ! set vertocal units
  if (iVertUnits == 0) then
    strVertUnits = 'Vertical units = US Survey Feet (=0)'
  else if (iVertUnits == 1) then
    strVertUnits = 'Vertical units = International Feet (=1)'
  else if (iVertUnits == 2) then
    strVertUnits = 'Vertical units = Meters (=2)'
  else
    strVertUnits = 'ERROR in reading the Vertical units'
  endif

  if (bHorizDatum >= 0) then
    SELECT CASE (iHorizDatum)
      CASE (HORIZ_DATUM_GEOGRAPHIC)
        call XF_GET_ELLIPSE(xCoordId, iEllipse, error)
        call XF_GET_LAT(xCoordId, iLat, error)
        call XF_GET_LON(xCoordId, iLon, error)
          ! Write Horizontal and Vertical Info
        WRITE(a_OutFile,*) 'Horizontal datum = Geographic'
        WRITE(a_OutFile,*) 'Horizontal units = ', strHorizUnits(1:LEN_TRIM(strHorizUnits))
        WRITE(a_OutFile,*) 'Vertical datum = ', strVertDatum(1:LEN_TRIM(strVertDatum))
        WRITE(a_OutFile,*) 'Vertical units = ', strVertUnits(1:LEN_TRIM(strVertUnits))
          ! Write Latitude data
        if (iLat == 0) then
          WRITE(a_OutFile,*) '  Latitude = North (=0)'
        else if (iLat == 1) then
          WRITE(a_OutFile,*) '  Latitude = South (=1)'
        else
          WRITE(a_OutFile,*) '  LATITUDE INFO INCORRECT'
        endif
          ! Write Longitude data
        if (iLon == 0) then
          WRITE(a_OutFile,*) '  Longitude = East (=0)'
        else if (iLon == 1) then
          WRITE(a_OutFile,*) '  Longitude = West (=1)'
        else
          WRITE(a_OutFile,*) '  LONGITUDE INFO INCORRECT'
        endif
          ! Ellipse Information
          ! User-defined Ellipse (==32)
        if (iEllipse == 32) then
          WRITE(a_OutFile,*) 'Ellipse = User-defined:'
          call XF_GET_MAJOR_R(xCoordId, dMajorR, error)
          call XF_GET_MINOR_R(xCoordId, dMinorR, error)
          WRITE(a_OutFile,*) '  MajorR = ', dMajorR
          WRITE(a_OutFile,*) '  MinorR = ', dMinorR
        else
          WRITE(a_OutFile,*) 'Ellipse = ', iEllipse
        endif
        WRITE(a_OutFile,*) ''
      CASE (HORIZ_DATUM_UTM)
        call XF_GET_UTM_ZONE(xCoordId, iUtmZone, error)
          ! output info to text file
        if (iHorizDatum == HORIZ_DATUM_UTM) then
          WRITE(a_OutFile,*) 'Horizontal datum = UTM'
        else if (iHorizDatum == HORIZ_DATUM_UTM_NAD27) then
          WRITE(a_OutFile,*) 'Horizontal datum = UTM NAD27 (US)'
        else
          WRITE(a_OutFile,*) 'Horizontal datum = UTM NAD83 (US)'
        endif
        WRITE(a_OutFile,*) 'Horizontal units = ', &
                        strHorizUnits(1:LEN_TRIM(strHorizUnits))
        WRITE(a_OutFile,*) 'Vertical datum = ', &
                        strVertDatum(1:LEN_TRIM(strVertDatum))
        WRITE(a_OutFile,*) 'Vertical units = ', &
                        strVertUnits(1:LEN_TRIM(strVertUnits))
        WRITE(a_OutFile,*) 'UTM Zone = ', iUtmZone
        WRITE(a_OutFile,*) ''

      CASE (HORIZ_DATUM_UTM_NAD27, HORIZ_DATUM_UTM_NAD83)
        call XF_GET_UTM_ZONE(xCoordId, iUtmZone, error)
          ! output info to text file
        if (iHorizDatum == HORIZ_DATUM_UTM) then
          WRITE(a_OutFile,*) 'Horizontal datum = UTM'
        else if (iHorizDatum == HORIZ_DATUM_UTM_NAD27) then
          WRITE(a_OutFile,*) 'Horizontal datum = UTM NAD27 (US)'
        else
          WRITE(a_OutFile,*) 'Horizontal datum = UTM NAD83 (US)'
        endif
        WRITE(a_OutFile,*) 'Horizontal units = ', &
                        strHorizUnits(1:LEN_TRIM(strHorizUnits))
        WRITE(a_OutFile,*) 'Vertical datum = ', &
                        strVertDatum(1:LEN_TRIM(strVertDatum))
        WRITE(a_OutFile,*) 'Vertical units = ', &
                        strVertUnits(1:LEN_TRIM(strVertUnits))
        WRITE(a_OutFile,*) 'UTM Zone = ', iUtmZone
        WRITE(a_OutFile,*) ''
      CASE (HORIZ_DATUM_STATE_PLANE_NAD27, HORIZ_DATUM_STATE_PLANE_NAD83)
        call XF_GET_SPC_ZONE(xCoordId, iSpcZone, error)
          ! output info to text file
        if (iHorizDatum == HORIZ_DATUM_STATE_PLANE_NAD27) then
          WRITE(a_OutFile,*) 'Horizontal datum = State Plane NAD27 (US)'
        else 
          WRITE(a_OutFile,*) 'Horizontal datum = State Plane NAD83 (US)'
        endif
        WRITE(a_OutFile,*) 'Horizontal units = ', &
                        strHorizUnits(1:LEN_TRIM(strHorizUnits))
        WRITE(a_OutFile,*) 'Vertical datum = ', &
                        strVertDatum(1:LEN_TRIM(strVertDatum))
        WRITE(a_OutFile,*) 'Vertical units = ', &
                        strVertUnits(1:LEN_TRIM(strVertUnits))
        WRITE(a_OutFile,*) 'SPC Zone = ', iSpcZone
        WRITE(a_OutFile,*) ''
      CASE (HORIZ_DATUM_UTM_HPGN, HORIZ_DATUM_STATE_PLANE_HPGN, &
                                  HORIZ_DATUM_GEOGRAPHIC_HPGN)
        call XF_GET_HPGN_AREA(xCoordId, iHpgnArea, error)
        if (iHorizDatum == HORIZ_DATUM_UTM_HPGN) then
          WRITE(a_OutFile,*) 'Horizontal datum = UTM HPGN (US)'
        else if (iHorizDatum == HORIZ_DATUM_STATE_PLANE_HPGN) then
          WRITE(a_OutFile,*) 'Horizontal datum = State Plane HPGN (US)'
        else
          WRITE(a_OutFile,*) 'Horizontal datum = Geographic HPGN (US)'
        endif
        WRITE(a_OutFile,*) 'Horizontal units = ', &
                        strHorizUnits(1:LEN_TRIM(strHorizUnits))
        WRITE(a_OutFile,*) 'Vertical datum = ', &
                        strVertDatum(1:LEN_TRIM(strVertDatum))
        WRITE(a_OutFile,*) 'Vertical units = ', &
                        strVertUnits(1:LEN_TRIM(strVertUnits))
        WRITE(a_OutFile,*) 'HPGN Area = ', iHpgnArea
        WRITE(a_OutFile,*) ''
      CASE (HORIZ_DATUM_CPP)
        call XF_GET_CPP_LAT(xCoordId, dCppLat, error)
        call XF_GET_CPP_LON(xCoordId, dCppLon, error)
        WRITE(a_OutFile,*) 'Horizontal datum = CPP (Carte Parallelo-Grammatique Projection)'
        WRITE(a_OutFile,*) 'Horizontal units = ', &
                        strHorizUnits(1:LEN_TRIM(strHorizUnits))
        WRITE(a_OutFile,*) 'Vertical datum = ', &
                        strVertDatum(1:LEN_TRIM(strVertDatum))
        WRITE(a_OutFile,*) 'Vertical units = ', &
                        strVertUnits(1:LEN_TRIM(strVertUnits))
        WRITE(a_OutFile,*) 'CPP Latitude = ', dCppLat
        WRITE(a_OutFile,*) 'CPP Longitude = ', dCppLon
        WRITE(a_OutFile,*) ''
    CASE (HORIZ_DATUM_LOCAL, HORIZ_DATUM_GEOGRAPHIC_NAD27, &
          HORIZ_DATUM_GEOGRAPHIC_NAD83)
          ! do other systems
        if (iHorizDatum == HORIZ_DATUM_LOCAL) then
          WRITE(a_OutFile,*) 'Horizontal datum = Local'
        else if (iHorizDatum == HORIZ_DATUM_GEOGRAPHIC_NAD27) then
          WRITE(a_OutFile,*) 'Horizontal datum = Geographic NAD27 (US)'
        else
          WRITE(a_OutFile,*) 'Horizontal datum = Geographic NAD83 (US)'
    endif
        WRITE(a_OutFile,*) 'Horizontal units = ', &
                        strHorizUnits(1:LEN_TRIM(strHorizUnits))
        WRITE(a_OutFile,*) 'Vertical datum = ', &
                        strVertDatum(1:LEN_TRIM(strVertDatum))
        WRITE(a_OutFile,*) 'Vertical units = ', &
                        strVertUnits(1:LEN_TRIM(strVertUnits))
        WRITE(a_OutFile,*) ''
      CASE DEFAULT
          WRITE(a_OutFile,*) 'ERROR: The coordinate information is not found in the .h5 file'
          error = -1
      return
    END SELECT
  else
    WRITE(a_OutFile,*) 'Coordinate information in HDF5 file is incomplete.'
    WRITE(a_OutFile,*) ''
  endif

  call XF_CLOSE_GROUP(xCoordId, error)
  xCoordId = 0

  return

END SUBROUTINE

SUBROUTINE TXI_TEST_GEOMETRIC_PATHS(error)
  INTEGER, INTENT(OUT) :: error
  INTEGER                 compression

  compression = -1

    ! test writing a geometric path file */
  WRITE(*,*)''
  WRITE(*,*)'Writing geometric path data'
  WRITE(*,*)''
 
  call TM_WRITE_TEST_PATHS(GEOMPATH_A_FILE_F, compression, error);
  if (error < 0) then
    WRITE(*,*) 'Error writing geometric path data A'
    return
  endif
  WRITE(*,*) 'Finished writing geometric path data A'

    ! test reading a geometric path file */
  call TM_READ_TEST_PATHS(GEOMPATH_A_FILE_F, GEOMPATH_A_FILE_F_OUT, error)

  return

END SUBROUTINE TXI_TEST_GEOMETRIC_PATHS

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

END MODULE TestsModule

PROGRAM TESTS

USE TestDatasets
USE Xmdf
USE TestMesh
USE TestDefs
USE TestsModule

INTEGER    error

call XF_INITIALIZE(error)

  ! test the dataset routines
call TXI_TEST_TIMESTEPS(error)
if (error < 0) then
  WRITE(*,*) 'Error in writing timesteps!'
endif

  ! test the dataset routines
call TXI_TEST_DATASETS(error)
if (error < 0) then
  WRITE(*,*) 'Error in writing datasets!'
endif

  ! test overwriting datasets
call TXI_TEST_OVERWRITE_DSETS(error)
if (error < 0) then
  WRITE(*,*) 'Error in overwriting datasets!'
else
  WRITE(*,*) 'Finished writing datasets...'
endif

  ! test mesh stuff
call TXI_TEST_MESHS(error)

  ! test grid stuff
call TXI_TEST_GRIDS(error)

  ! test c in fortran
call TXI_TEST_C(error)

  ! test calendar stuff
call TXI_TEST_CALENDAR(error)

  ! test version
call TXI_TEST_VERSION

   ! test geometric paths
call TXI_TEST_GEOMETRIC_PATHS(error)

call XF_CLOSE (error)

PAUSE 'Press ENTER to exit...'

END PROGRAM