Read_MultiSetModel3D_t Subroutine

  subroutine Read_MultiSetModel3D_t(this,filename)
    !! Read every set's solution from a file written by WriteModel, then copy it to the
    !! device. The file must hold the same element-set partition (setOfElem) and the
    !! same number of nodes per element; a variable with no dataset keeps its current
    !! value. The host copy is refreshed from the device first, so on GPU builds that
    !! value is the device's.
    !!
    !! Collective over the global mesh communicator.
    implicit none
    class(MultiSetModel3D_t),intent(inout) :: this
    character(*),intent(in) :: filename
    ! Local
    integer(HID_T) :: fileId
    character(4) :: setString
    character(LEN=:),allocatable :: dsetName
    integer,allocatable :: fileSetOfElem(:)
    integer :: k,ivar,elemOffset
    logical :: parallel

    parallel = this%globalMesh%decomp%mpiEnabled
    if(parallel) then
      call Open_HDF5(filename,H5F_ACC_RDONLY_F,fileId,this%globalMesh%decomp%mpiComm)
    else
      call Open_HDF5(filename,H5F_ACC_RDONLY_F,fileId)
    endif

    ! Check the stored extent first: ReadArray_HDF5 reads the whole dataset into the
    ! buffer, so a longer dataset would write past it.
    if(DatasetLength1D_HDF5(fileId,'/elementsets/setOfElem') /= size(this%setOfElem)) then
      print*,__FILE__," : ",trim(filename)," holds no element-set partition for ", &
        size(this%setOfElem)," elements"
      error stop "MultiSetModel3D : ReadModel partition mismatch"
    endif
    allocate(fileSetOfElem(1:size(this%setOfElem)))
    call ReadArray_HDF5(fileId,'/elementsets/setOfElem',fileSetOfElem)
    if(any(fileSetOfElem /= this%setOfElem)) then
      print*,__FILE__," : ",trim(filename)," holds a different element-set partition"
      error stop "MultiSetModel3D : ReadModel partition mismatch"
    endif
    deallocate(fileSetOfElem)

    do k = 1,this%nSets
      call this%sets(k)%model%solution%UpdateHost()
      write(setString,'(I4.4)') k
      elemOffset = this%sets(k)%mesh%decomp%offsetElem(this%globalMesh%decomp%rankId+1)
      do ivar = 1,this%sets(k)%model%solution%nVar
        dsetName = '/set_'//setString//'/controlgrid/solution/'// &
                   trim(this%sets(k)%model%solution%meta(ivar)%name)
        if(.not. DatasetExists_HDF5(fileId,dsetName)) then
          if(this%globalMesh%decomp%rankId == 0) then
            print*,__FILE__," : Pickup file holds no ",dsetName, &
              " - keeping the current value for this variable."
          endif
          cycle
        endif
        call ReadSetArray_HDF5(fileId,dsetName, &
                               this%sets(k)%model%solution%interior(:,:,:,:,ivar), &
                               elemOffset,this%sets(k)%mesh%nGlobalElem,parallel)
      enddo
      call this%sets(k)%model%solution%UpdateDevice()
    enddo

    call Close_HDF5(fileId)

  endsubroutine Read_MultiSetModel3D_t