From 0e1f9399238dc537bd6ac466e773414d77dc1108 Mon Sep 17 00:00:00 2001
From: Lionel GUEZ <guez@lmd.ens.fr>
Date: Tue, 31 Jan 2023 21:56:30 +0100
Subject: [PATCH] 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 | 74 ++++++++++++++++++++++++------------
1 file changed, 49 insertions(+), 25 deletions(-)
diff --git a/Inst_eddies/get_1_outerm.f90 b/Inst_eddies/get_1_outerm.f90
index b2a79f38..bcd01fc2 100644
--- a/Inst_eddies/get_1_outerm.f90
+++ b/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
--
GitLab