!>
!> @file pex4.f90
!>
!> @brief Create particle arrays using "putarr" (paralle version)
!>
!> @copyright
!> Copyright (©) 2021 EPFL (Ecole Polytechnique Fédérale de Lausanne)
!> SPC (Swiss Plasma Center)
!>
!> futils is free software: you can redistribute it and/or modify it under
!> the terms of the GNU Lesser General Public License as published by the Free
!> Software Foundation, either version 3 of the License, or (at your option)
!> any later version.
!>
!> futils is distributed in the hope that it will be useful, but WITHOUT ANY
!> WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
!> FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
!>
!> You should have received a copy of the GNU Lesser General Public License
!> along with this program. If not, see <https://www.gnu.org/licenses/>.
!>
!> @authors
!> (in alphabetical order)
!> @author Trach-Minh Tran <trach-minh.tran@epfl.ch>
!>
PROGRAM main
!
!   Particle array using fixed dims (parallel version)
!
  USE futils
  USE mpi
  IMPLICIT NONE
  CHARACTER(len=256) :: file='part.h5'
  CHARACTER(len=32) :: name
!!$  INTEGER, PARAMETER :: npart=1024*100, nattrs=8
  INTEGER, PARAMETER :: npart=1024, nattrs=2
  INTEGER :: cart
  INTEGER :: fid, istep, nrun=4
  INTEGER :: me, npes, n, ierr, nptot, nploc(1), s0(1)
  DOUBLE PRECISION :: r0=100.0, a=20.0, time, pi, dphi, s(npart), theta(npart)
  DOUBLE PRECISION, ALLOCATABLE :: part(:,:)   ! (r, z, phi) coordinates
  DOUBLE PRECISION :: t0, t1, twrite, mb, mbs_write
!===========================================================================
!                   1. Prologue
!
  CALL mpi_init(ierr)
  CALL mpi_comm_size(MPI_COMM_WORLD, npes, ierr)
  CALL mpi_comm_rank(MPI_COMM_WORLD, me, ierr)
!
  CALL mpi_cart_create(MPI_COMM_WORLD, 1, [npes], [.FALSE.], .FALSE., &
       &               cart, ierr)
!!$  PRINT*, 'cart created!'
!
  IF( command_argument_count() > 0 ) THEN
     CALL get_command_argument(1, file, n, ierr)
  END IF
!
  CALL creatf(file, fid, 'A simple simulation', &
       &      real_prec='d', mpicomm=cart)
  CALL creatg(fid, "/part", "Particles Coordinates")        ! Group
!===========================================================================
!                   2. Time loop
!
  nptot = npart
  CALL dist1d(1, nptot, s0, nploc(1))
  ALLOCATE(part(nattrs,nploc(1)))
  CALL initp(nattrs, nploc(1), part, s0(1))
  CALL mpi_barrier(MPI_COMM_WORLD, ierr)
  t0 = mpi_wtime()
!
  DO istep=1,nrun
     time = istep
     WRITE(name,'(a,i3.3)') "/part/", istep
     CALL creatg(fid, name)        ! Group
     CALL putarr(fid, TRIM(name)//'/s0',   s0, pardim=1)     ! Local scalar
     CALL putarr(fid, TRIM(name)//'/nploc', nploc, pardim=1) ! Local scalar
     CALL putarr(fid, TRIM(name)//'/part', part, pardim=2)   ! local 2d array
     CALL attach(fid, name, "time", time)   ! Attr on dataset
     CALL attach(fid, name, "step", istep)
  END DO
  nptot = npes*npart
  CALL attach(fid, "/part", "nptot", nptot)
  CALL attach(fid, "/part", "nattrs", nattrs)
  CALL attach(fid, "/part", "nsteps", nrun)
!
  CALL mpi_barrier(MPI_COMM_WORLD, ierr)
  twrite = mpi_wtime()-t0
  mb = 8.0*REAL(SIZE(part))/1024.0/1024*npes*nrun
  mbs_write = mb/twrite
  IF( me.EQ. 0) THEN
     WRITE(*,'(a,5(f8.3,a))') 'Write performance:', &
          &   twrite, ' s,',mb, ' MB,', mbs_write, ' MB/s'
  END IF
!===========================================================================
!                   9. Epilogue
!
  CALL mpi_barrier(MPI_COMM_WORLD, ierr)
  CALL closef(fid)
  CALL mpi_finalize(ierr)
END PROGRAM main

SUBROUTINE initp(natts, n, p, s)
  IMPLICIT NONE
  INTEGER :: natts, n,s
  DOUBLE PRECISION :: x, p(natts,n)
  INTEGER :: i, j
  x=s-1.0d0
  DO j=1,n
     x=x+1.0d0
     DO i=1,natts
        p(i,j) = 1.d0/x
     END DO
  END DO
END SUBROUTINE initp

SUBROUTINE dist1d(s0, ntot, s, nloc)
  USE mpi
  IMPLICIT NONE
  INTEGER, INTENT(in) :: s0, ntot
  INTEGER, INTENT(out) :: s, nloc
  INTEGER :: me, npes, ierr, naver, rem
!
  CALL MPI_COMM_SIZE(MPI_COMM_WORLD, npes, ierr)
  CALL MPI_COMM_RANK(MPI_COMM_WORLD, me, ierr)
  PRINT*, 'me, ntot, npes ', me, ntot, npes
  naver = ntot/npes
  rem = MODULO(ntot,npes)
  s = s0 + MIN(rem,me) + me*naver
  nloc = naver
  IF( me.LT.rem ) nloc = nloc+1
!
END SUBROUTINE dist1d