subroutine check(err, message)
          use mpi
          use pnetcdf
          implicit none
          integer err
          character(len=*) message

          ! It is a good idea to check returned value for possible error
          if (err .NE. NF90_NOERR) then
              write(6,*) trim(message), trim(nf90mpi_strerror(err))
              call MPI_Abort(MPI_COMM_WORLD, -1, err)
          end if
      end subroutine check

  program main
  use mpi
  use pnetcdf
  implicit none

  integer i, err, ncid
  integer ntracers, nx_glob, ny_glob, nz_glob
  integer(kind=MPI_OFFSET_KIND) :: len
  integer t_dimid, x_dimid, y_dimid, z_dimid, t_varid, x_varid, y_varid, z_varid
  integer hy_dens_varid,hy_dens_theta_varid,dens_varid,uwnd_varid,vwnd_varid
  integer wwnd_varid,theta_varid,windmag_varid,ke_varid,vort_xy_varid,vort_xz_varid
  integer tracer_varids(3)
  character(LEN=256) tname, tracer_names(3)

  call MPI_Init(err)

  ntracers = 3
  nx_glob = 400
  ny_glob = 400
  nz_glob = 200
  tracer_names(1) = "tracer1"
  tracer_names(2) = "tracer2"
  tracer_names(3) = "tracer3"

  err = nf90mpi_create( MPI_COMM_WORLD , 'testfile.nc' , ior( nf90_clobber , nf90_64bit_offset ) , MPI_INFO_NULL , ncid )
  call check(err, 'In nf90mpi_create: ')

  len=nf90_unlimited
  err = nf90mpi_def_dim( ncid , 't' , len , t_dimid )
  call check(err, 'In nf90mpi_def_dim: ')
  len=nx_glob
  err = nf90mpi_def_dim( ncid , 'x' , len , x_dimid )
  call check(err, 'In nf90mpi_def_dim: ')
  len=ny_glob
  err = nf90mpi_def_dim( ncid , 'y' , len , y_dimid )
  call check(err, 'In nf90mpi_def_dim: ')
  len=nz_glob
  err = nf90mpi_def_dim( ncid , 'z' , len , z_dimid )
  call check(err, 'In nf90mpi_def_dim: ')

  err = nf90mpi_def_var( ncid , 't' , nf90_double ,  (/ t_dimid /) , t_varid )
  call check(err, 'In nf90mpi_def_var: ')
  err = nf90mpi_def_var( ncid , 'x' , nf90_double ,  (/ x_dimid /) , x_varid )
  call check(err, 'In nf90mpi_def_var: ')
  err = nf90mpi_def_var( ncid , 'y' , nf90_double ,  (/ y_dimid /) , y_varid )
  call check(err, 'In nf90mpi_def_var: ')
  err = nf90mpi_def_var( ncid , 'z' , nf90_double ,  (/ z_dimid /) , z_varid )
  call check(err, 'In nf90mpi_def_var: ')

  err = nf90mpi_def_var( ncid , 'hy_dens'       , nf90_double ,  (/ z_dimid           /) ,       hy_dens_varid )
  call check(err, 'In nf90mpi_def_var: ')
  err = nf90mpi_def_var( ncid , 'hy_dens_theta' , nf90_double ,  (/ z_dimid           /) , hy_dens_theta_varid )
  call check(err, 'In nf90mpi_def_var: ')
  err = nf90mpi_def_var( ncid , 'dens'          , nf90_double ,  (/ x_dimid , y_dimid , z_dimid , t_dimid /) , dens_varid )
  call check(err, 'In nf90mpi_def_var: ')
  err = nf90mpi_def_var( ncid , 'uwnd'          , nf90_double ,  (/ x_dimid , y_dimid , z_dimid , t_dimid /) , uwnd_varid )
  call check(err, 'In nf90mpi_def_var: ')
  err = nf90mpi_def_var( ncid , 'vwnd'          , nf90_double ,  (/ x_dimid , y_dimid , z_dimid , t_dimid /) , vwnd_varid )
  call check(err, 'In nf90mpi_def_var: ')
  err = nf90mpi_def_var( ncid , 'wwnd'          , nf90_double ,  (/ x_dimid , y_dimid , z_dimid , t_dimid /) , wwnd_varid )
  call check(err, 'In nf90mpi_def_var: ')
  err = nf90mpi_def_var( ncid , 'theta'         , nf90_double ,  (/ x_dimid , y_dimid , z_dimid , t_dimid /) , theta_varid )
  call check(err, 'In nf90mpi_def_var: ')
  err = nf90mpi_def_var( ncid , 'windmag'       , nf90_double ,  (/ x_dimid , y_dimid , z_dimid , t_dimid /) , windmag_varid )
  call check(err, 'In nf90mpi_def_var: ')
  err = nf90mpi_def_var( ncid , 'ke'            , nf90_double ,  (/ x_dimid , y_dimid , z_dimid , t_dimid /) , ke_varid )
  call check(err, 'In nf90mpi_def_var: ')
  err = nf90mpi_def_var( ncid , 'vort_xy'       , nf90_double ,  (/ x_dimid , y_dimid , z_dimid , t_dimid /) , vort_xy_varid )
  call check(err, 'In nf90mpi_def_var: ')
  err = nf90mpi_def_var( ncid , 'vort_xz'       , nf90_double ,  (/ x_dimid , y_dimid , z_dimid , t_dimid /) , vort_xz_varid )
  call check(err, 'In nf90mpi_def_var: ')
  do i = 1 , ntracers
    tname = trim(tracer_names(i))
    err = nf90mpi_def_var( ncid , tname         , nf90_double ,  (/ x_dimid , y_dimid , z_dimid , t_dimid /) , tracer_varids(i) )
    call check(err, 'In nf90mpi_def_var: ')
  enddo
  err = nf90mpi_enddef( ncid )
  call check(err, 'In nf90mpi_enddef: ')
  err = nf90mpi_close( ncid )
  call check(err, 'In nf90mpi_close: ')
  call MPI_Finalize(err)
  end program main