! $Id$
module hdf5_read_module

  ! Error codes returned by routines in this module: 0 indicates success, non-zero 
  ! indicates a problem occurred. The meaning of an error code is specific to the 
  ! procedure that generated it.

  
  use hdf5 
  use iso_c_binding

  implicit none
   
  private  
  integer, parameter, private:: int1 = selected_int_kind(1)
  integer, parameter, private:: int2 = selected_int_kind(3)
  integer, parameter:: int4 = selected_int_kind(8)
  integer, parameter, private:: int8 = selected_int_kind(18)
  integer, parameter :: real4 = selected_real_kind(6,37)
  integer, parameter :: real8 = selected_real_kind(15,307)
  integer, parameter :: ipre = real4 
   
   
    ! specification of a 1d array subset
  type :: subset_1d
     integer(kind=int8) :: start
     integer(kind=int8) :: stride
     integer(kind=int8) :: edge
  end type subset_1d
  
    ! specification of a 2d array subset
  type :: subset_2d
     integer(kind=int8), dimension(2) :: start
     integer(kind=int8), dimension(2) :: stride
     integer(kind=int8), dimension(2) :: edge
  end type subset_2d 
   
  

  public :: hdf5_open_file, hdf5_close_file
  public :: hdf5_group_exists
  public :: hdf5_open_dset, hdf5_get_dims, hdf5_read_dset, hdf5_close_dset
  public :: hdf5_read_attr
  public :: HDF5_MAX_RANK

  interface hdf5_read_attr
     module procedure &
          hdf5_read_attr_string, &
          hdf5_read_attr_int1
  end interface hdf5_read_attr

  interface hdf5_read_dset
     module procedure &
          hdf5_read_dset_real_1d, &
          hdf5_read_dset_real_2d, &
          hdf5_read_dset_int_2d
  end interface hdf5_read_dset

  interface hdf5_get_dims
     module procedure &
          hdf5_get_dims_id, &
          hdf5_get_dims_name
  end interface hdf5_get_dims
         
  ! would like to use H5S_MAX_RANK from HDF5 library, but there doesn't seem to 
  ! be a Fortran equivalent. The value 32 is based on H5S_MAX_RANK defined in H5public.h
  ! in HDF5 v1.8.9.
  integer(kind=int4), parameter :: HDF5_MAX_RANK = 32

  logical, save :: hdf5_initialized = .false.

contains

  ! returns non-zero status on failure
  ! FUTURE: add optional argument flags to allow read or readwrite
  subroutine hdf5_open_file(filename, file_id, status)
    CHARACTER(LEN=*), INTENT(IN) :: filename
    INTEGER(KIND=INT4), intent(out) :: file_id
    INTEGER(KIND=INT4), intent(out) :: status

    ! initialize HDF5 library if needed
    call init_lib_if_needed(status)
    if(status < 0) return

    ! open file
    call h5fopen_f(filename, H5F_ACC_RDONLY_F, file_id, status)    
    if(status < 0) then
       status = -1
       return
    end if

    status = 0

  end subroutine hdf5_open_file


  subroutine hdf5_open_dset(file_id, dset_fullpath, dset_id, status)
    INTEGER(KIND=INT4), intent(in) :: file_id
    CHARACTER(LEN=*), INTENT(in) :: dset_fullpath
    INTEGER(KIND=INT4), intent(out) :: dset_id
    INTEGER(KIND=INT4), intent(out) :: status    

    ! initialize HDF5 library if needed
    call init_lib_if_needed(status)
    if(status < 0) return

    CALL h5dopen_f (file_id, dset_fullpath, dset_id, status)
    if(status < 0) then
       status = -1
       return
    end if

    status = 0

  end subroutine hdf5_open_dset

  ! read a "scalar" string attribute (actually a 1x1 array)
  subroutine hdf5_read_attr_string(obj_id, attr_name, attr_value, status)
    INTEGER(KIND=INT4), intent(in) :: obj_id
    CHARACTER(LEN=*), INTENT(in) :: attr_name
    CHARACTER(LEN=*), intent(out) :: attr_value
    INTEGER(KIND=INT4), intent(out) :: status

    INTEGER(SIZE_T) :: str_len
    INTEGER(HID_T)  :: attr
    INTEGER(HID_T)  :: filetype
    INTEGER(SIZE_T) :: size
    INTEGER(HID_T)  :: space
    INTEGER(HID_T)  :: memtype
    TYPE(C_PTR)     :: f_ptr

    ! viirs scalar string attributes are 1x1 arrays
    CHARACTER(LEN=LEN(attr_value)), dimension(1,1), target :: buf

    ! initialize HDF5 library if needed
    call init_lib_if_needed(status)
    if(status < 0) return

    str_len = len(attr_value)

    ! open attribute
    CALL h5aopen_f(obj_id, attr_name, attr, status)
    if(status < 0) then
       status = -1
       return
    end if

    ! get attribute type
    CALL H5Aget_type_f(attr, filetype, status)    
    if(status < 0) then
       CALL H5Aclose_f(attr, status)
       status = -2
       return
    end if

    ! get size
    CALL H5Tget_size_f(filetype, size, status)
    if(status < 0) then
       CALL H5Aclose_f(attr, status)
       CALL H5Tclose_f(filetype, status)
       status = -3
       return
    end if

    ! make sure string is big enough. assuming size includes null-terminator
    IF( size .GT. str_len + 1) THEN
       status = -4
       return
    ENDIF

    ! handle null-terminated string with length 0
    if( size <= 1) then
       attr_value(:) = ' '
       status = 0
       return
    end if

    ! get dataspace
    CALL H5Aget_space_f(attr, space, status)
    if(status < 0) then
       CALL H5Aclose_f(attr, status)
       CALL H5Tclose_f(filetype, status)
       status = -5
       return
    end if

    ! create the memory datatype.
    CALL H5Tcopy_f(H5T_FORTRAN_S1, memtype, status)
    if(status < 0) then
       CALL H5Aclose_f(attr, status)
       CALL H5Tclose_f(filetype, status)
       CALL H5Sclose_f(space, status)
       status = -6
       return
    end if

    CALL H5Tset_size_f(memtype, str_len, status)
    if(status < 0) then
       CALL H5Aclose_f(attr, status)
       CALL H5Tclose_f(filetype, status)
       CALL H5Sclose_f(space, status)
       CALL H5Tclose_f(memtype, status)
       status = -7
       return
    end if

    ! read the data.
    f_ptr = C_LOC(buf(1,1)(1:1))
    CALL H5Aread_f(attr, memtype, f_ptr, status)
    if(status < 0) then
       CALL H5Aclose_f(attr, status)
       CALL H5Tclose_f(filetype, status)
       CALL H5Sclose_f(space, status)
       CALL H5Tclose_f(memtype, status)
       status = -8
       return
    end if
    
    ! copy string to caller's buffer, padded with spaces
    attr_value(:) = ' '
    attr_value(1:size-1) = buf(1,1)(1:size-1)

    ! deallocate resources
    CALL H5Aclose_f(attr, status)
    CALL H5Tclose_f(filetype, status)
    CALL H5Sclose_f(space, status)
    CALL H5Tclose_f(memtype, status)

    status = 0

  end subroutine hdf5_read_attr_string

  ! read a "scalar" string attribute (actually a 1x1 array)
  subroutine hdf5_read_attr_int1(obj_id, attr_name, attr_value, status)
    INTEGER(KIND=INT4), intent(in) :: obj_id
    CHARACTER(LEN=*), INTENT(in) :: attr_name
    INTEGER(KIND=INT1), intent(out), target :: attr_value
    INTEGER(KIND=INT4), intent(out) :: status

    INTEGER(HID_T)  :: attr
    INTEGER(HID_T)  :: filetype
    INTEGER(SIZE_T) :: actual_size
    INTEGER(HID_T)  :: space
    INTEGER(HID_T)  :: memtype
    TYPE(C_PTR)     :: f_ptr

    ! initialize HDF5 library if needed
    call init_lib_if_needed(status)
    if(status < 0) return

    ! open attribute
    CALL h5aopen_f(obj_id, attr_name, attr, status)
    if(status < 0) then
       status = -1
       return
    end if

    ! get attribute type
    CALL H5Aget_type_f(attr, filetype, status)    
    if(status < 0) then
       CALL H5Aclose_f(attr, status)
       status = -2
       return
    end if

    ! get size
    CALL H5Tget_size_f(filetype, actual_size, status)
    if(status < 0) then
       CALL H5Aclose_f(attr, status)
       CALL H5Tclose_f(filetype, status)
       status = -3
       return
    end if

    ! get dataspace
    CALL H5Aget_space_f(attr, space, status)
    if(status < 0) then
       CALL H5Aclose_f(attr, status)
       CALL H5Tclose_f(filetype, status)
       status = -4
       return
    end if

    ! create the memory datatype. Specifying little endian, but endianness should not 
    ! matter for 8-bit unsigned.
    CALL H5Tcopy_f(H5T_STD_U8LE , memtype, status)         
    if(status < 0) then
       CALL H5Aclose_f(attr, status)
       CALL H5Tclose_f(filetype, status)
       CALL H5Sclose_f(space, status)
       status = -5
       return
    end if
 
    ! read the data.
    f_ptr = C_LOC(attr_value)
    CALL H5Aread_f(attr, memtype, f_ptr, status)
    if(status < 0) then
       CALL H5Aclose_f(attr, status)
       CALL H5Tclose_f(filetype, status)
       CALL H5Sclose_f(space, status)
       CALL H5Tclose_f(memtype, status)
       status = -6
       return
    end if
    
    ! deallocate resources
    CALL H5Aclose_f(attr, status)
    CALL H5Tclose_f(filetype, status)
    CALL H5Sclose_f(space, status)
    CALL H5Tclose_f(memtype, status)

    status = 0

  end subroutine hdf5_read_attr_int1

  ! FUTURE: there is a lot of duplicated code among the various read routines. Would like
  ! to approximate generic programming in some not-too-painful way... maybe via code generation?
  
  ! arr must be pre-allocated to at least the required size. 
  ! Warning: assumes Data has already been initialized to missing values
  ! assumes hyperslab (possible FUTURE: handle non-hyperslab by detecting from data).
  subroutine hdf5_read_dset_real_1d(file_id, dset_fullpath, arr, fill_value, fill_value_exists, status, subset)
    INTEGER(KIND=INT4), intent(in) :: file_id
    character(len=*), intent(in) :: dset_fullpath
    real(kind=real4), dimension(:), intent(inout) :: arr
    real(kind=real4), intent(out) :: fill_value
    logical, intent(out) :: fill_value_exists
    integer(kind=int4), intent(out) :: status
    type(subset_1d), intent(in), optional :: subset

    integer(kind=int4), parameter :: RANK = 1
    INTEGER(HSIZE_T), dimension(1), parameter :: ZERO_OFFSET = (/0/)

    INTEGER(HID_T)  :: dset_id
    INTEGER(HID_T)  :: fspace_id
    INTEGER(HSIZE_T), dimension(1) :: dims
    INTEGER(HSIZE_T), dimension(1) :: max_dims
    INTEGER(HID_T)  :: mspace_id
    INTEGER(HSIZE_T), dimension(1) :: buf_shape
    INTEGER(HID_T)  :: plist_id
    type(subset_1d) :: subset_

    ! initialize HDF5 library if needed
    call init_lib_if_needed(status)
    if(status < 0) return

    ! open dataset
    call h5dopen_f (file_id, dset_fullpath, dset_id, status)
    if(status < 0) then
       status = -1
       return
    end if

    ! get data space
    call h5dget_space_f(dset_id, fspace_id, status)
    if(status < 0) then
       call h5dclose_f(dset_id, status)
       status = -2
       return
    end if

    ! get dims
    call h5sget_simple_extent_dims_f(fspace_id, dims, max_dims, status)
    if(status < 0) then
       call h5sclose_f(fspace_id, status)
       call h5dclose_f(dset_id, status)
       status = -3
       return
    end if

    ! if subset was not specified, assume we want the whole thing
    if(present(subset)) then
       subset_ = subset
    else
       subset_%start = 1
       subset_%edge = dims(1)
       subset_%stride = 1
    end if

    ! check that the supplied array is large enough to hold the requested subset
    buf_shape = shape(arr) 
    if(buf_shape(1) < subset_%edge) then
       status = -4
       return
    end if
    
    ! check that dims are large enough for specified subset
    if(dims(1) < subset_%start + ((subset_%edge - 1) * subset_%stride)) then
       call h5sclose_f(fspace_id, status)
       call h5dclose_f(dset_id, status)
       status = -5
       return
    end if

    ! select hyperslab in the dataset. Note that start is specified as an offset (i.e. 0-indexed)
    CALL h5sselect_hyperslab_f(fspace_id, H5S_SELECT_SET_F, (/subset_%start - 1/), (/subset_%edge/), status, stride=(/subset_%stride/))
    if(status < 0) then
       call h5sclose_f(fspace_id, status)
       call h5dclose_f(dset_id, status)
       status = -6
       return
    endif

    ! create memory dataspace.
    CALL h5screate_simple_f(RANK, (/subset_%edge/), mspace_id, status)
    if(status < 0) then
       call h5sclose_f(fspace_id, status)
       call h5dclose_f(dset_id, status)
       status = -7
       return
    end if

    ! select hyperslab in memory.
    CALL h5sselect_hyperslab_f(mspace_id, H5S_SELECT_SET_F, ZERO_OFFSET, (/subset_%edge/), status)
    if(status < 0) then
       call h5sclose_f(mspace_id, status)
       call h5sclose_f(fspace_id, status)
       call h5dclose_f(dset_id, status)
       status = -8
       return
    end if

    ! read data
    call h5dread_f(dset_id, H5T_NATIVE_REAL, arr, buf_shape, status, mem_space_id=mspace_id, file_space_id=fspace_id)
    if(status < 0) then
       call h5sclose_f(mspace_id, status)
       call h5sclose_f(fspace_id, status)
       call h5dclose_f(dset_id, status)
       status = -9
       return
    end if

    ! get the dataset create property list 
    CALL h5dget_create_plist_f(dset_id, plist_id, status)
    if(status < 0) then
       call h5sclose_f(mspace_id, status)
       call h5sclose_f(fspace_id, status)
       call h5dclose_f(dset_id, status)
       status = -10
       return
    end if
    
    ! get the fill value. Negative status indicates there is no fill value defined,
    ! or some other error
    call h5pget_fill_value_f(plist_id, H5T_NATIVE_REAL, fill_value, status)
    if(status < 0) then
       fill_value_exists = .false.
    else
       fill_value_exists = .true.
!       print *, "fill_value = ", fill_value
    end if
    
    ! free resources
    call h5pclose_f(plist_id, status)
    call h5sclose_f(mspace_id, status)
    call h5sclose_f(fspace_id, status)
    call h5dclose_f(dset_id, status)

    status = 0
    
  end subroutine hdf5_read_dset_real_1d



  ! arr must be pre-allocated to at least the required size. 
  ! Warning: assumes Data has already been initialized to missing values
  ! assumes hyperslab (possible FUTURE: handle non-hyperslab by detecting from data)
  subroutine hdf5_read_dset_real_2d(file_id, dset_fullpath, arr, fill_value, fill_value_exists, status, subset)
    INTEGER(KIND=INT4), intent(in) :: file_id
    character(len=*), intent(in) :: dset_fullpath
    real(kind=real4), dimension(:,:), intent(inout) :: arr
    real(kind=real4), intent(out) :: fill_value
    logical, intent(out) :: fill_value_exists
    integer(kind=int4), intent(out) :: status
    type(subset_2d), intent(in), optional :: subset

    integer(kind=int4), parameter :: RANK = 2
    INTEGER(HSIZE_T), DIMENSION(2), parameter :: ZERO_OFFSET = (/0, 0/)

    INTEGER(HID_T)  :: dset_id
    INTEGER(HID_T)  :: fspace_id
    INTEGER(HSIZE_T), DIMENSION(2) :: dims
    INTEGER(HSIZE_T), DIMENSION(2) :: max_dims
    INTEGER(HID_T)  :: mspace_id
    INTEGER(HSIZE_T), DIMENSION(2) :: buf_shape
    INTEGER(HID_T)  :: plist_id
    type(subset_2d) :: subset_

    ! initialize HDF5 library if needed
    call init_lib_if_needed(status)
    if(status < 0) return

    ! open dataset
    call h5dopen_f (file_id, dset_fullpath, dset_id, status)
    if(status < 0) then
       status = -1
       return
    end if

    ! get data space
    call h5dget_space_f(dset_id, fspace_id, status)
    if(status < 0) then
       call h5dclose_f(dset_id, status)
       status = -2
       return
    end if

    ! get dims
    call h5sget_simple_extent_dims_f(fspace_id, dims, max_dims, status)
    if(status < 0) then
       call h5sclose_f(fspace_id, status)
       call h5dclose_f(dset_id, status)
       status = -3
       return
    end if

    ! if subset was not specified, assume we want the whole thing
    if(present(subset)) then
       subset_ = subset
    else
       subset_%start = 1
       subset_%edge = dims
       subset_%stride = 1
    end if

    ! check that the supplied array is large enough to hold the requested subset
    buf_shape = shape(arr)
    if(any(buf_shape < subset_%edge)) then
       status = -4
       return
    end if

    ! check that dims are large enough for specified subset
    if(any(dims < subset_%start + ((subset_%edge - 1) * subset_%stride))) then
       call h5sclose_f(fspace_id, status)
       call h5dclose_f(dset_id, status)
       status = -5
       return
    end if

    ! select hyperslab in the dataset. Note that start is specified as an offset (i.e. 0-indexed)
    CALL h5sselect_hyperslab_f(fspace_id, H5S_SELECT_SET_F, subset_%start - 1, subset_%edge, status, stride=subset_%stride)
    if(status < 0) then
       call h5sclose_f(fspace_id, status)
       call h5dclose_f(dset_id, status)
       status = -6
       return
    endif

    ! create memory dataspace.
    CALL h5screate_simple_f(RANK, subset_%edge, mspace_id, status)
    if(status < 0) then
       call h5sclose_f(fspace_id, status)
       call h5dclose_f(dset_id, status)
       status = -7
       return
    end if

    ! select hyperslab in memory.
    CALL h5sselect_hyperslab_f(mspace_id, H5S_SELECT_SET_F, ZERO_OFFSET, subset_%edge, status)
    if(status < 0) then
       call h5sclose_f(mspace_id, status)
       call h5sclose_f(fspace_id, status)
       call h5dclose_f(dset_id, status)
       status = -8
       return
    end if

    ! read data
    call h5dread_f(dset_id, H5T_NATIVE_REAL, arr, buf_shape, status, mem_space_id=mspace_id, file_space_id=fspace_id)
    if(status < 0) then
       call h5sclose_f(mspace_id, status)
       call h5sclose_f(fspace_id, status)
       call h5dclose_f(dset_id, status)
       status = -9
       return
    end if

    ! get the dataset create property list 
    CALL h5dget_create_plist_f(dset_id, plist_id, status)
    if(status < 0) then
       call h5sclose_f(mspace_id, status)
       call h5sclose_f(fspace_id, status)
       call h5dclose_f(dset_id, status)
       status = -10
       return
    end if
    
    ! get the fill value. Negative status indicates there is no fill value defined,
    ! or some other error
    call h5pget_fill_value_f(plist_id, H5T_NATIVE_REAL, fill_value, status)
    if(status < 0) then
       fill_value_exists = .false.
    else
       fill_value_exists = .true.
!       print *, "fill_value = ", fill_value
    end if
    
    ! free resources
    call h5pclose_f(plist_id, status)
    call h5sclose_f(mspace_id, status)
    call h5sclose_f(fspace_id, status)
    call h5dclose_f(dset_id, status)

    status = 0
    
  end subroutine hdf5_read_dset_real_2d

  ! arr must be pre-allocated to at least the required size. 
  ! Warning: assumes Data has already been initialized to missing values
  ! assumes hyperslab (possible FUTURE: handle non-hyperslab by detecting from data).
  ! Should be able to read any integer dataset, as long as HDF5 can map the file datatype
  ! to int4.
  subroutine hdf5_read_dset_int_2d(file_id, dset_fullpath, arr, fill_value, fill_value_exists, status, subset)
    INTEGER(KIND=INT4), intent(in) :: file_id
    character(len=*), intent(in) :: dset_fullpath
    integer(kind=int4), dimension(:,:), intent(inout) :: arr
    integer(kind=int4), intent(out) :: fill_value
    logical, intent(out) :: fill_value_exists
    integer(kind=int4), intent(out) :: status
    type(subset_2d), intent(in), optional :: subset


    integer(kind=int4), parameter :: RANK = 2
    INTEGER(HSIZE_T), DIMENSION(2), parameter :: ZERO_OFFSET = (/0, 0/)

    INTEGER(HID_T)  :: dset_id
    INTEGER(HID_T)  :: fspace_id
    INTEGER(HSIZE_T), DIMENSION(2) :: dims
    INTEGER(HSIZE_T), DIMENSION(2) :: max_dims
    INTEGER(HID_T)  :: mspace_id
    INTEGER(HSIZE_T), DIMENSION(2) :: buf_shape
    INTEGER(HID_T)  :: plist_id
    type(subset_2d) :: subset_

    ! initialize HDF5 library if needed
    call init_lib_if_needed(status)
    if(status < 0) return

    ! open dataset
    call h5dopen_f (file_id, dset_fullpath, dset_id, status)
    if(status < 0) then
       status = -1
       return
    end if

    ! get data space
    call h5dget_space_f(dset_id, fspace_id, status)
    if(status < 0) then
       call h5dclose_f(dset_id, status)
       status = -2
       return
    end if

    ! get dims
    call h5sget_simple_extent_dims_f(fspace_id, dims, max_dims, status)
    if(status < 0) then
       call h5sclose_f(fspace_id, status)
       call h5dclose_f(dset_id, status)
       status = -3
       return
    end if

    ! if subset was not specified, assume we want the whole thing
    if(present(subset)) then
       subset_ = subset
    else
       subset_%start = 1
       subset_%edge = dims
       subset_%stride = 1
    end if

    ! check that the supplied array is large enough to hold the requested subset
    buf_shape = shape(arr) 
    if(any(buf_shape < subset_%edge)) then
       status = -4
       return
    end if
    
    ! check that dims are large enough for specified subset
    if(any(dims < subset_%start + ((subset_%edge - 1) * subset_%stride))) then
       call h5sclose_f(fspace_id, status)
       call h5dclose_f(dset_id, status)
       status = -5
       return
    end if

    ! select hyperslab in the dataset. Note that start is specified as an offset (i.e. 0-indexed)
    CALL h5sselect_hyperslab_f(fspace_id, H5S_SELECT_SET_F, subset_%start - 1, subset_%edge, status, stride=subset_%stride)
    if(status < 0) then
       call h5sclose_f(fspace_id, status)
       call h5dclose_f(dset_id, status)
       status = -6
       return
    endif

    ! create memory dataspace.
    CALL h5screate_simple_f(RANK, subset_%edge, mspace_id, status)
    if(status < 0) then
       call h5sclose_f(fspace_id, status)
       call h5dclose_f(dset_id, status)
       status = -7
       return
    end if

    ! select hyperslab in memory.
    CALL h5sselect_hyperslab_f(mspace_id, H5S_SELECT_SET_F, ZERO_OFFSET, subset_%edge, status)
    if(status < 0) then
       call h5sclose_f(mspace_id, status)
       call h5sclose_f(fspace_id, status)
       call h5dclose_f(dset_id, status)
       status = -8
       return
    end if

    ! read data
    call h5dread_f(dset_id, H5T_NATIVE_INTEGER, arr, buf_shape, status, mem_space_id=mspace_id, file_space_id=fspace_id)
    if(status < 0) then
       call h5sclose_f(mspace_id, status)
       call h5sclose_f(fspace_id, status)
       call h5dclose_f(dset_id, status)
       status = -9
       return
    end if

    ! get the dataset create property list 
    CALL h5dget_create_plist_f(dset_id, plist_id, status)
    if(status < 0) then
       call h5sclose_f(mspace_id, status)
       call h5sclose_f(fspace_id, status)
       call h5dclose_f(dset_id, status)
       status = -10
       return
    end if
    
    ! get the fill value. Negative status indicates there is no fill value defined,
    ! or some other error
    call h5pget_fill_value_f(plist_id, H5T_NATIVE_INTEGER, fill_value, status)
    if(status < 0) then
       fill_value_exists = .false.
    else
       fill_value_exists = .true.
!       print *, "fill_value = ", fill_value
    end if
    
    ! free resources
    call h5pclose_f(plist_id, status)
    call h5sclose_f(mspace_id, status)
    call h5sclose_f(fspace_id, status)
    call h5dclose_f(dset_id, status)

    status = 0
    
  end subroutine hdf5_read_dset_int_2d


  subroutine hdf5_close_file(file_id, status)
    INTEGER(KIND=INT4), intent(in) :: file_id
    integer(kind=int4), intent(out) :: status

    ! initialize HDF5 library if needed
    call init_lib_if_needed(status)
    if(status < 0) return

    CALL h5fclose_f(file_id, status)
    if(status < 0) then
       status = -1
       return
    end if

    status = 0

  end subroutine hdf5_close_file

  subroutine hdf5_close_dset(dset_id, status)   
    INTEGER(KIND=INT4), intent(in) :: dset_id
    integer(kind=int4), intent(out) :: status
    
    ! initialize HDF5 library if needed
    call init_lib_if_needed(status)
    if(status < 0) return

    CALL h5dclose_f(dset_id, status)
    if(status < 0) then
       status = -1
       return
    end if

    status = 0
           
  end subroutine hdf5_close_dset

  subroutine hdf5_get_dims_name(file_id, dset_fullpath, dims, status)
    INTEGER(KIND=INT4), intent(in) :: file_id
    character(len=*), intent(in) :: dset_fullpath
    integer(kind=int8), dimension(:), intent(out) :: dims
    integer(kind=int4), intent(out) :: status

    INTEGER(HID_T)  :: dset_id

    ! initialize HDF5 library if needed
    call init_lib_if_needed(status)
    if(status < 0) return
   
    call hdf5_open_dset(file_id, dset_fullpath, dset_id, status)
    if(status < 0) then
       status = -1
       return    
    end if

    call hdf5_get_dims_id(dset_id, dims, status)
    if(status < 0) then
       call hdf5_close_dset(dset_id, status)
       status = -2
       return    
    end if

    call hdf5_close_dset(dset_id, status)    
    status = 0

  end subroutine hdf5_get_dims_name

  ! FUTURE: add support for attributes. Can check type of id via h5iget_type.
  ! Rank of 'dims' must be >= number of dimensions of the dataset. If number of dimensions 
  ! is not known, set rank to HDF5_MAX_RANK.
  subroutine hdf5_get_dims_id(dset_id, dims, status)
    INTEGER(KIND=INT4), intent(in) :: dset_id
    integer(kind=int8), dimension(:), intent(out) :: dims
    integer(kind=int4), intent(out) :: status

    INTEGER(HID_T)  :: dspace_id
    INTEGER(HSIZE_T), DIMENSION(HDF5_MAX_RANK) :: dims_hsize_t
    INTEGER(HSIZE_T), DIMENSION(HDF5_MAX_RANK) :: maxdims
    integer(kind=int4) :: rank

    ! initialize HDF5 library if needed
    call init_lib_if_needed(status)
    if(status < 0) return

    CALL H5Dget_space_f(dset_id, dspace_id, status)
    if(status < 0) then
       status = -1
       return    
    end if

    CALL H5Sget_simple_extent_dims_f(dspace_id, dims_hsize_t, maxdims, status)
    if(status < 0) then
       call h5sclose_f(dspace_id, status)
       status = -2
       return
    end if
    
    ! if hdf5 call succeeded, rank was returned via status variable. really. 
    rank = status

    ! rank cannot be greater than the max allowed rank. This would indicate an internal error.
    if(rank > HDF5_MAX_RANK) then
       call h5sclose_f(dspace_id, status)
       status = -3
       return
    end if

    ! check that rank of dims argument is sufficient
    if(rank > size(dims)) then
       call h5sclose_f(dspace_id, status)
       status = -4
       return
    end if

    dims(1:rank) = dims_hsize_t(1:rank)

    CALL H5Sclose_f(dspace_id, status)

    status = 0

  end subroutine hdf5_get_dims_id

  
  subroutine init_lib_if_needed(status)
    integer(kind=int4), intent(out) :: status

    ! initialize HDF5 library if needed
    if(.not. hdf5_initialized) then

       ! open library
       call h5open_f(status)
       if(status < 0) then 
          status = -50
          return
       end if

       ! disable error printing - this is needed so that existence checks do not generate errors
       call h5eset_auto_f(0, status)
       !if(status < 0) &
       !     print *, "Warning: failed to turn off HDF5 error printing. Ignore any error messages from HDF5"

       hdf5_initialized = .true.
    end if

    status = 0

  end subroutine init_lib_if_needed

  ! hdf5 lib does not offer a check for existence, so we try to open the group and check 
  ! the returned status. Error printing has already been disabled.
  function hdf5_group_exists(file_id, group_name) result(exists)
    INTEGER(KIND=INT4), intent(in) :: file_id
    character(len=*), intent(in) :: group_name
    logical :: exists

    integer(kind=int4) :: status
    integer(kind=int4) :: group_id

    exists = .false.

    ! initialize HDF5 library if needed
    call init_lib_if_needed(status)
    if(status < 0) return
    
    ! try to open group
    call h5gopen_f(file_id, trim(group_name), group_id, status)
    if(status == 0) then
       call h5gclose_f(group_id, status)
       exists = .true.
    end if

  end function hdf5_group_exists
  

end module hdf5_read_module