From cc13fd5ba5a0c0c2f2d26a19061c0798ba75a7d9 Mon Sep 17 00:00:00 2001
From: Lionel GUEZ <guez@lmd.ipsl.fr>
Date: Mon, 1 Jul 2024 18:53:28 +0200
Subject: [PATCH] Make `cont_list` a polyline
Make `cont_list` a polyline and `cont_list_proj` an SSH contour,
instead of vice-versa. The motivation is clarity in `get_1_outerm`: we
only define `cont_list` at the end.
---
Inst_eddies/Tests/test_get_1_outerm.f90 | 29 +++++------
Inst_eddies/Tests/test_set_max_speed.f90 | 33 +++++++------
Inst_eddies/get_1_outerm.f90 | 63 ++++++++++++------------
Inst_eddies/set_contours.f90 | 19 +++----
Inst_eddies/set_max_speed.f90 | 38 +++++++-------
5 files changed, 94 insertions(+), 88 deletions(-)
diff --git a/Inst_eddies/Tests/test_get_1_outerm.f90 b/Inst_eddies/Tests/test_get_1_outerm.f90
index f4752c89..89f60fa5 100644
--- a/Inst_eddies/Tests/test_get_1_outerm.f90
+++ b/Inst_eddies/Tests/test_get_1_outerm.f90
@@ -49,22 +49,23 @@ program test_get_1_outerm
logical:: grid_lon_lat = .true.
integer, parameter:: n_max_cont = 31
- type(ssh_contour) cont_list(n_max_cont)
+ type(polyline) cont_list(n_max_cont)
! Contour list, with points given by longitude and latitude, in
- ! radians. Defined only for subscripts 1:n_cont. If n_cont >= 1,
- ! the outermost contour is element with subscript n_cont. The
- ! contours are in monotonic order of ssh. For i <= n_cont,
- ! cont_list(i)%area has its default initialization value and the
- ! order of points in cont_list(i)%polyline (clockwise or
+ ! radians. Defined only for subscripts 1:n_cont. If n_cont >= 1, the
+ ! outermost contour is element with subscript n_cont. For i <=
+ ! n_cont, the order of points in cont_list(i) (clockwise or
! counter-clockwise) is not specified.
- type(polyline) cont_list_proj(n_max_cont)
- ! Contour list, with points given by projection
- ! coordinates. Defined only for subscripts 1:n_cont. The contour
- ! with a given subscript in cont_list_proj corresponds to the
- ! contour with the same subscript in cont_list. For i <= n_cont,
- ! the order of points in cont_list_proj(i) (clockwise or
- ! counter-clockwise) is not specified.
+ type(ssh_contour) cont_list_proj(n_max_cont)
+ ! Contour list, with points given by projection coordinates. Defined
+ ! only for subscripts 1:n_cont. If n_cont >= 1, the outermost
+ ! contour is element with subscript n_cont. The contours are in
+ ! monotonic order of ssh. The contour with a given subscript in
+ ! cont_list_proj corresponds to the contour with the same subscript
+ ! in cont_list. For i <= n_cont, cont_list_proj(i)%area has its
+ ! default initialization value and the order of points in
+ ! cont_list_proj(i) (clockwise or counter-clockwise) is not
+ ! specified.
integer n_cont ! number of good contours found and stored in cont_list
TYPE(shpfileobject) hshp_cont_list
@@ -117,7 +118,7 @@ program test_get_1_outerm
call shp_append_object_03(ishape, hshp_cont_list, shpt_polygon, &
cont_list(i)%points * rad_to_deg)
call dbf_write_attribute_03(hshp_cont_list, ishape, ifield, &
- cont_list(i)%ssh)
+ cont_list_proj(i)%ssh)
end do
CALL shpclose(hshp_cont_list)
diff --git a/Inst_eddies/Tests/test_set_max_speed.f90 b/Inst_eddies/Tests/test_set_max_speed.f90
index bce06775..8a3ccdc1 100644
--- a/Inst_eddies/Tests/test_set_max_speed.f90
+++ b/Inst_eddies/Tests/test_set_max_speed.f90
@@ -40,22 +40,23 @@ program test_set_max_speed
integer, parameter:: n_max_cont = 31
- type(ssh_contour) cont_list(n_max_cont)
+ type(polyline) cont_list(n_max_cont)
! Contour list, with points given by longitude and latitude, in
- ! radians. Defined only for subscripts 1:n_cont. If n_cont >= 1,
- ! the outermost contour is element with subscript n_cont. The
- ! contours are in monotonic order of ssh. For i <= n_cont,
- ! cont_list(i)%area has its default initialization value and the
- ! order of points in cont_list(i)%polyline (clockwise or
+ ! radians. Defined only for subscripts 1:n_cont. If n_cont >= 1, the
+ ! outermost contour is element with subscript n_cont. For i <=
+ ! n_cont, the order of points in cont_list(i) (clockwise or
! counter-clockwise) is not specified.
- type(polyline) cont_list_proj(n_max_cont)
- ! Contour list, with points given by projection
- ! coordinates. Defined only for subscripts 1:n_cont. The contour
- ! with a given subscript in cont_list_proj corresponds to the
- ! contour with the same subscript in cont_list. For i <= n_cont,
- ! the order of points in cont_list_proj(i) (clockwise or
- ! counter-clockwise) is not specified.
+ type(ssh_contour) cont_list_proj(n_max_cont)
+ ! Contour list, with points given by projection coordinates. Defined
+ ! only for subscripts 1:n_cont. If n_cont >= 1, the outermost
+ ! contour is element with subscript n_cont. The contours are in
+ ! monotonic order of ssh. The contour with a given subscript in
+ ! cont_list_proj corresponds to the contour with the same subscript
+ ! in cont_list. For i <= n_cont, cont_list_proj(i)%area has its
+ ! default initialization value and the order of points in
+ ! cont_list_proj(i) (clockwise or counter-clockwise) is not
+ ! specified.
TYPE(shpfileobject) hshp
type(polygon) p
@@ -102,13 +103,13 @@ program test_set_max_speed
call shp_read_pol(p, hshp, ishape)
if (p%nparts == 0) then
- cont_list(i)%polyline = null_polyline()
+ cont_list(i) = null_polyline()
else
- cont_list(i)%polyline = p%part(1)
+ cont_list(i) = p%part(1)
cont_list(i)%points = cont_list(i)%points * deg_to_rad
end if
- call dbf_read_attribute_03(cont_list(i)%ssh, hshp, ifield = 0, &
+ call dbf_read_attribute_03(cont_list_proj(i)%ssh, hshp, ifield = 0, &
ishape = ishape)
cont_list_proj(i)%n_points = cont_list(i)%n_points
cont_list_proj(i)%closed = cont_list(i)%closed
diff --git a/Inst_eddies/get_1_outerm.f90 b/Inst_eddies/get_1_outerm.f90
index 7cbe89c1..682cf39a 100644
--- a/Inst_eddies/get_1_outerm.f90
+++ b/Inst_eddies/get_1_outerm.f90
@@ -21,12 +21,12 @@ contains
! the outermost contour does not have sufficient area. If an
! outermost contour is found then its level is farther from the
! extremum than innermost_level_2. The outermost contour is
- ! duplicated in the array cont_list. In cont_list, the contours,
- ! including the outermost contour, are in monotonic order of ssh:
- ! ascending if the extremum is a minimum, descending if the
- ! extremum is a maximum. The area is not computed for contours in
- ! cont_list. On return, the array
- ! abs(ediff1d(cont_list(:n_cont)%ssh)) is in descending
+ ! duplicated in the array cont_list. In cont_list_proj, the
+ ! contours, including the outermost contour, are in monotonic
+ ! order of ssh: ascending if the extremum is a minimum, descending
+ ! if the extremum is a maximum. The area is not computed for
+ ! contours in cont_list_proj. On return, the array
+ ! abs(ediff1d(cont_list_proj(:n_cont)%ssh)) is in descending
! order. (See notes for a proof.) On return, n_cont == 1 means we
! have only found one good contour at innermost_level_2 and it has
! sufficient area.
@@ -55,22 +55,22 @@ contains
type(ssh_contour), intent(out):: out_cont
! points given by longitude, latitude, in rad
- type(ssh_contour), intent(out):: cont_list(:) ! (n_max_cont)
+ type(polyline), intent(out):: cont_list(:) ! (n_max_cont)
! Contour list, with points given by longitude and latitude, in
! radians. Defined only for subscripts 1:n_cont. If n_cont >= 1,
- ! the outermost contour is element with subscript n_cont. The
- ! contours are in monotonic order of ssh. For i <= n_cont,
- ! cont_list(i)%area has its default initialization value and the
- ! order of points in cont_list(i)%polyline (clockwise or
+ ! the outermost contour is element with subscript n_cont. For i <=
+ ! n_cont, the order of points in cont_list(i) (clockwise or
! counter-clockwise) is not specified.
- type(polyline), intent(out):: cont_list_proj(:) ! (n_max_cont)
+ type(ssh_contour), intent(out):: cont_list_proj(:) ! (n_max_cont)
! Contour list, with points given by projection
! coordinates. Defined only for subscripts 1:n_cont. If n_cont >=
! 1, the outermost contour is element with subscript n_cont. The
- ! contour with a given subscript in cont_list_proj corresponds to
- ! the contour with the same subscript in cont_list. For i <=
- ! n_cont, the order of points in cont_list_proj(i) (clockwise or
+ ! contours are in monotonic order of ssh. The contour with a given
+ ! subscript in cont_list_proj corresponds to the contour with the
+ ! same subscript in cont_list. For i <= n_cont,
+ ! cont_list_proj(i)%area has its default initialization value and
+ ! the order of points in cont_list_proj(i) (clockwise or
! counter-clockwise) is not specified.
integer, intent(out):: n_cont
@@ -111,23 +111,23 @@ contains
corner_proj = (corner - corner_whole) / step + 1.
coord_extr_proj = (coord_extr - corner_whole) / step + 1.
outside_points_proj = convert_to_ind(outside_points, corner_whole, step)
- n_max_cont = size(cont_list)
+ n_max_cont = size(cont_list_proj)
n_cont = 0
- cont_list_proj(1) = good_contour(corner_proj, ssh, innermost_level_2, &
- coord_extr_proj, outside_points_proj)
+ cont_list_proj(1)%polyline = good_contour(corner_proj, ssh, &
+ innermost_level_2, coord_extr_proj, outside_points_proj)
test_n_points: if (.not. cont_list_proj(1)%closed) then
! cont_list_proj(1)%n_points == 0
out_cont = null_ssh_contour()
else test_n_points
n_cont = 1
- cont_list(1)%ssh = innermost_level_2
+ cont_list_proj(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")
- cont_list_proj(2) = good_contour(corner_proj, ssh, level_try, &
+ cont_list_proj(2)%polyline = good_contour(corner_proj, ssh, level_try, &
coord_extr_proj, outside_points_proj)
if (cont_list_proj(2)%closed) then
@@ -136,24 +136,23 @@ contains
! everywhere around the path of contour. So it should
! probably never happen.
n_cont = 2
- cont_list(2)%ssh = level_try
+ cont_list_proj(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.
- cont_list_proj(n_cont + 1) = good_contour(corner_proj, ssh, &
- level_try, coord_extr_proj, outside_points_proj)
+ cont_list_proj(n_cont + 1)%polyline = good_contour(corner_proj, &
+ ssh, level_try, coord_extr_proj, outside_points_proj)
if (cont_list_proj(n_cont + 1)%closed) then
! cont_list_proj(n_cont + 1)%n_points /= 0
- cont_list(n_cont + 1)%ssh = level_try
+ cont_list_proj(n_cont + 1)%ssh = level_try
if (n_cont == n_max_cont - 1) then
! Replace the previous good contour found
cont_list_proj(n_cont) = cont_list_proj(n_max_cont)
- cont_list(n_cont)%ssh = cont_list(n_max_cont)%ssh
else
n_cont = n_cont + 1
end if
@@ -163,7 +162,7 @@ contains
level_bad = level_try
end if
- ! Assertion: cont_list_proj(n_cont) ==
+ ! Assertion: cont_list_proj(n_cont)%polyline ==
! good_contour(. . ., level_good) and
! cont_list_proj(n_cont)%closed and no good contour at
! level_bad
@@ -171,7 +170,8 @@ contains
! Assertion: n_cont <= n_max_cont - 1
end if
- out_cont%polyline = convert_to_reg_coord(cont_list_proj(n_cont), &
+ out_cont%polyline &
+ = convert_to_reg_coord(cont_list_proj(n_cont)%polyline, &
corner_whole, step)
out_cont%area = spher_polyline_area(out_cont%polyline)
@@ -180,10 +180,11 @@ contains
n_cont = 0
else
! n_cont >= 1
- forall (i = 1:n_cont - 1) cont_list(i)%polyline &
- = convert_to_reg_coord(cont_list_proj(i), corner_whole, step)
- cont_list(n_cont)%polyline = out_cont%polyline
- out_cont%ssh = cont_list(n_cont)%ssh
+ forall (i = 1:n_cont - 1) cont_list(i) &
+ = convert_to_reg_coord(cont_list_proj(i)%polyline, &
+ corner_whole, step)
+ cont_list(n_cont) = out_cont%polyline
+ out_cont%ssh = cont_list_proj(n_cont)%ssh
call ccw_orient(out_cont)
end if
end if test_n_points
diff --git a/Inst_eddies/set_contours.f90 b/Inst_eddies/set_contours.f90
index 5841b673..2c85c81c 100644
--- a/Inst_eddies/set_contours.f90
+++ b/Inst_eddies/set_contours.f90
@@ -64,20 +64,21 @@ contains
integer, parameter:: n_max_cont = 31 ! must be >= 3
- type(ssh_contour) cont_list(n_max_cont)
+ type(polyline) cont_list(n_max_cont)
! Contour list, with points given by longitude and latitude, in
! radians. Defined only for subscripts 1:n_cont. If n_cont >= 1,
- ! the outermost contour is element with subscript n_cont. The
- ! contours are in monotonic order of ssh. For i <= n_cont,
- ! cont_list(i)%area has its default initialization value and the
- ! order of points in cont_list(i)%polyline (clockwise or
+ ! the outermost contour is element with subscript n_cont. For i <=
+ ! n_cont, the order of points in cont_list(i) (clockwise or
! counter-clockwise) is not specified.
- type(polyline) cont_list_proj(n_max_cont)
+ type(ssh_contour) cont_list_proj(n_max_cont)
! Contour list, with points given by projection
- ! coordinates. Defined only for subscripts 1:n_cont. The contour
- ! with a given subscript in cont_list_proj corresponds to the
- ! contour with the same subscript in cont_list. For i <= n_cont,
+ ! coordinates. Defined only for subscripts 1:n_cont. If n_cont >=
+ ! 1, the outermost contour is element with subscript n_cont. The
+ ! contours are in monotonic order of ssh. The contour with a given
+ ! subscript in cont_list_proj corresponds to the contour with the
+ ! same subscript in cont_list. For i <= n_cont,
+ ! cont_list_proj(i)%area has its default initialization value and
! the order of points in cont_list_proj(i) (clockwise or
! counter-clockwise) is not specified.
diff --git a/Inst_eddies/set_max_speed.f90 b/Inst_eddies/set_max_speed.f90
index c02fec28..81b81275 100644
--- a/Inst_eddies/set_max_speed.f90
+++ b/Inst_eddies/set_max_speed.f90
@@ -42,20 +42,21 @@ contains
! speed_cont is not a null ssh contour then max_speed is the
! speed on speed_cont.
- type(ssh_contour), intent(inout):: cont_list(:) ! (n_max_cont)
+ type(polyline), intent(inout):: cont_list(:) ! (n_max_cont)
! Contour list, with points given by longitude and latitude, in
- ! radians. Defined only for subscripts 1:n_cont. On entry, the
- ! outermost contour is element with subscript n_cont. The contours
- ! between 1 and i_outer are in monotonic order of
- ! ssh. cont_list%area has its default initialization value and the
- ! order of points in the polyline components (clockwise or
+ ! radians. Defined only for subscripts 1:n_cont. If n_cont >= 1,
+ ! the outermost contour is element with subscript n_cont. For i <=
+ ! n_cont, the order of points in cont_list(i) (clockwise or
! counter-clockwise) is not specified.
- type(polyline), intent(inout):: cont_list_proj(:) ! (n_max_cont)
+ type(ssh_contour), intent(inout):: cont_list_proj(:) ! (n_max_cont)
! Contour list, with points given by projection
- ! coordinates. Defined only for subscripts 1:n_cont. The contour
- ! with a given subscript in cont_list_proj corresponds to the
- ! contour with the same subscript in cont_list. For i <= n_cont,
+ ! coordinates. Defined only for subscripts 1:n_cont. If n_cont >=
+ ! 1, the outermost contour is element with subscript n_cont. The
+ ! contours are in monotonic order of ssh. The contour with a given
+ ! subscript in cont_list_proj corresponds to the contour with the
+ ! same subscript in cont_list. For i <= n_cont,
+ ! cont_list_proj(i)%area has its default initialization value and
! the order of points in cont_list_proj(i) (clockwise or
! counter-clockwise) is not specified.
@@ -89,13 +90,13 @@ contains
corner_proj = (corner - corner_whole) / step + 1.
extr_coord_proj = (extr_coord - corner_whole) / step + 1.
i_outer = n_cont
- if (n_cont >= 2) call complete_ssh(cont_list%ssh, n_cont)
+ if (n_cont >= 2) call complete_ssh(cont_list_proj%ssh, n_cont)
! Now find the contours associated to the new values of SSH:
do i = i_outer + 1, n_cont
- cont_list_proj(i) = good_contour(corner_proj, ssh, cont_list(i)%ssh, &
- extr_coord_proj)
- cont_list(i)%polyline = convert_to_reg_coord(cont_list_proj(i), &
+ cont_list_proj(i)%polyline = good_contour(corner_proj, ssh, &
+ cont_list_proj(i)%ssh, extr_coord_proj)
+ cont_list(i) = convert_to_reg_coord(cont_list_proj(i)%polyline, &
corner_whole, step)
end do
@@ -103,8 +104,8 @@ contains
do i = 1, n_cont
if (cont_list(i)%closed) then
- speed(i) = mean_speed(u, v, cont_list(i)%polyline, &
- cont_list_proj(i)%points, extr_coord, corner, step)
+ speed(i) = mean_speed(u, v, cont_list(i), cont_list_proj(i)%points, &
+ extr_coord, corner, step)
else
! It is possible that good_contour returned a null polyline
! if there is a missing value of SSH inside the outermost
@@ -118,7 +119,7 @@ contains
call shp_append_object_03(ishape, hshp_cont_list, shpt_polygon, &
cont_list(i)%points * rad_to_deg)
call dbf_write_attribute_03(hshp_cont_list, ishape, ifield_ssh, &
- fieldvalue = cont_list(i)%ssh)
+ fieldvalue = cont_list_proj(i)%ssh)
if (ieee_is_nan(speed(i))) then
call dbf_write_attribute_03(hshp_cont_list, ishape, ifield_speed, &
@@ -145,7 +146,8 @@ contains
! Maximum speed on the outermost contour
speed_cont = null_ssh_contour()
else
- speed_cont = cont_list(i)
+ speed_cont%polyline = cont_list(i)
+ speed_cont%ssh = cont_list_proj(i)%ssh
speed_cont%area = spher_polyline_area(speed_cont%polyline)
call ccw_orient(speed_cont)
end if
--
GitLab