Use `read_snapshot` in program `test_nearby_extr`

Instead of reading `extr_map.nc`, `ishape_last.txt` and shapefile
extremum. The problem was the consistency between `extr_map.nc` and
the two other files. `extr_map.nc` contains eddies with the two
orientations while extremum now contains a single orientation, with
eddies renumbered.

Inst_eddies/Tests/CMakeLists.txt

@@ -79,10 +79,16 @@ target_include_directories(test_max_speed_contour_ssh PRIVATE
 
 add_executable(test_nearby_extr nearby_extr.f90
   ${CMAKE_SOURCE_DIR}/Common/derived_types.f90
-  ${CMAKE_CURRENT_LIST_DIR}/test_nearby_extr.f90)
+  ${CMAKE_CURRENT_LIST_DIR}/test_nearby_extr.f90
+  ${CMAKE_SOURCE_DIR}/Common/read_snapshot.f90
+  ${CMAKE_SOURCE_DIR}/Common/read_eddy.f90
+  ${CMAKE_SOURCE_DIR}/Common/read_field_indices.f90
+  ${CMAKE_SOURCE_DIR}/Common/shp_tr_open.f90
+  ${CMAKE_SOURCE_DIR}/Common/shp_tr_close.f90)
 
 target_link_libraries(test_nearby_extr NetCDF95::netcdf95 shapelib_03
-  contour_531 jumble nr_util NetCDF_Fortran::NetCDF_Fortran)
+  contour_531 jumble nr_util NetCDF_Fortran::NetCDF_Fortran ezmpi
+  MPI::MPI_Fortran gpc_f)
 
 target_include_directories(test_nearby_extr PRIVATE
   ${fortrangis_INCLUDE_DIR})

Inst_eddies/Tests/short_tests.json

@@ -325,13 +325,12 @@
         ]
     },
     {
-        "required": [
-            "$tests_old_dir/Inst_eddies_loop/SHPC_cyclo_all_dates/extremum.*",
-            "$tests_old_dir/Inst_eddies_loop/SHPC_cyclo_all_dates/ishape_last.txt",
-            "$PWD/Local_extrema/extr_map.nc"
-        ],
         "title": "Nearby_extr",
-        "command": "$build_dir/Inst_eddies/test_nearby_extr"
+        "command":
+	[
+	  "$build_dir/Inst_eddies/test_nearby_extr",
+	  "$tests_old_dir/Inst_eddies_loop/SHPC_cyclo_all_dates"
+	]
     },
     {
         "description": "test_set_all_outerm with periodicity.",

Inst_eddies/Tests/test_nearby_extr.f90

 program test_nearby_extr
 
   ! Libraries:
-  use jumble, only: new_unit
-  use netcdf, only: nf90_nowrite
-  use netcdf95, only: nf95_open, nf95_close, nf95_inq_varid, nf95_gw_var
-  use nr_util, only: deg_to_rad, rad_to_deg
-  use shapelib, only: shpfileobject, shpclose, shpobject, dbfreadattribute
-  use shapelib_03, only: shp_open_03, shp_read_object_03
-
-  use derived_types, only: eddy
+  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, rad_to_deg
+  use shapelib_03, only: dbf_read_attribute_03
+
+  use derived_types, only: snapshot, shp_tr
   use nearby_extr_m, only: nearby_extr
+  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
 
   implicit none
 
-  integer ncid, varid, ishape_last, ishape, i, l, unit
-  integer int_attr ! integer attribute
+  character(len = :), allocatable:: shpc_dir
+  type(snapshot) s
+  TYPE(shp_tr) hshp
+  integer k1, 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, allocatable:: extr_map(:, :)
-  ! At a point of extremum SSH: identification number or this
-  ! extremum. 0 at other points.
+  integer nlon, nlat
+  ! size of ssh array in input NetCDF, assuming no repeated point if
+  ! the grid is global
 
-  type(eddy), allocatable:: list_vis(:) ! visible eddies at a given date
-  TYPE(shpfileobject) hshp
-  TYPE(shpobject) psobject
+  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 nf95_open("extr_map.nc", nf90_nowrite, ncid)
-  call nf95_inq_varid(ncid, "extr_map", varid)
-  call nf95_gw_var(ncid, varid, extr_map)
-  call nf95_close(ncid)
+  call mpi_init
+  call get_command_arg_dyn(1, shpc_dir, "Required argument: SHPC-directory")
 
   call new_unit(unit)
-  open(unit, file = "ishape_last.txt", status = "old", action = "read", &
-       position = "rewind")
-  read(unit, fmt = *) ishape_last ! first date
+  open(unit, file = shpc_dir // "/grid_nml.txt", status = "old", &
+       action = "read", position  = "rewind")
+  read(unit, nml = grid_nml)
   close(unit)
 
-  call shp_open_03(hshp, "extremum", pszaccess = "rb")
-  allocate(list_vis(ishape_last + 1))
-
-  do ishape = 0, ishape_last
-     call dbfreadattribute(hshp, ishape, ifield = 2, attr = i)
-
-     call dbfreadattribute(hshp, ishape, ifield = 4, attr = int_attr)
-     list_vis(i)%cyclone = int_attr == 1
+  open(unit, file = shpc_dir // "/ishape_last.txt", status = "old", &
+       action = "read", position = "rewind")
+  read(unit, fmt = *) ishape_last ! first date
+  close(unit)
 
-     call dbfreadattribute(hshp, ishape, ifield = 5, attr = int_attr)
-     list_vis(i)%valid = int_attr == 1
+  call shp_tr_open(hshp, trim(shpc_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 = 0, k1 = k1, &
+       ishape_last = [ishape_last])
+  CALL shp_tr_close(hshp)
 
-     call shp_read_object_03(hshp, ishape, psobject)
-     list_vis(i)%coord_extr = [psobject%padfx(1), psobject%padfy(1)] &
-          * deg_to_rad
-  end do
-
-  CALL shpclose(hshp)
-  nearby = nearby_extr(extr_map, list_vis, i = 6)
+  nearby = nearby_extr(s%extr_map, s%list_vis, i = 4)
   print *, "nearby:"
 
   do l = 1, size(nearby, 2)
      print *, nearby(:, l) * rad_to_deg
   end do
 
+  call mpi_finalize
+
 end program test_nearby_extr