program testmpiio
  use mpi
!  use pnetcdf
  implicit none
#include 'pnetcdf.inc'

  integer :: ierr, ierr2, mype, npe

  integer(kind=mpi_offset_kind) :: start(3),kount(3), gdims(3)
  integer :: fh, dimid(4), vid, elemtype,request,status(MPI_STATUS_SIZE)
  double precision, allocatable :: bigfoo(:,:,:)
!  integer, allocatable :: bigfoo(:,:,:)


  !
  ! Initialize MPI
  !
  call MPI_Init(ierr)
  call CheckMPIreturn(__LINE__,ierr)
  call MPI_Comm_rank(MPI_COMM_WORLD, mype, ierr)
  call CheckMPIreturn(__LINE__,ierr)
  call MPI_Comm_size(MPI_COMM_WORLD, npe,  ierr)
  call CheckMPIreturn(__LINE__,ierr)

  if(npe/=8) then
     print *,' Test program for 8 tasks only, aborting'
     call mpi_abort(mpi_comm_world,ierr,ierr2)
  endif

  ierr = nfmpi_create(MPI_COMM_WORLD,'pnetcdf5.nc', &
       ior(NF_64BIT_DATA,NF_CLOBBER),MPI_INFO_NULL,fh)

  if(mype==0) then
     start = (/1,1,1/)
     kount = (/2,41,3932160/)
  else if(mype==1) then
     start = (/1,1,3932161/)
     kount = (/2,41,3932160/)
  else if(mype==2) then
     start = (/1,1,7864321/)
     kount = (/2,41,3932161/)
  else if(mype==3) then
     start = (/1,1,11796482/)
     kount = (/2,41,3932161/)
  else if(mype==4) then
     start = (/1,1,15728643/)
     kount = (/2,41,3932161/)
  else if(mype==5) then
     start = (/1,1,19660804/)
     kount = (/2,41,3932161/)
  else if(mype==6) then
     start = (/1,1,23592965/)
     kount = (/2,41,3932161/)
  else if(mype==7) then
     start = (/1,1,27525126/)
     kount = (/2,41,3932161/)
  endif

  call mpi_type_contiguous(int(product(kount)),MPI_DOUBLE,elemtype,ierr)
  call CheckMPIreturn(__LINE__,ierr)
  call mpi_type_commit(elemtype,ierr)
  call CheckMPIreturn(__LINE__,ierr)


  gdims = (/2,41,31457286/)

  ierr = nfmpi_def_dim(fh,'dim1',gdims(1),dimid(1))
  if(ierr /= nf_noerr) then
     print *,__LINE__,ierr
     call mpi_abort(mpi_comm_world,ierr,ierr2)
  endif
  ierr = nfmpi_def_dim(fh,'dim2',gdims(2),dimid(2))
  if(ierr /= nf_noerr) then
     print *,__LINE__,ierr
     call mpi_abort(mpi_comm_world,ierr,ierr2)
  endif
  ierr = nfmpi_def_dim(fh,'dim3',gdims(3),dimid(3))
  if(ierr /= nf_noerr) then
     print *,__LINE__,ierr
     call mpi_abort(mpi_comm_world,ierr,ierr2)
  endif
  
  ierr = nfmpi_def_var(fh,'var1',nf_double,3,dimid,vid)
  if(ierr /= nf_noerr) then
     print *,__LINE__,ierr
     call mpi_abort(mpi_comm_world,ierr,ierr2)
  endif

  ierr  = nfmpi_enddef(fh)
  if(ierr /= nf_noerr) then
     print *,__LINE__,ierr
     call mpi_abort(mpi_comm_world,ierr,ierr2)
  endif

  allocate(bigfoo(kount(1),kount(2),kount(3)))
  bigfoo = mype

  ierr = nfmpi_iput_vara(fh, vid, start,kount,bigfoo,1,elemtype,request)
  if(ierr /= nf_noerr) then
     print *,__LINE__,ierr
     call mpi_abort(mpi_comm_world,ierr,ierr2)
  endif

  ierr = nfmpi_wait_all(fh,1,(/request/),status)
  if(ierr /= nf_noerr) then
     print *,__LINE__,ierr
     call mpi_abort(mpi_comm_world,ierr,ierr2)
  endif

  ierr = nfmpi_close(fh)
  if(ierr /= nf_noerr) then
     print *,__LINE__,ierr
     call mpi_abort(mpi_comm_world,ierr,ierr2)
  endif

  deallocate(bigfoo)

  call MPI_Finalize(ierr)
  call CheckMPIreturn(__LINE__,ierr)
   
  stop
end program testmpiio
!=============================================
!  CheckMPIreturn:
!
!      Check and prints an error message
!  if an error occured in a MPI subroutine.
!=============================================
subroutine CheckMPIreturn(locmesg, errcode, file, line)
  use MPI
  implicit none
  character(len=*), intent(in) :: locmesg
  integer, intent(in) :: errcode
  character(len=*),optional :: file
  integer, intent(in),optional :: line
  character(len=MPI_MAX_ERROR_STRING) :: errorstring

  integer :: errorlen

  integer :: ierr

  if (errcode .ne. MPI_SUCCESS) then
     call MPI_Error_String(errcode,errorstring,errorlen,ierr)
     write(*,*) TRIM(ADJUSTL(locmesg))//errorstring(1:errorlen)
     call mpi_abort(mpi_comm_world,errcode,ierr)
  end if
  end subroutine CheckMPIreturn