program main
    use pnetcdf
    implicit none
    include 'mpif.h'

    integer status, ncid, varid, dimid(4)
    integer nprocs, rank, xrank, yrank, zrank, w, x, y, z
    double precision, dimension(:,:,:,:), allocatable :: val
    integer(kind=MPI_OFFSET_KIND) :: len, nv, count(4), offset(4)

    call MPI_Init(status)
    call MPI_Comm_rank(MPI_COMM_WORLD, rank, status)
    call MPI_Comm_size(MPI_COMM_WORLD, nprocs, status)
    if (nprocs .NE. 8) then
        print*,'Must run on 8 processes! ... exiting'
        call MPI_Finalize(status)
        stop
    endif

    ! create a new file
    status = nfmpi_create(MPI_COMM_WORLD, 'restart.nc', &
                          NF_CLOBBER, MPI_INFO_NULL, ncid)

    len = 10
    nv  = 2
    status = nfmpi_def_dim(ncid, "x", len, dimid(1))
    status = nfmpi_def_dim(ncid, "y", len, dimid(2))
    status = nfmpi_def_dim(ncid, "z", len, dimid(3))
    status = nfmpi_def_dim(ncid, "w", nv, dimid(4))
    status = nfmpi_def_var(ncid, "var", NF_DOUBLE, 4, dimid, varid)
    status = nfmpi_enddef(ncid)

    zrank = rank / 4
    yrank = mod(rank, 4) / 2
    xrank = mod(rank, 2)

    offset(1) = 5 * xrank + 1
    offset(2) = 5 * yrank + 1
    offset(3) = 5 * zrank + 1
    offset(4) = 1

    count(1:3) = 5
    count(4)   = nv

    123 format(A,i3,i3,i3,i3,A,i3,i3,i3,i3)
    print 123,'offset=',offset(:),' count=',count(:)

    ! allocate and initialize val buffer
    allocate( val(count(1), count(2), count(3), nv) )
    do w=1,nv
       do z=1,len/2
          do y=1,len/2
             do x=1,len/2
                val(x,y,z,w) = 1.0 * rank + 10.0 * w + z + 0.001 * y + 0.1 * x
             enddo
          enddo
       enddo
    enddo

    status = nfmpi_put_vara_double_all(ncid, varid, offset, count, val)

    status = nfmpi_close(ncid)

    !---- open for read
    status = nfmpi_open(MPI_COMM_WORLD, 'restart.nc', &
                      NF_NOWRITE, MPI_INFO_NULL, ncid)
    val(:,:,:,:) = 0.0

    status = nfmpi_inq_varid(ncid, "var", varid)
    status = nfmpi_get_vara_double_all(ncid, varid, offset, count, val)

    ! check the contents of read values
    456 format(i2,A,i2,i2,i2,i2,A,i3)
    do w=1,nv
       do z=1,len/2
          do y=1,len/2
             do x=1,len/2
                if (val(x,y,z,w) .NE. (1.0*rank+10.0*w+z+0.001*y+0.1*x)) then
                    print 456,rank,': read error at ',x,y,z,w,' read val=',val(x,y,z,w)
                    goto 999
                endif
             enddo
          enddo
       enddo
    enddo
    999 status = nfmpi_close(ncid)

    call MPI_Finalize(status)
end program