From 8b1785c05ccca92fadb1339255bf2d13a334feec Mon Sep 17 00:00:00 2001
From: Martin Diehl <m.diehl@mpie.de>
Date: Sun, 18 Nov 2018 09:47:50 +0100
Subject: [PATCH] parallel writing for pReal

---
 src/HDF5_utilities.f90 | 582 ++++++++++++++++++++++++++++++++---------
 1 file changed, 463 insertions(+), 119 deletions(-)

diff --git a/src/HDF5_utilities.f90 b/src/HDF5_utilities.f90
index 4c29fe980..e36e39e29 100644
--- a/src/HDF5_utilities.f90
+++ b/src/HDF5_utilities.f90
@@ -542,265 +542,609 @@ subroutine HDF5_read_pInt_7(dataset,loc_id,datasetName)
 
 end subroutine HDF5_read_pInt_7
 
+
 !--------------------------------------------------------------------------------------------------
-!> @brief subroutine for writing dataset of the type pReal with 1 dimensions
+!> @brief subroutine for writing dataset of type pReal with 1 dimensions
 !--------------------------------------------------------------------------------------------------
-subroutine HDF5_write_pReal1(dataset,loc_id,datasetName)
+subroutine HDF5_write_pReal1(dataset,loc_id,datasetName,parallel)
+ use numerics, only: &
+   worldrank, &
+   worldsize
 
  implicit none
  real(pReal),      intent(out), dimension(:) ::    dataset
  integer(HID_T),   intent(in) :: loc_id                                                             !< file or group handle
  character(len=*), intent(in) :: datasetName                                                        !< name of the dataset in the file
+ logical, intent(in), optional :: parallel
 
- integer(pInt), dimension(:), allocatable :: myShape                                                          !<shape of the dataset
- integer :: hdferr
- integer(HID_T) :: dset_id, space_id
 
- myShape = shape(dataset)
+ integer(pInt), dimension(:), allocatable :: &
+   globalShape, &                                                                                   !< shape of the dataset (all processes)
+   localShape, &                                                                                    !< shape of the dataset (this process)
+   outputSize                                                                                       !< contribution of all processes
+ integer :: hdferr,ierr
+ integer(HID_T)   :: dset_id, filespace_id, memspace_id, plist_id
+ integer(HSIZE_T), dimension(1) :: myStart
+
+ call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, ierr)
+
 !--------------------------------------------------------------------------------------------------
-! create dataspace
-call h5screate_simple_f(size(myShape), int(myShape,HSIZE_T), space_id, hdferr, &
-                           int(myShape,HSIZE_T))
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal1: h5screate_simple_f')
+! determine shape of dataset
+ localShape = shape(dataset)
+ allocate(outputSize(worldsize), source = 0_pInt)
+ outputSize(worldrank+1) = localShape(1)
+#ifdef PETSc
+ if (present(parallel)) then; if (parallel) then
+   call h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, hdferr)
+   if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal1: h5pset_dxpl_mpio_f')
+   call MPI_allreduce(MPI_IN_PLACE,outputSize,worldsize,MPI_LONG,MPI_SUM,PETSC_COMM_WORLD,ierr)       ! get total output size over each process
+   if (ierr /= 0) call IO_error(894_pInt,ext_msg='HDF5_write_pReal1: MPI_allreduce')
+ endif; endif
+#endif
+ myStart     = int([sum(outputSize(1:worldrank))],HSIZE_T)
+ globalShape = [localShape(1:0),sum(outputSize)]
+
+
+!--------------------------------------------------------------------------------------------------
+! create dataspace in memory (local shape)
+ call h5screate_simple_f(size(localShape), int(localShape,HSIZE_T), memspace_id, hdferr, &
+                         int(localShape,HSIZE_T))
+ if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal1: h5screate_simple_f/memspace_id')
+
+!--------------------------------------------------------------------------------------------------
+! create dataspace in file (global shape)
+ call h5screate_simple_f(size(globalShape), int(globalShape,HSIZE_T), filespace_id, hdferr, &
+                         int(globalShape,HSIZE_T))
+ if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal1: h5screate_simple_f/filespace_id')
+
 !--------------------------------------------------------------------------------------------------
 ! create dataset
-call h5dcreate_f(loc_id, trim(datasetName), H5T_NATIVE_DOUBLE, space_id, dset_id, hdferr)
-if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal1: h5dcreate_f')
+ call h5dcreate_f(loc_id, trim(datasetName), H5T_NATIVE_DOUBLE, filespace_id, dset_id, hdferr)
+ if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal1: h5dcreate_f')
 
-CALL h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE,dataset,int(myShape,HSIZE_T), hdferr)
-if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal1: h5dcreate_f')
+!--------------------------------------------------------------------------------------------------
+! select a hyperslab (the portion of the current process) in the file
+ call h5dget_space_f(dset_id, filespace_id, hdferr)
+ if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal1: h5dget_space_f')
+ call h5sselect_hyperslab_f(filespace_id, H5S_SELECT_SET_F, myStart, int(localShape,HSIZE_T), hdferr)
+ if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal1: h5sselect_hyperslab_f')
+
+!--------------------------------------------------------------------------------------------------
+! write
+ call h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE,dataset,int(globalShape,HSIZE_T), hdferr,&
+                file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id)
+ if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal1: h5dcreate_f')
 
 !--------------------------------------------------------------------------------------------------
 !close types, dataspaces
+ call h5pclose_f(plist_id, hdferr)
+ if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal1: plist_id')
  call h5dclose_f(dset_id, hdferr)
-if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_HDF5_write_pReal1: h5dclose_f')
-call h5sclose_f(space_id, hdferr)
-if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_HDF5_write_pReal1: h5sclose_f')
+ if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal1: h5dclose_f')
+ call h5sclose_f(filespace_id, hdferr)
+ if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal1: h5sclose_f/filespace_id')
+ call h5sclose_f(memspace_id, hdferr)
+ if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal1: h5sclose_f/memspace_id')
 
 end subroutine HDF5_write_pReal1
 
+
 !--------------------------------------------------------------------------------------------------
-!> @brief subroutine for writing dataset of the type pReal with 2 dimensions
+!> @brief subroutine for writing dataset of type pReal with 2 dimensions
 !--------------------------------------------------------------------------------------------------
-subroutine HDF5_write_pReal2(dataset,loc_id,datasetName)
+subroutine HDF5_write_pReal2(dataset,loc_id,datasetName,parallel)
+ use numerics, only: &
+   worldrank, &
+   worldsize
 
  implicit none
  real(pReal),      intent(out), dimension(:,:) ::    dataset
  integer(HID_T),   intent(in) :: loc_id                                                             !< file or group handle
  character(len=*), intent(in) :: datasetName                                                        !< name of the dataset in the file
+ logical, intent(in), optional :: parallel
 
- integer(pInt), dimension(:), allocatable :: myShape                                                          !<shape of the dataset
- integer :: hdferr
- integer(HID_T) :: dset_id, space_id
 
- myShape = shape(dataset)
+ integer(pInt), dimension(:), allocatable :: &
+   globalShape, &                                                                                   !< shape of the dataset (all processes)
+   localShape, &                                                                                    !< shape of the dataset (this process)
+   outputSize                                                                                       !< contribution of all processes
+ integer :: hdferr,ierr
+ integer(HID_T)   :: dset_id, filespace_id, memspace_id, plist_id
+ integer(HSIZE_T), dimension(2) :: myStart
+
+ call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, ierr)
+
 !--------------------------------------------------------------------------------------------------
-! create dataspace
-call h5screate_simple_f(size(myShape), int(myShape,HSIZE_T), space_id, hdferr, &
-                           int(myShape,HSIZE_T))
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal2: h5screate_simple_f')
+! determine shape of dataset
+ localShape = shape(dataset)
+ allocate(outputSize(worldsize), source = 0_pInt)
+ outputSize(worldrank+1) = localShape(2)
+#ifdef PETSc
+ if (present(parallel)) then; if (parallel) then
+   call h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, hdferr)
+   if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal2: h5pset_dxpl_mpio_f')
+   call MPI_allreduce(MPI_IN_PLACE,outputSize,worldsize,MPI_LONG,MPI_SUM,PETSC_COMM_WORLD,ierr)       ! get total output size over each process
+   if (ierr /= 0) call IO_error(894_pInt,ext_msg='HDF5_write_pReal2: MPI_allreduce')
+ endif; endif
+#endif
+ myStart     = int([0,sum(outputSize(1:worldrank))],HSIZE_T)
+ globalShape = [localShape(1:1),sum(outputSize)]
+
+
+!--------------------------------------------------------------------------------------------------
+! create dataspace in memory (local shape)
+ call h5screate_simple_f(size(localShape), int(localShape,HSIZE_T), memspace_id, hdferr, &
+                         int(localShape,HSIZE_T))
+ if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal2: h5screate_simple_f/memspace_id')
+
+!--------------------------------------------------------------------------------------------------
+! create dataspace in file (global shape)
+ call h5screate_simple_f(size(globalShape), int(globalShape,HSIZE_T), filespace_id, hdferr, &
+                         int(globalShape,HSIZE_T))
+ if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal2: h5screate_simple_f/filespace_id')
+
 !--------------------------------------------------------------------------------------------------
 ! create dataset
-call h5dcreate_f(loc_id, trim(datasetName), H5T_NATIVE_DOUBLE, space_id, dset_id, hdferr)
-if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal2: h5dcreate_f')
+ call h5dcreate_f(loc_id, trim(datasetName), H5T_NATIVE_DOUBLE, filespace_id, dset_id, hdferr)
+ if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal2: h5dcreate_f')
 
-CALL h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE,dataset,int(myShape,HSIZE_T), hdferr)
-if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal2: h5dcreate_f')
+!--------------------------------------------------------------------------------------------------
+! select a hyperslab (the portion of the current process) in the file
+ call h5dget_space_f(dset_id, filespace_id, hdferr)
+ if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal2: h5dget_space_f')
+ call h5sselect_hyperslab_f(filespace_id, H5S_SELECT_SET_F, myStart, int(localShape,HSIZE_T), hdferr)
+ if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal2: h5sselect_hyperslab_f')
+
+!--------------------------------------------------------------------------------------------------
+! write
+ call h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE,dataset,int(globalShape,HSIZE_T), hdferr,&
+                file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id)
+ if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal2: h5dcreate_f')
 
 !--------------------------------------------------------------------------------------------------
 !close types, dataspaces
+ call h5pclose_f(plist_id, hdferr)
+ if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal2: plist_id')
  call h5dclose_f(dset_id, hdferr)
-if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_HDF5_write_pReal2: h5dclose_f')
-call h5sclose_f(space_id, hdferr)
-if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_HDF5_write_pReal2: h5sclose_f')
+ if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal2: h5dclose_f')
+ call h5sclose_f(filespace_id, hdferr)
+ if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal2: h5sclose_f/filespace_id')
+ call h5sclose_f(memspace_id, hdferr)
+ if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal2: h5sclose_f/memspace_id')
 
 end subroutine HDF5_write_pReal2
 
+
 !--------------------------------------------------------------------------------------------------
-!> @brief subroutine for writing dataset of the type pReal with 3 dimensions
+!> @brief subroutine for writing dataset of type pReal with 3 dimensions
 !--------------------------------------------------------------------------------------------------
-subroutine HDF5_write_pReal3(dataset,loc_id,datasetName)
+subroutine HDF5_write_pReal3(dataset,loc_id,datasetName,parallel)
+ use numerics, only: &
+   worldrank, &
+   worldsize
 
  implicit none
  real(pReal),      intent(out), dimension(:,:,:) ::    dataset
  integer(HID_T),   intent(in) :: loc_id                                                             !< file or group handle
  character(len=*), intent(in) :: datasetName                                                        !< name of the dataset in the file
+ logical, intent(in), optional :: parallel
 
- integer(pInt), dimension(:), allocatable :: myShape                                                          !<shape of the dataset
- integer :: hdferr
- integer(HID_T) :: dset_id, space_id
 
- myShape = shape(dataset)
+ integer(pInt), dimension(:), allocatable :: &
+   globalShape, &                                                                                   !< shape of the dataset (all processes)
+   localShape, &                                                                                    !< shape of the dataset (this process)
+   outputSize                                                                                       !< contribution of all processes
+ integer :: hdferr,ierr
+ integer(HID_T)   :: dset_id, filespace_id, memspace_id, plist_id
+ integer(HSIZE_T), dimension(3) :: myStart
+
+ call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, ierr)
+
 !--------------------------------------------------------------------------------------------------
-! create dataspace
-call h5screate_simple_f(size(myShape), int(myShape,HSIZE_T), space_id, hdferr, &
-                           int(myShape,HSIZE_T))
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal3: h5screate_simple_f')
+! determine shape of dataset
+ localShape = shape(dataset)
+ allocate(outputSize(worldsize), source = 0_pInt)
+ outputSize(worldrank+1) = localShape(3)
+#ifdef PETSc
+ if (present(parallel)) then; if (parallel) then
+   call h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, hdferr)
+   if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal3: h5pset_dxpl_mpio_f')
+   call MPI_allreduce(MPI_IN_PLACE,outputSize,worldsize,MPI_LONG,MPI_SUM,PETSC_COMM_WORLD,ierr)       ! get total output size over each process
+   if (ierr /= 0) call IO_error(894_pInt,ext_msg='HDF5_write_pReal3: MPI_allreduce')
+ endif; endif
+#endif
+ myStart     = int([0,0,sum(outputSize(1:worldrank))],HSIZE_T)
+ globalShape = [localShape(1:2),sum(outputSize)]
+
+
+!--------------------------------------------------------------------------------------------------
+! create dataspace in memory (local shape)
+ call h5screate_simple_f(size(localShape), int(localShape,HSIZE_T), memspace_id, hdferr, &
+                         int(localShape,HSIZE_T))
+ if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal3: h5screate_simple_f/memspace_id')
+
+!--------------------------------------------------------------------------------------------------
+! create dataspace in file (global shape)
+ call h5screate_simple_f(size(globalShape), int(globalShape,HSIZE_T), filespace_id, hdferr, &
+                         int(globalShape,HSIZE_T))
+ if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal3: h5screate_simple_f/filespace_id')
+
 !--------------------------------------------------------------------------------------------------
 ! create dataset
-call h5dcreate_f(loc_id, trim(datasetName), H5T_NATIVE_DOUBLE, space_id, dset_id, hdferr)
-if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal3: h5dcreate_f')
+ call h5dcreate_f(loc_id, trim(datasetName), H5T_NATIVE_DOUBLE, filespace_id, dset_id, hdferr)
+ if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal3: h5dcreate_f')
 
-CALL h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE,dataset,int(myShape,HSIZE_T), hdferr)
-if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal3: h5dcreate_f')
+!--------------------------------------------------------------------------------------------------
+! select a hyperslab (the portion of the current process) in the file
+ call h5dget_space_f(dset_id, filespace_id, hdferr)
+ if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal3: h5dget_space_f')
+ call h5sselect_hyperslab_f(filespace_id, H5S_SELECT_SET_F, myStart, int(localShape,HSIZE_T), hdferr)
+ if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal3: h5sselect_hyperslab_f')
+
+!--------------------------------------------------------------------------------------------------
+! write
+ call h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE,dataset,int(globalShape,HSIZE_T), hdferr,&
+                file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id)
+ if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal3: h5dcreate_f')
 
 !--------------------------------------------------------------------------------------------------
 !close types, dataspaces
+ call h5pclose_f(plist_id, hdferr)
+ if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal3: plist_id')
  call h5dclose_f(dset_id, hdferr)
-if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_HDF5_write_pReal3: h5dclose_f')
-call h5sclose_f(space_id, hdferr)
-if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_HDF5_write_pReal3: h5sclose_f')
+ if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal3: h5dclose_f')
+ call h5sclose_f(filespace_id, hdferr)
+ if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal3: h5sclose_f/filespace_id')
+ call h5sclose_f(memspace_id, hdferr)
+ if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal3: h5sclose_f/memspace_id')
 
 end subroutine HDF5_write_pReal3
 
+
 !--------------------------------------------------------------------------------------------------
-!> @brief subroutine for writing dataset of the type pReal with 4 dimensions
+!> @brief subroutine for writing dataset of type pReal with 4 dimensions
 !--------------------------------------------------------------------------------------------------
-subroutine HDF5_write_pReal4(dataset,loc_id,datasetName)
+subroutine HDF5_write_pReal4(dataset,loc_id,datasetName,parallel)
+ use numerics, only: &
+   worldrank, &
+   worldsize
 
  implicit none
  real(pReal),      intent(out), dimension(:,:,:,:) ::    dataset
  integer(HID_T),   intent(in) :: loc_id                                                             !< file or group handle
  character(len=*), intent(in) :: datasetName                                                        !< name of the dataset in the file
+ logical, intent(in), optional :: parallel
 
- integer(pInt), dimension(:), allocatable :: myShape                                                          !<shape of the dataset
- integer :: hdferr
- integer(HID_T) :: dset_id, space_id
 
- myShape = shape(dataset)
+ integer(pInt), dimension(:), allocatable :: &
+   globalShape, &                                                                                   !< shape of the dataset (all processes)
+   localShape, &                                                                                    !< shape of the dataset (this process)
+   outputSize                                                                                       !< contribution of all processes
+ integer :: hdferr,ierr
+ integer(HID_T)   :: dset_id, filespace_id, memspace_id, plist_id
+ integer(HSIZE_T), dimension(4) :: myStart
+
+ call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, ierr)
+
 !--------------------------------------------------------------------------------------------------
-! create dataspace
-call h5screate_simple_f(size(myShape), int(myShape,HSIZE_T), space_id, hdferr, &
-                           int(myShape,HSIZE_T))
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal4: h5screate_simple_f')
+! determine shape of dataset
+ localShape = shape(dataset)
+ allocate(outputSize(worldsize), source = 0_pInt)
+ outputSize(worldrank+1) = localShape(4)
+#ifdef PETSc
+ if (present(parallel)) then; if (parallel) then
+   call h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, hdferr)
+   if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal4: h5pset_dxpl_mpio_f')
+   call MPI_allreduce(MPI_IN_PLACE,outputSize,worldsize,MPI_LONG,MPI_SUM,PETSC_COMM_WORLD,ierr)       ! get total output size over each process
+   if (ierr /= 0) call IO_error(894_pInt,ext_msg='HDF5_write_pReal4: MPI_allreduce')
+ endif; endif
+#endif
+ myStart     = int([0,0,0,sum(outputSize(1:worldrank))],HSIZE_T)
+ globalShape = [localShape(1:3),sum(outputSize)]
+
+
+!--------------------------------------------------------------------------------------------------
+! create dataspace in memory (local shape)
+ call h5screate_simple_f(size(localShape), int(localShape,HSIZE_T), memspace_id, hdferr, &
+                         int(localShape,HSIZE_T))
+ if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal4: h5screate_simple_f/memspace_id')
+
+!--------------------------------------------------------------------------------------------------
+! create dataspace in file (global shape)
+ call h5screate_simple_f(size(globalShape), int(globalShape,HSIZE_T), filespace_id, hdferr, &
+                         int(globalShape,HSIZE_T))
+ if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal4: h5screate_simple_f/filespace_id')
+
 !--------------------------------------------------------------------------------------------------
 ! create dataset
-call h5dcreate_f(loc_id, trim(datasetName), H5T_NATIVE_DOUBLE, space_id, dset_id, hdferr)
-if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal4: h5dcreate_f')
+ call h5dcreate_f(loc_id, trim(datasetName), H5T_NATIVE_DOUBLE, filespace_id, dset_id, hdferr)
+ if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal4: h5dcreate_f')
 
-CALL h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE,dataset,int(myShape,HSIZE_T), hdferr)
-if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal4: h5dcreate_f')
+!--------------------------------------------------------------------------------------------------
+! select a hyperslab (the portion of the current process) in the file
+ call h5dget_space_f(dset_id, filespace_id, hdferr)
+ if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal4: h5dget_space_f')
+ call h5sselect_hyperslab_f(filespace_id, H5S_SELECT_SET_F, myStart, int(localShape,HSIZE_T), hdferr)
+ if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal4: h5sselect_hyperslab_f')
+
+!--------------------------------------------------------------------------------------------------
+! write
+ call h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE,dataset,int(globalShape,HSIZE_T), hdferr,&
+                file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id)
+ if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal4: h5dcreate_f')
 
 !--------------------------------------------------------------------------------------------------
 !close types, dataspaces
+ call h5pclose_f(plist_id, hdferr)
+ if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal4: plist_id')
  call h5dclose_f(dset_id, hdferr)
-if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_HDF5_write_pReal4: h5dclose_f')
-call h5sclose_f(space_id, hdferr)
-if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_HDF5_write_pReal4: h5sclose_f')
+ if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal4: h5dclose_f')
+ call h5sclose_f(filespace_id, hdferr)
+ if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal4: h5sclose_f/filespace_id')
+ call h5sclose_f(memspace_id, hdferr)
+ if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal4: h5sclose_f/memspace_id')
 
 end subroutine HDF5_write_pReal4
 
+
 !--------------------------------------------------------------------------------------------------
-!> @brief subroutine for writing dataset of the type pReal with 5 dimensions
+!> @brief subroutine for writing dataset of type pReal with 5 dimensions
 !--------------------------------------------------------------------------------------------------
-subroutine HDF5_write_pReal5(dataset,loc_id,datasetName)
+subroutine HDF5_write_pReal5(dataset,loc_id,datasetName,parallel)
+ use numerics, only: &
+   worldrank, &
+   worldsize
 
  implicit none
  real(pReal),      intent(out), dimension(:,:,:,:,:) ::    dataset
  integer(HID_T),   intent(in) :: loc_id                                                             !< file or group handle
  character(len=*), intent(in) :: datasetName                                                        !< name of the dataset in the file
+ logical, intent(in), optional :: parallel
 
- integer(pInt), dimension(:), allocatable :: myShape                                                          !<shape of the dataset
- integer :: hdferr
- integer(HID_T) :: dset_id, space_id
 
- myShape = shape(dataset)
+ integer(pInt), dimension(:), allocatable :: &
+   globalShape, &                                                                                   !< shape of the dataset (all processes)
+   localShape, &                                                                                    !< shape of the dataset (this process)
+   outputSize                                                                                       !< contribution of all processes
+ integer :: hdferr,ierr
+ integer(HID_T)   :: dset_id, filespace_id, memspace_id, plist_id
+ integer(HSIZE_T), dimension(5) :: myStart
+
+ call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, ierr)
+
 !--------------------------------------------------------------------------------------------------
-! create dataspace
-call h5screate_simple_f(size(myShape), int(myShape,HSIZE_T), space_id, hdferr, &
-                           int(myShape,HSIZE_T))
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal5: h5screate_simple_f')
+! determine shape of dataset
+ localShape = shape(dataset)
+ allocate(outputSize(worldsize), source = 0_pInt)
+ outputSize(worldrank+1) = localShape(5)
+#ifdef PETSc
+ if (present(parallel)) then; if (parallel) then
+   call h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, hdferr)
+   if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal5: h5pset_dxpl_mpio_f')
+   call MPI_allreduce(MPI_IN_PLACE,outputSize,worldsize,MPI_LONG,MPI_SUM,PETSC_COMM_WORLD,ierr)       ! get total output size over each process
+   if (ierr /= 0) call IO_error(894_pInt,ext_msg='HDF5_write_pReal5: MPI_allreduce')
+ endif; endif
+#endif
+ myStart     = int([0,0,0,0,sum(outputSize(1:worldrank))],HSIZE_T)
+ globalShape = [localShape(1:4),sum(outputSize)]
+
+
+!--------------------------------------------------------------------------------------------------
+! create dataspace in memory (local shape)
+ call h5screate_simple_f(size(localShape), int(localShape,HSIZE_T), memspace_id, hdferr, &
+                         int(localShape,HSIZE_T))
+ if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal5: h5screate_simple_f/memspace_id')
+
+!--------------------------------------------------------------------------------------------------
+! create dataspace in file (global shape)
+ call h5screate_simple_f(size(globalShape), int(globalShape,HSIZE_T), filespace_id, hdferr, &
+                         int(globalShape,HSIZE_T))
+ if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal5: h5screate_simple_f/filespace_id')
+
 !--------------------------------------------------------------------------------------------------
 ! create dataset
-call h5dcreate_f(loc_id, trim(datasetName), H5T_NATIVE_DOUBLE, space_id, dset_id, hdferr)
-if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal5: h5dcreate_f')
+ call h5dcreate_f(loc_id, trim(datasetName), H5T_NATIVE_DOUBLE, filespace_id, dset_id, hdferr)
+ if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal5: h5dcreate_f')
 
-CALL h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE,dataset,int(myShape,HSIZE_T), hdferr)
-if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal5: h5dcreate_f')
+!--------------------------------------------------------------------------------------------------
+! select a hyperslab (the portion of the current process) in the file
+ call h5dget_space_f(dset_id, filespace_id, hdferr)
+ if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal5: h5dget_space_f')
+ call h5sselect_hyperslab_f(filespace_id, H5S_SELECT_SET_F, myStart, int(localShape,HSIZE_T), hdferr)
+ if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal5: h5sselect_hyperslab_f')
+
+!--------------------------------------------------------------------------------------------------
+! write
+ call h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE,dataset,int(globalShape,HSIZE_T), hdferr,&
+                file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id)
+ if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal5: h5dcreate_f')
 
 !--------------------------------------------------------------------------------------------------
 !close types, dataspaces
+ call h5pclose_f(plist_id, hdferr)
+ if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal5: plist_id')
  call h5dclose_f(dset_id, hdferr)
-if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_HDF5_write_pReal5: h5dclose_f')
-call h5sclose_f(space_id, hdferr)
-if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_HDF5_write_pReal5: h5sclose_f')
+ if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal5: h5dclose_f')
+ call h5sclose_f(filespace_id, hdferr)
+ if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal5: h5sclose_f/filespace_id')
+ call h5sclose_f(memspace_id, hdferr)
+ if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal5: h5sclose_f/memspace_id')
 
 end subroutine HDF5_write_pReal5
 
+
 !--------------------------------------------------------------------------------------------------
-!> @brief subroutine for writing dataset of the type pReal with 6 dimensions
+!> @brief subroutine for writing dataset of type pReal with 6 dimensions
 !--------------------------------------------------------------------------------------------------
-subroutine HDF5_write_pReal6(dataset,loc_id,datasetName)
+subroutine HDF5_write_pReal6(dataset,loc_id,datasetName,parallel)
+ use numerics, only: &
+   worldrank, &
+   worldsize
 
  implicit none
  real(pReal),      intent(out), dimension(:,:,:,:,:,:) ::    dataset
  integer(HID_T),   intent(in) :: loc_id                                                             !< file or group handle
  character(len=*), intent(in) :: datasetName                                                        !< name of the dataset in the file
+ logical, intent(in), optional :: parallel
 
- integer(pInt), dimension(:), allocatable :: myShape                                                          !<shape of the dataset
- integer :: hdferr
- integer(HID_T) :: dset_id, space_id
 
- myShape = shape(dataset)
+ integer(pInt), dimension(:), allocatable :: &
+   globalShape, &                                                                                   !< shape of the dataset (all processes)
+   localShape, &                                                                                    !< shape of the dataset (this process)
+   outputSize                                                                                       !< contribution of all processes
+ integer :: hdferr,ierr
+ integer(HID_T)   :: dset_id, filespace_id, memspace_id, plist_id
+ integer(HSIZE_T), dimension(6) :: myStart
+
+ call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, ierr)
+
 !--------------------------------------------------------------------------------------------------
-! create dataspace
-call h5screate_simple_f(size(myShape), int(myShape,HSIZE_T), space_id, hdferr, &
-                           int(myShape,HSIZE_T))
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal6: h5screate_simple_f')
+! determine shape of dataset
+ localShape = shape(dataset)
+ allocate(outputSize(worldsize), source = 0_pInt)
+ outputSize(worldrank+1) = localShape(6)
+#ifdef PETSc
+ if (present(parallel)) then; if (parallel) then
+   call h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, hdferr)
+   if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal6: h5pset_dxpl_mpio_f')
+   call MPI_allreduce(MPI_IN_PLACE,outputSize,worldsize,MPI_LONG,MPI_SUM,PETSC_COMM_WORLD,ierr)       ! get total output size over each process
+   if (ierr /= 0) call IO_error(894_pInt,ext_msg='HDF5_write_pReal6: MPI_allreduce')
+ endif; endif
+#endif
+ myStart     = int([0,0,0,0,0,sum(outputSize(1:worldrank))],HSIZE_T)
+ globalShape = [localShape(1:5),sum(outputSize)]
+
+
+!--------------------------------------------------------------------------------------------------
+! create dataspace in memory (local shape)
+ call h5screate_simple_f(size(localShape), int(localShape,HSIZE_T), memspace_id, hdferr, &
+                         int(localShape,HSIZE_T))
+ if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal6: h5screate_simple_f/memspace_id')
+
+!--------------------------------------------------------------------------------------------------
+! create dataspace in file (global shape)
+ call h5screate_simple_f(size(globalShape), int(globalShape,HSIZE_T), filespace_id, hdferr, &
+                         int(globalShape,HSIZE_T))
+ if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal6: h5screate_simple_f/filespace_id')
+
 !--------------------------------------------------------------------------------------------------
 ! create dataset
-call h5dcreate_f(loc_id, trim(datasetName), H5T_NATIVE_DOUBLE, space_id, dset_id, hdferr)
-if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal6: h5dcreate_f')
+ call h5dcreate_f(loc_id, trim(datasetName), H5T_NATIVE_DOUBLE, filespace_id, dset_id, hdferr)
+ if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal6: h5dcreate_f')
 
-CALL h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE,dataset,int(myShape,HSIZE_T), hdferr)
-if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal6: h5dcreate_f')
+!--------------------------------------------------------------------------------------------------
+! select a hyperslab (the portion of the current process) in the file
+ call h5dget_space_f(dset_id, filespace_id, hdferr)
+ if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal6: h5dget_space_f')
+ call h5sselect_hyperslab_f(filespace_id, H5S_SELECT_SET_F, myStart, int(localShape,HSIZE_T), hdferr)
+ if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal6: h5sselect_hyperslab_f')
+
+!--------------------------------------------------------------------------------------------------
+! write
+ call h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE,dataset,int(globalShape,HSIZE_T), hdferr,&
+                file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id)
+ if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal6: h5dcreate_f')
 
 !--------------------------------------------------------------------------------------------------
 !close types, dataspaces
+ call h5pclose_f(plist_id, hdferr)
+ if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal6: plist_id')
  call h5dclose_f(dset_id, hdferr)
-if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_HDF5_write_pReal6: h5dclose_f')
-call h5sclose_f(space_id, hdferr)
-if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_HDF5_write_pReal6: h5sclose_f')
+ if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal6: h5dclose_f')
+ call h5sclose_f(filespace_id, hdferr)
+ if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal6: h5sclose_f/filespace_id')
+ call h5sclose_f(memspace_id, hdferr)
+ if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal6: h5sclose_f/memspace_id')
 
 end subroutine HDF5_write_pReal6
 
+
 !--------------------------------------------------------------------------------------------------
-!> @brief subroutine for writing dataset of the type pReal with 7 dimensions
+!> @brief subroutine for writing dataset of type pReal with 7 dimensions
 !--------------------------------------------------------------------------------------------------
-subroutine HDF5_write_pReal7(dataset,loc_id,datasetName)
+subroutine HDF5_write_pReal7(dataset,loc_id,datasetName,parallel)
+ use numerics, only: &
+   worldrank, &
+   worldsize
 
  implicit none
  real(pReal),      intent(out), dimension(:,:,:,:,:,:,:) ::    dataset
  integer(HID_T),   intent(in) :: loc_id                                                             !< file or group handle
  character(len=*), intent(in) :: datasetName                                                        !< name of the dataset in the file
+ logical, intent(in), optional :: parallel
 
- integer(pInt), dimension(:), allocatable :: myShape                                                          !<shape of the dataset
- integer :: hdferr
- integer(HID_T) :: dset_id, space_id
 
- myShape = shape(dataset)
+ integer(pInt), dimension(:), allocatable :: &
+   globalShape, &                                                                                   !< shape of the dataset (all processes)
+   localShape, &                                                                                    !< shape of the dataset (this process)
+   outputSize                                                                                       !< contribution of all processes
+ integer :: hdferr,ierr
+ integer(HID_T)   :: dset_id, filespace_id, memspace_id, plist_id
+ integer(HSIZE_T), dimension(7) :: myStart
+
+ call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, ierr)
+
 !--------------------------------------------------------------------------------------------------
-! create dataspace
-call h5screate_simple_f(size(myShape), int(myShape,HSIZE_T), space_id, hdferr, &
-                           int(myShape,HSIZE_T))
- if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal7: h5screate_simple_f')
+! determine shape of dataset
+ localShape = shape(dataset)
+ allocate(outputSize(worldsize), source = 0_pInt)
+ outputSize(worldrank+1) = localShape(7)
+#ifdef PETSc
+ if (present(parallel)) then; if (parallel) then
+   call h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, hdferr)
+   if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal7: h5pset_dxpl_mpio_f')
+   call MPI_allreduce(MPI_IN_PLACE,outputSize,worldsize,MPI_LONG,MPI_SUM,PETSC_COMM_WORLD,ierr)       ! get total output size over each process
+   if (ierr /= 0) call IO_error(894_pInt,ext_msg='HDF5_write_pReal7: MPI_allreduce')
+ endif; endif
+#endif
+ myStart     = int([0,0,0,0,0,0,sum(outputSize(1:worldrank))],HSIZE_T)
+ globalShape = [localShape(1:6),sum(outputSize)]
+
+
+!--------------------------------------------------------------------------------------------------
+! create dataspace in memory (local shape)
+ call h5screate_simple_f(size(localShape), int(localShape,HSIZE_T), memspace_id, hdferr, &
+                         int(localShape,HSIZE_T))
+ if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal7: h5screate_simple_f/memspace_id')
+
+!--------------------------------------------------------------------------------------------------
+! create dataspace in file (global shape)
+ call h5screate_simple_f(size(globalShape), int(globalShape,HSIZE_T), filespace_id, hdferr, &
+                         int(globalShape,HSIZE_T))
+ if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal7: h5screate_simple_f/filespace_id')
+
 !--------------------------------------------------------------------------------------------------
 ! create dataset
-call h5dcreate_f(loc_id, trim(datasetName), H5T_NATIVE_DOUBLE, space_id, dset_id, hdferr)
-if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal7: h5dcreate_f')
+ call h5dcreate_f(loc_id, trim(datasetName), H5T_NATIVE_DOUBLE, filespace_id, dset_id, hdferr)
+ if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal7: h5dcreate_f')
 
-CALL h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE,dataset,int(myShape,HSIZE_T), hdferr)
-if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal7: h5dcreate_f')
+!--------------------------------------------------------------------------------------------------
+! select a hyperslab (the portion of the current process) in the file
+ call h5dget_space_f(dset_id, filespace_id, hdferr)
+ if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal7: h5dget_space_f')
+ call h5sselect_hyperslab_f(filespace_id, H5S_SELECT_SET_F, myStart, int(localShape,HSIZE_T), hdferr)
+ if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal7: h5sselect_hyperslab_f')
+
+!--------------------------------------------------------------------------------------------------
+! write
+ call h5dwrite_f(dset_id, H5T_NATIVE_DOUBLE,dataset,int(globalShape,HSIZE_T), hdferr,&
+                file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id)
+ if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal7: h5dcreate_f')
 
 !--------------------------------------------------------------------------------------------------
 !close types, dataspaces
+ call h5pclose_f(plist_id, hdferr)
+ if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal7: plist_id')
  call h5dclose_f(dset_id, hdferr)
-if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_HDF5_write_pReal7: h5dclose_f')
-call h5sclose_f(space_id, hdferr)
-if (hdferr < 0) call IO_error(1_pInt,ext_msg='HDF5_HDF5_write_pReal7: h5sclose_f')
+ if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal7: h5dclose_f')
+ call h5sclose_f(filespace_id, hdferr)
+ if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal7: h5sclose_f/filespace_id')
+ call h5sclose_f(memspace_id, hdferr)
+ if (hdferr /= 0) call IO_error(1_pInt,ext_msg='HDF5_write_pReal7: h5sclose_f/memspace_id')
 
 end subroutine HDF5_write_pReal7
 
+
 !--------------------------------------------------------------------------------------------------
 !> @brief subroutine for writing dataset of the type pInt with 1 dimensions
 !--------------------------------------------------------------------------------------------------