test_get_dispatch_snap.f90

program test_get_dispatch_snap

  use, intrinsic:: ISO_FORTRAN_ENV

  ! Libraries:
  use ezmpi, only: ezmpi_bcast
  use jumble, only: get_command_arg_dyn, read_column, new_unit
  use mpi_f08, only: mpi_init, mpi_finalize, MPI_Comm_rank, MPI_Comm_world, &
       MPI_Comm_size, mpi_abort
  use nr_util, only: assert, deg_to_rad
  use shapelib_03, only: dbf_read_attribute_03

  use derived_types, only: snapshot, shpc
  use dispatch_snapshot_m, only: dispatch_snapshot
  use get_snapshot_m, only: get_snapshot
  use send_snapshot_m, only: send_snapshot
  use shpc_close_m, only: shpc_close
  use shpc_open_m, only: shpc_open
  use write_snapshot_m, only: write_snapshot

  implicit none

  character(len = :), allocatable:: shpc_dir
  type(snapshot) s
  TYPE(shpc) hshp
  integer k_begin, d_init, copy, rank, n_proc, k_end, n_dates
  integer unit_isolated, unit_number_eddies

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

  real:: step(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.

  integer:: k = 0

  integer, allocatable:: ishape_last(:)
  ! shape index (0-based) in the collection of shapefiles of the last
  ! shape at a given date index

  namelist /main_nml/ corner, step, nlon, nlat, dist_lim, k

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

  call mpi_init
  call MPI_Comm_rank(MPI_Comm_world, rank)
  call MPI_Comm_size(MPI_Comm_world, n_proc)

  if (n_proc > 2) then
     if (rank == 0) print *, "test_get_dispatch_snap: 1 or 2 processes only"
     call mpi_abort(MPI_Comm_world, errorcode = 1)
  end if

  call get_command_arg_dyn(1, shpc_dir, "Required argument: SHPC-directory")

  if (rank == 0) then
     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)

     ! As we are requiring the grid spacing to be uniform, the value of
     ! "periodic" may be deduced from the values of step(1) and nlon:
     periodic = nint(360. / step(1)) == nlon
     print *, "periodic = ", periodic
     if (periodic) call assert(2 * dist_lim * step(1) < 180., &
          "test_get_dispatch_snap dist_lim")
     copy = merge(dist_lim, 0, periodic)
     call read_column(ishape_last, &
          file = trim(shpc_dir) // "/ishape_last.txt")
     n_dates = size(ishape_last)
     call new_unit(unit_isolated)
     open(unit_isolated, file = "isolated_nodes.txt", status = "replace", &
          action = "write")
     call new_unit(unit_number_eddies)
     open(unit_number_eddies, file = "number_eddies.csv", status = "replace", &
          action = "write")
  end if

  call ezmpi_bcast(corner, root = 0)
  call ezmpi_bcast(step, root = 0)
  call ezmpi_bcast(nlon, root = 0)
  call ezmpi_bcast(nlat, root = 0)
  call ezmpi_bcast(k, root = 0)
  call ezmpi_bcast(copy, root = 0)
  call ezmpi_bcast(n_dates, root = 0)
  if (rank /= 0) allocate(ishape_last(n_dates))
  call ezmpi_bcast(ishape_last, root = 0)
  call shpc_open(hshp, trim(shpc_dir), rank)
  if (rank == 0) call dbf_read_attribute_03(d_init, hshp%extremum, &
       hshp%extr_date, ishape = 0)
  call ezmpi_bcast(d_init, root = 0)
  k_begin = (rank * n_dates) / n_proc
  
  if (rank < n_proc - 1) then
     k_end = ((rank + 1) * n_dates) / n_proc
  else
     k_end = n_dates - 1
  end if
  
  call get_snapshot(s, nlon, nlat, ishape_last, corner * deg_to_rad, &
       step * deg_to_rad, copy, hshp, d_init, k, k_end, rank, n_proc, &
       max_delta = 1)
  CALL shpc_close(hshp)
  call dispatch_snapshot(s, unit_isolated, unit_number_eddies, rank, k_begin, &
       max_delta = 1, d_init = d_init, k = k)

  if (rank == 0) then
     call write_snapshot(s, corner, step, nlon, nlat, copy, d = d_init + k)
     close(unit_isolated)
     close(unit_number_eddies)
  end if

  call mpi_finalize

end program test_get_dispatch_snap