diff --git a/Documentation_texfol/documentation.tex b/Documentation_texfol/documentation.tex
index b2e8580fdc6cc7c3582b8d911c0095ef97dcb108..de818de56f2632ce0914cc52c746a6bfc31cfa84 100644
--- a/Documentation_texfol/documentation.tex
+++ b/Documentation_texfol/documentation.tex
@@ -505,7 +505,7 @@ soit environ \np{0.5} KiB.
 \item[snapshot\%list\_vis] On a besoin d'encapsuler le vecteur de type
   eddy dans un scalaire de type snapshot pour pouvoir considérer un
   vecteur de snapshots. Notons $n_{ve}$ le nombre de tourbillons à une
-  date donnée (number\_vis\_eddies). Mémoire utilisée par ce champ :
+  date donnée (number\_vis\_extr). Mémoire utilisée par ce champ :
   \begin{equation*}
     n_{ve} \langle \mathrm{taille}(\mathtt{eddy}) \rangle
     = n_{ve}
@@ -1106,7 +1106,7 @@ flow(j)\%list\_vis\%delta\_in, flow\%number\_eddies. Cf. algorithme
   \begin{algorithmic}
     \STATE \COMMENT{recouvrements à dates non successives}
 
-    \FOR{i1 = 1 \TO flow(j - delta)\%number\_vis\_eddies}
+    \FOR{i1 = 1 \TO flow(j - delta)\%number\_vis\_extr}
 
     \IF{flow(j - delta)\%list\_vis(i1)\%valid}
 
diff --git a/Tests/long_tests.json b/Tests/long_tests.json
new file mode 100644
index 0000000000000000000000000000000000000000..29e95939682c01dfd484bc330c31f33cd6a2c9bf
--- /dev/null
+++ b/Tests/long_tests.json
@@ -0,0 +1,32 @@
+[
+    {
+	"args": ["$compil_prod_dir/test_set_all_outerm",
+		 "$input_dir/h_region_3.nc"],
+	"title": "Set_all_outerm",
+	"input": "&main_nml min_amp = 0.001/\n"
+    },
+    {
+	"args": "$compil_prod_dir/test_get_snapshot",
+	"title": "Get_snapshot_region_3",
+	"required": [["$input_dir/h_region_3.nc", "h.nc"],
+		     ["$input_dir/uv_region_3.nc", "uv.nc"]],
+	"input": "&main_nml /\n",
+	"description": "Larger region, 120 x 120. Includes degenerate extrema."
+    },
+    {
+	"args": "$compil_prod_dir/test_get_snapshot",
+	"title": "Get_snapshot_region_3_min",
+	"required": [["$input_dir/h_region_3.nc", "h.nc"],
+		     ["$input_dir/uv_region_3.nc", "uv.nc"]],
+	"input": "&main_nml min_amp = 1e-3/\n",
+	"description": "Same as Get_snapshot_region_3 except with 1 mm minimum amplitude."
+    },
+    {
+	"args": "$compil_prod_dir/test_get_snapshot",
+	"title": "Get_snapshot_region_5",
+	"required": [["$input_dir/h_region_5.nc", "h.nc"],
+		     ["$input_dir/uv_region_5.nc", "uv.nc"]],
+	"input": "&main_nml min_amp = 1e-3/\n",
+	"description": "Same as Get_snapshot_region_4 with larger domain."
+    }
+]
diff --git a/Tests/plot_snapshot.py b/Tests/plot_snapshot.py
index 488e1a0b6bbaaa93b98594f789e1c6164e42e8a1..e2781e4baa5d7476952f592d2ab69efb499044de 100755
--- a/Tests/plot_snapshot.py
+++ b/Tests/plot_snapshot.py
@@ -55,8 +55,8 @@ if args.window:
     lon_mask = longitude <= urcrnrlon
     lat_mask = np.logical_and(latitude >= llcrnrlat, latitude <= urcrnrlat)
 else:
-    lon_mask = True
-    lat_mask = True
+    lon_mask = np.ones(len(longitude), dtype = bool)
+    lat_mask = np.ones(len(latitude), dtype = bool)
 
 plt.figure()
 src_crs = ccrs.PlateCarree()
@@ -147,8 +147,8 @@ for shape_rec_extr, shape_outer, shape_m_s \
 if args.velocity:
     with netCDF4.Dataset("uv.nc") as f:
         quiver_return = ax.quiver(longitude[lon_mask], latitude[lat_mask],
-                                  f["u"][0, lat_mask, lon_mask],
-                                  f["v"][0, lat_mask, lon_mask],
+                                  f["u"][0, lat_mask][:, lon_mask],
+                                  f["v"][0, lat_mask][:, lon_mask],
                                   scale = args.scale, scale_units = "width",
                                   transform = src_crs)
     plt.quiverkey(quiver_return, 0.9, 0.9, 1, r"1 m s$^{-1}$",
diff --git a/Tests/tests.json b/Tests/short_tests.json
similarity index 91%
rename from Tests/tests.json
rename to Tests/short_tests.json
index bc9b88d3038c81dc06176e50a9bded5763a3c070..e0552b964bf8056316b7471cf81a7dba4358b3e1 100644
--- a/Tests/tests.json
+++ b/Tests/short_tests.json
@@ -129,12 +129,6 @@
 	],
 	"input": "&MAIN_NML IND_TARG_EXTR=19,11 /\n"
     },
-    {
-	"args": ["$compil_prod_dir/test_set_all_outerm",
-		 "$input_dir/h_region_3.nc"],
-	"title": "Set_all_outerm",
-	"input": "&main_nml min_amp = 0.001/\n"
-    },
     {
 	"args": "$compil_prod_dir/test_get_snapshot",
 	"title": "Get_snapshot_region_1",
@@ -176,22 +170,6 @@
 	"title": "Inside_4_false",
 	"stdin_filename": "$input_dir/inside_4_false_nml.txt"
     },
-    {
-	"args": "$compil_prod_dir/test_get_snapshot",
-	"title": "Get_snapshot_region_3",
-	"required": [["$input_dir/h_region_3.nc", "h.nc"],
-		     ["$input_dir/uv_region_3.nc", "uv.nc"]],
-	"input": "&main_nml /\n",
-	"description": "Larger region, 120 x 120. Includes degenerate extrema."
-    },
-    {
-	"args": "$compil_prod_dir/test_get_snapshot",
-	"title": "Get_snapshot_region_3_min",
-	"required": [["$input_dir/h_region_3.nc", "h.nc"],
-		     ["$input_dir/uv_region_3.nc", "uv.nc"]],
-	"input": "&main_nml min_amp = 1e-3/\n",
-	"description": "Same as Get_snapshot_region_3 except with 1 mm minimum amplitude."
-    },
     {
 	"args": "$compil_prod_dir/test_get_snapshot",
 	"title": "Get_snapshot_region_4",
@@ -199,14 +177,6 @@
 	"input": "&main_nml min_amp = 1e-3/\n",
 	"description": "Part of the domain has missing values."
     },
-    {
-	"args": "$compil_prod_dir/test_get_snapshot",
-	"title": "Get_snapshot_region_5",
-	"required": [["$input_dir/h_region_5.nc", "h.nc"],
-		     ["$input_dir/uv_region_5.nc", "uv.nc"]],
-	"input": "&main_nml min_amp = 1e-3/\n",
-	"description": "Same as Get_snapshot_region_4 with larger domain."
-    },
     {"args": "$compil_prod_dir/test_weight", "title": "Weight"},
     {
 	"args": "$compil_prod_dir/test_spherical_polyline_area",
diff --git a/Tests/stat.py b/Tests/stat.py
new file mode 100755
index 0000000000000000000000000000000000000000..f12b46d7a14cc377ad07a3f00f1a5a5545cb9a4a
--- /dev/null
+++ b/Tests/stat.py
@@ -0,0 +1,74 @@
+#!/usr/bin/env python3
+
+import shapefile
+import numpy as np
+
+n_valid_cycl = 0
+n_valid = 0
+n_minima = 0
+n_radius4 = np.zeros(3, dtype = int)
+n_outer = 0
+n_points_valid = 0
+n_points = 0
+
+with shapefile.Reader("extremum_1") as extremum, \
+     shapefile.Reader("outermost_contour_1") as outermost_contour:
+    n_extr = len(extremum)
+    
+    for rec_extr, shape_rec_outer in zip(extremum.iterRecords(),
+                                         outermost_contour):
+        if rec_extr.cyclone == 1: n_minima += 1
+        l = len(shape_rec_outer.shape.points)
+        n_points += l
+
+        if rec_extr.valid == 1:
+            n_valid += 1
+            n_points_valid += l
+            
+        if rec_extr.cyclone == 1 and rec_extr.valid == 1: n_valid_cycl += 1
+        if shape_rec_outer.record.r_eq_area >= 0: n_outer += 1
+        
+        if shape_rec_outer.record.radius4 == 0:
+            n_radius4[0] += 1
+        elif shape_rec_outer.record.radius4 == 1:
+            n_radius4[1] += 1
+        else:
+            n_radius4[2] += 1
+
+print("number of extrema of SSH:", len(extremum))
+print("number of minima of SSH:", n_minima)
+print("number of valid eddies (outermost contour found with sufficient area):",
+      n_valid)
+print("number of valid cyclones:", n_valid_cycl)
+print("number of extrema with radius4 == 0 (no outermost contour or "
+      "insufficient area):", n_radius4[0])
+print("number of extrema with radius4 == 1 (valid eddies with radius4 == 1):",
+      n_radius4[1])
+print("number of extrema with radius4 >= 2 (valid eddies with radius4 >= 2):",
+      n_radius4[2])
+print("number of outermost contours (outermost contour found, but may be with "
+      "sufficient area or not):", n_outer)
+print("average number of points of outermost contour per valid eddy:",
+      n_points_valid / n_valid)
+print("average number of points of outermost contour per extremum (valid eddy "
+      "or not):", n_points / n_extr)
+
+n_max_speed = 0
+n_points_valid = 0
+n_points = 0
+
+with shapefile.Reader("max_speed_contour_1") as max_speed_contour:
+    for shape_rec in max_speed_contour:
+        l = len(shape_rec.shape.points)
+        n_points += l
+        if shape_rec.record.r_eq_area >= 0:
+            n_max_speed += 1
+            n_points_valid += l
+
+print("number of non-null maximum-speed contours (distinct from the "
+      "outermost contour):",
+      n_max_speed)
+print("average number of points of non-null maximum-speed contours:",
+      n_points_valid / n_max_speed)
+print("average number of points of maximum-speed contours, null or not, per "
+      "extremum:", n_points / n_extr)
diff --git a/Tests/test_get_snapshot.f b/Tests/test_get_snapshot.f
index 6f7a3801b8866cf4c96939fc511a0c1efc2d2c91..2d5c5d2dae81bc0a70914e337653112dbd5678ac 100644
--- a/Tests/test_get_snapshot.f
+++ b/Tests/test_get_snapshot.f
@@ -110,13 +110,4 @@ program test_get_snapshot
   close(unit_number_eddies)
   print *, 'Created "number_eddies_1.csv".'
 
-  if (s%number_vis_extr /= 0) then
-     print *, "Average number of points per outermost contour: ", &
-          sum(s%list_vis%out_cont%n_points) / real(s%number_vis_extr)
-     print *, "Average number of points per max-speed contour: ", &
-          sum(s%list_vis%speed_cont%n_points) / real(s%number_vis_extr)
-  else
-     print *, "No visible eddy."
-  end if
-
 end program test_get_snapshot
diff --git a/set_max_speed.f b/set_max_speed.f
index 29bfe7bf491587f712dbbda49121508d61fd2982..47b435eb26e7d6187bb6b08367f33dff1a69f11a 100644
--- a/set_max_speed.f
+++ b/set_max_speed.f
@@ -9,7 +9,8 @@ contains
 
     ! Defines the components speed_cont, max_speed and radius4 of
     ! argument e. On return, speed_cont may be a null ssh contour but
-    ! max_speed is never set to missing_speed.
+    ! max_speed is never set to missing_speed. The value of radius4
+    ! set here is always >= 1.
 
     ! The length of the interval of longitudes of the domain, step(1)
     ! * (size(ssh, 1) - 1), should be < 180 degrees, so that the