From 0ce6b1b49e114def56714d85a766cf5f7bca7992 Mon Sep 17 00:00:00 2001
From: Lionel GUEZ <guez@lmd.ens.fr>
Date: Mon, 4 Dec 2017 19:10:26 +0100
Subject: [PATCH] If an outermost contour has exactly the minimum amplitude, it
 is acceptable. So define flat_extr by a strict inequality in procedure
 get_snapshot.

Procedure set_outermost_contour no longer has a dummy argument
noise_around. We no longer call outermost_possible_level so we do not
have an artificial discontinuity of results with and without maximum
amplitude. We abandon the idea of reducing the amplitude of initial
level_good compared to innermost_level: problem of consistence with
get_snapshot, and it does not seem worth the trouble, just abandon
those problematic extrema. Allow for null outermost_contour instead of
aborting.

In procedure get_snapshot, s%list_vis(i)%suff_amp is first defined
only if flat_extr(i). Also we are no longer sure that
set_outermost_contour finds an outermost contour so we have to test
this to define s%list_vis(i)%suff_amp. noise_around is now defined
only for flat_extr(i) and s%list_vis(i)%suff_amp. Also, since we may
not find an outermost contour even if not flat_extr(i), we update
s%extr_map after the second call to set_outermost_contour.

In procedure local_extrema, we use the mask diff_center instead of
mask_center. So we no longer need procedure construct_mask_center. In
procedure local_extrema, we ignore one in two adjacent degenerate
extrema instead of aborting.

In plot_snapshot.py, color extrema and allow for null outermost
contour.

Synthesize output in test_local_extrema.py.
---
 Sources/derived_types.f                   | 12 +--
 Sources/get_snapshot.f                    | 69 ++++++++++------
 Sources/local_extrema.f                   | 78 ++++++++++--------
 Sources/set_outermost_contour.f           | 96 +++++++++++------------
 Tests/Stdin/set_outermost_contour_nml.txt |  1 -
 Tests/examine_ssh_values.py               | 27 +++++++
 Tests/plot_snapshot.py                    | 38 ++++-----
 Tests/test_get_snapshot.f                 |  2 -
 Tests/test_local_extrema.f                |  3 +-
 Tests/test_local_extrema.py               | 21 +++--
 Tests/test_set_outermost_contour.f        |  6 +-
 Tests/tests.json                          |  4 +-
 12 files changed, 203 insertions(+), 154 deletions(-)
 delete mode 100644 Tests/Stdin/set_outermost_contour_nml.txt
 create mode 100755 Tests/examine_ssh_values.py

diff --git a/Sources/derived_types.f b/Sources/derived_types.f
index f818a221..1d096114 100644
--- a/Sources/derived_types.f
+++ b/Sources/derived_types.f
@@ -17,8 +17,9 @@ module derived_types
      logical cyclone
      type(ssh_contour) outermost_contour, max_speed_contour
      real max_speed ! average on max_speed_contour, in m s-1
-     logical suff_amp ! sufficient amplitude of ssh, between extremum
-                      ! and outermost contour
+     logical suff_amp ! outermost_contour found and sufficient
+                      ! amplitude of ssh between extremum and
+                      ! outermost contour
      
      integer:: delta_in = huge(0)
      ! Maximum difference in time index where there is a direct
@@ -36,9 +37,10 @@ module derived_types
      integer number_vis_eddies ! number of visible eddies
 
      integer, allocatable:: extr_map(:, :) ! (nlon, nlat) At a point
-     ! of extremum SSH with sufficient amplitude: identification
-     ! number or this extremum. At a point of extremum SSH with
-     ! insufficient amplitude: the opposite of the identification
+     ! of extremum SSH with outermost_contour found, of sufficient
+     ! amplitude: identification number or this extremum. At a point
+     ! of extremum SSH with no outermost_contour or outermost_contour
+     ! of insufficient amplitude: the opposite of the identification
      ! number or this extremum. 0 at other points.
 
      integer, allocatable:: ind_extr(:, :) ! (2, number_vis_eddies)
diff --git a/Sources/get_snapshot.f b/Sources/get_snapshot.f
index af910cb9..f7abcf38 100644
--- a/Sources/get_snapshot.f
+++ b/Sources/get_snapshot.f
@@ -93,7 +93,7 @@ contains
             ind_targ_extr(2, s%number_vis_eddies), &
             flat_extr(s%number_vis_eddies), noise_around(s%number_vis_eddies))
 
-       do i = 1, s%number_vis_eddies
+       forall (i = 1:s%number_vis_eddies)
           s%list_vis(i)%coord_extr = corner + (s%ind_extr(:, i) - 1) * step
           s%list_vis(i)%ssh_extremum = ssh(s%ind_extr(1, i), s%ind_extr(2, i))
           s%list_vis(i)%cyclone = cyclone(i)
@@ -103,12 +103,12 @@ contains
           urc(:, i) = min(s%ind_extr(:, i) + max_radius, [nlon, nlat])
           corner_window(:, i) = corner + (llc(:, i) - 1) * step
           ind_targ_extr(:, i) = s%ind_extr(:, i) - llc(:, i) + 1
-       end do
+       end forall
 
        if (min_amp == 0.) then
           flat_extr = .false.
        else
-          flat_extr = abs(innermost_level - s%list_vis%ssh_extremum) <= min_amp
+          flat_extr = abs(innermost_level - s%list_vis%ssh_extremum) < min_amp
        end if
 
        do i = 1, s%number_vis_eddies
@@ -117,41 +117,58 @@ contains
                   innermost_level(i), &
                   s%extr_map(llc(1, i):urc(1, i), llc(2, i):urc(2, i)), &
                   ssh(llc(1, i):urc(1, i), llc(2, i):urc(2, i)), &
-                  corner_window(:, i), step, noise_around = .false.)
-             s%list_vis(i)%suff_amp = abs(s%list_vis(i)%outermost_contour%ssh &
-                  - s%list_vis(i)%ssh_extremum) >= min_amp
-          else
-             s%list_vis(i)%suff_amp = .true.
+                  corner_window(:, i), step)
+             if (s%list_vis(i)%outermost_contour%n_points == 0) then
+                s%list_vis(i)%suff_amp = .false.
+             else
+                s%list_vis(i)%suff_amp &
+                     = abs(s%list_vis(i)%outermost_contour%ssh &
+                     - s%list_vis(i)%ssh_extremum) >= min_amp
+             end if
           end if
        end do
 
+       ! We must modify s%extr_map in this loop, separate from the
+       ! previous loop, so we do not influence the batch of
+       ! set_outermost_contour.
        do i = 1, s%number_vis_eddies
-          if (.not. s%list_vis(i)%suff_amp) s%extr_map(s%ind_extr(1, i), &
-               s%ind_extr(2, i)) = - s%extr_map(s%ind_extr(1, i), &
-               s%ind_extr(2, i))
+          if (flat_extr(i) .and. .not. s%list_vis(i)%suff_amp) &
+               s%extr_map(s%ind_extr(1, i), s%ind_extr(2, i)) &
+               = - s%extr_map(s%ind_extr(1, i), s%ind_extr(2, i))
        end do
+       ! (no forall because of gfortran internal compiler error)
 
        ! Define noise_around:
-       if (min_amp == 0.) then
-          noise_around = .false.
-       else
-          do i = 1, s%number_vis_eddies
-             if (s%list_vis(i)%suff_amp) noise_around(i) &
-                  = any(s%extr_map(llc(1, i):urc(1, i), llc(2, i):urc(2, i)) &
-                  < 0)
-          end do
-       end if
+       if (min_amp /= 0.) &
+            forall &
+            (i = 1:s%number_vis_eddies, &
+            flat_extr(i) .and. s%list_vis(i)%suff_amp) &
+            noise_around(i) &
+            = any(s%extr_map(llc(1, i):urc(1, i), llc(2, i):urc(2, i)) < 0)
 
        do i = 1, s%number_vis_eddies
           if (s%list_vis(i)%suff_amp .and. noise_around(i) &
-               .or. .not. flat_extr(i)) &
-               call set_outermost_contour(s%list_vis(i), ind_targ_extr(:, i), &
-               innermost_level(i), &
-               s%extr_map(llc(1, i):urc(1, i), llc(2, i):urc(2, i)), &
-               ssh(llc(1, i):urc(1, i), llc(2, i):urc(2, i)), &
-               corner_window(:, i), step, noise_around(i))
+               .or. .not. flat_extr(i)) then
+             call set_outermost_contour(s%list_vis(i), ind_targ_extr(:, i), &
+                  innermost_level(i), &
+                  s%extr_map(llc(1, i):urc(1, i), llc(2, i):urc(2, i)), &
+                  ssh(llc(1, i):urc(1, i), llc(2, i):urc(2, i)), &
+                  corner_window(:, i), step)
+             s%list_vis(i)%suff_amp &
+                  = s%list_vis(i)%outermost_contour%n_points /= 0
+          end if
        end do
        
+       ! We must modify s%extr_map in this loop, separate from the
+       ! previous loop, so we do not influence the batch of
+       ! set_outermost_contour.
+       do i = 1, s%number_vis_eddies
+          if (.not. flat_extr(i) .and. .not. s%list_vis(i)%suff_amp) &
+               s%extr_map(s%ind_extr(1, i), s%ind_extr(2, i)) &
+               = - s%extr_map(s%ind_extr(1, i), s%ind_extr(2, i))
+       end do
+       ! (no forall because of gfortran internal compiler error)
+
        ! Done with outermost contours, now let us take care of the
        ! max-speed contours.
        
diff --git a/Sources/local_extrema.f b/Sources/local_extrema.f
index 3ebd3b1d..bcb80210 100644
--- a/Sources/local_extrema.f
+++ b/Sources/local_extrema.f
@@ -2,8 +2,6 @@ module local_extrema_m
 
   implicit none
 
-  logical, save, private::  mask_center(3, 3)
-
 contains
 
   subroutine local_extrema(field, extr_map, ind_extr, innermost_level, &
@@ -33,12 +31,13 @@ contains
     logical, intent(out), allocatable:: local_min(:) ! (n_extr)
 
     ! Local:
-    integer n_extr, i, j, k, m, n
+    integer n_extr, i, j, k, m, n, count_diff
 
     real innermost_map(size(field, 1), size(field, 2)) ! (m, n) map of
     ! levels for innermost contours, defined only for extremum points
 
     real max_around, min_around
+    logical diff_center(3, 3) ! values different from the center value
 
     !--------------------------------------------------------------------
 
@@ -55,27 +54,51 @@ contains
 
     do j = 2, n - 1
        do i = 2, m - 1
-          max_around = maxval(field(i - 1:i + 1, j - 1:j + 1), mask_center)
+          diff_center = field(i - 1:i + 1, j - 1:j + 1) /= field(i, j)
+          count_diff = count(diff_center)
 
-          if (field(i, j) > max_around) then
-             n_extr = n_extr + 1
-             extr_map(i, j) = n_extr
-             innermost_map(i, j) = max_around
+          if (count_diff == 0) then
+             extr_map(i, j) = 0
           else
-             min_around = minval(field(i - 1:i + 1, j - 1:j + 1), mask_center)
-
-             if (field(i, j) == max_around .or. field(i, j) == min_around) then
-                print *, "local_extrema: we require strict extrema"
-                print *, "field(i = ", i, ", j = ", j, ") = ", field(i, j)
-                print *, "max_around = ", max_around
-                print *, "min_around = ", min_around
-                stop 1
-             else if (field(i, j) < min_around) then
-                n_extr = n_extr + 1
-                extr_map(i, j) = - n_extr
-                innermost_map(i, j) = min_around
+             ! {count_diff >= 1}
+             max_around = maxval(field(i - 1:i + 1, j - 1:j + 1), diff_center)
+
+             if (field(i, j) > max_around) then
+                ! {local maximum, possibly degenerate}
+                if (count_diff <= 7 .and. any([extr_map(i - 1:i + 1, j - 1), &
+                     extr_map(i - 1, j)] > 0)) then
+                   ! {degenerate maximum and one of previously examined
+                   ! values is a maximum, so field(i, j) equals one of
+                   ! those previously examined values}
+                   extr_map(i, j) = 0
+                else
+                   n_extr = n_extr + 1
+                   extr_map(i, j) = n_extr
+                   innermost_map(i, j) = max_around
+                end if
              else
-                extr_map(i, j) = 0
+                ! {field(i, j) < max_around}
+                min_around = minval(field(i - 1:i + 1, j - 1:j + 1), &
+                     diff_center)
+
+                if (field(i, j) < min_around) then
+                   ! {local minimum, possibly degenerate}
+                   if (count_diff <= 7 .and. &
+                        any([extr_map(i - 1:i + 1, j - 1), &
+                        extr_map(i - 1, j)] < 0)) then
+                      ! {degenerate minimum and one of previously examined
+                      ! values is a minimum, so field(i, j) equals one of
+                      ! those previously examined values}
+                      extr_map(i, j) = 0
+                   else
+                      n_extr = n_extr + 1
+                      extr_map(i, j) = - n_extr
+                      innermost_map(i, j) = min_around
+                   end if
+                else
+                   ! {min_around < field(i, j) < max_around}
+                   extr_map(i, j) = 0
+                end if
              end if
           end if
        end do
@@ -89,12 +112,12 @@ contains
 
     i = 1
     j = 2
-    
+
     do k = 1, n_extr
        do
           ! next (i, j)
           i = i + 1
-          
+
           if (i == m) then
              ! Go to next column
              i = 2
@@ -122,13 +145,4 @@ contains
 
   end subroutine local_extrema
 
-  !***********************************************************************
-
-  subroutine construct_mask_center
-
-    mask_center = .true.
-    mask_center(2, 2) = .false.
-
-  end subroutine construct_mask_center
-
 end module local_extrema_m
diff --git a/Sources/set_outermost_contour.f b/Sources/set_outermost_contour.f
index 716a65ff..5b376455 100644
--- a/Sources/set_outermost_contour.f
+++ b/Sources/set_outermost_contour.f
@@ -5,7 +5,7 @@ module set_outermost_contour_m
 contains
 
   subroutine set_outermost_contour(e, ind_targ_extr, innermost_level, &
-       extr_map, ssh, corner, step, noise_around)
+       extr_map, ssh, corner, step)
 
     ! Defines e%outermost_contour. Method: dichotomy on level of ssh.
 
@@ -36,12 +36,10 @@ contains
     real, intent(in):: step(:) ! (2)
     ! longitude and latitude steps, in rad
 
-    logical, intent(in):: noise_around
-
     ! Local:
 
     integer i ! identifying number of the eddy
-    
+
     real, allocatable:: nearby_extr(:, :) ! (2, :) longitude and
     ! latitude, in rad, of all the significant extrema, except the
     ! target extremum
@@ -53,56 +51,52 @@ contains
     !-----------------------------------------------------------------
 
     i = extr_map(ind_targ_extr(1), ind_targ_extr(2))
-
-    if (noise_around) then
-       nearby_extr = convert_to_reg_coord(argwhere(extr_map > 0 &
-            .and. extr_map /= i), corner, step)
-       level_try = merge(maxval(ssh), minval(ssh), e%cyclone)
+    nearby_extr = convert_to_reg_coord(argwhere(extr_map > 0 &
+         .and. extr_map /= i), corner, step)
+    e%outermost_contour%polyline = good_contour(corner, step, ssh, &
+         innermost_level, e%coord_extr, nearby_extr)
+
+    if (e%outermost_contour%n_points == 0) then
+       e%outermost_contour%ssh = e%ssh_extremum
+       e%outermost_contour%area = 0.
     else
-       nearby_extr = convert_to_reg_coord(pack_indices(extr_map, &
-            excluded = [0, i]), corner, step)
-       level_try = outermost_possible_level(ssh, ind_targ_extr, e%cyclone)
-    end if
-
-    level_good = e%ssh_extremum &
-         + (innermost_level - e%ssh_extremum) * (1. - 100. * epsilon(0.))
-    e%outermost_contour%polyline = good_contour(corner, step, ssh, level_good, &
-         e%coord_extr, nearby_extr)
-    call assert(e%outermost_contour%n_points /= 0, &
-         "set_outermost_contour innermost_level")
-    call assert((level_try - e%ssh_extremum) / (level_good - e%ssh_extremum) &
-         > 1., "set_outermost_contour level_try")
-
-    tentative_contour = good_contour(corner, step, ssh, level_try, &
-         e%coord_extr, nearby_extr)
-
-    if (tentative_contour%n_points /= 0) then
-       e%outermost_contour%polyline = tentative_contour
-       level_good = level_try
-    else
-       ! {no good contour at level_try}
-       level_bad = level_try
-
-       do while (abs(level_bad - level_good) > precision)
-          level_try = (level_good + level_bad) / 2.
-          tentative_contour = good_contour(corner, step, ssh, level_try, &
-               e%coord_extr, nearby_extr)
-          if (tentative_contour%n_points /= 0) then
-             level_good = level_try
-             e%outermost_contour%polyline = tentative_contour
-          else
-             level_bad = level_try
-          end if
-
-          ! {e%outermost_contour%polyline == good_contour(. . ., level_good) and
-          ! e%outermost_contour%%n_points /= 0 and no good contour at level_bad}
-       end do
+       level_good = innermost_level
+       level_try = merge(maxval(ssh), minval(ssh), e%cyclone)
+       call assert((level_try - e%ssh_extremum) &
+            / (level_good - e%ssh_extremum) > 1., &
+            "set_outermost_contour level_try")
+       tentative_contour = good_contour(corner, step, ssh, level_try, &
+            e%coord_extr, nearby_extr)
+
+       if (tentative_contour%n_points /= 0) then
+          e%outermost_contour%polyline = tentative_contour
+          level_good = level_try
+       else
+          ! {no good contour at level_try}
+          level_bad = level_try
+
+          do while (abs(level_bad - level_good) > precision)
+             level_try = (level_good + level_bad) / 2.
+             tentative_contour = good_contour(corner, step, ssh, level_try, &
+                  e%coord_extr, nearby_extr)
+             if (tentative_contour%n_points /= 0) then
+                level_good = level_try
+                e%outermost_contour%polyline = tentative_contour
+             else
+                level_bad = level_try
+             end if
+
+             ! {e%outermost_contour%polyline == good_contour(. . .,
+             ! level_good) and e%outermost_contour%%n_points /= 0 and
+             ! no good contour at level_bad}
+          end do
+       end if
+
+       e%outermost_contour%ssh = level_good
+       e%outermost_contour%area &
+            = spherical_polygon_area(e%outermost_contour%polyline)
     end if
 
-    e%outermost_contour%ssh = level_good
-    e%outermost_contour%area &
-         = spherical_polygon_area(e%outermost_contour%polyline)
-
   end subroutine set_outermost_contour
 
 end module set_outermost_contour_m
diff --git a/Tests/Stdin/set_outermost_contour_nml.txt b/Tests/Stdin/set_outermost_contour_nml.txt
deleted file mode 100644
index 3933fc28..00000000
--- a/Tests/Stdin/set_outermost_contour_nml.txt
+++ /dev/null
@@ -1 +0,0 @@
-&main_nml NOISE_AROUND=t /
diff --git a/Tests/examine_ssh_values.py b/Tests/examine_ssh_values.py
new file mode 100755
index 00000000..48b50ef6
--- /dev/null
+++ b/Tests/examine_ssh_values.py
@@ -0,0 +1,27 @@
+#!/usr/bin/env python3
+
+import netCDF4
+import matplotlib.pyplot as plt
+import numpy as np
+
+with netCDF4.Dataset("h.nc") as f:
+    ssh = f.variables["adt"][:].squeeze()
+
+i_min = int(input("imin = ? "))
+i_max = int(input("imax = ? "))
+j_min = int(input("jmin = ? "))
+j_max = int(input("jmax = ? "))
+level = float(input("level = ? "))
+
+plt.figure()
+plt.contour(ssh, (level,))
+i_array = np.arange(i_min,i_max + 1)
+j_array = np.arange(j_min,j_max + 1)
+i_array_2d, j_array_2d = np.meshgrid(i_array, j_array)
+plt.plot(i_array_2d.reshape(-1), j_array_2d.reshape(-1), "+")
+
+for i in i_array:
+    for j in j_array:
+        plt.annotate(ssh[j,i], (i, j))
+
+plt.show()
diff --git a/Tests/plot_snapshot.py b/Tests/plot_snapshot.py
index 3ab4b701..826913ae 100755
--- a/Tests/plot_snapshot.py
+++ b/Tests/plot_snapshot.py
@@ -9,6 +9,7 @@ import shapefile
 import numpy as np
 import matplotlib.pyplot as plt
 import netCDF4
+import matplotlib.markers
 
 parser = argparse.ArgumentParser()
 parser.add_argument("-v", "--velocity", help = "plot velocity field",
@@ -37,14 +38,14 @@ for shape_rec in reader.iterShapeRecords():
     suff_amp_list.append(shape_rec.record[5])
 
 points = np.array(points)
-plt.plot(points[:, 0], points[:, 1], "o", markersize = 10, fillstyle = "none",
-         color = "black")
+colors = np.choose(cyclone, ("red", "blue"))
+
+plt.scatter(points[:, 0], points[:, 1], s = 100, edgecolors = colors,
+            facecolors = "none")
 
 for i, e in enumerate(points, start = 1):
     plt.annotate(str(i), e, xytext = (3, 3), textcoords = "offset points")
 
-colors = np.choose(cyclone, ("red", "blue"))
-
 # Outermost and max-speed contours:
 
 reader_outer = shapefile.Reader("outermost_contour_1")
@@ -53,20 +54,21 @@ reader_m_s = shapefile.Reader("max_speed_contour_1")
 for shape_outer, shape_m_s, my_color, suff_amp in zip(reader_outer.iterShapes(),
                                                       reader_m_s.iterShapes(),
                                                       colors, suff_amp_list):
-    points = np.array(shape_outer.points)
-    lines = plt.plot(points[:, 0], points[:, 1], color = my_color)
-
-    if suff_amp == 0:
-        lines[0].set_marker("s")
-        lines[0].set_fillstyle("none")
-    elif shape_m_s.shapeType == shapefile.NULL:
-        lines[0].set_marker("x")
-    else:
-        lines[0].set_marker("o")
-            
-    if shape_m_s.shapeType != shapefile.NULL:
-        points = np.array(shape_m_s.points)
-        plt.plot(points[:, 0], points[:, 1], color = my_color, marker = "o")
+    if shape_outer.shapeType != shapefile.NULL:
+        points = np.array(shape_outer.points)
+        lines = plt.plot(points[:, 0], points[:, 1], color = my_color)
+
+        if suff_amp == 0:
+            lines[0].set_marker("s")
+            lines[0].set_fillstyle("none")
+        elif shape_m_s.shapeType == shapefile.NULL:
+            lines[0].set_marker("x")
+        else:
+            lines[0].set_marker("o")
+
+        if shape_m_s.shapeType != shapefile.NULL:
+            points = np.array(shape_m_s.points)
+            plt.plot(points[:, 0], points[:, 1], color = my_color, marker = "o")
 
 if args.velocity:
     with netCDF4.Dataset("uv.nc") as f:
diff --git a/Tests/test_get_snapshot.f b/Tests/test_get_snapshot.f
index 6d544b1d..44a2ff35 100644
--- a/Tests/test_get_snapshot.f
+++ b/Tests/test_get_snapshot.f
@@ -3,7 +3,6 @@ program test_get_snapshot
   use derived_types, only: snapshot
   use dispatch_snapshot_m, only: dispatch_snapshot
   use get_snapshot_m, only: get_snapshot
-  use local_extrema_m, only: construct_mask_center
   use jumble, only: new_unit
   use netcdf, only: nf90_nowrite
   use netcdf95, only: nf95_open, find_coord, nf95_inquire_dimension, &
@@ -57,7 +56,6 @@ program test_get_snapshot
 
   call nf95_close(ncid)
 
-  call construct_mask_center
   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, &
diff --git a/Tests/test_local_extrema.f b/Tests/test_local_extrema.f
index e528add3..5cc2a1af 100644
--- a/Tests/test_local_extrema.f
+++ b/Tests/test_local_extrema.f
@@ -1,7 +1,7 @@
 program test_local_extrema
 
   use jumble, only: new_unit
-  use local_extrema_m, only: local_extrema, construct_mask_center
+  use local_extrema_m, only: 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, &
@@ -54,7 +54,6 @@ program test_local_extrema
 
   call nf95_close(ncid)
 
-  call construct_mask_center
   call local_extrema(ssh, extr_map, ind_extr, innermost_level, cyclone)
 
   call new_unit(unit)
diff --git a/Tests/test_local_extrema.py b/Tests/test_local_extrema.py
index 00332c6f..6bfd67fd 100755
--- a/Tests/test_local_extrema.py
+++ b/Tests/test_local_extrema.py
@@ -30,23 +30,22 @@ my_min[:,0] = False
 my_min[:,-1] = False
 
 ilat_max, ilon_max = np.where(my_max)
-print("Maxima (longitude index, latitude index) (1-based):")
-print(np.column_stack((ilon_max + 1, ilat_max + 1)))
-print("Maxima (longitude, latitude):")
-print(np.column_stack((longitude[ilon_max], latitude[ilat_max])))
-
 ilat_min, ilon_min = np.where(my_min)
-print("Minima (longitude index, latitude index) (1-based):")
-print(np.column_stack((ilon_min + 1, ilat_min + 1)))
-print("Minima:")
-print(np.column_stack((longitude[ilon_min], latitude[ilat_min])))
 
-print("Values around maxima:")
+print("Maxima:")
+print("(longitude index, latitude index) (1-based) -- (longitude, latitude):")
+print("Values around maxima\n")
 for j,i in np.argwhere(my_max):
+    print("({}, {}) -- ({}, {})".format(i + 1, j + 1, longitude[i],
+                                        latitude[j]))
     print(ssh[j-1:j+2, i-1:i+2], "\n")
 
-print("Values around minima:")
+print("Minima:")
+print("(longitude index, latitude index) (1-based) -- (longitude, latitude):")
+print("Values around minima\n")
 for j,i in np.argwhere(my_min):
+    print("({}, {}) -- ({}, {})".format(i + 1, j + 1, longitude[i],
+                                        latitude[j]))
     print(ssh[j-1:j+2, i-1:i+2], "\n")
 
 f.close()
diff --git a/Tests/test_set_outermost_contour.f b/Tests/test_set_outermost_contour.f
index 7785b7f9..db387dcc 100644
--- a/Tests/test_set_outermost_contour.f
+++ b/Tests/test_set_outermost_contour.f
@@ -34,10 +34,9 @@ program test_set_outermost_contour
   real, parameter:: step = 0.25 / 180. * pi ! in rad
   TYPE(shpfileobject) shphandle
   integer field_number, shape_number
-  logical:: noise_around = .false.
 
   namelist /main_nml/ ilon_llc, ilat_llc, ilon_urc, ilat_urc, ind_targ_extr, &
-       innermost_level, cyclone, noise_around
+       innermost_level, cyclone
 
   !----------------------------------------------------------------
 
@@ -82,8 +81,7 @@ program test_set_outermost_contour
   e%coord_extr = [longitude(ind_targ_extr(1)), latitude(ind_targ_extr(2))]
   e%cyclone = cyclone
   call set_outermost_contour(e, ind_targ_extr, innermost_level, extr_map, &
-       ssh, corner = [longitude(1), latitude(1)], step = [step, step], &
-       noise_around = noise_around)
+       ssh, corner = [longitude(1), latitude(1)], step = [step, step])
 
   print *, "e%outermost_contour%closed = ", e%outermost_contour%closed
   print *, "Radius of disk of equal area: ", &
diff --git a/Tests/tests.json b/Tests/tests.json
index f46c1bd5..d920bd6c 100644
--- a/Tests/tests.json
+++ b/Tests/tests.json
@@ -46,7 +46,7 @@
 		  "$input_dir/extr_map_negative_2_8.nc"],
 	"stdout" : "test_set_outermost_contour_stdout.txt",
 	"directory" : "Tests_new/Set_outermost_contour_noise_2_8",
-	"stdin" : "$stdin_dir/set_outermost_contour_nml.txt"
+	"stdin" : "$stdin_dir/default_main_nml.txt"
     },
     {
 	"args" : ["$compil_prod_dir/test_set_outermost_contour",
@@ -54,7 +54,7 @@
 		  "$input_dir/extr_map_negative_2.nc"],
 	"stdout" : "test_set_outermost_contour_stdout.txt",
 	"directory" : "Tests_new/Set_outermost_contour_noise_2",
-	"stdin" : "$stdin_dir/set_outermost_contour_nml.txt"
+	"stdin" : "$stdin_dir/default_main_nml.txt"
     },
     {
 	"args" : ["$compil_prod_dir/test_max_speed_contour_ssh",
-- 
GitLab