Allow dist_lim to be chosen at run-time. So make it an argument of

read_snapshot and successive_overlap instead of a named constant of
module successive_overlap_m.

Bug fix in main programs test_read_snapshot and
test_successive_overlap: take into account that step is in degrees.

Tests/CMakeLists.txt

@@ -49,8 +49,7 @@ target_link_libraries(test_successive_overlap ${contour_531_LIBRARY}
 
 add_executable(test_read_snapshot derived_types.f90
   init_shapefiles.f90 read_snapshot.f90 write_eddy.f90 read_eddy.f90
-  read_field_indices.f90 successive_overlap.f90 spher_polygon_area.f90
-  spher_polyline_area.f90 weight.f90 candidate_overlap.f90
+  read_field_indices.f90
   ${CMAKE_CURRENT_LIST_DIR}/test_read_snapshot.f90
   ${CMAKE_CURRENT_LIST_DIR}/write_extr_map.f90)
 
@@ -61,7 +60,7 @@ target_include_directories(test_read_snapshot PRIVATE
 target_link_libraries(test_read_snapshot ${contour_531_LIBRARY}
   ${shapelib_03_LIBRARY} ${nr_util_LIBRARY} ${fortrangis_LIBRARY}
   ${fortranc_LIBRARY} ${shapelib_LIBRARY} ${GPC_F_LIBRARY} ${GPC_LIBRARY}
-  ${netcdf95_LIBRARY} ${netcdff_LIBRARY} ${NetCDF_LIBRARY} ${geometry_LIBRARY})
+  ${netcdf95_LIBRARY} ${netcdff_LIBRARY} ${NetCDF_LIBRARY})
 
 # test_spher_polygon_area
 

Tests/test_read_snapshot.f90

@@ -3,14 +3,13 @@ program test_read_snapshot
   use, intrinsic:: ISO_FORTRAN_ENV
 
   ! Libraries:
-  use nr_util, only: deg_to_rad, twopi, pi, arth, assert
+  use nr_util, only: deg_to_rad, arth, assert
   use shapelib, only: shpfileobject, shpclose
   use shapelib_03, only: shp_open_03
 
   use derived_types, only: snapshot
   use init_shapefiles_m, only: init_shapefiles
   use read_snapshot_m, only: read_snapshot
-  use successive_overlap_m, only: dist_lim
   use write_eddy_m, only: write_eddy
   use write_extr_map_m, only: write_extr_map
 
@@ -28,7 +27,12 @@ program test_read_snapshot
   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
-  namelist /main_nml/ corner, step, nlon, nlat
+
+  integer:: dist_lim = 12
+  ! We look for an overlapping eddy at dist_lim (in grid points) of
+  ! the first extremum.
+
+  namelist /main_nml/ corner, step, nlon, nlat, dist_lim
 
   !-------------------------------------------------------------------------
 
@@ -38,10 +42,10 @@ program test_read_snapshot
   write(unit = *, nml = main_nml)
 
   ! As we are requiring the grid spacing to be uniform, the value of
-  ! "periodic" must be consistent with the values of step(1) and nlon:
-  periodic = nint(twopi / step(1)) == nlon
+  ! "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) < pi, &
+  if (periodic) call assert(2 * dist_lim * step(1) < 180., &
        "test_read_snapshot dist_lim")
 
   call shp_open_03("extremum_old", pszaccess = "rb", hshp = hshp_extremum)
@@ -50,7 +54,7 @@ program test_read_snapshot
   call shp_open_03("max_speed_contour_old", pszaccess = "rb", &
        hshp = hshp_max_speed)
   call read_snapshot(s, k, hshp_extremum, hshp_outermost, hshp_max_speed, &
-       corner * deg_to_rad, step * deg_to_rad, nlon, nlat, periodic)
+       corner * deg_to_rad, step * deg_to_rad, nlon, nlat, periodic, dist_lim)
   CALL shpclose(hshp_extremum)
   CALL shpclose(hshp_outermost)
   CALL shpclose(hshp_max_speed)

Tests/test_set_all_outerm.f90

@@ -79,7 +79,7 @@ program test_set_all_outerm
        * deg_to_rad
 
   ! As we are requiring the grid spacing to be uniform, the value of
-  ! "periodic" must be consistent with the values of step(1) and nlon:
+  ! "periodic" may be deduced from the values of step(1) and nlon:
   periodic = nint(twopi / step(1)) == nlon
   print *, "periodic = ", periodic
   call assert(2 * max_radius(1) * step(1) < pi, &

Tests/test_successive_overlap.f90

@@ -4,13 +4,13 @@ program test_successive_overlap
 
   ! Libraries:
   use jumble, only: new_unit
-  use nr_util, only: deg_to_rad, twopi, assert, pi
+  use nr_util, only: deg_to_rad, assert
   use shapelib, only: shpfileobject, shpclose
   use shapelib_03, only: shp_open_03
 
   use derived_types, only: snapshot
   use read_snapshot_m, only: read_snapshot
-  use successive_overlap_m, only: successive_overlap, dist_lim
+  use successive_overlap_m, only: successive_overlap
 
   implicit none
 
@@ -26,7 +26,12 @@ program test_successive_overlap
   real:: step(2) = [0.25, 0.25] ! longitude and latitude steps, in degrees
   logical periodic ! grid is periodic in longitude
   integer:: nlon = - 1, nlat = - 1
-  namelist /main_nml/ corner, step, nlon, nlat
+
+  integer:: dist_lim = 12
+  ! We look for an overlapping eddy at dist_lim (in grid points) of
+  ! the first extremum.
+
+  namelist /main_nml/ corner, step, nlon, nlat, dist_lim
 
   !-------------------------------------------------------------------------
 
@@ -36,10 +41,10 @@ program test_successive_overlap
   write(unit = *, nml = main_nml)
 
   ! As we are requiring the grid spacing to be uniform, the value of
-  ! "periodic" must be consistent with the values of step(1) and nlon:
-  periodic = nint(twopi / step(1)) == nlon
+  ! "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) < pi, &
+  if (periodic) call assert(2 * dist_lim * step(1) < 180., &
        "test_successive_overlap dist_lim")
 
   call shp_open_03("extremum", pszaccess = "rb", hshp = hshp_extremum)
@@ -49,14 +54,14 @@ program test_successive_overlap
        hshp = hshp_max_speed)
   call read_snapshot(flow(1), k, hshp_extremum, hshp_outermost, &
        hshp_max_speed, corner * deg_to_rad, step * deg_to_rad, nlon, nlat, &
-       periodic)
+       periodic, dist_lim)
   CALL shpclose(hshp_extremum)
   CALL shpclose(hshp_outermost)
   CALL shpclose(hshp_max_speed)
 
   flow(2) = flow(1)
   call new_unit(unit_edgelist)
-  open(unit_edgelist, file = "edgelist_1.csv", status = "replace", &
+  open(unit_edgelist, file = "edgelist.csv", status = "replace", &
        action  = "write")
 
   ! Title lines:
@@ -64,11 +69,11 @@ program test_successive_overlap
        // '"predecessor eddy subscript" "successor date subscript" ' &
        // '"successor eddy subscript"'
   write(unit_edgelist, fmt = *) "k1 i1 k2 i2 weight"
-  
-  call successive_overlap(flow, unit_edgelist, nlon, nlat, periodic, j = 2, &
-       k = k + 1)
+
+  call successive_overlap(flow, unit_edgelist, nlon, nlat, periodic, dist_lim, &
+       j = 2, k = k + 1)
   close(unit_edgelist)
-  print *, 'Created file "edgelist_1.csv".'
+  print *, 'Created file "edgelist.csv".'
 
   do j = 1, 2
      print *, "i, flow(", j, ")%list_vis%delta_in:"
@@ -79,7 +84,7 @@ program test_successive_overlap
      end do
 
      print *, "i, flow(", j, ")%list_vis%delta_out:"
-     
+
      do i = 1, flow(j)%number_vis_extr
         if (flow(j)%list_vis(i)%delta_out /= huge(0)) &
              print *, i, flow(j)%list_vis(i)%delta_out

extraction_eddies.f90

@@ -109,7 +109,7 @@ program extraction_eddies
        * deg_to_rad
 
   ! As we are requiring the grid spacing to be uniform, the value of
-  ! "periodic" must be consistent with the values of step(1) and nlon:
+  ! "periodic" may be deduced from the values of step(1) and nlon:
   periodic = nint(twopi / step(1)) == nlon
   print *, "periodic = ", periodic
   

read_snapshot.f90

@@ -5,7 +5,7 @@ module read_snapshot_m
 contains
 
   subroutine read_snapshot(s, k, hshp_extremum, hshp_outermost, &
-       hshp_max_speed, corner, step, nlon, nlat, periodic)
+       hshp_max_speed, corner, step, nlon, nlat, periodic, dist_lim)
 
     ! Libraries:
     use contour_531, only: convert_to_ind 
@@ -16,7 +16,6 @@ contains
     use derived_types, only: snapshot, eddy
     use read_eddy_m, only: read_eddy
     use read_field_indices_m, only: read_field_indices
-    use successive_overlap_m, only: dist_lim
 
     type(snapshot), intent(out):: s ! completely defined
     integer, intent(out):: k ! date index
@@ -35,6 +34,10 @@ contains
     integer, intent(in):: nlon, nlat
     logical, intent(in):: periodic ! grid is periodic in longitude
 
+    integer, intent(in):: dist_lim
+    ! We look for an overlapping eddy at dist_lim (in grid points) of
+    ! the first extremum.
+
     ! Local:
     integer ishape, copy
     integer k2 ! date index

successive_overlap.f90

@@ -2,13 +2,10 @@ module successive_overlap_m
 
   implicit none
 
-    integer, parameter:: dist_lim = 12
-    ! We look for an overlapping eddy at dist_lim (in grid points) of
-    ! the first extremum.
-
 contains
 
-  subroutine successive_overlap(flow, unit_edgelist, nlon, nlat, periodic, j, k)
+  subroutine successive_overlap(flow, unit_edgelist, nlon, nlat, periodic, &
+       dist_lim, j, k)
 
     ! Finds edges between flow(j - 1) and flow(j), corresponding to
     ! dates k - 1 and k. Writes these edges to unit_edgelist. Updates
@@ -29,6 +26,11 @@ contains
     integer, intent(in):: unit_edgelist ! logical unit for edgelist file
     integer, intent(in):: nlon, nlat
     logical, intent(in):: periodic ! grid is periodic in longitude
+
+    integer, intent(in):: dist_lim
+    ! We look for an overlapping eddy at dist_lim (in grid points) of
+    ! the first extremum.
+
     integer, intent(in):: j, k
 
     ! Local:
@@ -81,7 +83,7 @@ contains
              m = floor((flow(j - 1)%list_vis(i1)%coord_extr(1) &
                   - flow(j)%list_vis(i2)%coord_extr(1)) / twopi + 0.5)
              polyline_2%points(1, :) = polyline_2%points(1, :) + m * twopi
-             
+
              call gpc_polygon_clip_f(GPC_INT, polygon(nparts = 1, &
                   part = [polyline_1], hole = [.false.]), &
                   polygon(nparts = 1, part = [polyline_2], hole = [.false.]), &