Add arguments corner, nlon and nlat to procedure read_snapshot. In

read_snapshot, assert that there is no duplicate tuple in
s%ind_extr. Since we are not sure of that, we must use a do construct
and not a forall in order to define s%extr_map.

Bug fix in program test_local_extrema: support missing values. Add
corresponding test.

Output s%ind_extr in program test_read_snapshot.

Tests/Stdin/read_snapshot_nml.txt

+ &MAIN_NML
+ CORNER  = 16.12500       , -38.87500      ,
+ NLON    = 20,
+ NLAT    = 20
+ /

Tests/read_snapshot.f

@@ -4,13 +4,14 @@ module read_snapshot_m
 
 contains
 
-  subroutine read_snapshot(s, hshp_extremum, hshp_outermost, hshp_max_speed)
+  subroutine read_snapshot(s, hshp_extremum, hshp_outermost, hshp_max_speed, &
+       corner, nlon, nlat)
 
     use, intrinsic:: ISO_C_BINDING
 
     use contour_531, only: convert_to_ind
     use derived_types, only: snapshot, eddy
-    use nr_util, only: pi
+    use nr_util, only: pi, assert, deg_to_rad
     use read_eddy_m, only: read_eddy
     use shapelib, only: shpfileobject, shpgetinfo
 
@@ -23,12 +24,16 @@ contains
     TYPE(shpfileobject), intent(inout):: hshp_max_speed
     ! shapefile max_speed_contour_1
 
+    real, intent(in):: corner(:) ! (2) longitude and latitude of the
+                                 ! corner of the whole grid, in rad
+
+    integer, intent(in):: nlon, nlat
+
     ! Local:
     integer ishape, k, i
     integer shapetype, dbffieldcount, dbfrecordcount
     real(c_double) minbound(4), maxbound(4)
     type(eddy) e
-    real, parameter:: deg_to_rad = pi / 180.
 
     !---------------------------------------------------------------------
 
@@ -43,17 +48,20 @@ contains
     end do
 
     allocate(s%ind_extr(2, s%number_vis_eddies))
-
     forall (i = 1:s%number_vis_eddies) s%ind_extr(:, i) &
-         = nint(convert_to_ind(s%list_vis(i)%coord_extr, &
-         corner = [0.125, -59.875] * deg_to_rad, &
+         = nint(convert_to_ind(s%list_vis(i)%coord_extr, corner, &
          step = [0.25, 0.25] * deg_to_rad))
 
-    allocate(s%extr_map(120, 180))
+    allocate(s%extr_map(nlon, nlat))
     s%extr_map = 0
-    forall (i = 1:s%number_vis_eddies) &
-         s%extr_map(s%ind_extr(1, i), s%ind_extr(2, i)) &
-         = merge(i, -i, s%list_vis(i)%suff_amp)
+    
+    do i = 1, s%number_vis_eddies
+       s%extr_map(s%ind_extr(1, i), s%ind_extr(2, i)) &
+            = merge(i, -i, s%list_vis(i)%suff_amp)
+    end do
+
+    call assert(count(s%extr_map /= 0) == s%number_vis_eddies, &
+         "read_snapshot: one eddy per grid point")
 
     s%number_eddies = s%number_vis_eddies
 

Tests/test_get_snapshot.f

@@ -11,14 +11,12 @@ program test_get_snapshot
   use netcdf, only: nf90_nowrite
   use netcdf95, only: nf95_open, find_coord, nf95_inquire_dimension, &
        nf95_get_var, nf95_close
-  use nr_util, only: pi, assert
+  use nr_util, only: pi, assert, deg_to_rad
   use shapelib, only: shpfileobject, shpclose
 
   implicit none
 
   type(snapshot) s
-  real, parameter:: deg_over_rad = pi / 180.
-
   TYPE(shpfileobject) hshp_extremum ! shapefile extremum_$m
   TYPE(shpfileobject) hshp_outermost ! shapefile outermost_contour_$m
   TYPE(shpfileobject) hshp_max_speed ! shapefile max_speed_contour_$m
@@ -57,9 +55,9 @@ program test_get_snapshot
 
   call get_snapshot(s, min_amp, m = 1, n_proc = 1, k_end = 1, max_delta = 4, &
        nlon = nlon, nlat = nlat, k = 1, max_radius = [20, 20], &
-       corner = [longitude(1), latitude(1)] * deg_over_rad, &
+       corner = [longitude(1), latitude(1)] * deg_to_rad, &
        step = &
-       [longitude(2) - longitude(1), latitude(2) - latitude(1)]  * deg_over_rad)
+       [longitude(2) - longitude(1), latitude(2) - latitude(1)]  * deg_to_rad)
 
   call init_shapefiles(hshp_extremum, hshp_outermost, hshp_max_speed)
   

Tests/test_local_extrema.f

@@ -5,7 +5,7 @@ program test_local_extrema
   use netcdf, only: nf90_nowrite, nf90_clobber, NF90_FLOAT, NF90_INT
   use netcdf95, only: nf95_open, nf95_close, nf95_inq_varid, nf95_get_var, &
        nf95_put_var, nf95_create, nf95_def_dim, nf95_def_var, nf95_put_att, &
-       nf95_enddef, nf95_gw_var
+       nf95_enddef, nf95_gw_var, nf95_get_att
   use nr_util, only: assert
 
   implicit none
@@ -28,6 +28,7 @@ program test_local_extrema
   ! level of innermost contour, for each extremum
 
   logical, allocatable:: cyclone(:) ! (n_extr)
+  real Fill_Value
 
   !---------------------------------------------------------------------
 
@@ -51,6 +52,8 @@ program test_local_extrema
 
   call nf95_inq_varid(ncid, "adt", varid)
   call nf95_get_var(ncid, varid, ssh)
+  call nf95_get_att(ncid, varid, "_FillValue", Fill_Value)
+  where (ssh == Fill_Value) ssh = huge(0.)
 
   call nf95_close(ncid)
 
@@ -59,10 +62,13 @@ program test_local_extrema
   call new_unit(unit)
   open(unit, file = "test_local_extrema.csv", status = "replace", &
        action = "write")
-  write(unit, fmt = *) '"longitude index" "latitude index" cyclone ', &
-       '"longitude (degrees east)" ', &
-       '"latitude (degrees north)" "ssh (m)" "innermost_level (m)"'
-  ! (title line)
+
+  ! Title lines:
+  write(unit, fmt = *) '"longitude index" "latitude index"'
+  write(unit, fmt = *) '"" "" "" "degrees east" "degrees north" "m" "m"'
+  write(unit, fmt = *) 'ilon ilat cyclone longitude latitude ssh ', &
+       'innermost_level'
+
   do i = 1, size(ind_extr, 2)
      write(unit, fmt = *) ind_extr(1, i), ind_extr(2, i), cyclone(i), &
           longitude(ind_extr(1, i)), latitude(ind_extr(2, i)), &

Tests/test_read_snapshot.f

 program test_read_snapshot
 
+  use, intrinsic:: ISO_FORTRAN_ENV
+
   use derived_types, only: snapshot
   use dispatch_snapshot_m, only: dispatch_snapshot
   use init_shapefiles_m, only: init_shapefiles
   use jumble, only: new_unit
+  use nr_util, only: deg_to_rad
   use read_snapshot_m, only: read_snapshot
   use shapelib, only: shpfileobject, shpclose
   use shapelib_03, only: shp_open_03
@@ -14,22 +17,34 @@ program test_read_snapshot
   TYPE(shpfileobject) hshp_extremum ! shapefile extremum_1
   TYPE(shpfileobject) hshp_outermost ! shapefile outermost_contour_1
   TYPE(shpfileobject) hshp_max_speed ! shapefile max_speed_contour_1
-  integer unit_isolated, unit_number_eddies
+  integer unit_isolated, unit_number_eddies, i
+
+  real:: corner(2) = [0.125, - 59.875]
+  ! longitude and latitude of the corner of the whole grid, in degrees
+
+  integer:: nlon = 120, nlat = 120
+  namelist /main_nml/ corner, nlon, nlat
 
   !-------------------------------------------------------------------------
 
+  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 shp_open_03("extremum_1_old", pszaccess = "rb", hshp = hshp_extremum)
   call shp_open_03("outermost_contour_1_old", pszaccess = "rb", &
        hshp = hshp_outermost)
   call shp_open_03("max_speed_contour_1_old", pszaccess = "rb", &
        hshp = hshp_max_speed)
-  call read_snapshot(s, hshp_extremum, hshp_outermost, hshp_max_speed)
+  call read_snapshot(s, hshp_extremum, hshp_outermost, hshp_max_speed, &
+       corner * deg_to_rad, nlon, nlat)
   CALL shpclose(hshp_extremum)
   CALL shpclose(hshp_outermost)
   CALL shpclose(hshp_max_speed)
 
   call init_shapefiles(hshp_extremum, hshp_outermost, hshp_max_speed)
-  
+
   call new_unit(unit_isolated)
   open(unit_isolated, file = "isolated_nodes_1.csv", status = "replace", &
        action = "write")
@@ -57,4 +72,10 @@ program test_read_snapshot
   close(unit_number_eddies)
   print *, 'Created "number_eddies_1.csv".'
 
+  print *, "s%ind_extr:"
+
+  do i = 1, s%number_vis_eddies
+     print *, s%ind_extr(:, i)
+  end do
+
 end program test_read_snapshot

Tests/tests.json

@@ -211,6 +211,13 @@
     {
 	"args": "$compil_prod_dir/test_read_snapshot",
 	"required": ["$input_dir/Test_read_eddy/*"],
-	"title": "Read_snapshot"
+	"title": "Read_snapshot",
+	"stdin": "$stdin_dir/read_snapshot_nml.txt"
+    },
+    {
+	"args" : ["$compil_prod_dir/test_local_extrema",
+		  "$input_dir/2006_01/h_region_4.nc"],
+	"title" : "Local_extrema_missing",
+	"description": "test_local_extrema with input file containing missing values."
     }
 ]

get_snapshot.f

@@ -36,7 +36,7 @@ contains
     ! eddy in longitude and latitude, in number of grid points
 
     real, intent(in):: corner(:) ! (2) longitude and latitude of the
-                                 ! corner of the global grid, in rad
+                                 ! corner of the whole grid, in rad
     real, intent(in):: step(:) ! (2) longitude and latitude steps, in rad
 
     ! Local: