pfb_read.f90

 subroutine pfb_read(value,fname,nx,ny,nz)
  implicit none
  real*8 value(nx,ny,nz)
  real*8 dx, dy, dz, x1, y1, z1
  integer*4 i,j,k, nni, nnj, nnk, ix, iy, iz,&
            ns,  rx, ry, rz,nx,ny,nz, nnx, nny, nnz,is
  character*200 fname

!print*, nx, ny, nz, fname

open(100,file=trim(adjustl(fname)),form='unformatted',   &
access='stream',convert='BIG_ENDIAN')         !binary outputfile of Parflow
! Read in header info

!print*, fname

! Start: reading of domain spatial information
  read(100) x1 !X
  read(100) y1 !Y
  read(100) z1 !Z

!print*, x1, y1, z1

  read(100) nx !NX
  read(100) ny !NY
  read(100) nz !NZ

!print*, nx, ny, nz

  read(100) dx !DX
  read(100) dy !DY
  read(100) dz !DZ

!print*, dx, dy, dz

  read(100) ns !num_subgrids
! End: reading of domain spatial information

!print*, ns

! Start: loop over number of sub grids
  do is = 0, (ns-1)

! Start: reading of sub-grid spatial information
   read(100) ix
   read(100) iy
   read(100) iz

   read(100) nnx
   read(100) nny
   read(100) nnz
   read(100) rx
   read(100) ry
   read(100) rz

!print*, ix, iy, iz
!print*, nnx, nny, nnz
!print*, rx, ry, rz



! End: reading of sub-grid spatial information

! Start: read in saturation data from each individual subgrid
!  do  k=iz +1 , iz + nnz
!   do  j=iy +1 , iy + nny
!    do  i=ix +1 , ix + nnx
read(100) value((ix+1):(ix+nnx),(iy+1):(iy+nny),(iz+1):(iz+nnz))
!     read(100) value(i,j,k)
!    end do
!   end do
!  end do
! End: read in saturation data from each individual subgrid

! End: read in saturation data from each individual subgrid

  end do
! End: loop over number of sub grids



  close(100)
  end subroutine