hdf5_file.F90

Go to the documentation of this file.

! Copyright (c) 2024-2026, The Neko Authors
! All rights reserved.
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions
! are met:
!
!   * Redistributions of source code must retain the above copyright
!     notice, this list of conditions and the following disclaimer.
!
!   * Redistributions in binary form must reproduce the above
!     copyright notice, this list of conditions and the following
!     disclaimer in the documentation and/or other materials provided
!     with the distribution.
!
!   * Neither the name of the authors nor the names of its
!     contributors may be used to endorse or promote products derived
!     from this software without specific prior written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
! FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
! COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
! INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
! BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
! LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
! ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
! POSSIBILITY OF SUCH DAMAGE.
!
module hdf5_file
  use num_types, only : rp, dp, sp, i8
  use generic_file, only : generic_file_t
  use checkpoint, only : chkp_t
  use checkpoint_payload, only : checkpoint_payload_t, checkpoint_array_t, &
       checkpoint_mesh_array_t
  use utils, only : neko_error, neko_warning, filename_suffix_pos, &
       filename_split
  use mesh, only : mesh_t
  use field, only : field_t, field_ptr_t
  use field_list, only : field_list_t
  use field_series, only : field_series_t, field_series_ptr_t
  use dofmap, only : dofmap_t
  use space, only : neko_space_t => space_t, gll
  use interpolation, only : interpolator_t
  use logger, only : neko_log, log_size, neko_log_debug
  use vector, only : vector_t
  use matrix, only : matrix_t
  use datadist, only : linear_dist_t
  use global_interpolation, only : global_interpolation_t
  use math, only : rzero
  use comm, only : pe_rank, pe_size, neko_comm, mpi_real_precision
  use mpi_f08, only : mpi_info_null, mpi_allreduce, mpi_allgather, &
       mpi_in_place, mpi_integer, mpi_sum, mpi_max, mpi_comm_size, mpi_exscan, &
       mpi_barrier, mpi_integer8, mpi_scan, mpi_bcast, &
       mpi_double_precision
  use hdf5_session, only : hdf5_session_init, hdf5_session_finalize
#ifdef HAVE_HDF5
  use hdf5
#endif
  implicit none
  private
 
  type, public, extends(generic_file_t) :: hdf5_file_t
 
     ! HDF5 members
#ifdef HAVE_HDF5
     integer(hid_t) :: file_id = -1_hid_t
     integer(hid_t) :: active_group_id = -1_hid_t
     integer(hid_t) :: plist_id = -1_hid_t
#endif
     character(len=1) :: mode
     integer :: precision = -1
     integer :: offset = 0
     integer :: count = 0
 
   contains
     ! General methods for reading/writing HDF5 files
     procedure :: read => hdf5_file_read
     procedure :: write => hdf5_file_write
     procedure :: get_next_output_fname => hdf5_file_get_next_output_fname
     procedure :: set_overwrite => hdf5_file_set_overwrite
     ! Granular methods for dealing with HDF5 files
     procedure :: open => hdf5_file_open
     procedure :: close => hdf5_file_close
     procedure :: set_active_group => hdf5_file_set_group
     procedure :: set_precision => hdf5_file_set_precision
     procedure, pass(this) :: write_vector => hdf5_file_write_vector
     procedure, pass(this) :: write_matrix => hdf5_file_write_matrix
     procedure, pass(this) :: write_field => hdf5_file_write_field
     procedure, pass(this) :: write_int_attribute => &
          hdf5_file_write_int_attribute
     procedure, pass(this) :: write_rp_attribute => hdf5_file_write_rp_attribute
     procedure, pass(this) :: read_vector => hdf5_file_read_vector
     procedure, pass(this) :: read_matrix => hdf5_file_read_matrix
     procedure, pass(this) :: read_int_attribute => hdf5_file_read_int_attribute
     procedure, pass(this) :: read_rp_attribute => hdf5_file_read_rp_attribute
     procedure :: write_dataset => hdf5_file_write_dataset
     procedure :: read_dataset => hdf5_file_read_dataset
     procedure :: write_attribute => hdf5_file_write_attribute
     procedure :: read_attribute => hdf5_file_read_attribute
  end type hdf5_file_t
 
  type :: hdf5_checkpoint_layout_t
     type(mesh_t), pointer :: msh => null()
     type(neko_space_t) :: xh
     logical :: mesh2mesh = .false.
     type(dofmap_t) :: src_dof
     type(global_interpolation_t) :: global_interp
   contains
     procedure, pass(this) :: free => hdf5_checkpoint_layout_free
  end type hdf5_checkpoint_layout_t
 
contains
 
  function hdf5_file_get_next_output_fname(this) result(fname)
    class(hdf5_file_t), intent(in) :: this
    character(len=1024) :: fname
    character(len=5) :: id_str
    integer :: counter, suffix_pos
 
    if (this%overwrite) then
       fname = this%get_fname()
    else
       counter = this%get_counter()
       if (counter .eq. -1) then
          counter = this%get_start_counter()
       else
          counter = counter + 1
       end if
 
       fname = this%get_base_fname()
       suffix_pos = filename_suffix_pos(fname)
       write(id_str, '(i5.5)') counter
       fname = trim(fname(1:suffix_pos-1)) // id_str // fname(suffix_pos:)
    end if
 
  end function hdf5_file_get_next_output_fname
 
  subroutine hdf5_file_set_overwrite(this, overwrite)
    class(hdf5_file_t), intent(inout) :: this
    logical, intent(in) :: overwrite
    this%overwrite = overwrite
  end subroutine hdf5_file_set_overwrite
 
  function file_get_fname(this) result(base_fname)
    class(hdf5_file_t), intent(in) :: this
    character(len=1024) :: base_fname
    character(len=1024) :: fname
    character(len=1024) :: path, name, suffix
 
    fname = trim(this%get_base_fname())
    call filename_split(fname, path, name, suffix)
 
    ! Append a counter
    !write(base_fname, '(A,A,"_",I0,A)') &
    !     trim(path), trim(name), this%get_start_counter(), trim(suffix)
 
    ! Do not append anything
    base_fname = trim(fname)
 
  end function file_get_fname
 
  subroutine hdf5_file_set_precision(this, precision)
    class(hdf5_file_t), intent(inout) :: this
    integer, intent(in) :: precision
    this%precision = precision
  end subroutine hdf5_file_set_precision
 
  subroutine hdf5_checkpoint_layout_free(this)
    class(hdf5_checkpoint_layout_t), intent(inout) :: this
 
    if (this%mesh2mesh) then
       call this%global_interp%free()
       call this%src_dof%free()
    end if
    call this%Xh%free()
    nullify(this%msh)
    this%mesh2mesh = .false.
 
  end subroutine hdf5_checkpoint_layout_free
 
 
#ifdef HAVE_HDF5
 
  ! ===============
  ! General methods
  ! ===============
 
  subroutine hdf5_file_write(this, data, t)
    class(hdf5_file_t), intent(inout) :: this
    class(*), target, intent(in) :: data
    real(kind=dp), intent(in), optional :: t
    type(mesh_t), pointer :: msh
    type(dofmap_t), pointer :: dof
    type(field_ptr_t), allocatable :: fp(:)
    type(field_series_ptr_t), allocatable :: fsp(:)
    real(kind=dp), pointer :: dtlag(:)
    real(kind=dp), pointer :: tlag(:)
    integer :: ierr, info, drank, i, j
    integer(hid_t) :: plist_id, fapl_id
    integer(hid_t) :: file_id, dset_id, grp_id, attr_id
    integer(hid_t) :: filespace, dspace_id, memspace
    integer(hid_t) :: H5T_NEKO_REAL
    integer(hsize_t), dimension(1) :: ddim, dcount, doffset
    integer :: suffix_pos
    character(len=5) :: id_str
    character(len=1024) :: fname
    logical :: checkpoint_data
 
    call hdf5_file_determine_data(data, msh, dof, fp, fsp, dtlag, tlag)
    checkpoint_data = .false.
    select type (data)
    type is (chkp_t)
       checkpoint_data = .true.
    end select
 
    if (.not. this%overwrite) call this%increment_counter()
    fname = trim(this%get_fname())
 
    call hdf5_session_init()
 
    call hdf5_file_determine_real(h5t_neko_real)
 
    ! The file-access and dataset-transfer property lists are distinct
    ! objects; reusing one identifier for both would leak the first, and
    ! HDF5 then refuses to close the library at the end of the session.
    call h5pcreate_f(h5p_file_access_f, fapl_id, ierr)
    info = mpi_info_null%mpi_val
    call h5pset_fapl_mpio_f(fapl_id, neko_comm%mpi_val, info, ierr)
 
    call h5fcreate_f(fname, h5f_acc_trunc_f, &
         file_id, ierr, access_prp = fapl_id)
    call h5pclose_f(fapl_id, ierr)
 
    call h5pcreate_f(h5p_dataset_xfer_f, plist_id, ierr)
    call h5pset_dxpl_mpio_f(plist_id, h5fd_mpio_collective_f, ierr)
 
    call h5screate_f(h5s_scalar_f, filespace, ierr)
    ddim = 1
 
    if (present(t)) then
       call h5acreate_f(file_id, "Time", h5t_native_double, filespace, attr_id, &
            ierr, h5p_default_f, h5p_default_f)
       call h5awrite_f(attr_id, h5t_native_double, t, ddim, ierr)
       call h5aclose_f(attr_id, ierr)
    end if
 
    if (associated(dof)) then
       call h5acreate_f(file_id, "Lx", h5t_native_integer, filespace, attr_id, &
            ierr, h5p_default_f, h5p_default_f)
       call h5awrite_f(attr_id, h5t_native_integer, dof%Xh%lx, ddim, ierr)
       call h5aclose_f(attr_id, ierr)
    end if
 
    if (associated(msh)) then
       call h5gcreate_f(file_id, "Mesh", grp_id, ierr, &
            lcpl_id = h5p_default_f, gcpl_id = h5p_default_f, &
            gapl_id = h5p_default_f)
 
       call h5acreate_f(grp_id, "Elements", h5t_native_integer, filespace, &
            attr_id, ierr, h5p_default_f, h5p_default_f)
       call h5awrite_f(attr_id, h5t_native_integer, msh%glb_nelv, ddim, ierr)
       call h5aclose_f(attr_id, ierr)
 
       call h5acreate_f(grp_id, "Dimension", h5t_native_integer, filespace, &
            attr_id, ierr, h5p_default_f, h5p_default_f)
       call h5awrite_f(attr_id, h5t_native_integer, msh%gdim, ddim, ierr)
       call h5aclose_f(attr_id, ierr)
 
       call h5gclose_f(grp_id, ierr)
    end if
 
 
    call h5sclose_f(filespace, ierr)
 
    !
    ! Write restart group (tlag, dtlag)
    !
    if (associated(tlag) .and. associated(dtlag)) then
       call h5gcreate_f(file_id, "Restart", grp_id, ierr, &
            lcpl_id = h5p_default_f, gcpl_id = h5p_default_f, &
            gapl_id = h5p_default_f)
 
       drank = 1
       ddim = size(tlag)
       doffset(1) = 0
       if (pe_rank .eq. 0) then
          dcount = size(tlag)
       else
          dcount = 0
       end if
 
       call h5screate_simple_f(drank, ddim, filespace, ierr)
 
       call h5dcreate_f(grp_id, 'tlag', h5t_native_double, &
            filespace, dset_id, ierr)
       call h5dget_space_f(dset_id, dspace_id, ierr)
       call h5sselect_hyperslab_f (dspace_id, h5s_select_set_f, &
            doffset, dcount, ierr)
       call h5dwrite_f(dset_id, h5t_native_double, tlag, &
            ddim, ierr, xfer_prp = plist_id)
       call h5sclose_f(dspace_id, ierr)
       call h5dclose_f(dset_id, ierr)
 
       call h5dcreate_f(grp_id, 'dtlag', h5t_native_double, &
            filespace, dset_id, ierr)
       call h5dget_space_f(dset_id, dspace_id, ierr)
       call h5sselect_hyperslab_f (dspace_id, h5s_select_set_f, &
            doffset, dcount, ierr)
       call h5dwrite_f(dset_id, h5t_native_double, dtlag, &
            ddim, ierr, xfer_prp = plist_id)
       call h5sclose_f(dspace_id, ierr)
       call h5dclose_f(dset_id, ierr)
 
       call h5sclose_f(filespace, ierr)
       call h5gclose_f(grp_id, ierr)
 
    end if
 
 
    !
    ! Write fields group
    !
    if (checkpoint_data) then
       select type (data)
       type is (chkp_t)
          call hdf5_checkpoint_write_payloads( &
               file_id, plist_id, h5t_neko_real, data)
       end select
       if (allocated(fp)) deallocate(fp)
       if (allocated(fsp)) deallocate(fsp)
    else if (allocated(fp) .or. allocated(fsp)) then
       call h5gcreate_f(file_id, "Fields", grp_id, ierr, &
            lcpl_id = h5p_default_f, gcpl_id = h5p_default_f, &
            gapl_id = h5p_default_f)
 
       ! `lxyz` rather than `lx**3`: they agree in 3D, but a 2D case has
       ! `lz = 1`, and the reader derives its offsets the same way.
       dcount(1) = int(dof%size(), 8)
       doffset(1) = int(msh%offset_el, 8) * int(dof%Xh%lxyz, 8)
       ddim = int(dof%size(), 8)
       drank = 1
       call mpi_allreduce(mpi_in_place, ddim(1), 1, &
            mpi_integer8, mpi_sum, neko_comm, ierr)
 
       call h5screate_simple_f(drank, ddim, filespace, ierr)
       call h5screate_simple_f(drank, dcount, memspace, ierr)
 
 
       if (allocated(fp)) then
          do i = 1, size(fp)
             call h5dcreate_f(grp_id, fp(i)%ptr%name, h5t_neko_real, &
                  filespace, dset_id, ierr)
             call h5dget_space_f(dset_id, dspace_id, ierr)
             call h5sselect_hyperslab_f(dspace_id, h5s_select_set_f, &
                  doffset, dcount, ierr)
             call h5dwrite_f(dset_id, h5t_neko_real, &
                  fp(i)%ptr%x(1,1,1,1), &
                  ddim, ierr, file_space_id = dspace_id, &
                  mem_space_id = memspace, xfer_prp = plist_id)
             call h5sclose_f(dspace_id, ierr)
             call h5dclose_f(dset_id, ierr)
          end do
          deallocate(fp)
       end if
 
       if (allocated(fsp)) then
          do i = 1, size(fsp)
             do j = 1, fsp(i)%ptr%size()
                call h5dcreate_f(grp_id, fsp(i)%ptr%lf(j)%name, &
                     h5t_neko_real, filespace, dset_id, ierr)
                call h5dget_space_f(dset_id, dspace_id, ierr)
                call h5sselect_hyperslab_f(dspace_id, h5s_select_set_f, &
                     doffset, dcount, ierr)
                call h5dwrite_f(dset_id, h5t_neko_real, &
                     fsp(i)%ptr%lf(j)%x(1,1,1,1), &
                     ddim, ierr, file_space_id = dspace_id, &
                     mem_space_id = memspace, xfer_prp = plist_id)
                call h5sclose_f(dspace_id, ierr)
                call h5dclose_f(dset_id, ierr)
             end do
          end do
          deallocate(fsp)
       end if
 
       call h5gclose_f(grp_id, ierr)
       call h5sclose_f(filespace, ierr)
       call h5sclose_f(memspace, ierr)
    end if
 
    call h5pclose_f(plist_id, ierr)
    call h5fclose_f(file_id, ierr)
    call hdf5_session_finalize()
 
  end subroutine hdf5_file_write
 
  subroutine hdf5_file_read(this, data)
    class(hdf5_file_t) :: this
    class(*), target, intent(inout) :: data
    integer(hid_t) :: plist_id, fapl_id
    integer(hid_t) :: file_id, dset_id, grp_id, attr_id
    integer(hid_t) :: filespace, memspace
    integer(hid_t) :: H5T_NEKO_REAL
    integer(hsize_t), dimension(1) :: ddim, dcount, doffset
    integer :: i,j, ierr, info, glb_nelv, gdim, lx, drank
    type(mesh_t), pointer :: msh
    type(dofmap_t), pointer :: dof
    type(field_ptr_t), allocatable :: fp(:)
    type(field_series_ptr_t), allocatable :: fsp(:)
    type(hdf5_checkpoint_layout_t), target :: layout
    real(kind=dp), pointer :: dtlag(:)
    real(kind=dp), pointer :: tlag(:)
    real(kind=dp) :: t
    character(len=1024) :: fname
    ! Not LOG_SIZE: these messages are longer than 79 characters, and an
    ! overflowing internal write is a runtime error, not a truncation.
    character(len=256) :: err_msg
    logical :: payloads_exist
 
    fname = trim(this%get_fname())
 
    call hdf5_session_init()
 
    call hdf5_file_determine_data(data, msh, dof, fp, fsp, dtlag, tlag)
    call hdf5_file_determine_real(h5t_neko_real)
 
    ! As in the write path, keep the file-access and dataset-transfer
    ! property lists in separate identifiers so neither is leaked.
    call h5pcreate_f(h5p_file_access_f, fapl_id, ierr)
    info = mpi_info_null%mpi_val
    call h5pset_fapl_mpio_f(fapl_id, neko_comm%mpi_val, info, ierr)
 
    call h5fopen_f(fname, h5f_acc_rdonly_f, &
         file_id, ierr, access_prp = fapl_id)
    call h5pclose_f(fapl_id, ierr)
 
    call h5lexists_f(file_id, 'Payloads', payloads_exist, ierr)
 
    call h5pcreate_f(h5p_dataset_xfer_f, plist_id, ierr)
    call h5pset_dxpl_mpio_f(plist_id, h5fd_mpio_collective_f, ierr)
 
    ddim = 1
    call h5aopen_name_f(file_id, 'Time', attr_id, ierr)
    call h5aread_f(attr_id, h5t_native_double, t, ddim, ierr)
    call h5aclose_f(attr_id, ierr)
 
    select type (data)
    type is (chkp_t)
       data%t = t
    end select
 
    call h5aopen_name_f(file_id, 'Lx', attr_id, ierr)
    call h5aread_f(attr_id, h5t_native_integer, lx, ddim, ierr)
    call h5aclose_f(attr_id, ierr)
 
    call h5gopen_f(file_id, 'Mesh', grp_id, ierr, gapl_id = h5p_default_f)
 
    call h5aopen_name_f(grp_id, 'Elements', attr_id, ierr)
    call h5aread_f(attr_id, h5t_native_integer, glb_nelv, ddim, ierr)
    call h5aclose_f(attr_id, ierr)
 
    call h5aopen_name_f(grp_id, 'Dimension', attr_id, ierr)
    call h5aread_f(attr_id, h5t_native_integer, gdim, ddim, ierr)
    call h5aclose_f(attr_id, ierr)
    call h5gclose_f(grp_id, ierr)
 
    ! A mesh supplied through `case.restart_mesh_file` is the one the
    ! checkpoint was written on, and the file's layout follows it.
    layout%msh => msh
    layout%mesh2mesh = .false.
    select type (data)
    type is (chkp_t)
       if (allocated(data%previous_mesh%elements)) then
          layout%msh => data%previous_mesh
          layout%mesh2mesh = .true.
       end if
    end select
 
    if (gdim .ne. layout%msh%gdim) then
       write(err_msg, '(A,I0,A,I0,A)') 'HDF5 checkpoint is for a ', gdim, &
            'D mesh but this case is ', layout%msh%gdim, 'D'
       call neko_error(trim(err_msg))
    end if
 
    if (glb_nelv .ne. layout%msh%glb_nelv) then
       write(err_msg, '(A,I0,A,I0,A)') 'HDF5 checkpoint has ', glb_nelv, &
            ' elements but the mesh it is read on has ', &
            layout%msh%glb_nelv, '; set case.restart_mesh_file to the ' // &
            'mesh the checkpoint was written on'
       call neko_error(trim(err_msg))
    end if
 
    if (gdim .eq. 3) then
       call layout%Xh%init(gll, lx, lx, lx)
    else
       call layout%Xh%init(gll, lx, lx)
    end if
 
    select type (data)
    type is (chkp_t)
       if (gdim .eq. 3) then
          call data%previous_Xh%init(gll, lx, lx, lx)
       else
          call data%previous_Xh%init(gll, lx, lx)
       end if
 
       if (layout%mesh2mesh) then
          call layout%src_dof%init(layout%msh, layout%Xh)
          call layout%global_interp%init(layout%src_dof, neko_comm, &
               tol = data%mesh2mesh_tol)
          call layout%global_interp%find_points(dof%x%x, dof%y%x, dof%z%x, &
               dof%size())
       end if
    end select
 
    if (associated(tlag) .and. associated(dtlag)) then
       drank = 1
       ddim = size(tlag)
       doffset(1) = 0
       if (pe_rank .eq. 0) then
          dcount = size(tlag)
       else
          dcount = 0
       end if
 
       call h5gopen_f(file_id, 'Restart', grp_id, ierr, &
            gapl_id = h5p_default_f)
       call h5dopen_f(grp_id, 'tlag', dset_id, ierr)
       call h5dget_space_f(dset_id, filespace, ierr)
       call h5sselect_hyperslab_f (filespace, h5s_select_set_f, &
            doffset, dcount, ierr)
       call h5dread_f(dset_id, h5t_native_double, tlag, ddim, ierr, &
            xfer_prp = plist_id)
       call h5dclose_f(dset_id, ierr)
       call h5sclose_f(filespace, ierr)
 
       call h5dopen_f(grp_id, 'dtlag', dset_id, ierr)
       call h5dget_space_f(dset_id, filespace, ierr)
       call h5sselect_hyperslab_f (filespace, h5s_select_set_f, &
            doffset, dcount, ierr)
       call h5dread_f(dset_id, h5t_native_double, dtlag, ddim, ierr, &
            xfer_prp = plist_id)
       call h5dclose_f(dset_id, ierr)
       call h5sclose_f(filespace, ierr)
 
       call h5gclose_f(grp_id, ierr)
    end if
 
    if (payloads_exist) then
       select type (data)
       type is (chkp_t)
          call hdf5_checkpoint_read_payloads( &
               file_id, plist_id, h5t_neko_real, data, layout)
       class default
          call neko_error( &
               "HDF5 payload checkpoints require a checkpoint object")
       end select
       if (allocated(fp)) deallocate(fp)
       if (allocated(fsp)) deallocate(fsp)
    else if (allocated(fp) .or. allocated(fsp)) then
       call h5gopen_f(file_id, 'Fields', grp_id, ierr, gapl_id = h5p_default_f)
 
       if (allocated(fp)) then
          do i = 1, size(fp)
             call hdf5_checkpoint_read_field(grp_id, plist_id, &
                  h5t_neko_real, fp(i)%ptr, layout)
          end do
       end if
 
       if (allocated(fsp)) then
          do i = 1, size(fsp)
             do j = 1, fsp(i)%ptr%size()
                call hdf5_checkpoint_read_field(grp_id, plist_id, &
                     h5t_neko_real, fsp(i)%ptr%lf(j), layout)
             end do
          end do
       end if
       call h5gclose_f(grp_id, ierr)
    end if
 
    call h5pclose_f(plist_id, ierr)
    call h5fclose_f(file_id, ierr)
 
    call layout%free()
 
    call hdf5_session_finalize()
 
  end subroutine hdf5_file_read
 
  subroutine hdf5_checkpoint_write_payloads(file_id, plist_id, &
       h5_neko_real, chkp)
    integer(hid_t), intent(in) :: file_id, plist_id, h5_neko_real
    type(chkp_t), intent(in) :: chkp
    integer(hid_t) :: payloads_id, payload_id
    integer :: i, j, k, ierr
 
    call h5gcreate_f(file_id, "Payloads", payloads_id, ierr)
 
    do i = 1, chkp%payload_count()
       call hdf5_checkpoint_open_group(payloads_id, &
            chkp%payloads(i)%ptr%name, .true., payload_id)
 
       do j = 1, chkp%payloads(i)%ptr%field_count()
          call hdf5_checkpoint_write_field(payload_id, plist_id, &
               h5_neko_real, chkp%payloads(i)%ptr%fields(j)%ptr)
       end do
 
       do j = 1, chkp%payloads(i)%ptr%series_count()
          do k = 1, chkp%payloads(i)%ptr%series(j)%ptr%size()
             call hdf5_checkpoint_write_field(payload_id, plist_id, &
                  h5_neko_real, &
                  chkp%payloads(i)%ptr%series(j)%ptr%lf(k))
          end do
       end do
 
       do j = 1, chkp%payloads(i)%ptr%mesh_array_count()
          call hdf5_checkpoint_write_mesh_array(payload_id, plist_id, &
               h5_neko_real, chkp%payloads(i)%ptr%mesh_arrays(j)%ptr)
       end do
 
       do j = 1, chkp%payloads(i)%ptr%array_count()
          call hdf5_checkpoint_write_array(payload_id, plist_id, &
               h5_neko_real, chkp%payloads(i)%ptr%arrays(j)%ptr)
       end do
 
       call h5gclose_f(payload_id, ierr)
    end do
 
    call h5gclose_f(payloads_id, ierr)
 
  end subroutine hdf5_checkpoint_write_payloads
 
  subroutine hdf5_checkpoint_read_payloads(file_id, plist_id, &
       h5_neko_real, chkp, layout)
    integer(hid_t), intent(in) :: file_id, plist_id, h5_neko_real
    type(chkp_t), intent(inout) :: chkp
    type(hdf5_checkpoint_layout_t), target, intent(inout) :: layout
    integer(hid_t) :: payloads_id, payload_id
    integer :: i, j, k, ierr
 
    call h5gopen_f(file_id, "Payloads", payloads_id, ierr)
 
    do i = 1, chkp%payload_count()
       call hdf5_checkpoint_open_group(payloads_id, &
            chkp%payloads(i)%ptr%name, .false., payload_id)
 
       do j = 1, chkp%payloads(i)%ptr%field_count()
          call hdf5_checkpoint_read_field(payload_id, plist_id, &
               h5_neko_real, chkp%payloads(i)%ptr%fields(j)%ptr, &
               layout)
       end do
 
       do j = 1, chkp%payloads(i)%ptr%series_count()
          do k = 1, chkp%payloads(i)%ptr%series(j)%ptr%size()
             call hdf5_checkpoint_read_field(payload_id, plist_id, &
                  h5_neko_real, &
                  chkp%payloads(i)%ptr%series(j)%ptr%lf(k), layout)
          end do
       end do
 
       do j = 1, chkp%payloads(i)%ptr%mesh_array_count()
          call hdf5_checkpoint_read_mesh_array(payload_id, plist_id, &
               h5_neko_real, chkp%payloads(i)%ptr%mesh_arrays(j)%ptr, &
               layout)
       end do
 
       do j = 1, chkp%payloads(i)%ptr%array_count()
          call hdf5_checkpoint_read_array(payload_id, plist_id, &
               h5_neko_real, chkp%payloads(i)%ptr%arrays(j)%ptr)
       end do
 
       call h5gclose_f(payload_id, ierr)
    end do
 
    call h5gclose_f(payloads_id, ierr)
 
  end subroutine hdf5_checkpoint_read_payloads
 
  subroutine hdf5_checkpoint_open_group(root_id, path, create, group_id)
    integer(hid_t), intent(in) :: root_id
    character(len=*), intent(in) :: path
    logical, intent(in) :: create
    integer(hid_t), intent(out) :: group_id
    integer(hid_t) :: current_id, next_id
    integer :: first, last, slash, path_len, ierr
    logical :: group_exists
    character(len=:), allocatable :: group_name
 
    current_id = root_id
    first = 1
    path_len = len_trim(path)
 
    do
       slash = index(path(first:path_len), "/")
       if (slash .eq. 0) then
          last = path_len
       else
          last = first + slash - 2
       end if
 
       group_name = path(first:last)
       call h5lexists_f(current_id, group_name, group_exists, ierr)
       if (group_exists) then
          call h5gopen_f(current_id, group_name, next_id, ierr)
       else if (create) then
          call h5gcreate_f(current_id, group_name, next_id, ierr)
       else
          call neko_error("Checkpoint payload group '" // trim(path) // &
               "' is missing")
       end if
 
       if (current_id .ne. root_id) call h5gclose_f(current_id, ierr)
       current_id = next_id
 
       if (slash .eq. 0) exit
       first = last + 2
    end do
 
    group_id = current_id
 
  end subroutine hdf5_checkpoint_open_group
 
  subroutine hdf5_checkpoint_write_field(group_id, plist_id, &
       h5_neko_real, fld)
    integer(hid_t), intent(in) :: group_id, plist_id, h5_neko_real
    type(field_t), intent(in) :: fld
    integer(hid_t) :: dset_id, filespace, memspace
    integer(hsize_t), dimension(1) :: ddim, dcount, doffset
    integer :: ierr
 
    dcount(1) = int(fld%dof%size(), 8)
    doffset(1) = int(fld%msh%offset_el, 8) * int(fld%Xh%lxyz, 8)
    ddim = dcount
    call mpi_allreduce(mpi_in_place, ddim(1), 1, &
         mpi_integer8, mpi_sum, neko_comm, ierr)
 
    call h5screate_simple_f(1, ddim, filespace, ierr)
    call h5screate_simple_f(1, dcount, memspace, ierr)
    call h5sselect_hyperslab_f(filespace, h5s_select_set_f, &
         doffset, dcount, ierr)
    call h5dcreate_f(group_id, trim(fld%name), h5_neko_real, &
         filespace, dset_id, ierr)
    call h5dwrite_f(dset_id, h5_neko_real, fld%x(1,1,1,1), &
         dcount, ierr, file_space_id = filespace, mem_space_id = memspace, &
         xfer_prp = plist_id)
 
    call h5dclose_f(dset_id, ierr)
    call h5sclose_f(filespace, ierr)
    call h5sclose_f(memspace, ierr)
 
  end subroutine hdf5_checkpoint_write_field
 
  subroutine hdf5_checkpoint_read_field(group_id, plist_id, &
       h5_neko_real, fld, layout)
    integer(hid_t), intent(in) :: group_id, plist_id, h5_neko_real
    type(field_t), intent(inout) :: fld
    type(hdf5_checkpoint_layout_t), target, intent(inout) :: layout
    integer(hid_t) :: dset_id, filespace, memspace
    integer(hsize_t), dimension(1) :: dcount, doffset
    real(kind=rp), allocatable :: checkpoint_data(:)
    type(interpolator_t) :: space_interp
    integer :: ierr
    logical :: dataset_exists
 
    call h5lexists_f(group_id, trim(fld%name), dataset_exists, ierr)
    if (.not. dataset_exists) then
       call neko_error("Checkpoint field '" // trim(fld%name) // &
            "' is missing")
    end if
 
    call h5dopen_f(group_id, trim(fld%name), dset_id, ierr)
    call h5dget_space_f(dset_id, filespace, ierr)
 
    ! The hyperslab follows the file's own layout, which coincides with
    ! the running case's only when neither the mesh nor the order changed.
    dcount(1) = int(layout%msh%nelv, hsize_t) * int(layout%Xh%lxyz, hsize_t)
    doffset(1) = int(layout%msh%offset_el, hsize_t) * &
         int(layout%Xh%lxyz, hsize_t)
    call h5screate_simple_f(1, dcount, memspace, ierr)
    call h5sselect_hyperslab_f(filespace, h5s_select_set_f, &
         doffset, dcount, ierr)
 
    if (layout%mesh2mesh) then
       allocate(checkpoint_data(int(dcount(1))))
       call h5dread_f(dset_id, h5_neko_real, checkpoint_data, dcount, ierr, &
            file_space_id = filespace, mem_space_id = memspace, &
            xfer_prp = plist_id)
       call rzero(fld%x, fld%dof%size())
       call layout%global_interp%evaluate(fld%x, checkpoint_data, .true.)
       deallocate(checkpoint_data)
    else if (layout%Xh%lxyz .ne. fld%Xh%lxyz) then
       allocate(checkpoint_data(int(dcount(1))))
       call h5dread_f(dset_id, h5_neko_real, checkpoint_data, dcount, ierr, &
            file_space_id = filespace, mem_space_id = memspace, &
            xfer_prp = plist_id)
       call space_interp%init(fld%Xh, layout%Xh)
       call space_interp%map_host(fld%x, checkpoint_data, fld%msh%nelv, &
            fld%Xh)
       call space_interp%free()
       deallocate(checkpoint_data)
    else
       call h5dread_f(dset_id, h5_neko_real, fld%x(1,1,1,1), &
            dcount, ierr, file_space_id = filespace, &
            mem_space_id = memspace, xfer_prp = plist_id)
    end if
 
    call h5dclose_f(dset_id, ierr)
    call h5sclose_f(filespace, ierr)
    call h5sclose_f(memspace, ierr)
 
  end subroutine hdf5_checkpoint_read_field
 
  subroutine hdf5_checkpoint_write_mesh_array(group_id, plist_id, &
       h5_neko_real, array)
    integer(hid_t), intent(in) :: group_id, plist_id, h5_neko_real
    type(checkpoint_mesh_array_t), intent(in) :: array
    integer(hid_t) :: dset_id, filespace, memspace, attr_id, attr_space
    integer(hsize_t), dimension(1) :: ddim, dcount, doffset, attr_dims
    integer :: nodal_shape(3)
    integer :: ierr
 
    dcount(1) = int(size(array%x), hsize_t)
    doffset(1) = int(array%msh%offset_el, hsize_t) * &
         int(array%Xh%lxyz, hsize_t)
    ddim(1) = int(array%msh%glb_nelv, hsize_t) * &
         int(array%Xh%lxyz, hsize_t)
 
    call h5screate_simple_f(1, ddim, filespace, ierr)
    call h5screate_simple_f(1, dcount, memspace, ierr)
    call h5sselect_hyperslab_f(filespace, h5s_select_set_f, &
         doffset, dcount, ierr)
    call h5dcreate_f(group_id, trim(array%name), h5_neko_real, &
         filespace, dset_id, ierr)
 
    nodal_shape = [array%Xh%lx, array%Xh%ly, array%Xh%lz]
    attr_dims(1) = size(nodal_shape)
    call h5screate_simple_f(1, attr_dims, attr_space, ierr)
    call h5acreate_f(dset_id, "NodalShape", h5t_native_integer, &
         attr_space, attr_id, ierr)
    call h5awrite_f(attr_id, h5t_native_integer, nodal_shape, &
         attr_dims, ierr)
    call h5aclose_f(attr_id, ierr)
    call h5sclose_f(attr_space, ierr)
 
    call h5dwrite_f(dset_id, h5_neko_real, array%x, dcount, ierr, &
         file_space_id = filespace, mem_space_id = memspace, &
         xfer_prp = plist_id)
 
    call h5dclose_f(dset_id, ierr)
    call h5sclose_f(filespace, ierr)
    call h5sclose_f(memspace, ierr)
 
  end subroutine hdf5_checkpoint_write_mesh_array
 
  subroutine hdf5_checkpoint_read_mesh_array(group_id, plist_id, &
       h5_neko_real, array, layout)
    integer(hid_t), intent(in) :: group_id, plist_id, h5_neko_real
    type(checkpoint_mesh_array_t), intent(inout) :: array
    type(hdf5_checkpoint_layout_t), target, intent(inout) :: layout
    integer(hid_t) :: dset_id, filespace, memspace, attr_id, attr_space
    integer(hsize_t), dimension(1) :: dcount, doffset
    integer(hsize_t), dimension(1) :: dataset_dims, dataset_maxdims
    integer(hsize_t), dimension(1) :: attr_dims, attr_maxdims
    integer :: stored_shape(3)
    real(kind=rp), allocatable :: stored_data(:)
    type(neko_space_t), target :: stored_xh
    type(interpolator_t) :: space_interp
    integer :: ierr
    logical :: dataset_exists, shape_exists
 
    ! Mesh arrays hold ALE state (mesh coordinates, lagged mass matrices),
    ! which has no meaning on a different mesh. Same check and message as
    ! chkp_file.
    if (layout%mesh2mesh) then
       call neko_error('ALE does not yet support mesh2mesh ' // &
            'interpolation for restart!')
    end if
 
    call h5lexists_f(group_id, trim(array%name), dataset_exists, ierr)
    if (.not. dataset_exists) then
       call neko_error("Checkpoint mesh array '" // trim(array%name) // &
            "' is missing")
    end if
 
    call h5dopen_f(group_id, trim(array%name), dset_id, ierr)
    call h5dget_space_f(dset_id, filespace, ierr)
    call h5sget_simple_extent_dims_f(filespace, dataset_dims, &
         dataset_maxdims, ierr)
 
    stored_shape = [layout%Xh%lx, layout%Xh%ly, layout%Xh%lz]
    call h5aexists_f(dset_id, "NodalShape", shape_exists, ierr)
    if (shape_exists) then
       call h5aopen_f(dset_id, "NodalShape", attr_id, ierr)
    else
       ! Compatibility with development checkpoints using generic blocks.
       call h5aexists_f(dset_id, "BlockShape", shape_exists, ierr)
       if (shape_exists) then
          call h5aopen_f(dset_id, "BlockShape", attr_id, ierr)
       end if
    end if
    if (shape_exists) then
       call h5aget_space_f(attr_id, attr_space, ierr)
       call h5sget_simple_extent_dims_f(attr_space, attr_dims, &
            attr_maxdims, ierr)
       if (attr_dims(1) .ne. 3) then
          call neko_error("Checkpoint mesh-array nodal shape must have " // &
               "three dimensions")
       end if
       call h5aread_f(attr_id, h5t_native_integer, stored_shape, &
            attr_dims, ierr)
       call h5sclose_f(attr_space, ierr)
       call h5aclose_f(attr_id, ierr)
    end if
 
    if (any(stored_shape .le. 0)) then
       call neko_error("Checkpoint mesh-array nodal shape must be positive")
    end if
    if (array%msh%gdim .eq. 3) then
       call stored_xh%init(gll, stored_shape(1), stored_shape(2), &
            stored_shape(3))
    else
       if (stored_shape(3) .ne. 1) then
          call neko_error("Two-dimensional checkpoint mesh arrays must " // &
               "have one nodal plane")
       end if
       call stored_xh%init(gll, stored_shape(1), stored_shape(2))
    end if
 
    if (int(dataset_dims(1), i8) .ne. &
         int(layout%msh%glb_nelv, i8) * int(stored_xh%lxyz, i8)) then
       call neko_error("Checkpoint mesh array '" // trim(array%name) // &
            "' does not match the current mesh")
    end if
 
    dcount(1) = int(layout%msh%nelv, hsize_t) * &
         int(stored_xh%lxyz, hsize_t)
    doffset(1) = int(layout%msh%offset_el, hsize_t) * &
         int(stored_xh%lxyz, hsize_t)
    call h5screate_simple_f(1, dcount, memspace, ierr)
    call h5sselect_hyperslab_f(filespace, h5s_select_set_f, &
         doffset, dcount, ierr)
 
    if (stored_xh%lxyz .eq. array%Xh%lxyz) then
       call h5dread_f(dset_id, h5_neko_real, array%x, dcount, ierr, &
            file_space_id = filespace, mem_space_id = memspace, &
            xfer_prp = plist_id)
    else
       allocate(stored_data(int(dcount(1))))
       call h5dread_f(dset_id, h5_neko_real, stored_data, dcount, ierr, &
            file_space_id = filespace, mem_space_id = memspace, &
            xfer_prp = plist_id)
       call space_interp%init(array%Xh, stored_xh)
       call space_interp%map_host(array%x, stored_data, array%msh%nelv, &
            array%Xh)
       call space_interp%free()
       deallocate(stored_data)
    end if
 
    call stored_xh%free()
    call h5dclose_f(dset_id, ierr)
    call h5sclose_f(filespace, ierr)
    call h5sclose_f(memspace, ierr)
 
  end subroutine hdf5_checkpoint_read_mesh_array
 
  subroutine hdf5_checkpoint_write_array(group_id, plist_id, &
       h5_neko_real, array)
    integer(hid_t), intent(in) :: group_id, plist_id, h5_neko_real
    type(checkpoint_array_t), intent(in) :: array
    integer(hid_t) :: dset_id, filespace, memspace, h5_type
    integer(hsize_t), dimension(1) :: ddim, dcount, doffset
    integer :: ierr
    logical :: stored_dp
 
    ! Arrays registered through `add_array_dp` keep double precision in the
    ! file whatever `rp` is, so that the time history survives a restart into
    ! a build of the other precision.
    stored_dp = associated(array%x_dp)
    if (stored_dp) then
       h5_type = h5t_native_double
       dcount(1) = int(size(array%x_dp), hsize_t)
    else
       h5_type = h5_neko_real
       dcount(1) = int(size(array%x), hsize_t)
    end if
    ddim(1) = int(array%global_count, hsize_t)
    doffset(1) = int(array%offset, hsize_t)
 
    call h5screate_simple_f(1, ddim, filespace, ierr)
    call h5screate_simple_f(1, dcount, memspace, ierr)
    call h5dcreate_f(group_id, trim(array%name), h5_type, &
         filespace, dset_id, ierr)
 
    if (array%replicated .and. pe_rank .ne. 0) then
       call h5sselect_none_f(filespace, ierr)
       call h5sselect_none_f(memspace, ierr)
    else
       call h5sselect_hyperslab_f(filespace, h5s_select_set_f, &
            doffset, dcount, ierr)
    end if
 
    if (stored_dp) then
       call h5dwrite_f(dset_id, h5_type, array%x_dp, dcount, ierr, &
            file_space_id = filespace, mem_space_id = memspace, &
            xfer_prp = plist_id)
    else
       call h5dwrite_f(dset_id, h5_type, array%x, dcount, ierr, &
            file_space_id = filespace, mem_space_id = memspace, &
            xfer_prp = plist_id)
    end if
 
    call h5dclose_f(dset_id, ierr)
    call h5sclose_f(filespace, ierr)
    call h5sclose_f(memspace, ierr)
 
  end subroutine hdf5_checkpoint_write_array
 
  subroutine hdf5_checkpoint_read_array(group_id, plist_id, &
       h5_neko_real, array)
    integer(hid_t), intent(in) :: group_id, plist_id, h5_neko_real
    type(checkpoint_array_t), intent(inout) :: array
    integer(hid_t) :: dset_id, filespace, memspace, h5_type
    integer(hsize_t), dimension(1) :: dcount, doffset
    integer(hsize_t), dimension(1) :: dataset_dims, dataset_maxdims
    integer :: ierr
    logical :: dataset_exists, stored_dp
 
    call h5lexists_f(group_id, trim(array%name), dataset_exists, ierr)
    if (.not. dataset_exists) then
       call neko_error("Checkpoint array '" // trim(array%name) // &
            "' is missing")
    end if
 
    call h5dopen_f(group_id, trim(array%name), dset_id, ierr)
    call h5dget_space_f(dset_id, filespace, ierr)
    call h5sget_simple_extent_dims_f(filespace, dataset_dims, &
         dataset_maxdims, ierr)
 
    if (int(dataset_dims(1), i8) .ne. array%global_count) then
       call neko_error("Checkpoint array '" // trim(array%name) // &
            "' has an incompatible global extent")
    end if
 
    stored_dp = associated(array%x_dp)
    if (stored_dp) then
       h5_type = h5t_native_double
       dcount(1) = int(size(array%x_dp), hsize_t)
    else
       h5_type = h5_neko_real
       dcount(1) = int(size(array%x), hsize_t)
    end if
    doffset(1) = int(array%offset, hsize_t)
    call h5screate_simple_f(1, dcount, memspace, ierr)
 
    if (array%replicated .and. pe_rank .ne. 0) then
       call h5sselect_none_f(filespace, ierr)
       call h5sselect_none_f(memspace, ierr)
    else
       call h5sselect_hyperslab_f(filespace, h5s_select_set_f, &
            doffset, dcount, ierr)
    end if
    if (stored_dp) then
       call h5dread_f(dset_id, h5_type, array%x_dp, dcount, ierr, &
            file_space_id = filespace, mem_space_id = memspace, &
            xfer_prp = plist_id)
    else
       call h5dread_f(dset_id, h5_type, array%x, dcount, ierr, &
            file_space_id = filespace, mem_space_id = memspace, &
            xfer_prp = plist_id)
    end if
 
    call h5dclose_f(dset_id, ierr)
    call h5sclose_f(filespace, ierr)
    call h5sclose_f(memspace, ierr)
 
    ! Only rank zero selected the dataset above, so the copy every rank holds
    ! has to come from there.
    if (array%replicated) then
       if (stored_dp) then
          call mpi_bcast(array%x_dp, size(array%x_dp), &
               mpi_double_precision, 0, neko_comm, ierr)
       else
          call mpi_bcast(array%x, size(array%x), mpi_real_precision, 0, &
               neko_comm, ierr)
       end if
    end if
 
  end subroutine hdf5_checkpoint_read_array
 
  subroutine hdf5_file_determine_data(data, msh, dof, fp, fsp, dtlag, tlag)
    class(*), target, intent(in) :: data
    type(mesh_t), pointer, intent(inout) :: msh
    type(dofmap_t), pointer, intent(inout) :: dof
    type(field_ptr_t), allocatable, intent(inout) :: fp(:)
    type(field_series_ptr_t), allocatable, intent(inout) :: fsp(:)
    real(kind=dp), pointer, intent(inout) :: dtlag(:)
    real(kind=dp), pointer, intent(inout) :: tlag(:)
    integer :: i, j, fp_size, fp_cur, fsp_size, fsp_cur
 
    select type (data)
    type is (field_t)
       dof => data%dof
       msh => data%msh
       fp_size = 1
       allocate(fp(fp_size))
       fp(1)%ptr => data
 
       nullify(dtlag)
       nullify(tlag)
 
    type is (field_list_t)
 
       if (data%size() .gt. 0) then
          allocate(fp(data%size()))
 
          dof => data%dof(1)
          msh => data%msh(1)
 
          do i = 1, data%size()
             fp(i)%ptr => data%items(i)%ptr
          end do
       else
          call neko_error('Empty field list')
       end if
 
       nullify(dtlag)
       nullify(tlag)
 
    type is (chkp_t)
       block
         type(checkpoint_payload_t), pointer :: fluid
         type(field_t), pointer :: u
 
         fluid => data%get_payload("fluid")
         u => fluid%find_field("u")
         if (.not. associated(u) .or. &
              .not. associated(fluid%find_field("v")) .or. &
              .not. associated(fluid%find_field("w")) .or. &
              .not. associated(fluid%find_field("p"))) then
            call neko_error("Checkpoint not initialized")
         end if
         dof => u%dof
         msh => u%msh
 
         fp_size = 0
         fsp_size = 0
         do i = 1, data%payload_count()
            fp_size = fp_size + data%payloads(i)%ptr%field_count()
            fsp_size = fsp_size + data%payloads(i)%ptr%series_count()
         end do
 
         if (fp_size .gt. 0) allocate(fp(fp_size))
         if (fsp_size .gt. 0) allocate(fsp(fsp_size))
 
         fp_cur = 1
         fsp_cur = 1
         do i = 1, data%payload_count()
            do j = 1, data%payloads(i)%ptr%field_count()
               fp(fp_cur)%ptr => data%payloads(i)%ptr%fields(j)%ptr
               fp_cur = fp_cur + 1
            end do
            do j = 1, data%payloads(i)%ptr%series_count()
               fsp(fsp_cur)%ptr => data%payloads(i)%ptr%series(j)%ptr
               fsp_cur = fsp_cur + 1
            end do
         end do
 
         call data%get_time_history(tlag, dtlag)
       end block
 
    class default
       call neko_log%error('Invalid data')
    end select
 
  end subroutine hdf5_file_determine_data
 
  subroutine hdf5_file_determine_real(H5T_NEKO_REAL)
    integer(hid_t), intent(inout) :: H5T_NEKO_REAL
    select case (rp)
    case (dp)
       h5t_neko_real = h5t_native_double
    case (sp)
       h5t_neko_real = h5t_native_real
    case default
       call neko_error("Unsupported real type")
    end select
  end subroutine hdf5_file_determine_real
 
  ! ================
  ! Granular methods
  ! ================
 
  subroutine hdf5_file_open(this, mode)
    class(hdf5_file_t), intent(inout) :: this
    character(len=1), intent(in) :: mode
    integer :: ierr, mpi_info, mpi_comm, i, n_fields, counter
    logical :: file_exists
    character(len=1024) :: fname
    character(len=LOG_SIZE) :: log_buf
 
    ! Set the mode for the file
    this%mode = mode
 
    ! Ensure precision is set and are valid.
    if (this%precision .gt. rp) then
       this%precision = rp
       call neko_warning('Requested precision is higher than working precision')
    else if (this%precision .eq. -1) then
       this%precision = rp
    end if
 
    fname = trim(file_get_fname(this))
    counter = this%get_counter() - this%get_start_counter()
 
    ! Set the configuration for MPI IO
    call hdf5_session_init()
 
    mpi_info = mpi_info_null%mpi_val
    mpi_comm = neko_comm%mpi_val
    call h5pcreate_f(h5p_file_access_f, this%plist_id, ierr)
    call h5pset_fapl_mpio_f(this%plist_id, mpi_comm, mpi_info, ierr)
 
    ! Open the file
    inquire(file = fname, exist = file_exists)
    if (file_exists) then
       call h5fopen_f(fname, h5f_acc_rdwr_f, this%file_id, ierr, &
            access_prp = this%plist_id)
    else
       call h5fcreate_f(fname, h5f_acc_trunc_f, &
            this%file_id, ierr, access_prp = this%plist_id)
    end if
 
    ! Set the active group to the root of the file
    call this%set_active_group()
 
    write (log_buf, *) "Opened HDF5 file: ", trim(fname), " with counter: ", &
         counter
    call neko_log%message(log_buf, lvl = neko_log_debug)
 
  end subroutine hdf5_file_open
 
  subroutine hdf5_file_close(this)
    class(hdf5_file_t), intent(inout) :: this
    integer :: ierr
 
    if (this%active_group_id .ne. -1_hid_t .and. &
         this%active_group_id .ne. this%file_id) then
       call h5gclose_f(this%active_group_id, ierr)
    end if
    this%active_group_id = -1_hid_t
 
    call h5pclose_f(this%plist_id, ierr)
    this%plist_id = -1_hid_t
    call h5fclose_f(this%file_id, ierr)
    this%file_id = -1_hid_t
    call hdf5_session_finalize()
 
    call neko_log%message("Closed HDF5 file: " // trim(this%get_fname()), &
         lvl = neko_log_debug)
 
  end subroutine hdf5_file_close
 
 
  subroutine hdf5_file_set_group(this, group_name_path)
    class(hdf5_file_t), intent(inout) :: this
    character(len=*), intent(in), optional :: group_name_path
    character(len=1000), allocatable :: group_name(:)
 
    integer(hid_t) :: current_id, group_id
    integer :: ierr, i, j, num_groups, name_len, group_loc
    logical :: group_exists
 
 
    ! Close previous active group if one is open
    if (this%active_group_id .ne. -1_hid_t .and. this%active_group_id .ne. &
         this%file_id) then
       call h5gclose_f(this%active_group_id, ierr)
    end if
    this%active_group_id = -1_hid_t
 
    ! Start from root location = file
    current_id = this%file_id
    ! Return the root directory if no group name is given
    if (.not. present(group_name_path)) then
       this%active_group_id = current_id
       return
    end if
 
    ! Split the input string into group names using "/" as a delimiter
    name_len = len(trim(group_name_path))
    ! Count how many groups
    num_groups = 1 ! There is at least one group if this was passed
    do i = 1, name_len
       if (group_name_path .eq. "/") then
          num_groups = num_groups + 1
       end if
    end do
 
    ! Allocate the group array and populate it
    allocate(group_name(num_groups))
    j = 1
    group_loc = 1
    do i = 1, name_len
       if (group_name_path .eq. "/") then
          group_name(group_loc) = group_name_path(j:i-1)
          group_loc = group_loc + 1
          j = i + 1
       end if
    end do
    if (j .ne. name_len) then
       group_name(group_loc) = group_name_path(j:name_len)
    end if
 
    ! Iterate over the groups in the path
    do i = 1, num_groups
       call h5lexists_f(current_id, trim(group_name(i)), group_exists, ierr)
 
       ! Only create groups if they dont exist and we are in write mode "w"
       if (group_exists) then
          call h5gopen_f(current_id, trim(group_name(i)), group_id, ierr)
       else
          if (this%mode == "r") then
             call neko_error("Group " // trim(group_name(i)) // &
                  " does not exist in file " // trim(file_get_fname(this)))
          end if
          call h5gcreate_f(current_id, trim(group_name(i)), group_id, ierr)
       end if
 
       ! Close previous location only if it was an opened group, not the file
       if (i > 1) then
          call h5gclose_f(current_id, ierr)
       end if
 
       current_id = group_id
    end do
 
    this%active_group_id = current_id
  end subroutine hdf5_file_set_group
 
 
  subroutine hdf5_file_write_dataset(this, data)
    class(hdf5_file_t), intent(inout) :: this
    class(*), intent(inout) :: data
 
    select type (d => data)
    type is (vector_t)
       call this%write_vector(d)
    type is (matrix_t)
       call this%write_matrix(d)
    type is (field_t)
       call this%write_field(d)
    class default
       call neko_error("write_dataset not implemented for this data type")
    end select
  end subroutine hdf5_file_write_dataset
 
  subroutine hdf5_file_read_dataset(this, data_name, data, strategy)
    class(hdf5_file_t), intent(inout) :: this
    character(len=*), intent(in) :: data_name
    class(*), intent(inout) :: data
    character(len=*), intent(in), optional :: strategy
 
    select type (d => data)
    type is (vector_t)
       call this%read_vector(data_name, d, strategy)
    type is (matrix_t)
       call this%read_matrix(data_name, d, strategy)
    type is (field_t)
       call neko_error("Reading a field_t is not supported yet")
    class default
       call neko_error("read_dataset not implemented for this data type")
    end select
  end subroutine hdf5_file_read_dataset
 
  subroutine hdf5_file_write_attribute(this, data_name, data)
    class(hdf5_file_t), intent(inout) :: this
    character(len=*), intent(in) :: data_name
    class(*), intent(inout) :: data
 
    select type (d => data)
    type is (integer)
       call this%write_int_attribute(data_name, d)
    type is (real(kind=rp))
       call this%write_rp_attribute(data_name, d)
    class default
       call neko_error("write_attribute not implemented for this data type")
    end select
  end subroutine hdf5_file_write_attribute
 
  subroutine hdf5_file_read_attribute(this, data_name, data, exist)
    class(hdf5_file_t), intent(inout) :: this
    character(len=*), intent(in) :: data_name
    class(*), intent(inout) :: data
    logical, intent(inout) :: exist
 
    select type (d => data)
    type is (integer)
       call this%read_int_attribute(data_name, d, exist)
    type is (real(kind=rp))
       call this%read_rp_attribute(data_name, d, exist)
    class default
       call neko_error("read_attribute not implemented for this data type")
    end select
  end subroutine hdf5_file_read_attribute
 
 
  subroutine hdf5_file_write_vector(this, vec)
    class(hdf5_file_t), intent(inout) :: this
    type(vector_t), intent(inout) :: vec
    integer :: ierr, counts, offset, total_count, dset_rank, max_count
    integer(hsize_t) :: append_offset
    integer(hid_t) :: precision_hdf
    integer(hid_t) :: xf_id, filespace, dset_id, memspace, dcpl_id
    integer(hsize_t), dimension(1) :: dcount, doffset
    integer(hsize_t), dimension(1) :: ddims, ddims_max, chunkdims
    integer(hsize_t), dimension(1) :: tempddims, tempmaxddims
    logical :: dset_exists
    real(kind=sp), allocatable :: write_buffer_sp(:) ! Write buffer single
    real(kind=dp), allocatable :: write_buffer_dp(:) ! Write buffer double
 
    ! ===============
    ! Get vector info
    ! ===============
    counts = vec%size()
    append_offset = 0_hsize_t
    offset = 0
    total_count = 0
    max_count = 0
    call mpi_scan(counts, offset, 1, mpi_integer, &
         mpi_sum, neko_comm, ierr)
    offset = offset - counts ! Not using exclusive scan
    call mpi_allreduce(counts, total_count, 1, mpi_integer, &
         mpi_sum, neko_comm, ierr)
    call mpi_allreduce(counts, max_count, 1, mpi_integer, &
         mpi_max, neko_comm, ierr)
 
    ! ===============
    ! Configure MPIIO
    ! ===============
    call h5pcreate_f(h5p_dataset_xfer_f, xf_id, ierr)
    call h5pset_dxpl_mpio_f(xf_id, h5fd_mpio_collective_f, ierr)
    precision_hdf = h5kind_to_type(this%precision, h5_real_kind)
 
    ! ===================
    ! Create the data set
    ! ===================
    dset_rank = 1 ! rank 1 array, i.e. a vector
    ddims = [int(total_count, hsize_t)] ! global size of the vector
 
    ! Enable chunking to be able to append
    chunkdims = [max(int(max_count, hsize_t), 1_hsize_t)]
    ddims_max = [h5s_unlimited_f] ! allow unlimited size for appending
    call h5lexists_f(this%active_group_id, trim(vec%name), dset_exists, ierr)
    if (dset_exists) then
       if (this%overwrite) then
          ! retrieve the dset id for the existing data set
          call h5dopen_f(this%active_group_id, trim(vec%name), dset_id, ierr)
       else
          ! Retreive the existing data set
          call h5dopen_f(this%active_group_id, trim(vec%name), dset_id, ierr)
          ! Retrieve the current filespace (shape space)
          call h5dget_space_f(dset_id, filespace, ierr)
          ! Get the current shape
          call h5sget_simple_extent_dims_f(filespace, tempddims, tempmaxddims, &
               ierr)
          ! Clean up the opened file space
          call h5sclose_f(filespace, ierr)
          ! Overwrite the new full shape
          ddims(1) = ddims(1) + tempddims(1) ! New size
          append_offset = tempddims(1) ! current size which is the offset
          ! Extend the data set to the new shape
          call h5dset_extent_f(dset_id, ddims, ierr)
       end if
    else
       ! create file space of this shape
       call h5screate_simple_f(dset_rank, ddims, filespace, ierr, ddims_max)
       ! Create chunk property list (needed to be able to append)
       call h5pcreate_f(h5p_dataset_create_f, dcpl_id, ierr)
       call h5pset_chunk_f(dcpl_id, dset_rank, chunkdims, ierr)
       ! create the data set with the given shape
       call h5dcreate_f(this%active_group_id, trim(vec%name), precision_hdf, &
            filespace, dset_id, ierr, dcpl_id = dcpl_id)
       ! clean opened ids
       call h5sclose_f(filespace, ierr)
       call h5pclose_f(dcpl_id, ierr)
    end if
 
    ! ===========================
    ! Set up writing the data set
    ! ===========================
    dcount = [int(counts, hsize_t)] ! local size of the vector
 
    ! offset for this rank in the global vector
    doffset = [int(offset, hsize_t) + append_offset]
    ! Get the total file space (shape) of the data set
    call h5dget_space_f(dset_id, filespace, ierr)
    ! Get only the slice where my rank writes
    call h5sselect_hyperslab_f(filespace, h5s_select_set_f, doffset, dcount, &
         ierr)
    ! Create the corresponding memory space (buffer) for my local data
    call h5screate_simple_f(dset_rank, dcount, memspace, ierr)
 
 
    ! =======================
    ! Cast and write the data
    ! =======================
    if (this%precision == sp) then
       allocate(write_buffer_sp(vec%size()))
       if (vec%size() > 0) write_buffer_sp = real(vec%x, kind=sp)
       ! Write the data
       call h5dwrite_f(dset_id, precision_hdf, write_buffer_sp, dcount, ierr, &
            file_space_id = filespace, mem_space_id = memspace, &
            xfer_prp = xf_id)
       deallocate(write_buffer_sp)
    else if (this%precision == dp) then
       allocate(write_buffer_dp(vec%size()))
       if (vec%size() > 0) write_buffer_dp = real(vec%x, kind=dp)
       ! Write the data
       call h5dwrite_f(dset_id, precision_hdf, write_buffer_dp, dcount, ierr, &
            file_space_id = filespace, mem_space_id = memspace, &
            xfer_prp = xf_id)
       deallocate(write_buffer_dp)
    else
       call neko_error("Unsupported precision")
    end if
 
    ! =======================
    ! Clean up
    ! =======================
    call h5pclose_f(xf_id, ierr)
    call h5sclose_f(memspace, ierr)
    call h5sclose_f(filespace, ierr)
    call h5dclose_f(dset_id, ierr)
 
  end subroutine hdf5_file_write_vector
 
  subroutine hdf5_file_write_matrix(this, mat)
    class(hdf5_file_t), intent(inout) :: this
    type(matrix_t), intent(inout) :: mat
    integer :: ierr, counts, offset, total_count, dset_rank, strides, max_count
    integer(hsize_t) :: append_offset
    integer(hid_t) :: precision_hdf
    integer(hid_t) :: xf_id, filespace, dset_id, memspace, dcpl_id
    integer(hsize_t), dimension(2) :: dcount, doffset
    integer(hsize_t), dimension(2) :: ddims, ddims_max, chunkdims
    integer(hsize_t), dimension(2) :: tempddims, tempmaxddims
    logical :: dset_exists
    real(kind=sp), allocatable :: write_buffer_sp(:,:) ! Write buffer single
    real(kind=dp), allocatable :: write_buffer_dp(:,:) ! Write buffer double
 
    ! ===============
    ! Get Matrix info
    ! ===============
    strides = mat%get_nrows()
    counts = mat%get_ncols()
    append_offset = 0_hsize_t
    total_count = 0
    max_count = 0
    offset = 0
    call mpi_scan(counts, offset, 1, mpi_integer, &
         mpi_sum, neko_comm, ierr)
    offset = offset - counts ! Not using exclusive scan
    call mpi_allreduce(counts, total_count, 1, mpi_integer, &
         mpi_sum, neko_comm, ierr)
    call mpi_allreduce(counts, max_count, 1, mpi_integer, &
         mpi_max, neko_comm, ierr)
 
    ! ===============
    ! Configure MPIIO
    ! ===============
    call h5pcreate_f(h5p_dataset_xfer_f, xf_id, ierr)
    call h5pset_dxpl_mpio_f(xf_id, h5fd_mpio_collective_f, ierr)
    precision_hdf = h5kind_to_type(this%precision, h5_real_kind)
 
    ! ===================
    ! Create the data set
    ! ===================
    dset_rank = 2 ! rank 2 array, i.e. a matrix
    ! global size of the matrix
    ddims = [int(strides, hsize_t), int(total_count, hsize_t)]
    chunkdims = [int(strides, hsize_t), max(int(max_count, hsize_t), 1_hsize_t)]
    ddims_max = [int(strides, hsize_t), h5s_unlimited_f]
    call h5lexists_f(this%active_group_id, trim(mat%name), dset_exists, ierr)
    if (dset_exists) then
       if (this%overwrite) then
 
          if (pe_rank .eq. 0) then
             call neko_warning("Dataset " // trim(mat%name) // &
                  " already exists and wil be overwritten")
          end if
          ! retrieve the dset id for the existing data set
          call h5dopen_f(this%active_group_id, trim(mat%name), dset_id, ierr)
       else
          call h5dopen_f(this%active_group_id, trim(mat%name), dset_id, ierr)
          call h5dget_space_f(dset_id, filespace, ierr)
          call h5sget_simple_extent_dims_f(filespace, tempddims, tempmaxddims, &
               ierr)
          call h5sclose_f(filespace, ierr)
          ddims(2) = ddims(2) + tempddims(2)
          append_offset = tempddims(2)
          call h5dset_extent_f(dset_id, ddims, ierr)
       end if
    else
       ! create file space of this shape
       call h5screate_simple_f(dset_rank, ddims, filespace, ierr, ddims_max)
       call h5pcreate_f(h5p_dataset_create_f, dcpl_id, ierr)
       call h5pset_chunk_f(dcpl_id, dset_rank, chunkdims, ierr)
       ! create the data set with the given shape
       call h5dcreate_f(this%active_group_id, trim(mat%name), precision_hdf, &
            filespace, dset_id, ierr, dcpl_id = dcpl_id)
       call h5sclose_f(filespace, ierr)
       call h5pclose_f(dcpl_id, ierr)
    end if
 
    ! ===========================
    ! Set up writing the data set
    ! ===========================
    ! local size of the matrix
    dcount = [int(strides, hsize_t), int(counts, hsize_t)]
    ! offset for this rank in the global matrix
    doffset = [0_hsize_t, int(offset, hsize_t) + append_offset]
    ! Get the total file space (shape) of the data set
    call h5dget_space_f(dset_id, filespace, ierr)
    ! Get only the slice where my rank writes
    call h5sselect_hyperslab_f(filespace, h5s_select_set_f, doffset, dcount, &
         ierr)
    ! Create the corresponding memory space (buffer) for my local data
    call h5screate_simple_f(dset_rank, dcount, memspace, ierr)
 
    ! =======================
    ! Cast and write the data
    ! =======================
    if (this%precision == sp) then
       allocate(write_buffer_sp(mat%get_nrows(), mat%get_ncols()))
       if (mat%size() > 0) write_buffer_sp = real(mat%x, kind=sp)
       ! Write the data
       call h5dwrite_f(dset_id, precision_hdf, write_buffer_sp, dcount, ierr, &
            file_space_id = filespace, mem_space_id = memspace, &
            xfer_prp = xf_id)
       deallocate(write_buffer_sp)
    else if (this%precision == dp) then
       allocate(write_buffer_dp(mat%get_nrows(), mat%get_ncols()))
       if (mat%size() > 0) write_buffer_dp = real(mat%x, kind=dp)
       ! Write the data
       call h5dwrite_f(dset_id, precision_hdf, write_buffer_dp, dcount, ierr, &
            file_space_id = filespace, mem_space_id = memspace, &
            xfer_prp = xf_id)
       deallocate(write_buffer_dp)
    else
       call neko_error("Unsupported precision")
    end if
 
    ! =======================
    ! Clean up
    ! =======================
    call h5pclose_f(xf_id, ierr)
    call h5sclose_f(memspace, ierr)
    call h5sclose_f(filespace, ierr)
    call h5dclose_f(dset_id, ierr)
 
  end subroutine hdf5_file_write_matrix
 
  subroutine hdf5_file_write_field(this, field)
    class(hdf5_file_t), intent(inout) :: this
    type(field_t), intent(inout) :: field
    integer :: ierr, counts, offset, total_count, dset_rank, max_count
    integer :: stride_ax_1, stride_ax_2, stride_ax_3
    integer(hsize_t) :: append_offset
    integer(hid_t) :: precision_hdf
    integer(hid_t) :: xf_id, filespace, dset_id, memspace, dcpl_id
    integer(hsize_t), dimension(4) :: dcount, doffset
    integer(hsize_t), dimension(4) :: ddims, ddims_max, chunkdims
    integer(hsize_t), dimension(4) :: tempddims, tempmaxddims
    logical :: dset_exists
    real(kind=sp), allocatable :: write_buffer_sp(:,:,:,:) ! Write buffer single
    real(kind=dp), allocatable :: write_buffer_dp(:,:,:,:) ! Write buffer double
 
    ! ==============
    ! Get Field info
    ! ==============
    stride_ax_1 = field%Xh%lx
    stride_ax_2 = field%Xh%ly
    stride_ax_3 = field%Xh%lz
    counts = field%msh%nelv
    append_offset = 0_hsize_t
    total_count = field%msh%glb_nelv
    max_count = 0
    offset = field%msh%offset_el
    call mpi_allreduce(counts, max_count, 1, mpi_integer, &
         mpi_max, neko_comm, ierr)
 
    ! ===============
    ! Configure MPIIO
    ! ===============
    call h5pcreate_f(h5p_dataset_xfer_f, xf_id, ierr)
    call h5pset_dxpl_mpio_f(xf_id, h5fd_mpio_collective_f, ierr)
    precision_hdf = h5kind_to_type(this%precision, h5_real_kind)
 
    ! ===================
    ! Create the data set
    ! ===================
    dset_rank = 4 ! rank 4 array, i.e. a 4D tensor
    ddims = [int(stride_ax_1, hsize_t), &
         int(stride_ax_2, hsize_t), &
         int(stride_ax_3, hsize_t), &
         int(total_count, hsize_t)] ! global size of the tensor
    chunkdims = [int(stride_ax_1, hsize_t), &
         int(stride_ax_2, hsize_t), &
         int(stride_ax_3, hsize_t), &
         max(int(max_count, hsize_t), 1_hsize_t)]
    ddims_max = [int(stride_ax_1, hsize_t), &
         int(stride_ax_2, hsize_t), &
         int(stride_ax_3, hsize_t), &
         h5s_unlimited_f]
    call h5lexists_f(this%active_group_id, trim(field%name), dset_exists, ierr)
    if (dset_exists) then
       if (this%overwrite) then
          ! retrieve the dset id for the existing data set
          if (pe_rank .eq. 0) then
             call neko_warning("Overwriting dataset: " // trim(field%name))
          end if
          call h5dopen_f(this%active_group_id, trim(field%name), dset_id, ierr)
       else
          call h5dopen_f(this%active_group_id, trim(field%name), dset_id, ierr)
          call h5dget_space_f(dset_id, filespace, ierr)
          call h5sget_simple_extent_dims_f(filespace, tempddims, tempmaxddims, &
               ierr)
          call h5sclose_f(filespace, ierr)
          ddims(4) = ddims(4) + tempddims(4)
          append_offset = tempddims(4)
          call h5dset_extent_f(dset_id, ddims, ierr)
       end if
    else
       ! create file space of this shape
       call h5screate_simple_f(dset_rank, ddims, filespace, ierr, ddims_max)
       call h5pcreate_f(h5p_dataset_create_f, dcpl_id, ierr)
       call h5pset_chunk_f(dcpl_id, dset_rank, chunkdims, ierr)
       ! create the data set with the given shape
       call h5dcreate_f(this%active_group_id, trim(field%name), precision_hdf, &
            filespace, dset_id, ierr, dcpl_id = dcpl_id)
       call h5sclose_f(filespace, ierr)
       call h5pclose_f(dcpl_id, ierr)
    end if
 
    ! ===========================
    ! Set up writing the data set
    ! ===========================
    dcount = [int(stride_ax_1, hsize_t), &
         int(stride_ax_2, hsize_t), &
         int(stride_ax_3, hsize_t), &
         int(counts, hsize_t)] ! local size of the tensor
    doffset = [0_hsize_t, 0_hsize_t, 0_hsize_t, &
         int(offset, hsize_t) + append_offset] ! offset in the global tensor
    ! Get the total file space (shape) of the data set
    call h5dget_space_f(dset_id, filespace, ierr)
    ! Get only the slice where my rank writes
    call h5sselect_hyperslab_f(filespace, h5s_select_set_f, doffset, dcount, &
         ierr)
    ! Create the corresponding memory space (buffer) for my local data
    call h5screate_simple_f(dset_rank, dcount, memspace, ierr)
 
    ! =======================
    ! Cast and write the data
    ! =======================
    if (this%precision == sp) then
       allocate(write_buffer_sp(field%Xh%lx, field%Xh%ly, field%Xh%lz, &
            field%msh%nelv))
       if (field%msh%nelv > 0) write_buffer_sp = real(field%x, kind=sp)
       ! Write the data
       call h5dwrite_f(dset_id, precision_hdf, write_buffer_sp, dcount, ierr, &
            file_space_id = filespace, mem_space_id = memspace, &
            xfer_prp = xf_id)
       deallocate(write_buffer_sp)
    else if (this%precision == dp) then
       allocate(write_buffer_dp(field%Xh%lx, field%Xh%ly, field%Xh%lz, &
            field%msh%nelv))
       if (field%msh%nelv > 0) write_buffer_dp = real(field%x, kind=dp)
       ! Write the data
       call h5dwrite_f(dset_id, precision_hdf, write_buffer_dp, dcount, ierr, &
            file_space_id = filespace, mem_space_id = memspace, &
            xfer_prp = xf_id)
       deallocate(write_buffer_dp)
    else
       call neko_error("Unsupported precision")
    end if
 
    ! =======================
    ! Clean up
    ! =======================
    call h5pclose_f(xf_id, ierr)
    call h5sclose_f(memspace, ierr)
    call h5sclose_f(filespace, ierr)
    call h5dclose_f(dset_id, ierr)
 
  end subroutine hdf5_file_write_field
 
  subroutine hdf5_file_read_vector(this, data_name, vec, strategy)
    class(hdf5_file_t) :: this
    character(len=*), intent(in) :: data_name
    type(vector_t), intent(inout) :: vec
    character(len=*), intent(in), optional :: strategy
    character(len=1000) :: strategy_
    integer :: ierr, counts, offset, total_count, dset_rank
    integer(hid_t) :: precision_hdf
    integer(hid_t) :: xf_id, filespace, dset_id, memspace
    integer(hsize_t), dimension(1) :: dcount, doffset
    integer(hsize_t), dimension(1) :: tempddims, tempmaxddims
    integer :: temprank
    logical :: dset_exists
    type(linear_dist_t) :: dist
 
    ! Set up strategy
    if (present(strategy)) then
       if (trim(strategy) .eq. "linear" .or. &
            trim(strategy) .eq. "rank_0") then
          strategy_ = strategy
       else
          call neko_error("Unsupported strategy: " // trim(strategy))
       end if
    else
       strategy_ = "linear"
    end if
 
    ! Free the input
    call vec%free()
 
    ! ===============
    ! Configure MPIIO
    ! ===============
    call h5pcreate_f(h5p_dataset_xfer_f, xf_id, ierr)
    call h5pset_dxpl_mpio_f(xf_id, h5fd_mpio_collective_f, ierr)
    precision_hdf = h5kind_to_type(rp, h5_real_kind)
 
    ! ===================
    ! Get the data set info
    ! ===================
    call h5lexists_f(this%active_group_id, trim(data_name), dset_exists, ierr)
    if (dset_exists) then
       ! Open the data set
       call h5dopen_f(this%active_group_id, trim(data_name), dset_id, ierr)
       ! Get the current rank of the dataset
       call h5dget_space_f(dset_id, filespace, ierr)
       call h5sget_simple_extent_ndims_f(filespace, temprank, ierr)
       if (temprank .ne. 1) then
          call neko_error("Dataset " // trim(data_name) // &
               " is not a rank 1 vector in file " // trim(file_get_fname(this)))
       end if
       ! Get the current shape and close the filespace
       call h5sget_simple_extent_dims_f(filespace, tempddims, tempmaxddims, &
            ierr)
       call h5sclose_f(filespace, ierr)
    else
       call neko_error("Dataset " // trim(data_name) // &
            " does not exist in current group " // trim(file_get_fname(this)))
    end if
 
    ! =============================
    ! Perform the data distribution
    ! =============================
    total_count = int(tempddims(1))
    if (strategy_ .eq. "linear") then
       dist = linear_dist_t(total_count, pe_rank, pe_size, neko_comm)
       counts = dist%num_local()
       offset = 0
    else if (strategy_ .eq. "rank_0") then
       if (pe_rank .eq. 0) then
          counts = total_count
       else
          counts = 0
       end if
       offset = 0
    end if
    call mpi_exscan(counts, offset, 1, mpi_integer, &
         mpi_sum, neko_comm, ierr)
 
    ! ===========================
    ! Set up reading the data set
    ! ===========================
    dset_rank = 1 ! rank 1 array, i.e. a vector
    dcount = [int(counts, hsize_t)] ! local size of the vector
    doffset = [int(offset, hsize_t)] ! offset for this rank in the global vector
    ! Get the total file space (shape) of the data set
    call h5dget_space_f(dset_id, filespace, ierr)
    ! Get only the slice where my rank reads
    call h5sselect_hyperslab_f(filespace, h5s_select_set_f, doffset, dcount, &
         ierr)
    ! Create the corresponding memory space (buffer) for my local data
    call h5screate_simple_f(dset_rank, dcount, memspace, ierr)
 
    ! =============================
    ! Allocate data. HDF5 will cast
    ! =============================
    call vec%init(counts, trim(data_name)) ! this is rp
    call h5dread_f(dset_id, precision_hdf, vec%x, dcount, ierr, &
         file_space_id = filespace, mem_space_id = memspace, &
         xfer_prp = xf_id)
 
    ! =======================
    ! Clean up
    ! =======================
    call h5pclose_f(xf_id, ierr)
    call h5sclose_f(memspace, ierr)
    call h5sclose_f(filespace, ierr)
    call h5dclose_f(dset_id, ierr)
 
  end subroutine hdf5_file_read_vector
 
  subroutine hdf5_file_read_matrix(this, data_name, mat, strategy)
    class(hdf5_file_t) :: this
    character(len=*), intent(in) :: data_name
    type(matrix_t), intent(inout) :: mat
    character(len=*), intent(in), optional :: strategy
    character(len=1000) :: strategy_
    integer :: ierr, counts, offset, total_count, dset_rank
    integer(hid_t) :: precision_hdf
    integer(hid_t) :: xf_id, filespace, dset_id, memspace
    integer(hsize_t), dimension(2) :: dcount, doffset
    integer(hsize_t), dimension(2) :: tempddims, tempmaxddims
    integer :: temprank
    logical :: dset_exists
    type(linear_dist_t) :: dist
 
    ! Set up strategy
    if (present(strategy)) then
       if (trim(strategy) .eq. "linear" .or. &
            trim(strategy) .eq. "rank_0") then
          strategy_ = strategy
       else
          call neko_error("Unsupported strategy: " // trim(strategy))
       end if
    else
       strategy_ = "linear"
    end if
 
    ! Free the input
    call mat%free()
 
    ! ===============
    ! Configure MPIIO
    ! ===============
    call h5pcreate_f(h5p_dataset_xfer_f, xf_id, ierr)
    call h5pset_dxpl_mpio_f(xf_id, h5fd_mpio_collective_f, ierr)
    precision_hdf = h5kind_to_type(rp, h5_real_kind)
 
    ! ===================
    ! Get the data set info
    ! ===================
    call h5lexists_f(this%active_group_id, trim(data_name), dset_exists, ierr)
    if (dset_exists) then
       ! Openr the data set
       call h5dopen_f(this%active_group_id, trim(data_name), dset_id, ierr)
       ! Get the current rank of the of the dataset
       call h5dget_space_f(dset_id, filespace, ierr)
       call h5sget_simple_extent_ndims_f(filespace, temprank, ierr)
       if (temprank .ne. 2) then
          call neko_error("Dataset " // trim(data_name) // &
               " is not a rank 2 matrix in file " // trim(file_get_fname(this)))
       end if
       ! Get the current shape and close the filespace
       call h5sget_simple_extent_dims_f(filespace, tempddims, tempmaxddims, &
            ierr)
       call h5sclose_f(filespace, ierr)
    else
       call neko_error("Dataset " // trim(data_name) &
            // " does not exist in current group " // &
            trim(file_get_fname(this)))
    end if
 
    ! =============================
    ! Perform the data distribution
    ! =============================
    total_count = int(tempddims(2))
    if (strategy_ .eq. "linear") then
       dist = linear_dist_t(total_count, pe_rank, pe_size, neko_comm)
       counts = dist%num_local()
       offset = 0
    else if (strategy_ .eq. "rank_0") then
       if (pe_rank .eq. 0) then
          counts = total_count
       else
          counts = 0
       end if
       offset = 0
    end if
    call mpi_scan(counts, offset, 1, mpi_integer, &
         mpi_sum, neko_comm, ierr)
    offset = offset - counts ! Not using exclusive scan
 
    ! ===========================
    ! Set up reading the data set
    ! ===========================
    dset_rank = 2 ! rank 2 array, i.e. a matrix
    dcount = [int(tempddims(1), hsize_t), int(counts, hsize_t)] ! local size
    doffset = [0_hsize_t, int(offset, hsize_t)] ! offset in the global matrix
    ! Get the total file space (shape) of the data set
    call h5dget_space_f(dset_id, filespace, ierr)
    ! Get only the slice where my rank reads
    call h5sselect_hyperslab_f(filespace, h5s_select_set_f, doffset, dcount, &
         ierr)
    ! Create the corresponding memory space (buffer) for my local data
    call h5screate_simple_f(dset_rank, dcount, memspace, ierr)
 
    ! =============================
    ! Allocate data. HDF5 will cast
    ! =============================
    call mat%init(int(tempddims(1)), counts, trim(data_name)) ! this is rp
    call h5dread_f(dset_id, precision_hdf, mat%x, dcount, ierr, &
         file_space_id = filespace, mem_space_id = memspace, &
         xfer_prp = xf_id)
 
    ! =======================
    ! Clean up
    ! =======================
    call h5pclose_f(xf_id, ierr)
    call h5sclose_f(memspace, ierr)
    call h5sclose_f(filespace, ierr)
    call h5dclose_f(dset_id, ierr)
 
 
  end subroutine hdf5_file_read_matrix
 
  subroutine hdf5_file_write_int_attribute(this, attr_name, attr)
    class(hdf5_file_t), intent(inout) :: this
    character(len=*), intent(in) :: attr_name
    integer, intent(in) :: attr
    integer :: ierr
    integer(hid_t) :: filespace, attr_id
    integer(hsize_t), dimension(1) :: dcount
    logical :: attr_exists
 
    ! ====================
    ! Create the attribute
    ! ====================
    dcount = [int(1, hsize_t)]
    call h5aexists_f(this%active_group_id, trim(attr_name), attr_exists, ierr)
    if (attr_exists) then
       ! retrieve the attr id for the existing attribute
       call h5aopen_f(this%active_group_id, trim(attr_name), attr_id, ierr)
    else
       ! create file space of this shape
       call h5screate_f(h5s_scalar_f, filespace, ierr)
       ! create the data set with the given shape
       call h5acreate_f(this%active_group_id, trim(attr_name), &
            h5t_native_integer, &
            filespace, attr_id, ierr, h5p_default_f, h5p_default_f)
       call h5sclose_f(filespace, ierr)
    end if
 
    ! ===========================
    ! Set up writing the data set
    ! ===========================
    call h5awrite_f(attr_id, h5t_native_integer, attr, dcount, ierr)
 
    ! =======================
    ! Clean up
    ! =======================
    call h5aclose_f(attr_id, ierr)
 
  end subroutine hdf5_file_write_int_attribute
 
  subroutine hdf5_file_write_rp_attribute(this, attr_name, attr)
    class(hdf5_file_t), intent(inout) :: this
    character(len=*), intent(in) :: attr_name
    real(kind=rp), intent(in) :: attr
    integer :: ierr
    integer(hid_t) :: precision_hdf
    integer(hid_t) :: filespace, attr_id
    integer(hsize_t), dimension(1) :: dcount
    logical :: attr_exists
 
    ! Get the precision
    precision_hdf = h5kind_to_type(rp, h5_real_kind)
 
    ! ====================
    ! Create the attribute
    ! ====================
    dcount = [int(1, hsize_t)]
    call h5aexists_f(this%active_group_id, trim(attr_name), attr_exists, ierr)
    if (attr_exists) then
       ! retrieve the attr id for the existing attribute
       call h5aopen_f(this%active_group_id, trim(attr_name), attr_id, ierr)
    else
       ! create file space of this shape
       call h5screate_f(h5s_scalar_f, filespace, ierr)
       ! create the data set with the given shape
       call h5acreate_f(this%active_group_id, trim(attr_name), precision_hdf, &
            filespace, attr_id, ierr, h5p_default_f, h5p_default_f)
       call h5sclose_f(filespace, ierr)
    end if
 
    ! ===========================
    ! Set up writing the data set
    ! ===========================
    call h5awrite_f(attr_id, precision_hdf, attr, dcount, ierr)
 
    ! =======================
    ! Clean up
    ! =======================
    call h5aclose_f(attr_id, ierr)
 
  end subroutine hdf5_file_write_rp_attribute
 
  subroutine hdf5_file_read_int_attribute(this, attr_name, attr, attr_exists)
    class(hdf5_file_t), intent(inout) :: this
    character(len=*), intent(in) :: attr_name
    integer, intent(inout) :: attr
    logical, intent(inout) :: attr_exists
    integer :: ierr
    integer(hid_t) :: filespace, attr_id
    integer(hsize_t), dimension(1) :: dcount
 
    ! ====================
    ! Create the attribute
    ! ====================
    dcount = [int(1, hsize_t)]
    call h5aexists_f(this%active_group_id, trim(attr_name), attr_exists, ierr)
    if (attr_exists) then
       ! retrieve the attr id for the existing attribute
       call h5aopen_f(this%active_group_id, trim(attr_name), attr_id, ierr)
    else
       return
    end if
 
    ! ===========================
    ! Set up writing the data set
    ! ===========================
    call h5aread_f(attr_id, h5t_native_integer, attr, dcount, ierr)
 
    ! =======================
    ! Clean up
    ! =======================
    call h5aclose_f(attr_id, ierr)
 
  end subroutine hdf5_file_read_int_attribute
 
  subroutine hdf5_file_read_rp_attribute(this, attr_name, attr, attr_exists)
    class(hdf5_file_t), intent(inout) :: this
    character(len=*), intent(in) :: attr_name
    real(kind=rp), intent(inout) :: attr
    logical, intent(inout) :: attr_exists
    integer :: ierr
    integer(hid_t) :: precision_hdf
    integer(hid_t) :: filespace, attr_id
    integer(hsize_t), dimension(1) :: dcount
 
    ! Get the precision
    precision_hdf = h5kind_to_type(rp, h5_real_kind)
 
    ! ====================
    ! Create the attribute
    ! ====================
    dcount = [int(1, hsize_t)]
    call h5aexists_f(this%active_group_id, trim(attr_name), attr_exists, ierr)
    if (attr_exists) then
       ! retrieve the attr id for the existing attribute
       call h5aopen_f(this%active_group_id, trim(attr_name), attr_id, ierr)
    else
       return
    end if
 
    ! ===========================
    ! Set up writing the data set
    ! ===========================
    call h5aread_f(attr_id, precision_hdf, attr, dcount, ierr)
 
    ! =======================
    ! Clean up
    ! =======================
    call h5aclose_f(attr_id, ierr)
 
  end subroutine hdf5_file_read_rp_attribute
 
#else
 
  subroutine hdf5_file_open(this, mode)
    class(hdf5_file_t), intent(inout) :: this
    character(len=1), intent(in) :: mode
    call neko_error('Neko needs to be built with HDF5 support')
  end subroutine hdf5_file_open
 
  subroutine hdf5_file_close(this)
    class(hdf5_file_t), intent(inout) :: this
    call neko_error('Neko needs to be built with HDF5 support')
  end subroutine hdf5_file_close
 
  subroutine hdf5_file_set_group(this, group_name_path)
    class(hdf5_file_t), intent(inout) :: this
    character(len=*), intent(in), optional :: group_name_path
    call neko_error('Neko needs to be built with HDF5 support')
  end subroutine hdf5_file_set_group
 
  subroutine hdf5_file_write(this, data, t)
    class(hdf5_file_t), intent(inout) :: this
    class(*), target, intent(in) :: data
    real(kind=dp), intent(in), optional :: t
    call neko_error('Neko needs to be built with HDF5 support')
  end subroutine hdf5_file_write
 
  subroutine hdf5_file_read(this, data)
    class(hdf5_file_t) :: this
    class(*), target, intent(inout) :: data
    call neko_error('Neko needs to be built with HDF5 support')
  end subroutine hdf5_file_read
 
  subroutine hdf5_file_write_dataset(this, data)
    class(hdf5_file_t), intent(inout) :: this
    class(*), intent(inout) :: data
    call neko_error('Neko needs to be built with HDF5 support')
  end subroutine hdf5_file_write_dataset
 
  subroutine hdf5_file_read_dataset(this, data_name, data, strategy)
    class(hdf5_file_t), intent(inout) :: this
    character(len=*), intent(in) :: data_name
    class(*), intent(inout) :: data
    character(len=*), intent(in), optional :: strategy
    call neko_error('Neko needs to be built with HDF5 support')
  end subroutine hdf5_file_read_dataset
 
  subroutine hdf5_file_write_attribute(this, data_name, data)
    class(hdf5_file_t), intent(inout) :: this
    character(len=*), intent(in) :: data_name
    class(*), intent(inout) :: data
    call neko_error('Neko needs to be built with HDF5 support')
  end subroutine hdf5_file_write_attribute
 
  subroutine hdf5_file_read_attribute(this, data_name, data, exist)
    class(hdf5_file_t), intent(inout) :: this
    character(len=*), intent(in) :: data_name
    class(*), intent(inout) :: data
    logical, intent(inout) :: exist
    call neko_error('Neko needs to be built with HDF5 support')
  end subroutine hdf5_file_read_attribute
 
  subroutine hdf5_file_write_vector(this, vec)
    class(hdf5_file_t), intent(inout) :: this
    type(vector_t), intent(inout) :: vec
    call neko_error('Neko needs to be built with HDF5 support')
  end subroutine hdf5_file_write_vector
 
  subroutine hdf5_file_write_matrix(this, mat)
    class(hdf5_file_t), intent(inout) :: this
    type(matrix_t), intent(inout) :: mat
    call neko_error('Neko needs to be built with HDF5 support')
  end subroutine hdf5_file_write_matrix
 
  subroutine hdf5_file_write_field(this, fld)
    class(hdf5_file_t), intent(inout) :: this
    type(field_t), intent(inout) :: fld
    call neko_error('Neko needs to be built with HDF5 support')
  end subroutine hdf5_file_write_field
 
  subroutine hdf5_file_read_vector(this, data_name, vec, strategy)
    class(hdf5_file_t) :: this
    character(len=*), intent(in) :: data_name
    type(vector_t), intent(inout) :: vec
    character(len=*), intent(in), optional :: strategy
    call neko_error('Neko needs to be built with HDF5 support')
  end subroutine hdf5_file_read_vector
 
  subroutine hdf5_file_read_matrix(this, data_name, mat, strategy)
    class(hdf5_file_t) :: this
    character(len=*), intent(in) :: data_name
    type(matrix_t), intent(inout) :: mat
    character(len=*), intent(in), optional :: strategy
    call neko_error('Neko needs to be built with HDF5 support')
  end subroutine hdf5_file_read_matrix
 
  subroutine hdf5_file_write_int_attribute(this, attr_name, attr)
    class(hdf5_file_t), intent(inout) :: this
    character(len=*), intent(in) :: attr_name
    integer, intent(in) :: attr
    call neko_error('Neko needs to be built with HDF5 support')
  end subroutine hdf5_file_write_int_attribute
 
  subroutine hdf5_file_write_rp_attribute(this, attr_name, attr)
    class(hdf5_file_t), intent(inout) :: this
    character(len=*), intent(in) :: attr_name
    real(kind=rp), intent(in) :: attr
    call neko_error('Neko needs to be built with HDF5 support')
  end subroutine hdf5_file_write_rp_attribute
 
  subroutine hdf5_file_read_int_attribute(this, attr_name, attr, attr_exists)
    class(hdf5_file_t), intent(inout) :: this
    character(len=*), intent(in) :: attr_name
    integer, intent(inout) :: attr
    logical, intent(inout) :: attr_exists
    call neko_error('Neko needs to be built with HDF5 support')
  end subroutine hdf5_file_read_int_attribute
 
  subroutine hdf5_file_read_rp_attribute(this, attr_name, attr, attr_exists)
    class(hdf5_file_t), intent(inout) :: this
    character(len=*), intent(in) :: attr_name
    real(kind=rp), intent(inout) :: attr
    logical, intent(inout) :: attr_exists
    call neko_error('Neko needs to be built with HDF5 support')
  end subroutine hdf5_file_read_rp_attribute
 
#endif
 
end module hdf5_file