PROGRAM posix_io

  IMPLICIT NONE
  INCLUDE 'mpif.h'

  INTEGER :: i, myid, nproc, bufsize, ierr
  INTEGER, PARAMETER :: nx = 32
  INTEGER, ALLOCATABLE :: buf(:,:)
  INTEGER :: fp, intsize

! Variables for cartesian topology
  INTEGER, DIMENSION(2) :: mpi_dims, coords
  LOGICAL, DIMENSION(2) :: period = (/.FALSE., .FALSE./)
  INTEGER :: cart_comm

! Variables for domain decomposition and subarray view 
  INTEGER, DIMENSION(2) :: dimsf = (/nx, nx/)
  INTEGER, DIMENSION(2) :: count, offset

  INTEGER(KIND=MPI_OFFSET_KIND) :: disp
  INTEGER :: filetype


  CALL MPI_INIT(ierr)
  CALL MPI_COMM_RANK(MPI_COMM_WORLD, myid, ierr);
  CALL MPI_COMM_SIZE(MPI_COMM_WORLD, nproc, ierr);

  mpi_dims(1) = INT(SQRT(REAL(nproc)))
  mpi_dims(2) = mpi_dims(1)
  IF (mpi_dims(1) * mpi_dims(2) /= nproc ) STOP 'Wrong number of processes' 
  CALL MPI_CART_CREATE(MPI_COMM_WORLD, 2, mpi_dims, period, 0, cart_comm, ierr)
  CALL MPI_COMM_RANK(cart_comm, myid, ierr)
  CALL MPI_CART_COORDS(cart_comm, myid, 2, coords, ierr)

  DO i=1,2
    count(i) = nx / mpi_dims(i)
    offset(i) = coords(i) * count(i)
  END DO
  bufsize = count(1) * count(2)

  ALLOCATE(buf(count(1), count(2)))

  buf(:,:) = myid

  CALL MPI_FILE_OPEN(MPI_COMM_WORLD, "test.dat", &
       MPI_MODE_CREATE+MPI_MODE_WRONLY, &
       MPI_INFO_NULL, fp, ierr);
  
! Create the file view
  CALL MPI_TYPE_CREATE_SUBARRAY(2, dimsf, count, offset, MPI_ORDER_FORTRAN, &
                                MPI_INTEGER, filetype, ierr)
  CALL MPI_TYPE_COMMIT(filetype, ierr)
  CALL MPI_FILE_SET_VIEW(fp, 0, MPI_INTEGER, filetype, "native", &
                         MPI_INFO_NULL, ierr)

! Write the data
  CALL MPI_FILE_WRITE_ALL(fp, buf, bufsize, &
                          MPI_INTEGER, MPI_STATUS_IGNORE, ierr)

  CALL MPI_FILE_CLOSE(fp, ierr)
  CALL MPI_TYPE_FREE(filetype, ierr)
  CALL MPI_FINALIZE(ierr)

END PROGRAM
