module hdf5io
  use hdf5
  use mpi
  use parallel
  use mesh
  implicit none

  private
  public :: init_output, fin_output, output

  character(len=100), parameter :: filename = "data/iotest2.h5"
  character(len=100), dimension(nd), parameter ::               &
       & varname = (/"Rho", " Vx", " Vy", " Vz", " Bx", " By", " Bz", " Pr" /)

  integer(HID_T) :: file_id  ! File identifier
  integer(HID_T) :: plist_id ! Property list identifier

  integer :: ierr

contains


  subroutine init_output

    ! Initialize HDF5 library and Fortran interfaces.

    call h5open_f (ierr)

    ! Setup file access property list with parallel I/O access.
 
    call h5pcreate_f (H5P_FILE_ACCESS_F, plist_id, ierr)
    call h5pset_fapl_mpio_f (plist_id, cart3d%comm, MPI_INFO_NULL, ierr)

    ! Create the file collectively.

    call h5fcreate_f (filename, H5F_ACC_TRUNC_F, file_id, ierr, access_prp = plist_id)
    call h5pclose_f (plist_id, ierr)
    
    return
  end subroutine init_output


  subroutine fin_output

    ! Close the file.

    call h5fclose_f (file_id, ierr)

    ! Close HDF5 library and Fortran interfaces.

    call h5close_f (ierr)

    return
  end subroutine fin_output


  subroutine output (v, t, nout, tout, dtio)

    ! v             :   dumped data in each processor
    ! nout          :   I/O step
    ! t             :   current simulation time
    ! tout          :   I/O simulation time
    ! dtio          :   I/O time interval

    real(8), intent(in)    :: v(nd,nx,ny,nz), t, dtio
    real(8), intent(inout) :: tout
    integer, intent(inout) :: nout

    character(len=100) :: stepname, dsetname
    integer(HID_T)     :: step_id
    integer(HID_T)     :: dset_id
    integer(HID_T)     :: filespace
    integer(HID_T)     :: memspace

    integer(HSIZE_T), dimension(ndim) :: dimsf
    integer(HSIZE_T), dimension(ndim) :: chunk_dims
    integer(HSIZE_T), dimension(ndim) :: offset
    integer(HSIZE_T), dimension(ndim) :: counts
    integer(HSIZE_T), dimension(ndim) :: stride
    integer(HSIZE_T), dimension(ndim) :: blocks

    integer :: d, is, ie, js, je, ks, ke
    
    ! Create the timestep group for the specific I/O step.

    write(stepname, '("TimeStep", i4.4)') nout
    call h5gcreate_f (file_id, stepname, step_id, ierr)

    ! Each process defines dataset in memory and file.

    dimsf(1) = nxg
    dimsf(2) = nyg
    dimsf(3) = nzg
    chunk_dims(1) = dnx(cart3d%coords(1)+1)
    chunk_dims(2) = dny(cart3d%coords(2)+1)
    chunk_dims(3) = dnz(cart3d%coords(3)+1)

    stride(:) = 1
    counts(:) = 1
    blocks(:) = chunk_dims(:)
    offset(1) = sum(dnx(1:cart3d%coords(1)))
    offset(2) = sum(dny(1:cart3d%coords(2)))
    offset(3) = sum(dnz(1:cart3d%coords(3)))

    is = ng + 1
    ie = nx - ng
    js = ng + 1
    je = ny - ng
    ks = ng + 1
    ke = nz - ng

    do d = 1, nd

       ! Create the data space for the dataset.

       call h5screate_simple_f (ndim, dimsf, filespace, ierr)
       call h5screate_simple_f (ndim, chunk_dims, memspace, ierr)

       ! Create chunked dataset.

       dsetname = trim(varname(d))
       call h5pcreate_f (H5P_DATASET_CREATE_F, plist_id, ierr)
       call h5pset_chunk_f (plist_id, ndim, chunk_dims, ierr)
       call h5dcreate_f (step_id, dsetname, H5T_NATIVE_DOUBLE, filespace, &
            &            dset_id, ierr, plist_id)
       call h5sclose_f (filespace, ierr)

       ! Select hyperslab in the file.

       call h5dget_space_f (dset_id, filespace, ierr)
       call h5sselect_hyperslab_f (filespace, H5S_SELECT_SET_F, offset, counts, ierr, &
            &                      stride, blocks)

       ! 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)

       ! Write the dataset collectively.

       call h5dwrite_f (dset_id, H5T_NATIVE_DOUBLE, v(d,is:ie,js:je,ks:ke), &
            &           chunk_dims, ierr, file_space_id = filespace,        &
            &           mem_space_id = memspace, xfer_prp = plist_id        )

       ! Close dataspaces.

       call h5sclose_f (filespace, ierr)
       call h5sclose_f (memspace, ierr)

       ! Close the dataset.

       call h5dclose_f (dset_id, ierr)

       ! Close the property list.

       call h5pclose_f (plist_id, ierr)

    end do

    ! Close the timestep group.

    call h5gclose_f (step_id, ierr)


    if (ranks%me .eq. master_process) then
       write(*,*) "I/O number", nout, ",  simulation time", t
    end if

    nout = nout + 1
    tout = tout + dtio

    return
  end subroutine output


end module hdf5io