From 94515cc5b5fb96d761132c92f060704d0921818f Mon Sep 17 00:00:00 2001
From: Lionel GUEZ <guez@lmd.ens.fr>
Date: Fri, 30 Mar 2018 00:41:43 +0200
Subject: [PATCH] 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 |  5 ++---
 derived_types.f           | 16 ++++++++++++++++
 get_1_outerm.f            | 20 +++++++++++---------
 get_snapshot.f            |  6 ++----
 set_all_outerm.f          | 14 +++++++-------
 set_max_speed.f           | 10 +++-------
 6 files changed, 41 insertions(+), 30 deletions(-)

diff --git a/Tests/test_get_1_outerm.f b/Tests/test_get_1_outerm.f
index 2adec744..7fdde880 100644
--- a/Tests/test_get_1_outerm.f
+++ b/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, &
diff --git a/derived_types.f b/derived_types.f
index dc3f47f3..073deec2 100644
--- a/derived_types.f
+++ b/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
diff --git a/get_1_outerm.f b/get_1_outerm.f
index cbb2bd2a..09e66bcc 100644
--- a/get_1_outerm.f
+++ b/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)
 
diff --git a/get_snapshot.f b/get_snapshot.f
index a4798165..032775bb 100644
--- a/get_snapshot.f
+++ b/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
diff --git a/set_all_outerm.f b/set_all_outerm.f
index 6a3aa92d..c65035bc 100644
--- a/set_all_outerm.f
+++ b/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)
diff --git a/set_max_speed.f b/set_max_speed.f
index e4c10ba0..135362cf 100644
--- a/set_max_speed.f
+++ b/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.
-- 
GitLab