In procedure get_snapshot, define nearby_extr in two steps. First,

define a selection of extrema. Then get their coordinates. This
prepares the way for a more complex selection: taking cylonicity into account.

get_snapshot.f

@@ -10,8 +10,7 @@ contains
     use, intrinsic:: ieee_arithmetic, only: ieee_value, IEEE_QUIET_NAN
 
     ! Libraries:
-    use contour_531, only: convert_to_ind, null_polyline, convert_to_reg_coord
-    use jumble, only: argwhere
+    use contour_531, only: convert_to_ind, null_polyline
     use netcdf, only: nf90_nowrite
     use netcdf95, only: nf95_open, nf95_close, nf95_inq_varid, nf95_get_var, &
          nf95_get_att
@@ -48,10 +47,13 @@ contains
     real ssh(nlon, nlat) ! sea-surface height, in m
     real u(nlon, nlat), v(nlon, nlat) ! wind, in m s-1
     real Fill_Value
-    integer i
+    integer i, n_select, l
+
+    integer, allocatable:: selection(:)
+    ! identifying numbers of a selection of eddies
 
     real, allocatable:: nearby_extr(:, :) ! (2, :) longitude and
-    ! latitude, in rad, of all the extrema except the target extremum
+    ! latitude, in rad, of extrema near the target extremum
 
     ! Window around each extremum:
     integer llc(2) ! indices in global grid of lower left corner
@@ -107,15 +109,19 @@ contains
                   [nlon, nlat])
              ! (min should have no effect except because of roundup error)
              
-             nearby_extr &
-                  = convert_to_reg_coord(argwhere(s%extr_map(llc(1):urc(1), &
-                  llc(2):urc(2)) > 0 .and. s%extr_map(llc(1):urc(1), &
-                  llc(2):urc(2)) /= i), corner + (llc - 1) * step, step)
+             selection = abs(pack(s%extr_map(llc(1):urc(1), llc(2):urc(2)), &
+                  s%extr_map(llc(1):urc(1), llc(2):urc(2)) > 0 &
+                  .and. s%extr_map(llc(1):urc(1), llc(2):urc(2)) /= i))
+             n_select = size(selection)
+             allocate(nearby_extr(2, n_select))
+             forall (l = 1:n_select) &
+                  nearby_extr(:, l) = s%list_vis(selection(l))%coord_extr
              call set_max_speed(s%list_vis(i), s%ind_extr(:, i) - llc + 1, &
                   nearby_extr, ssh(llc(1):urc(1), llc(2):urc(2)), &
                   u(llc(1):urc(1), llc(2):urc(2)), &
                   v(llc(1):urc(1), llc(2):urc(2)), corner + (llc - 1) * step, &
                   step)
+             deallocate(nearby_extr)
           else
              s%list_vis(i)%speed_cont = null_ssh_contour()
              s%list_vis(i)%max_speed = missing_speed