! This file is largely based on the example by the HDF Group
! http://www.hdfgroup.org/ftp/HDF5/examples/parallel/dataset.f90
! Modifications Pierre de Buyl 2014

PROGRAM DATASET

  USE HDF5 ! This module contains all necessary modules 
  use mpi
  IMPLICIT NONE

  CHARACTER(LEN=10), PARAMETER :: filename = "sds.h5"  ! File name
  CHARACTER(LEN=8), PARAMETER :: dsetname = "IntArray" ! Dataset name
  CHARACTER(LEN=14), PARAMETER :: author = "Pierre de Buyl" ! Dataset name

  INTEGER(HID_T) :: file_id       ! File identifier 
  INTEGER(HID_T) :: dset_id       ! Dataset identifier 
  INTEGER(HID_T) :: filespace     ! Dataspace identifier in file 
  INTEGER(HID_T) :: memspace     ! Dataspace identifier in file 
  INTEGER(HID_T) :: plist_id      ! Property list identifier 
  INTEGER(HID_T) :: aplist_id      ! Property list identifier 

  INTEGER(HSIZE_T), DIMENSION(2) :: dimsf = (/5,8/) ! Dataset dimensions.
  INTEGER(HSIZE_T), DIMENSION(2) :: dimsfi
  INTEGER(HSIZE_T), DIMENSION(2) :: count, offset
  integer(hsize_t), dimension(1) :: adims

  INTEGER(HID_T)  :: filetype, attr_id, aspace_id

  INTEGER, ALLOCATABLE :: data(:,:)   ! Data to write
  INTEGER :: rank = 2 ! Dataset rank 

  INTEGER :: error, error_n  ! Error flags
  INTEGER :: i, j
  !
  ! MPI definitions and calls.
  !
  INTEGER :: mpierror       ! MPI error flag
  INTEGER :: comm, info
  INTEGER :: mpi_size, mpi_rank
  comm = MPI_COMM_WORLD
  info = MPI_INFO_NULL
  CALL MPI_INIT(mpierror)
  CALL MPI_COMM_SIZE(comm, mpi_size, mpierror)
  CALL MPI_COMM_RANK(comm, mpi_rank, mpierror) 
  ! 
  ! Initialize data buffer with trivial data.
  !
  ALLOCATE ( data(dimsf(1),dimsf(2)))
  do i = 1, dimsf(2)
     do j = 1, dimsf(1)
        data(j,i) = j - 1 + (i-1)*dimsf(1) + 100*mpi_rank
     enddo
  enddo
  !
  ! Initialize FORTRAN interface
  !
  CALL h5open_f(error) 

  ! 
  ! Setup file access property list with parallel I/O access.
  !
  CALL h5pcreate_f(H5P_FILE_ACCESS_F, plist_id, error)
  CALL h5pset_fapl_mpio_f(plist_id, comm, info, error)

  !
  ! Create the file collectively.
  ! 
  CALL h5fcreate_f(filename, H5F_ACC_TRUNC_F, file_id, error, access_prp = plist_id)
  CALL h5pclose_f(plist_id, error)

  !
  ! Create the data space for the  dataset. 
  !
  dimsf(2) = 8*mpi_size
  CALL h5screate_simple_f(rank, dimsf, filespace, error)

  !
  ! Create the dataset with default properties.
  !
  CALL h5dcreate_f(file_id, dsetname, H5T_NATIVE_INTEGER, filespace, &
       dset_id, error)
  !
  ! Create property list for collective dataset write
  !
  CALL h5pcreate_f(H5P_DATASET_XFER_F, plist_id, error) 
  CALL h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, error)


  ! write author attribute
  call h5screate_f(H5S_SCALAR_F, aspace_id, error)
  adims(1) = len(author)
  call h5tcopy_f(H5T_NATIVE_CHARACTER, filetype, error)
  call h5tset_size_f(filetype, adims(1), error)
  call h5acreate_f(dset_id, 'author', filetype, aspace_id, attr_id, error)
  call h5awrite_f(attr_id, filetype, author, adims, error)
  call h5aclose_f(attr_id, error)
  call h5sclose_f(aspace_id, error)
  call h5tclose_f(filetype, error)

  count = (/ 5, 8 /)
  offset = (/ 0, mpi_rank*8 /)

  CALL h5screate_simple_f(rank, count, memspace, error)

  call H5Sselect_hyperslab_f(filespace, H5S_SELECT_SET_F, offset, count, error)
  !
  ! Write the dataset collectively. 
  !
  CALL h5dwrite_f(dset_id, H5T_NATIVE_INTEGER, data, dimsfi, error, &
       mem_space_id=memspace, file_space_id=filespace, xfer_prp = plist_id)
  !
  ! Deallocate data buffer.
  !
  DEALLOCATE(data)

  !
  ! Close resources.
  !
  CALL h5sclose_f(memspace, error)
  CALL h5sclose_f(filespace, error)
  CALL h5dclose_f(dset_id, error)
  CALL h5pclose_f(plist_id, error)
  CALL h5fclose_f(file_id, error)

  !
  ! Close FORTRAN interface
  !
  CALL h5close_f(error)

  CALL MPI_FINALIZE(mpierror)

END PROGRAM DATASET