Keep the good contours

Keep the good contours found when looking for the outermost contour,
subject to a maximum of `n_max_cont` contours. The idea is to use them
later when looking for the maximum-speed contour. Inside function
`get_1_outerm`, we temporarily store the good contours in a local
explicit-shape array, `cont_list`.

Inst_eddies/get_1_outerm.f90

@@ -9,10 +9,20 @@ contains
 
     ! This procedure gets one outermost contour, if it can find
     ! one. If the procedure cannot find an outermost contour, it
-    ! returns a null ssh_contour. An outermost contour is not found
-    ! if, and only if, there is no good contour at innermost level. If
-    ! a contour is found then its level is farther from the extremum
-    ! than innermost level.
+    ! returns a null ssh_contour. The procedure also finds some inner
+    ! contours. The number of contours stored, including the outermost
+    ! contour, is lower than or equal to n_max_cont. If we find more
+    ! than n_max_cont -1 good contours excluding the outermost
+    ! contour, then we store only the first n_max_cont -1 found,
+    ! because they provide a better sample of the ssh range. If the
+    ! procedure cannot find an outermost contour, it does not store
+    ! inner contours. An outermost contour is not found if, and only
+    ! if, there is no good contour at innermost level. If an outermost
+    ! contour is found then its level is farther from the extremum
+    ! than innermost level. In the array of stored contours, the
+    ! contours are in monotonic order of ssh: ascending if the
+    ! extremum is a minimum, descending if the extremum is a
+    ! maximum. The area is only computed for the outermost contour.
 
     ! Method: dichotomy on level of ssh.
 
@@ -29,7 +39,7 @@ contains
     use jumble, only: assert
 
     use ccw_orient_m, only: ccw_orient
-    use derived_types, only: ssh_contour, missing_ssh
+    use derived_types, only: ssh_contour, null_ssh_contour
     use good_contour_m, only: good_contour
     use spher_polyline_area_m, only: spher_polyline_area
 
@@ -55,60 +65,74 @@ contains
 
     ! Local:
     real level_try, level_good, level_bad ! in m
-    type(polyline) tentative_contour
     real, parameter:: precision = 5e-4 ! in m
     logical mask(size(ssh, 1), size(ssh, 2))
+    integer, parameter:: n_max_cont = 31
+
+    type(ssh_contour) cont_list(n_max_cont + 1)
+    ! Contour list. The outermost contour is element n_cont. The
+    ! contours are in monotonic order of ssh.
+
+    integer n_cont ! number of good contours found and stored
 
     !-----------------------------------------------------------------
 
-    get_1_outerm%polyline = good_contour(corner, step, ssh, innermost_level_2, &
+    n_cont = 0
+    cont_list(1)%polyline = good_contour(corner, step, ssh, innermost_level_2, &
          coord_extr, outside_points)
 
-    if (get_1_outerm%n_points == 0) then
-       ! null ssh contour
-       get_1_outerm%ssh = missing_ssh
-       get_1_outerm%area = - 1e6
+    if (cont_list(1)%n_points == 0) then
+       get_1_outerm = null_ssh_contour()
     else
+       n_cont = 1
+       cont_list(1)%ssh = innermost_level_2
        level_good = innermost_level_2
        mask = ssh /= huge(0.)
        level_try = merge(maxval(ssh, mask), minval(ssh, mask), cyclone)
        call assert(merge(level_try > level_good, level_try < level_good, &
             cyclone), "get_1_outerm level_try")
-       tentative_contour = good_contour(corner, step, ssh, level_try, &
+       cont_list(2)%polyline = good_contour(corner, step, ssh, level_try, &
             coord_extr, outside_points)
 
-       if (tentative_contour%n_points /= 0) then
+       if (cont_list(2)%n_points /= 0) then
           ! This can only happen if there are degenerate extrema
           ! everywhere around the path of contour. So it should
           ! probably never happen.
-          get_1_outerm%polyline = tentative_contour
-          get_1_outerm%ssh = level_try
+          n_cont = 2
+          cont_list(2)%ssh = level_try
        else
           ! Assertion: {no good contour at maxval or minval}
           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, &
-                  coord_extr, outside_points)
+             cont_list(n_cont + 1)%polyline = good_contour(corner, step, ssh, &
+                  level_try, coord_extr, outside_points)
+
+             if (cont_list(n_cont + 1)%n_points /= 0) then
+                cont_list(n_cont + 1)%ssh = level_try
+
+                if (n_cont == n_max_cont) then
+                   ! Replace the previous good contour found
+                   cont_list(n_cont) = cont_list(n_cont + 1)
+                else
+                   n_cont = n_cont + 1
+                end if
 
-             if (tentative_contour%n_points /= 0) then
                 level_good = level_try
-                get_1_outerm%polyline = tentative_contour
              else
                 level_bad = level_try
              end if
 
-             ! Assertion: get_1_outerm%polyline == good_contour(. . .,
-             ! level_good) and get_1_outerm%%n_points /= 0 and
+             ! Assertion: cont_list(n_cont)%polyline == good_contour(. . .,
+             ! level_good) and cont_list(n_cont)%%n_points /= 0 and
              ! no good contour at level_bad
           end do
-
-          get_1_outerm%ssh = level_good
        end if
 
-       get_1_outerm%area = spher_polyline_area(get_1_outerm%polyline)
-       call ccw_orient(get_1_outerm)
+       cont_list(n_cont)%area = spher_polyline_area(cont_list(n_cont)%polyline)
+       call ccw_orient(cont_list(n_cont))
+       get_1_outerm = cont_list(n_cont)
     end if
 
   end function get_1_outerm