program test_read_snapshot

  use, intrinsic:: ISO_FORTRAN_ENV

  ! Libraries:
  use jumble, only: get_command_arg_dyn, new_unit
  use mpi_f08, only: mpi_init, mpi_finalize
  use nr_util, only: deg_to_rad, assert
  use shapelib_03, only: dbf_read_attribute_03

  use derived_types, only: snapshot, shp_tr
  use read_snapshot_m, only: read_snapshot
  use shp_tr_close_m, only: shp_tr_close
  use shp_tr_open_m, only: shp_tr_open
  use write_snapshot_m, only: write_snapshot

  implicit none

  character(len = :), allocatable:: shp_tr_dir
  type(snapshot) s
  TYPE(shp_tr) hshp
  integer k1, copy, unit

  real:: corner_deg(2) = [0.125, - 59.875]
  ! longitude and latitude of the corner of the whole grid in input
  ! NetCDF, in degrees

  real:: step_deg(2) = [0.25, 0.25] ! longitude and latitude steps, in degrees
  logical periodic ! grid is periodic in longitude

  integer:: nlon = 120, nlat = 120
  ! size of ssh array in input NetCDF, assuming no repeated point if
  ! the grid is global

  integer:: dist_lim = 12
  ! We look for an overlapping eddy at dist_lim (in grid points) of
  ! the first extremum.

  namelist /main_nml/ dist_lim
  namelist /grid_nml/ corner_deg, step_deg, nlon, nlat

  !-------------------------------------------------------------------------

  call mpi_init
  call get_command_arg_dyn(1, shp_tr_dir, &
       "Required argument: SHP-triplet-directory")

  write(unit = error_unit, nml = main_nml)
  write(unit = error_unit, fmt = *) "Enter namelist main_nml."
  read(unit = *, nml = main_nml)
  write(unit = *, nml = main_nml)

  call new_unit(unit)
  open(unit, file = shp_tr_dir // "/grid_nml.txt", status = "old", &
       action = "read", position  = "rewind")
  read(unit, nml = grid_nml)
  close(unit)

  ! As we are requiring the grid spacing to be uniform, the value of
  ! "periodic" may be deduced from the values of step_deg(1) and nlon:
  periodic = nint(360. / step_deg(1)) == nlon
  print *, "periodic = ", periodic
  if (periodic) call assert(2 * dist_lim * step_deg(1) < 180., &
       "test_read_snapshot dist_lim")
  copy = merge(dist_lim, 0, periodic)

  call shp_tr_open(hshp, trim(shp_tr_dir), rank = 0)
  call dbf_read_attribute_03(k1, hshp%extremum, hshp%extr_date, ishape = 0)
  call read_snapshot(s, k1, hshp, corner_deg * deg_to_rad, &
       step_deg * deg_to_rad, nlon, nlat, copy, k1, ishape_last = [3])
  CALL shp_tr_close(hshp)

  call write_snapshot(s, corner_deg, step_deg, nlon, nlat, copy, k1)
  call mpi_finalize

end program test_read_snapshot