program output
  !============================================================================
  ! PROGRAM output
  !============================================================================
  ! Author: Marc B.R. Joos
  !
  ! Created/last modified: jan 25, 2013/jan 28, 2013
  !
  ! This file is distributed under GNU/GPL licence, 
  ! see <http://www.gnu.org/licenses/>.
  !============================================================================
  ! This program tests the ParallelNetCDF I/O functions in Fortran 90.
  !
  ! One thread creates and write metadata in a file, which is re-open by all
  ! the threads for a collective I/O.
  ! 
  ! It writes a 20x20x20 array with 8 MPI process, distributed in a 2x2x2 grid.
  ! Each process writes its rank in a 10x10x10 subarray.
  ! 
  ! parallelio.nc {
  ! dimensions:
  !       boxdim = 3 ;
  !       x = 20 ;
  !       y = 20 ;
  !       z = 20 :
  ! variables:
  !       int boxsize(boxdim) ;
  !       double var(x, y, z) ;
  ! data:
  !
  ! var = 
  ! 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1
  ! 0 0 ...             1 1 ...
  ! ...
  ! 4 4 4 4 4 4 4 4 4 4 5 5 5 5 5 5 5 5 5 5 
  ! 4 4 ...             5 5 ...
  ! ...
  ! }
  !============================================================================

  use pnetcdf

  implicit none
  include 'mpif.h'

  integer, parameter :: xdim=10, ydim=10, zdim=10
  integer, parameter :: nx=2, ny=2, nz=2
  real(8),dimension(xdim,ydim,zdim) :: data
  real(8), dimension(xdim) :: x
  real(8), dimension(ydim) :: y
  real(8), dimension(zdim) :: z
  integer, dimension(3) :: boxsize
  integer :: ierr, myrank
  integer :: xpos, ypos, zpos
  character(LEN=13) :: filename

  ! PnetCDF variables
  integer(kind=MPI_OFFSET_KIND) :: nxtot, nytot, nztot
  integer :: nout, ncid, bdid, xdimid, ydimid, zdimid
  integer, dimension(3) :: sdimid
  integer :: bsid, varid
  integer(kind=MPI_OFFSET_KIND), dimension(3) :: dims, start, count

  call MPI_Init(ierr)
  call MPI_Comm_rank(MPI_COMM_WORLD, myrank, ierr)

  ! Position of the thread in a 3D cartesian grid
  xpos = mod(myrank, nx)
  ypos = mod(myrank/(nx*nz), ny)
  zpos = mod(myrank/nx, nz)

  ! Fill data
  data = myrank

  ! Define the global dimension and the local dimension
  nxtot = nx*xdim
  nytot = ny*ydim
  nztot = nz*zdim

  dims = (/ xdim, ydim, zdim /)
 
  ! Create file & write metadata (by one thread only)
  filename = 'parallelio.nc'
  if(myrank.eq.0) then
     nout = nfmpi_create(MPI_COMM_SELF, filename, NF_CLOBBER, MPI_INFO_NULL &
          , ncid)

     ! Define and write some metadata
     boxsize = (/ nx, ny, nz /)
     
     nout = nfmpi_def_dim(ncid, "boxdim", 3_8, bdid)
     nout = nfmpi_def_var(ncid, "boxsize", NF_INT, 1, (/bdid/), bsid)
     nout = nfmpi_enddef(ncid)
     
     ! Write the metadata & close the file
     nout = nfmpi_put_vara_int_all(ncid, bsid, (/1_8/), (/3_8/), boxsize)
     nout = nfmpi_close(ncid)
  end if

  ! Synchronize the threads
  call MPI_Barrier(MPI_COMM_WORLD, ierr)

  ! Reopen the file and allow the definition of new variables
  nout = nfmpi_open(MPI_COMM_WORLD, filename, NF_WRITE, MPI_INFO_NULL, ncid)
  nout = nfmpi_redef(ncid)

  ! Define dimensions
  nout = nfmpi_def_dim(ncid, "x", nxtot, xdimid)
  nout = nfmpi_def_dim(ncid, "y", nytot, ydimid)
  nout = nfmpi_def_dim(ncid, "z", nztot, zdimid)
  sdimid = (/ xdimid, ydimid, zdimid /)

  ! Create variable
  nout = nfmpi_def_var(ncid, "var", NF_DOUBLE, 3, sdimid, varid)

  ! End of definitions
  nout = nfmpi_enddef(ncid)

  ! Indices start at 1!
  start = (/ xpos, ypos, zpos /)*dims+1
  count = dims

  ! Write data
  nout = nfmpi_put_vara_double_all(ncid, varid, start, count, data)

  ! Close file
  nout = nfmpi_close(ncid)

  call MPI_Finalize(ierr)

end program output