test_nearby_extr.f90

program test_nearby_extr

  ! Libraries:
  use jumble, only: get_command_arg_dyn, new_unit
  use nr_util, only: deg_to_rad, rad_to_deg
  use shapelib_03, only: dbf_read_attribute_03

  use derived_types, only: snapshot, shpc
  use nearby_extr_m, only: nearby_extr
  use read_snapshot_m, only: read_snapshot
  use shpc_close_m, only: shpc_close
  use shpc_open_m, only: shpc_open

  implicit none

  character(len = :), allocatable:: shpc_dir
  type(snapshot) s
  TYPE(shpc) hshp
  integer d_init, unit, ishape_last, l

  real corner_deg(2)
  ! longitude and latitude of the corner of the whole grid in input
  ! NetCDF, in degrees

  real step_deg(2) ! longitude and latitude steps, in degrees

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

  namelist /grid_nml/ corner_deg, step_deg, nlon, nlat

  real, allocatable:: nearby(:, :) ! (2, :) longitude and
  ! latitude, in rad, of extrema near the target extremum   

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

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

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

  open(unit, file = shpc_dir // "/ishape_last.txt", status = "old", &
       action = "read", position = "rewind")
  read(unit, fmt = *) ishape_last ! first date
  close(unit)

  call shpc_open(hshp, trim(shpc_dir), pszaccess = "rb")
  call dbf_read_attribute_03(d_init, hshp%extremum, hshp%extr_date, ishape = 0)
  call read_snapshot(s, hshp, nlon, nlat, d_init, k = d_init, &
       corner = corner_deg * deg_to_rad, step = step_deg * deg_to_rad, &
       copy = 0, ishape_last = [ishape_last])
  CALL shpc_close(hshp)

  nearby = nearby_extr(s%extr_map, s%list, i = 4)
  print *, "nearby:"

  do l = 1, size(nearby, 2)
     print *, nearby(:, l) * rad_to_deg
  end do

end program test_nearby_extr