program example2G_bug

   use mpi
   use hdf5
   implicit none

   character(80)    :: fn
   integer :: myProc, nProcs, comm, istat
   integer :: ierr, i, lSz, m, iseed, info, nDim

   integer, parameter :: LocalSZ = 646
   integer, parameter :: Numvar = 1


   integer, parameter :: MAXDIM = 3

   type grid_t
      integer    :: nd
      integer(8) :: is(MAXDIM), ie(MAXDIM)
      integer(8) :: num(MAXDIM)
   end type grid_t

   integer(hid_t)   :: file_id      ! File    identifier
   integer(hid_t)   :: group_id     ! Group   identifier
   integer(hid_t)   :: dset_id      ! Dataset identifier
   integer(hid_t)   :: filespace    ! Dataspace identifier in file
   integer(hid_t)   :: memspace     ! Dataspace identifier in memory
   integer(hid_t)   :: plist_id     ! Property list identifier
   integer(HSIZE_t) :: sz(3),gsz(3), starts(3)
   integer, allocatable :: newSeed(:)
   real(8), allocatable :: u(:)

   type (grid_t) global, local

   !------------------------------------------------------------
   ! Start parallel

   comm = MPI_COMM_WORLD

   call MPI_init(ierr)
   call MPI_Comm_size(comm, nProcs, ierr)
   call MPI_Comm_rank(comm, myProc, ierr)

   nDim        = MAXDIM
   global % nd = MAXDIM
   do i = 1, global % nd 
      global % is(i)  = 1
      global % ie(i)  = LocalSZ
      global % num(i) = LocalSZ
   end do

   local = global

   !------------------------------------------------------------
   ! Form a random value using system random_seed
   
   call random_seed(size=m)
   allocate(newSeed(m))
   
   iseed = 12345678
   newSeed = iseed
   call random_seed(put=newSeed)

   !------------------------------------------------------------
   ! Compute global size of solution vector(s)

   lSz = 1
   do i = 1, local % nd
      starts(i)   = local % is(i) - 1
      sz(i)       = local % num(i)
      gsz(i)      = global % num(i)
      lSz         = lSz * local % num(i)
   end do
   

   info = MPI_INFO_NULL
   call MPI_Info_create(info, ierr)

   fn = "myfile.h5"
   if ( myProc == 0) then
      call MPI_File_delete(fn,MPI_INFO_NULL,ierr)
   end if

   
   !
   ! (1) Initialize FORTRAN predefined datatypes

   call H5open_f(ierr)
   
   !
   ! (2) Setup file access property list w/ parallel I/O access.

   call H5Pcreate_f(H5P_FILE_ACCESS_F,plist_id,ierr)
   call H5Pset_fapl_mpio_f(plist_id,  comm, info, ierr);

   !
   ! (3.0) Create the file collectively
   call H5Fcreate_f(fn, H5F_ACC_TRUNC_F, file_id, ierr, access_prp = plist_id)
   call H5Pclose_f(plist_id, ierr)
   

   ! (3.1) Create group
   call H5Gcreate_f(file_id,"Solution", group_id, ierr)

   allocate(u(lSz), stat = istat)

   if (istat /= 0) then
      print *, "unable to allocate soln vector u, istat: ",istat
   end if

         
   call init3d(local, u)

   !
   ! (4) Create the data space for the dataset: filespace, memspace
   call H5Screate_simple_f(nDim, gsz, filespace, ierr)
   
   !
   ! Each process defines dataset in memory and writes it to the hyperslab
   ! in the file.
   call H5Screate_simple_f(nDim, sz, memspace, ierr)

   !
   ! (5) Create the dataset with default properties.
   !
   CALL H5Dcreate_f(group_id, "T", H5T_NATIVE_DOUBLE, filespace, &
        dset_id, ierr)
   CALL H5Sclose_f(filespace, ierr)
   !
   !
   ! (6) Select hyperslab in the file.
   !
   CALL H5Dget_space_f(dset_id, filespace, ierr)
   CALL H5Sselect_hyperslab_f (filespace, H5S_SELECT_SET_F, starts, sz, ierr)


   !
   ! (7a) Create property list for collective dataset write
   call H5Pcreate_f(H5P_DATASET_XFER_F, plist_id, ierr)
   call H5Pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, ierr)

   
   call H5Dwrite_f(dset_id, H5T_NATIVE_DOUBLE, u, gsz, ierr, &
        file_space_id = filespace, mem_space_id = memspace, xfer_prp = plist_id)

   !
   ! (9) Close dataspaces.
   !
   call H5Sclose_f(filespace, ierr)
   call H5Sclose_f(memspace, ierr)

   !
   ! (10) Close the dataset and property list.
   !
   call H5Dclose_f(dset_id, ierr)
   call H5Pclose_f(plist_id, ierr)

   !
   ! (12) Close the group and file.
   !
   call H5Gclose_f(group_id, ierr)
   call H5Fclose_f(file_id,  ierr)

   !
   ! (12) Close FORTRAN predefined datatypes.
   !
   call H5Close_f(ierr)

   deallocate(u)

      
   !------------------------------------------------------------
   ! End parallel
   call MPI_Finalize(ierr)
end program example2G_bug

subroutine init3d(local, u)
   implicit none
   integer, parameter :: MAXDIM = 3
   type grid_t
      integer    :: nd
      integer(8) :: is(MAXDIM), ie(MAXDIM)
      integer(8) :: num(MAXDIM)
   end type grid_t

   type(grid_t) local
   integer :: kk, k, jj, j, ii, i, ja
   real(8) :: u(local % num(1)* local % num(2) * local % num(3))
   real(8) :: rval
   
   ja = 0
   kk = local % is(3) - 1
   do k = 1, local % num(3)
      kk = kk + 1
      jj = local % is(2) - 1
      do j = 1, local % num(2)
         jj = jj + 1
         ii = local % is(1) - 1
         do i = 1, local % num(1)
            call random_number(rval)
            ii = ii + 1
            ja = ja + 1
            u(ja) = dble(kk*10000 + jj*100 + ii)+rval
         end do
      end do
   end do

end subroutine init3d
