From 4399e353c21746c7aba07950e531b71cefc71de1 Mon Sep 17 00:00:00 2001
From: Lionel GUEZ <guez@lmd.ens.fr>
Date: Wed, 19 Jun 2024 16:11:45 +0200
Subject: [PATCH] Compute good contour in projection coordinates
Procedure `good_contour` now takes input in projection coordinates and
returns a contour in projection coordinates. The dummy argument step
of `good_contour` is thus removed.
---
Inst_eddies/Tests/Input/good_contour.txt | 5 ++--
Inst_eddies/Tests/Input/good_contour_2.txt | 4 ++--
Inst_eddies/Tests/Input/no_good_contour.txt | 2 +-
Inst_eddies/Tests/Input/outside_points_1.csv | 4 ++--
Inst_eddies/Tests/Input/outside_points_2.csv | 4 ++--
Inst_eddies/Tests/test_good_contour.f90 | 5 ++--
Inst_eddies/get_1_outerm.f90 | 25 ++++++++++++++------
Inst_eddies/good_contour.f90 | 6 ++---
Inst_eddies/set_max_speed.f90 | 9 +++++--
9 files changed, 38 insertions(+), 26 deletions(-)
diff --git a/Inst_eddies/Tests/Input/good_contour.txt b/Inst_eddies/Tests/Input/good_contour.txt
index edde4556..9f0246b1 100644
--- a/Inst_eddies/Tests/Input/good_contour.txt
+++ b/Inst_eddies/Tests/Input/good_contour.txt
@@ -1,5 +1,4 @@
&MAIN_NML
- INSIDE_POINT=-0.7,0
+ INSIDE_POINT=-8.75,0.
LEVEL=5
-STEP=2*0.08
-CORNER=-2,-2 /
+CORNER=-25.,-25. /
diff --git a/Inst_eddies/Tests/Input/good_contour_2.txt b/Inst_eddies/Tests/Input/good_contour_2.txt
index cbb59b81..6d7bd203 100644
--- a/Inst_eddies/Tests/Input/good_contour_2.txt
+++ b/Inst_eddies/Tests/Input/good_contour_2.txt
@@ -1,3 +1,3 @@
-&MAIN_NML INSIDE_POINT=0.3,0,
+&MAIN_NML INSIDE_POINT=3.75,0,
LEVEL=5
-corner = -2., -2., step = 0.08, 0.08 /
+corner = -25., -25. /
diff --git a/Inst_eddies/Tests/Input/no_good_contour.txt b/Inst_eddies/Tests/Input/no_good_contour.txt
index e255617d..5ae3c521 100644
--- a/Inst_eddies/Tests/Input/no_good_contour.txt
+++ b/Inst_eddies/Tests/Input/no_good_contour.txt
@@ -1 +1 @@
-&MAIN_NML LEVEL= 5, corner = -2., -2., step = 0.08, 0.08 /
+&MAIN_NML LEVEL= 5, corner = -25., -25. /
diff --git a/Inst_eddies/Tests/Input/outside_points_1.csv b/Inst_eddies/Tests/Input/outside_points_1.csv
index b6fb92e0..7b7aeaa4 100644
--- a/Inst_eddies/Tests/Input/outside_points_1.csv
+++ b/Inst_eddies/Tests/Input/outside_points_1.csv
@@ -1,2 +1,2 @@
--0.5, -1.
-0.3, 0.
+-6.25, -12.5
+3.75, 0.
diff --git a/Inst_eddies/Tests/Input/outside_points_2.csv b/Inst_eddies/Tests/Input/outside_points_2.csv
index 3260f5c4..396d0790 100644
--- a/Inst_eddies/Tests/Input/outside_points_2.csv
+++ b/Inst_eddies/Tests/Input/outside_points_2.csv
@@ -1,2 +1,2 @@
--0.5, -1.
--0.7, 0.
+-6.25, -12.5
+-8.75, 0.
diff --git a/Inst_eddies/Tests/test_good_contour.f90 b/Inst_eddies/Tests/test_good_contour.f90
index 9099ce40..c0fe3888 100644
--- a/Inst_eddies/Tests/test_good_contour.f90
+++ b/Inst_eddies/Tests/test_good_contour.f90
@@ -21,7 +21,6 @@ program test_good_contour
real, allocatable:: Z(:, :)
real:: inside_point(2) = [0., 0.]
real, allocatable:: outside_points(:, :) ! (2, n_out)
- real:: step(2) = [1., 1.]
real:: corner(2) = [1., 1.]
TYPE(shpfileobject) shphandle
integer field_number, shape_number
@@ -29,7 +28,7 @@ program test_good_contour
integer, allocatable:: dimids(:)
character(len = 10):: varname = "z"
- namelist /main_nml/ inside_point, level, step, corner, varname
+ namelist /main_nml/ inside_point, level, corner, varname
!-------------------------------------------------------------------
@@ -70,7 +69,7 @@ program test_good_contour
print *, outside_points(:, i)
end do
- c = good_contour(corner, step, z, level, inside_point, outside_points)
+ c = good_contour(corner, z, level, inside_point, outside_points)
if (c%n_points /= 0) then
call shp_create_03("test_good_contour", shpt_polygon, shphandle)
diff --git a/Inst_eddies/get_1_outerm.f90 b/Inst_eddies/get_1_outerm.f90
index ff651eb0..008afc2f 100644
--- a/Inst_eddies/get_1_outerm.f90
+++ b/Inst_eddies/get_1_outerm.f90
@@ -41,13 +41,14 @@ contains
! 360 degrees.
! Libraries:
- use contour_531, only: polyline
+ use contour_531, only: polyline, convert_to_reg_coord, convert_to_ind
use jumble, only: assert
use ccw_orient_m, only: ccw_orient
use config_m, only: min_area
use derived_types, only: ssh_contour, null_ssh_contour
use good_contour_m, only: good_contour
+ use input_ssh_m, only: corner_whole => corner
use spher_polyline_area_m, only: spher_polyline_area
type(ssh_contour), intent(out):: out_cont
@@ -88,13 +89,20 @@ contains
real, parameter:: precision = 5e-4 ! in m
logical mask(size(ssh, 1), size(ssh, 2))
integer n_max_cont ! >= 3
+ integer i
+ real corner_proj(2), coord_extr_proj(2)
+ real outside_points_proj(2, size(outside_points, 2))
!-----------------------------------------------------------------
+ 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_cont = 0
- cont_list(1)%polyline = good_contour(corner, step, ssh, innermost_level_2, &
- coord_extr, outside_points)
+ cont_list(1)%polyline = convert_to_reg_coord(good_contour(corner_proj, &
+ ssh, innermost_level_2, coord_extr_proj, outside_points_proj), &
+ corner_whole, step)
test_n_points: if (cont_list(1)%n_points == 0) then
out_cont = null_ssh_contour()
@@ -106,8 +114,9 @@ contains
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(2)%polyline = good_contour(corner, step, ssh, level_try, &
- coord_extr, outside_points)
+ cont_list(2)%polyline = convert_to_reg_coord(good_contour(corner_proj, &
+ ssh, level_try, coord_extr_proj, outside_points_proj), &
+ corner_whole, step)
if (cont_list(2)%n_points /= 0) then
! This can only happen if there are degenerate extrema
@@ -121,8 +130,10 @@ contains
do while (abs(level_bad - level_good) > precision)
level_try = (level_good + level_bad) / 2.
- cont_list(n_cont + 1)%polyline = good_contour(corner, step, ssh, &
- level_try, coord_extr, outside_points)
+ cont_list(n_cont + 1)%polyline &
+ = convert_to_reg_coord(good_contour(corner_proj, ssh, &
+ level_try, coord_extr_proj, outside_points_proj), &
+ corner_whole, step)
if (cont_list(n_cont + 1)%n_points /= 0) then
cont_list(n_cont + 1)%ssh = level_try
diff --git a/Inst_eddies/good_contour.f90 b/Inst_eddies/good_contour.f90
index 856107a1..9eb3e0b1 100644
--- a/Inst_eddies/good_contour.f90
+++ b/Inst_eddies/good_contour.f90
@@ -4,7 +4,7 @@ module good_contour_m
contains
- type(polyline) function good_contour(corner, step, z, level, inside_point, &
+ type(polyline) function good_contour(corner, z, level, inside_point, &
outside_points)
! Finds a contour in the plane, if it exists, which is closed,
@@ -20,8 +20,6 @@ contains
use geometry, only: polygon_contains_point
real, intent(in):: corner(:) ! (2) [x, y] corresponding to z(1, 1)
- real, intent(in):: step(:) ! (2)
-
real, intent(in):: z(:, :)
! We are assuming that significant values of z cannot be greater
! than zmax and that the missing value flag in z is greater than
@@ -48,7 +46,7 @@ contains
n_out = 0
end if
- call find_contours_reg_grid(corner, step, z, level, contours, zmax)
+ call find_contours_reg_grid(corner, [1., 1.], z, level, contours, zmax)
n_cont = size(contours)
found_good_contour = .false.
i = 1
diff --git a/Inst_eddies/set_max_speed.f90 b/Inst_eddies/set_max_speed.f90
index a9594985..27a763b0 100644
--- a/Inst_eddies/set_max_speed.f90
+++ b/Inst_eddies/set_max_speed.f90
@@ -22,11 +22,13 @@ contains
use jumble, only: rad_to_deg
use shapelib, only: shpt_polygon
use shapelib_03, only: shp_append_object_03, dbf_write_attribute_03
+ use contour_531, only: convert_to_reg_coord
use ccw_orient_m, only: ccw_orient
use complete_ssh_m, only: complete_ssh
use derived_types, only: eddy, null_ssh_contour, ssh_contour, missing_speed
use good_contour_m, only: good_contour
+ use input_ssh_m, only: corner_whole => corner
use mean_speed_m, only: mean_speed
use spher_polyline_area_m, only: spher_polyline_area
use cont_list_m, only: write_cont_list, hshp_cont_list, ifield_ssh, &
@@ -71,16 +73,19 @@ contains
real, allocatable:: speed(:) ! (n_cont) speed on the contour
integer i, i_outer, ishape
+ real corner_proj(2), extr_coord_proj(2)
!---------------------------------------------------------------
+ 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)
! Now find the contours associated to the new values of SSH:
do i = i_outer + 1, n_cont
- cont_list(i)%polyline = good_contour(corner, step, ssh, &
- cont_list(i)%ssh, extr_coord)
+ cont_list(i)%polyline = convert_to_reg_coord(good_contour(corner_proj, &
+ ssh, cont_list(i)%ssh, extr_coord_proj), corner_whole, step)
end do
allocate(speed(n_cont))
--
GitLab