marker_set_ReadH5HUT Subroutine

    impure subroutine marker_set_ReadH5HUT(this,iter,time,step)
      !> Reads marker data from file in parallel using H5HUT.
      implicit none
      class(marker_set), intent(inout) :: this                                 !! A collection of tessellation elements
      integer,           intent(out)   :: iter                                 !! Iteration at write
      real(wp),          intent(out)   :: time                                 !! Time at write
      integer,           intent(in),    &
                              optional :: step                                 !! User supplied step to open
      ! Work variables
      type(h5hut_obj)               :: h5                                      !! H5hut structure
      integer,          allocatable :: buffi(:)                                !! Integer buffer
      integer(leapI8),  allocatable :: buffi8(:)                               !! Integer (kind=8) buffer
      real(wp),         allocatable :: buffr(:)                                !! Real buffer
      integer,          allocatable :: size_per_rank(:)                        !! Number of markers per rank
      integer :: max_chunk                                                     !! Max number of chunks on any MPI rank
      integer :: chunk_size                                                    !! Number of markers in chunk
      integer :: goff                                                          !! Global offset
      integer :: loff                                                          !! Local offset
      integer :: n                                                             !! Iterator

      ! Open the file
      call h5%Initialize(trim(adjustl(this%read_file)),"R",this%parallel)

      if (.not.present(step)) then
        ! Jump to last step
        call h5%LastTimeStep(iter,time)
      else
        ! Jump to desired step
        call h5%JumpToStep(step,iter,time)
      end if

      ! Get total number of point objects
      call h5%GetNPoints(this%count)

      ! Distribute equally the data on each mpi rank
      associate(mpi=>this%parallel)
        allocate(size_per_rank(mpi%nproc))

        do n=1,mpi%nproc
          size_per_rank(n) = int(this%count/mpi%nproc)
          if (n.le.mod(this%count,mpi%nproc)) size_per_rank(n)=size_per_rank(n)+1
        end do
      end associate

      ! Allocate buffers
      call this%parallel%Max(ceiling(size_per_rank(this%parallel%rank%mine)/real(this%R_chunk_size,wp)),max_chunk)

      loff  = 0
      chunk_size = 0
      do n=1,max_chunk
        ! Local offset
        loff    = loff + chunk_size

        ! Determine next chunk size
        if (loff.lt.size_per_rank(this%parallel%rank%mine)) then
          chunk_size = min(size_per_rank(this%parallel%rank%mine) - loff,this%R_chunk_size)

          ! Resize the data array to contain read Lagrangian objects
          call this%Resize(this%count_+chunk_size)
        else
          chunk_size = 0
        end if

        ! Global offset
        goff = sum(size_per_rank(1:this%parallel%rank%mine-1)) + loff

        allocate(buffi8(chunk_size),buffi(chunk_size),buffr(chunk_size))

        ! Read the base Lagrangian object structure
        if (chunk_size.eq.0) then
          call h5%Read("id",buffi8,offset=goff)
          call h5%Read("x", buffr, offset=goff)
          call h5%Read("y", buffr, offset=goff)
          call h5%Read("z", buffr, offset=goff)
          call h5%Read("i", buffi, offset=goff)
          call h5%Read("j", buffi, offset=goff)
          call h5%Read("k", buffi, offset=goff)
          call h5%Read("s", buffi, offset=goff)
          call h5%Read("SA",buffr, offset=goff)
          call h5%Read("nx",buffr, offset=goff)
          call h5%Read("ny",buffr, offset=goff)
          call h5%Read("nz",buffr, offset=goff)
          call h5%Read("vx",buffr, offset=goff)
          call h5%Read("vy",buffr, offset=goff)
          call h5%Read("vz",buffr, offset=goff)
          call h5%Read("fx",buffr, offset=goff)
          call h5%Read("fy",buffr, offset=goff)
          call h5%Read("fz",buffr, offset=goff)
        else
          select type (markers=>this%p)
          type is (marker_obj)
            call h5%Read("id",buffi8,offset=goff); markers(this%count_-chunk_size+1:this%count_)%id   = buffi8
            call h5%Read("x", buffr, offset=goff); markers(this%count_-chunk_size+1:this%count_)%p(1) = buffr
            call h5%Read("y", buffr, offset=goff); markers(this%count_-chunk_size+1:this%count_)%p(2) = buffr
            call h5%Read("z", buffr, offset=goff); markers(this%count_-chunk_size+1:this%count_)%p(3) = buffr
            call h5%Read("i", buffi, offset=goff); markers(this%count_-chunk_size+1:this%count_)%c(1) = buffi
            call h5%Read("j", buffi, offset=goff); markers(this%count_-chunk_size+1:this%count_)%c(2) = buffi
            call h5%Read("k", buffi, offset=goff); markers(this%count_-chunk_size+1:this%count_)%c(3) = buffi
            call h5%Read("s", buffi, offset=goff); markers(this%count_-chunk_size+1:this%count_)%s    = buffi
            call h5%Read("SA",buffr, offset=goff); markers(this%count_-chunk_size+1:this%count_)%SA   = buffr
            call h5%Read("nx",buffr, offset=goff); markers(this%count_-chunk_size+1:this%count_)%n(1) = buffr
            call h5%Read("ny",buffr, offset=goff); markers(this%count_-chunk_size+1:this%count_)%n(2) = buffr
            call h5%Read("nz",buffr, offset=goff); markers(this%count_-chunk_size+1:this%count_)%n(3) = buffr
            call h5%Read("vx",buffr, offset=goff); markers(this%count_-chunk_size+1:this%count_)%v(1) = buffr
            call h5%Read("vy",buffr, offset=goff); markers(this%count_-chunk_size+1:this%count_)%v(2) = buffr
            call h5%Read("vz",buffr, offset=goff); markers(this%count_-chunk_size+1:this%count_)%v(3) = buffr
            call h5%Read("fx",buffr, offset=goff); markers(this%count_-chunk_size+1:this%count_)%f(1) = buffr
            call h5%Read("fy",buffr, offset=goff); markers(this%count_-chunk_size+1:this%count_)%f(2) = buffr
            call h5%Read("fz",buffr, offset=goff); markers(this%count_-chunk_size+1:this%count_)%f(3) = buffr
          end select
        end if

        ! Deallocate buffers
        deallocate(buffi8,buffi,buffr)

        ! Send markers to their mpiranks
        call this%Communicate()
      end do

      ! Free data
      deallocate(size_per_rank)
      call h5%Finalize

      ! Set unread variables to zero
      select type (markers=>this%p)
      class is (marker_obj)
        do n=1,this%count_
          markers(n)%pold  = 0.0_wp
          markers(n)%vold  = 0.0_wp
        end do
      end select

      ! Localize markers on the grid
      call this%Localize()

      return
    end subroutine marker_set_ReadH5HUT