!>
!> @file ex13.f90
!>
!> @brief Read slices (in time) of 3d array
!>
!> @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
!
!   Create 3D (+time as 3rd extendible dim) datasets with "append"
!   Read slices (in time) of 3d array
!
  USE futils
  IMPLICIT NONE
  CHARACTER(len=32)  :: file='prof_slice.h5'
  INTEGER, PARAMETER :: nx=256, ny=256, nrun=20, slice=2
  INTEGER            :: fid, rank, dims(7), offsets(3)
  INTEGER            :: i, j, istep, is, nerrs
  DOUBLE PRECISION   :: xg(nx), yg(ny), phi3(nx,ny,nrun), phi_read(nx,ny,slice)
  DOUBLE PRECISION   :: dx, dy, pi, time, x0, s
!===========================================================================
!
!  Create file using DB reals
  CALL creatf(file, fid, real_prec='d')
  WRITE(*,'(a)') file(1:LEN_TRIM(file))//' open'
!
!  Create group
  CALL creatg(fid, "/profile_2d", "2D profiles")
  WRITE(*,'(a)') "group profile_12d created"
!
!  Create extendible dataset
  CALL creatd(fid, 2, (/nx, ny/), "/profile_2d/phi")
!
!  Mesh
  pi = 4*ATAN(1.0d0)
  dx = 2.*pi/REAL(nx-1)
  dy = pi/real(ny-1)
  DO i=1,nx
     xg(i) = (i-1)*dx
  END DO
  DO j=1,ny
     yg(j) = (j-1)*dy
  END DO
  CALL putarr(fid, "/profile_2d/xg", xg, "x mesh")
  CALL putarr(fid, "/profile_2d/yg", yg, "y mesh")
!
!  Time loop
  x0 = 0.5
  DO istep=1,nrun
     time = istep-1
     DO i=1,nx
        DO j=1,ny
           phi3(i,j,istep) = SIN(xg(i)) * COS(yg(j)) * COS(0.04*pi*time)
        END DO
     END DO
     CALL append(fid, "/profile_2d/phi", phi3(:,:,istep))
  END DO
  CALL closef(fid)
!
!   Reopen for read
  CALL openf(file, fid)
  WRITE(*,'(a)') file(1:LEN_TRIM(file))//' reopen'
!
!   Read 2d profiles and check
  offsets = 0
  DO istep=1,nrun,slice
     offsets(3) = istep-1
     CALL getarr(fid, "/profile_2d/phi", phi_read, offsets=offsets)
     nerrs=0
     DO i=1,nx
        DO j=1,ny
           DO is=1,slice
              IF(phi_read(i,j,is) .NE. phi3(i,j,is+istep-1) ) nerrs=nerrs+1
           END DO
        END DO
     END DO
     WRITE(*,'(a,i3,i6)') 'istep, nerrs', istep, nerrs
  END DO
!
!  Clean up
  CALL closef(fid)
END PROGRAM main