When we do not find an outermost contour or a max-speed contour,

instead of defining the component ssh to the ssh of the extremum or to
the ssh of the outermost contour, use a missing value flag
(1e4). Also, when radius4 >= 2 and the max speed contour is the
outermost contour, use missing value flag in the ssh component of max
speed contour instead of ssh on outermost contour. This makes clearer
output.

Add function null_ssh_contour.

In procedure get_1_outerm, simplify the comparison of level_try and
level_good. So we do not need any longer the dummy argument
ssh_extremum.

Tests/test_get_1_outerm.f

@@ -29,9 +29,8 @@ program test_get_1_outerm
   type(ssh_contour) outermost_contour
   TYPE(shpfileobject) shphandle
   integer field_number, shape_number
-  real:: ssh_extremum = 0.2759
 
-  namelist /main_nml/ ind_targ_extr, ssh_extremum, innermost_level, cyclone
+  namelist /main_nml/ ind_targ_extr, innermost_level, cyclone
 
   !----------------------------------------------------------------
 
@@ -70,7 +69,7 @@ program test_get_1_outerm
   call nf95_close(ncid)
 
   step = [longitude(2) - longitude(1), latitude(2) - latitude(1)]
-  outermost_contour = get_1_outerm(ssh_extremum, cyclone, &
+  outermost_contour = get_1_outerm(cyclone, &
        coord_extr = [longitude(1), latitude(1)] + (ind_targ_extr - 1) * step, &
        i = extr_map(ind_targ_extr(1), ind_targ_extr(2)), &
        innermost_level = innermost_level, extr_map = extr_map, ssh = ssh, &

derived_types.f

@@ -58,4 +58,20 @@ module derived_types
      integer number_eddies ! number of eddies, including interpolated eddies
   end type snapshot
 
+  real, parameter:: missing_ssh = 1e4 ! flag for undefined contour
+
+contains
+
+  pure type(ssh_contour) function null_ssh_contour()
+
+    use contour_531, only: null_polyline
+
+    !--------------------------------------------------------------
+
+    null_ssh_contour%polyline = null_polyline()
+    null_ssh_contour%ssh = missing_ssh
+    null_ssh_contour%area = 0.
+
+  end function null_ssh_contour
+
 end module derived_types

get_1_outerm.f

@@ -4,24 +4,23 @@ module get_1_outerm_m
 
 contains
 
-  type(ssh_contour) function get_1_outerm(ssh_extremum, cyclone, coord_extr, &
-       i, innermost_level, extr_map, ssh, corner, step)
+  type(ssh_contour) function get_1_outerm(cyclone, coord_extr, i, &
+       innermost_level, extr_map, ssh, corner, step)
 
     ! Gets one outermost contour, if it can find one. Method:
     ! dichotomy on level of ssh. If the procedure cannot find an
     ! outermost contour, it return a null polyline, zero area and ssh
-    ! equal to ssh extremum. An outermost contour is not found if, and
+    ! equal to missing_ssh. An outermost contour is not found if, and
     ! only if, there is no good contour at innermost level.
 
     use contour_531, only: polyline, convert_to_reg_coord
-    use derived_types, only: ssh_contour
+    use derived_types, only: ssh_contour, missing_ssh
     use good_contour_m, only: good_contour
     use jumble, only: argwhere, pack_indices
     use nr_util, only: assert
     use outermost_possible_level_m, only: outermost_possible_level
     use spherical_polygon_area_m, only: spherical_polygon_area
 
-    real, intent(in):: ssh_extremum
     logical, intent(in):: cyclone
     real, intent(in):: coord_extr(:) ! (2)
 
@@ -30,7 +29,9 @@ contains
     ! at coord_extr
 
     real, intent(in):: innermost_level
-    ! ssh level of the innermost contour around the target extremum, in m
+    ! ssh level of the innermost contour around the target extremum,
+    ! in m. By construction, innermost_level < extremum for a maximum,
+    ! > extremum for a minimum.
 
     integer, intent(in):: extr_map(:, :)
     real, intent(in):: ssh(:, :) ! in m
@@ -60,14 +61,15 @@ contains
          coord_extr, nearby_extr)
 
     if (get_1_outerm%n_points == 0) then
-       get_1_outerm%ssh = ssh_extremum
+       ! null ssh contour
+       get_1_outerm%ssh = missing_ssh
        get_1_outerm%area = 0.
     else
        level_good = innermost_level
        mask = ssh /= huge(0.)
        level_try = merge(maxval(ssh, mask), minval(ssh, mask), cyclone)
-       call assert((level_try - ssh_extremum) &
-            / (level_good - ssh_extremum) > 1., "get_1_outerm level_try")
+       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, &
             coord_extr, nearby_extr)
 

get_snapshot.f

@@ -10,7 +10,7 @@ contains
     use, intrinsic:: ieee_arithmetic, only: ieee_value, IEEE_QUIET_NAN
 
     use contour_531, only: convert_to_ind, null_polyline
-    use derived_types, only: snapshot
+    use derived_types, only: snapshot, null_ssh_contour
     use netcdf, only: nf90_nowrite
     use netcdf95, only: nf95_open, nf95_close, nf95_inq_varid, nf95_get_var, &
          nf95_get_att
@@ -102,9 +102,7 @@ contains
                   v(llc(1):urc(1), llc(2):urc(2)), corner + (llc - 1) * step, &
                   step)
           else
-             s%list_vis(i)%speed_cont%polyline = null_polyline()
-             s%list_vis(i)%speed_cont%ssh = s%list_vis(i)%ssh_extr
-             s%list_vis(i)%speed_cont%area = 0.
+             s%list_vis(i)%speed_cont = null_ssh_contour()
              s%list_vis(i)%max_speed = 0.
              s%list_vis(i)%radius4 = 0
           end if

set_all_outerm.f

@@ -34,7 +34,9 @@ contains
     ! Local:
 
     real, allocatable:: innermost_level(:) ! (s%number_vis_eddies)
-    ! level of innermost contour, for each extremum
+    ! ssh level of the innermost contour around each extremum, in
+    ! m. By construction, innermost_level < extremum for a maximum, >
+    ! extremum for a minimum.
 
     logical, allocatable:: cyclone(:) ! (s%number_vis_eddies)
     integer i
@@ -84,9 +86,8 @@ contains
 
     do i = 1, s%number_vis_eddies
        if (flat_extr(i)) then
-          s%list_vis(i)%out_cont = get_1_outerm(s%list_vis(i)%ssh_extr, &
-               s%list_vis(i)%cyclone, s%list_vis(i)%coord_extr, i, &
-               innermost_level(i), &
+          s%list_vis(i)%out_cont = get_1_outerm(s%list_vis(i)%cyclone, &
+               s%list_vis(i)%coord_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)
@@ -122,9 +123,8 @@ contains
     do i = 1, s%number_vis_eddies
        if (s%list_vis(i)%suff_amp .and. noise_around(i) &
             .or. .not. flat_extr(i)) then
-          s%list_vis(i)%out_cont = get_1_outerm(s%list_vis(i)%ssh_extr, &
-               s%list_vis(i)%cyclone, s%list_vis(i)%coord_extr, i, &
-               innermost_level(i), &
+          s%list_vis(i)%out_cont = get_1_outerm(s%list_vis(i)%cyclone, &
+               s%list_vis(i)%coord_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)

set_max_speed.f

@@ -11,7 +11,7 @@ contains
     use, intrinsic:: IEEE_ARITHMETIC, only: IEEE_IS_NAN
 
     use contour_531, only: convert_to_reg_coord, convert_to_ind, null_polyline
-    use derived_types, only: eddy
+    use derived_types, only: eddy, null_ssh_contour
     use geometry, only: polygon_point_dist_2d
     use good_contour_m, only: good_contour
     use inside_4_m, only: inside_4
@@ -65,9 +65,7 @@ contains
 
     if (e%radius4 == 1) then
        ! We cannot look for a maximum speed contour other than out_cont.
-       e%speed_cont%ssh = e%ssh_extr
-       e%speed_cont%polyline = null_polyline()
-       e%speed_cont%area = 0.
+       e%speed_cont = null_ssh_contour()
 
        e%max_speed = speed_outerm
        ! (We will maybe need this to compute the weight of edges in
@@ -97,9 +95,7 @@ contains
 
        if (abs(speed_outerm) > abs(e%max_speed)) then
           ! Abandon the contour coming from max_speed_contour_ssh:
-          e%speed_cont%ssh = e%out_cont%ssh
-          e%speed_cont%polyline = null_polyline()
-          e%speed_cont%area = e%out_cont%area
+          e%speed_cont = null_ssh_contour()
           e%max_speed = speed_outerm
        else
           ! Stick to the contour coming from max_speed_contour_ssh.