In procedure local_extrema, add possibility of periodicity in the

first dimension of input field. So there may be extrema on the borders
and we have to create a private procedure local_extremum.

Added scale option to script plot_snapshot.py.

Documentation_texfol/Graphiques/.gitignore

@@ -5,3 +5,4 @@ m3.pdf
 processes.pdf
 regions.pdf
 window.pdf
+periodicity.pdf

Documentation_texfol/Graphiques/GNUmakefile

-sources = window processes m3 call_graph input_output degeneracy
+sources = window processes m3 call_graph input_output degeneracy periodicity
 
 objects := $(addsuffix .pdf, ${sources})
 

Documentation_texfol/documentation.tex

@@ -1095,6 +1095,20 @@ de départ. \verb+scipy.ndimage.maximum_filter+ applique
 \verb+scipy.ndimage.maximum_filter1d+ dans chaque dimension. Cf. aussi
 Douglas (1996 k0968).
 
+Cf. figure (\ref{fig:periodicity}).
+\begin{figure}[htbp]
+  \centering
+  \includegraphics{periodicity}
+  \caption{Cas de périodicité dans local\_extrema. Le tableau
+    représenté est extr\_map. Pour définir la case (1, j), il faut
+    avoir copié la case (m, j - 1) (qui peut être non nulle) dans
+    local\_map. Pour définir la case (m, j), il faut avoir copié les
+    cases (1, j - 1) et (1, j) (dont l'une peut être non nulle) dans
+    local\_map.}
+  \label{fig:periodicity}
+\end{figure}
+Cas de périodicité : field(i + m, :) = field(i, :).
+
 \subsection{get\_1\_outerm}
 
 Au lieu d'initialiser l'amplitude maximale à 4 m plus du bruit, comme
@@ -1268,4 +1282,13 @@ comparant les fichiers outermost\_contour\_1.dbf) : dans le cas avec
 minimum, ils englobent un extremum d'amplitude insuffisante. Dans le
 cas avec minimum, 48 contours extérieurs sont calculés deux fois.
 
+Sur la région 3, avec un mm de minimum d'amplitude, le contour
+extérieur 168 a un diamètre de 27 points, environ 7 degrés en
+longitude. \`A 40° S, donc un diamètre de 600 km environ. On ne peut
+donc pas prendre max\_radius = 12. Ayant fait le test avec max\_radius
+= 12, je constate aussi que le temps d'exécution est très sensible à
+max\_radius : un facteur \np{2.7} environ entre max\_radius = 12 et
+max\_radius = 20. Le temps d'exécution semble donc proportionnel à
+max\_radius$^2$, ce qui peut se comprendre.
+
 \end{document}

Tests/plot_snapshot.py

@@ -21,6 +21,8 @@ import itertools
 parser = argparse.ArgumentParser()
 parser.add_argument("-v", "--velocity", help = "plot velocity field",
                     action = "store_true")
+parser.add_argument("-s", "--scale", default = 20,
+                    help = "scale of arrows for the velocity field")
 args = parser.parse_args()
 
 with netCDF4.Dataset("h.nc") as f:
@@ -93,8 +95,9 @@ if args.velocity:
                                                latitude)
         
     quiver_return = m.quiver(lon_2d, lat_2d, u_rotated, v_rotated,
-                             latlon = True, scale = 5, scale_units = "width")
-    plt.quiverkey(quiver_return, 0.75, 0.9, 0.1, r"0.1 m s$^{-1}$",
+                             latlon = True, scale = args.scale,
+                             scale_units = "width")
+    plt.quiverkey(quiver_return, 0.75, 0.9, 1, r"1 m s$^{-1}$",
               coordinates = "figure")
 
 m.drawmeridians(np.linspace(llcrnrlon, urcrnrlon, 6), labels = [0, 0, 0, 1])

Tests/test_local_extrema.f

@@ -3,21 +3,21 @@ program test_local_extrema
   ! Reads ADT from a NetCDF file and writes the list of extrema to a
   ! CSV file. Also creates a NetCDF file containing the map of extrema.
 
-  use jumble, only: new_unit
-  use local_extrema_m, only: local_extrema
+  ! Libraries:
+  use jumble, only: new_unit, get_command_arg_dyn
   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_get_att
-  use nr_util, only: assert
+
+  use local_extrema_m, only: local_extrema
 
   implicit none
 
   integer nlon, nlat
-
   integer ncid, varid, dimid_lat, dimid_lon
   integer varid_lat, varid_lon, varid_extr_map
-  integer i, length, status, unit
+  integer i, unit
   character(len = :), allocatable:: filename
   real, allocatable:: longitude(:) ! (nlon) in degrees
   real, allocatable:: latitude(:) ! (nlat) in degrees
@@ -32,14 +32,14 @@ program test_local_extrema
 
   logical, allocatable:: cyclone(:) ! (n_extr)
   real Fill_Value
+  logical periodic ! in longitude
 
   !---------------------------------------------------------------------
 
-  call get_command_argument(number = 1, length = length, status = status)
-  call assert(status == 0, "Required argument: ADT file")
-  allocate(character(len = length):: filename)
-  call get_command_argument(1, filename)
-
+  call get_command_arg_dyn(1, filename, message = "Required argument: ADT file")
+  print *, "periodic = ? "
+  read *, periodic
+  print *, "periodic = ", periodic
   print *, "Reading from ", filename, "..."
   call nf95_open(filename, nf90_nowrite, ncid)
 
@@ -60,7 +60,8 @@ program test_local_extrema
 
   call nf95_close(ncid)
 
-  call local_extrema(ssh, extr_map, ind_extr, innermost_level, cyclone)
+  call local_extrema(extr_map, ind_extr, innermost_level, cyclone, ssh, &
+       periodic)
 
   call new_unit(unit)
   open(unit, file = "test_local_extrema.csv", status = "replace", &

Tests/tests.json

@@ -25,12 +25,14 @@
     {
 	"args" : ["$compil_prod_dir/test_local_extrema",
 		  "$input_dir/h_region_1.nc"],
-	"title" : "Local_extrema"
+	"title" : "Local_extrema",
+	"input" : "f"
     },
     {
 	"args" : ["$compil_prod_dir/test_local_extrema",
 		  "$input_dir/h_region_2.nc"],
-	"title" : "Local_extrema_larger"
+	"title" : "Local_extrema_larger",
+	"input" : "f"
     },
     {
 	"args" : "$compil_prod_dir/test_get_1_outerm",
@@ -289,7 +291,8 @@
 	"args" : ["$compil_prod_dir/test_local_extrema",
 		  "$input_dir/Region_4/h.nc"],
 	"title" : "Local_extrema_missing",
-	"description": "test_local_extrema with input file containing missing values."
+	"description": "test_local_extrema with input file containing missing values.",
+	"input" : "f"
     },
     {
 	"args" : "$compil_prod_dir/test_successive_overlap",
@@ -337,5 +340,12 @@
 	    "$tests_old_dir/Get_snapshot_region_5/outermost_contour_1.*"
 	],
 	"stdin": "$input_dir/successive_overlap_region_5_nml.txt"
+    },
+    {
+	"args" : ["$compil_prod_dir/test_local_extrema",
+		  "$input_dir/h_region_1.nc"],
+	"title" : "Local_extrema_periodic",
+	"input" : "t",
+	"description": "Same as Local_extrema but with periodicity."
     }
 ]

local_extrema.f

@@ -2,22 +2,23 @@ module local_extrema_m
 
   implicit none
 
+  private local_extremum
+
 contains
 
-  subroutine local_extrema(field, extr_map, ind_extr, innermost_level, &
-       local_min)
+  subroutine local_extrema(extr_map, ind_extr, innermost_level, local_min, &
+       field, periodic)
 
     use nr_util, only: assert_eq
 
-    real, intent(in):: field(:, :) ! (m, n) Missing value is huge(0.).
-
     integer, intent(out):: extr_map(:, :) ! (m, n) Map of
     ! identification numbers of extrema. Identification numbers are
     ! indices in the second dimension of ind_extr. During a call to
     ! local_extrema, we use the array extr_map to temporarily record
-    ! whether the extremum is a maximum or a minimum: extr_map_map is
-    ! positive for a maximum, negative for a minimum. But, at the end
-    ! of the call, all values of extr_map are set positive.
+    ! whether the extremum is a maximum or a minimum: extr_map is
+    ! positive for a maximum, negative for a minimum, 0
+    ! elsewhere. But, at the end of the call, all values of extr_map
+    ! are set positive.
 
     integer, intent(out), allocatable:: ind_extr(:, :) ! (2, n_extr)
     ! indices in the two dimensions of extrema of field
@@ -28,94 +29,79 @@ contains
     ! minimum.
 
     logical, intent(out), allocatable:: local_min(:) ! (n_extr)
+    real, intent(in):: field(:, :) ! (m, n) Missing value is huge(0.).
+    logical, intent(in):: periodic ! field is periodic in the first dimension
 
     ! Local:
-    integer n_extr, i, j, k, m, n, count_diff
+
+    integer n_extr, i, j, k, m, n
 
     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
+    integer local_map(- 1:1, - 1:1)
+    real local_field(- 1:1, - 1:1)
 
     !--------------------------------------------------------------------
 
     m = assert_eq(size(field, 1), size(extr_map, 1), "local_extrema m")
     n = assert_eq(size(field, 2), size(extr_map, 2), "local_extrema n")
 
-    ! No extremum on the borders of field:
-    extr_map(1, 2:n - 1) = 0
-    extr_map(m, 2:n - 1) = 0
-    extr_map(:, 1) = 0
-    extr_map(:, n) = 0
-
+    extr_map = 0
     n_extr = 0
 
-    do j = 2, n - 1
-       do i = 2, m - 1
-          if (any(field(i - 1:i + 1, j - 1:j + 1) == huge(0.))) then
-             extr_map(i, j) = 0
-          else
-             diff_center = field(i - 1:i + 1, j - 1:j + 1) /= field(i, j)
-             count_diff = count(diff_center)
-
-             if (count_diff == 0) then
-                extr_map(i, j) = 0
-             else
-                ! {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
-                   ! {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 if
+    ! No extremum on the borders of the second dimension of field so
+    ! loop on columns 2 to n - 1.
+    
+    if (periodic) then
+       do j = 2, n - 1
+          ! Row 1:
+          
+          local_map(- 1, :) = extr_map(m, j - 1:j + 1)
+          local_map(0:, :) = extr_map(:2, j - 1:j + 1)
+
+          local_field(- 1, :) = field(m, j - 1:j + 1)
+          local_field(0:, :) = field(:2, j - 1:j + 1)
+
+          call local_extremum(innermost_map(1, j), local_map, n_extr, &
+               local_field)
+          extr_map(1, j) = local_map(0, 0)
+
+          do i = 2, m - 1
+             call local_extremum(innermost_map(i, j), &
+                  extr_map(i - 1:i + 1, j - 1:j + 1), n_extr, &
+                  field(i - 1:i + 1, j - 1:j + 1))
+          end do
+
+          ! Row m:
+          
+          local_map(:0, :) = extr_map(m - 1:, j - 1:j + 1)
+          local_map(1, :) = extr_map(1, j - 1:j + 1)
+
+          local_field(:0, :) = field(m - 1:, j - 1:j + 1)
+          local_field(1, :) = field(1, j - 1:j + 1)
+
+          call local_extremum(innermost_map(m, j), local_map, n_extr, &
+               local_field)
+          extr_map(m, j) = local_map(0, 0)
        end do
-    end do
+    else
+       do j = 2, n - 1
+          do i = 2, m - 1
+             call local_extremum(innermost_map(i, j), &
+                  extr_map(i - 1:i + 1, j - 1:j + 1), n_extr, &
+                  field(i - 1:i + 1, j - 1:j + 1))
+          end do
+       end do
+    end if
 
-    ! We do not call pack_indices to define ind_extr because we
-    ! already know the number of extrema so the following algorithm is
-    ! faster (stops at n_extr) and uses less memory.
+    ! Define ind_extr. We do not call pack_indices because we already
+    ! know the number of extrema so the following algorithm is faster
+    ! (stops at n_extr) and uses less memory.
 
     allocate(ind_extr(2, n_extr))
 
-    i = 1
+    i = 0
     j = 2
 
     do k = 1, n_extr
@@ -123,9 +109,9 @@ contains
           ! next (i, j)
           i = i + 1
 
-          if (i == m) then
+          if (i == m + 1) then
              ! Go to next column
-             i = 2
+             i = 1
              j = j + 1
           end if
 
@@ -150,4 +136,64 @@ contains
 
   end subroutine local_extrema
 
+  !*********************************************************************
+
+  subroutine local_extremum(extr_around, local_map, n_extr, local_field)
+
+    ! Look for an extremum at the center of a 3x3 domain.
+
+    real, intent(out):: extr_around
+    ! Value of extremum different from the center value. May not be
+    ! defined, but guaranteed to be defined if the center is a local
+    ! extremum.
+
+    integer, intent(inout):: local_map(- 1:, - 1:) ! (3, 3) Map of
+    ! identification numbers of extrema. "local_map" is positive for a
+    ! maximum, negative for a minimum.  Only local_map(0, 0) is
+    ! possibly updated.
+
+    integer, intent(inout):: n_extr
+
+    real, intent(in):: local_field(- 1:, - 1:) ! (3, 3) Missing value
+    ! is huge(0.).
+
+    ! Local:
+    logical diff_center(3, 3) ! values different from the center value
+    integer count_diff
+
+    !--------------------------------------------------------------------
+
+    if (all(local_field /= huge(0.))) then
+       diff_center = local_field /= local_field(0, 0)
+       count_diff = count(diff_center) ! {between 0 and 8}
+
+       if (count_diff /= 0) then
+          ! {count_diff >= 1}
+          extr_around = maxval(local_field, diff_center)
+
+          if (local_field(0, 0) > extr_around) then
+             ! {local maximum, possibly degenerate}
+             if (count_diff == 8 .or. all(local_map <= 0)) then
+                ! {not a degenerate maximum or none of the previously
+                ! examined values is a maximum}
+                n_extr = n_extr + 1
+                local_map(0, 0) = n_extr
+             end if
+          else
+             ! {local_field(0, 0) < maxval(local_field, diff_center)}
+             extr_around = minval(local_field, diff_center)
+
+             if (local_field(0, 0) < extr_around &
+                  .and. (count_diff == 8 .or. all(local_map >= 0))) then
+                ! {local minimum, not degenerate or none of previously
+                ! examined values is a minimum}
+                n_extr = n_extr + 1
+                local_map(0, 0) = - n_extr
+             end if
+          end if
+       end if
+    end if
+
+  end subroutine local_extremum
+
 end module local_extrema_m

set_all_outerm.f

@@ -68,7 +68,8 @@ contains
     nlon = size(ssh, 1)
     nlat = size(ssh, 2)
     allocate(s%extr_map(nlon, nlat))
-    call local_extrema(ssh, s%extr_map, s%ind_extr, innermost_level, cyclone)
+    call local_extrema(s%extr_map, s%ind_extr, innermost_level, cyclone, ssh, &
+         periodic = .false.)
     s%number_vis_eddies = size(s%ind_extr, 2)
     allocate(s%list_vis(s%number_vis_eddies), sorted_extr(s%number_vis_eddies))