From f79616fa61a4db348ddddfc8e9f1b36250c36f94 Mon Sep 17 00:00:00 2001
From: Lionel GUEZ <guez@lmd.ipsl.fr>
Date: Fri, 5 Jul 2024 17:32:19 +0200
Subject: [PATCH] Redefine `outside_points` in projection space
In main program units `inst_eddies`, `examine_eddy` and
`test_get_1_outerm`, redefine `outside_points` as coordinates in
projection space, avoiding a round trip to longitude and latitude. In
`examine_eddy`, we therefore store `coord_proj` instead of coord.
---
.../Tests/Input/Region_1/outside_points.csv | 13 ++--
Inst_eddies/Tests/examine_eddy.f90 | 61 +++++++++----------
Inst_eddies/Tests/test_get_1_outerm.f90 | 18 +++---
Inst_eddies/inst_eddies.f90 | 30 ++++-----
Inst_eddies/set_contours.f90 | 6 +-
5 files changed, 57 insertions(+), 71 deletions(-)
diff --git a/Inst_eddies/Tests/Input/Region_1/outside_points.csv b/Inst_eddies/Tests/Input/Region_1/outside_points.csv
index 6b40dd18..752a7b9f 100644
--- a/Inst_eddies/Tests/Input/Region_1/outside_points.csv
+++ b/Inst_eddies/Tests/Input/Region_1/outside_points.csv
@@ -1,7 +1,6 @@
-"degrees east" "degrees north"
-"longitude" "latitude"
-8.375 -36.125
-6.125 -35.375
-8.125 -34.625
-10.875 -34.625
-5.875 -33.125
+ x y
+ 13 8
+ 24 8
+ 4 14
+ 14 2
+ 5 5
diff --git a/Inst_eddies/Tests/examine_eddy.f90 b/Inst_eddies/Tests/examine_eddy.f90
index a1b56731..501f7049 100644
--- a/Inst_eddies/Tests/examine_eddy.f90
+++ b/Inst_eddies/Tests/examine_eddy.f90
@@ -20,15 +20,13 @@ program examine_eddy
! The main difficulty here was to get outside_points. We need in
! outside_points all the extrema of opposite sign (including those
! not associated to an eddy) plus the extrema of same sign
- ! associated to an eddy. Note that coord and outside_points may not
- ! be exactly equal to what they were in inst_eddies because there is
- ! a multiplication of coord by rad_to_deg when writing the shapefile
- ! and a multiplication by deg_to_rad when reading it.
+ ! associated to an eddy. Note that e%extr%coord may not be exactly
+ ! equal to what it was in inst_eddies because there is a
+ ! multiplication by rad_to_deg when writing the shapefile and a
+ ! multiplication by deg_to_rad when reading it.
! Libraries:
- use contour_531, only: convert_to_ind
- use jumble, only: get_command_arg_dyn, read_column, deg_to_rad, rad_to_deg, &
- twopi, new_unit
+ use jumble, only: get_command_arg_dyn, read_column, deg_to_rad, new_unit
use shapelib, only: shpfileobject, shpclose
use shapelib_03, only: dbf_read_attribute_03, dbf_get_field_index_03, &
shp_read_point, shp_open_03
@@ -47,7 +45,7 @@ program examine_eddy
implicit none
integer:: date = huge(0), eddy_i_target = huge(0)
- integer ishape_first, d0, ishape, i, l, n1, n2, number_extr, unit
+ integer ishape_first, d0, ishape, i, l, n1, n2, number_extr, unit, nlon
logical:: cyclone = .true. ! orientation of target eddy
namelist /main_nml/ cyclone, date, eddy_i_target
TYPE(shpfileobject) hshp
@@ -71,26 +69,24 @@ program examine_eddy
! (nlon, nlat). Wind, in m s-1.
type(snapshot) s
+ real coord(2) ! longitude and latitude, in rad
- real, allocatable:: coord(:, :)
- ! (2, number_extr) longitude and latitude, in rad
+ integer, allocatable:: coord_proj(:, :)
+ ! (2, number_extr) coordinates in projection space
integer, allocatable:: selection(:)
! identifying numbers of a selection of eddies
logical, allocatable:: in_window(:)
- real, allocatable:: outside_points(:, :) ! (2, :) longitude and
- ! latitude, in rad, of all the significant extrema, except the
- ! target extremum
+ integer, allocatable:: outside_points(:, :) ! (2, :)
+ ! coordinates in projection space of all the significant extrema,
+ ! except the target extremum
TYPE(shpc_slice_handler) hshpc
integer llc(2) ! indices in global grid of lower left corner
integer urc(2) ! indices in global grid of upper right corner
- real corner_window(2) ! longitude and latitude of the window around each
- ! extremum, in rad
-
!------------------------------------------------------------------------
call config
@@ -122,15 +118,16 @@ program examine_eddy
end if
number_extr = ishape_last(date) - ishape_first + 1
- allocate(coord(2, number_extr), in_window(number_extr))
+ allocate(coord_proj(2, number_extr), in_window(number_extr))
do ishape = ishape_first, ishape_last(date)
- call shp_read_point(coord(:, ishape - ishape_first + 1), hshp, ishape)
+ call shp_read_point(coord, hshp, ishape)
+ coord_proj(:, ishape - ishape_first + 1) = nint((coord * deg_to_rad &
+ - corner) / step + 1.)
end do
- coord = coord * deg_to_rad
- e%extr%coord = coord(:, eddy_i_target)
- e%extr%coord_proj = nint(convert_to_ind(e%extr%coord, corner, step))
+ e%extr%coord = (coord_proj(:, eddy_i_target) - 1) * step + corner
+ e%extr%coord_proj = coord_proj(:, eddy_i_target)
call dbf_get_field_index_03(hshp, "ssh", ifield)
call dbf_read_attribute_03(e%extr%ssh, hshp, ifield, &
ishape = ishape_first + eddy_i_target - 1)
@@ -149,18 +146,18 @@ program examine_eddy
urc(1) = min(urc(1), size(ssh, 1))
end if
- corner_window = corner + (llc - 1) * step
-
! Done defining geographical window
+ if (periodic) nlon = size(ssh, 1) - 2 * max_radius(1)
+
do i = 1, number_extr
if (i /= eddy_i_target) then
! Shift the longitude to a value close to the longitude of the
! target extremum:
- if (periodic) coord(1, i) = coord(1, i) &
- + ceiling((corner_window(1) - coord(1, i)) / twopi) * twopi
+ if (periodic) coord_proj(1, i) = coord_proj(1, i) &
+ + ceiling((llc(1) - coord_proj(1, i)) / real(nlon)) * nlon
- in_window(i) = all(abs(nint((coord(:, i) - e%extr%coord) / step)) &
+ in_window(i) = all(abs(coord_proj(:, i) - e%extr%coord_proj) &
< max_radius)
else
in_window(i) = .false.
@@ -174,12 +171,12 @@ program examine_eddy
selection = pack(selection, s%list(selection)%cyclone .neqv. cyclone)
n1 = size(selection)
allocate(outside_points(2, n1 + n2))
- forall (i = 1:n1) outside_points(:, i) = s%list(selection(i))%extr%coord
+ forall (i = 1:n1) outside_points(:, i) = s%list(selection(i))%extr%coord_proj
l = n1 + 1
do i = 1, number_extr
if (in_window(i)) then
- outside_points(:, l) = coord(:, i)
+ outside_points(:, l) = coord_proj(:, i)
l = l + 1
end if
end do
@@ -190,7 +187,7 @@ program examine_eddy
call set_contours(e%out_cont, e%speed_cont, e%max_speed, cyclone, e%extr, &
e%innermost_level, step, ssh(llc(1):urc(1), llc(2):urc(2)), &
u(llc(1):urc(1), llc(2):urc(2)), v(llc(1):urc(1), llc(2):urc(2)), llc, &
- convert_to_ind(outside_points, corner, step))
+ outside_points)
call close_cont_list
call shpc_create(hshpc, shpc_dir = "SHPC", cyclone = cyclone, slice = 0, &
grid_lon_lat = .true.)
@@ -199,7 +196,6 @@ program examine_eddy
CALL shpc_close(hshpc)
print *, "llc = ", llc
print *, "urc = ", urc
- print *, "corner_window (in degrees):", corner_window * rad_to_deg
print *, "Number of extrema of opposite orientation which must be outside ", &
"outermost contour:", n1
print *, "Number of (valid) eddies of same orientation which must be ", &
@@ -207,11 +203,10 @@ program examine_eddy
call new_unit(unit)
open(unit, file = "outside_points.csv", status = "replace", action = "write")
- write(unit, fmt = *) '"degrees east" "degrees north"'
- write(unit, fmt = *) "longitude latitude"
+ write(unit, fmt = *) "x y"
do i = 1, n1 + n2
- write(unit, fmt = *) outside_points(:, i) * rad_to_deg
+ write(unit, fmt = *) outside_points(:, i)
end do
close(unit)
diff --git a/Inst_eddies/Tests/test_get_1_outerm.f90 b/Inst_eddies/Tests/test_get_1_outerm.f90
index 86b68a2c..201f50f6 100644
--- a/Inst_eddies/Tests/test_get_1_outerm.f90
+++ b/Inst_eddies/Tests/test_get_1_outerm.f90
@@ -6,8 +6,8 @@ program test_get_1_outerm
! cont_list and cont_list_proj.
! Libraries:
- use contour_531, only: polyline, convert_to_ind
- use jumble, only: csvread, deg_to_rad, rad_to_deg
+ use contour_531, only: polyline
+ use jumble, only: csvread, rad_to_deg
use shapelib, only: shpfileobject, ftdouble, shpt_polygon, shpclose
use shapelib_03, only: shp_create_03, dbf_add_field_03, &
shp_append_object_03, dbf_write_attribute_03
@@ -47,9 +47,9 @@ program test_get_1_outerm
TYPE(shpc_slice_handler) hshpc
type(eddy) e
- real, allocatable:: outside_points(:, :) ! (2, :) longitude and
- ! latitude, in rad, of all the significant extrema, except the
- ! target extremum
+ integer, allocatable:: outside_points(:, :) ! (2, :)
+ ! coordinates in projection space of all the significant extrema,
+ ! except the target extremum
logical:: grid_lon_lat = .true.
integer, parameter:: n_max_cont = 31
@@ -89,8 +89,8 @@ program test_get_1_outerm
inquire(file = "outside_points.csv", exist = exist)
if (exist) then
- call csvread("outside_points.csv", outside_points, first_r = 3)
- outside_points = transpose(outside_points) * deg_to_rad
+ call csvread("outside_points.csv", outside_points, first_r = 2)
+ outside_points = transpose(outside_points)
else
print *, '"outside_points.csv" not found. Assuming no outside points...'
allocate(outside_points(2, 0))
@@ -99,8 +99,8 @@ program test_get_1_outerm
e%extr%coord_proj = ind_targ_extr
e%extr%coord = corner + (ind_targ_extr - 1) * step
call get_1_outerm(e%out_cont, cont_list, cont_list_proj, n_cont, cyclone, &
- real(ind_targ_extr), innermost_level, convert_to_ind(outside_points, &
- corner, step), ssh, corner = [1., 1.], step = step)
+ real(ind_targ_extr), innermost_level, real(outside_points), ssh, &
+ corner = [1., 1.], step = step)
if (e%out_cont%closed) then
e%extr%ssh = ssh(ind_targ_extr(1), ind_targ_extr(2))
diff --git a/Inst_eddies/inst_eddies.f90 b/Inst_eddies/inst_eddies.f90
index cc173a18..cb146bee 100644
--- a/Inst_eddies/inst_eddies.f90
+++ b/Inst_eddies/inst_eddies.f90
@@ -7,13 +7,12 @@ program inst_eddies
use, intrinsic:: ISO_FORTRAN_ENV, only: ERROR_UNIT
! Libraries:
- use contour_531, only: convert_to_ind, convert_to_reg_coord
- use jumble, only: assert, new_unit, argwhere, twopi
+ use jumble, only: assert, new_unit, argwhere
use numer_rec_95, only: indexx
use config_m, only: config
use derived_types, only: snapshot, shpc_slice_handler
- use input_ssh_m, only: input_ssh, max_radius, corner
+ use input_ssh_m, only: input_ssh, max_radius
use nearby_extr_m, only: nearby_extr
use set_all_extr_m, only: set_all_extr
use set_contours_m, only: set_contours
@@ -59,12 +58,9 @@ program inst_eddies
integer llc(2) ! indices in global grid of lower left corner
integer urc(2) ! indices in global grid of upper right corner
- real corner_window(2) ! longitude and latitude of the window around each
- ! extremum, in rad
-
- real, allocatable:: outside_points(:, :) ! (2, :) longitude and
- ! latitude, in rad, of all the significant extrema, except the
- ! target extremum
+ integer, allocatable:: outside_points(:, :) ! (2, :)
+ ! coordinates in projection space of all the significant extrema,
+ ! except the target extremum
!--------------------------------------------------------------
@@ -137,7 +133,7 @@ program inst_eddies
! Done sorting
- if (.not. periodic) nlon = size(ssh, 1)
+ nlon = size(ssh, 1) - 2 * merge(max_radius(1), 0, periodic)
nlat = size(ssh, 2)
do l = 1, s%number_extr
@@ -156,25 +152,21 @@ program inst_eddies
urc(1) = min(urc(1), nlon)
end if
- corner_window = corner + (llc - 1) * step
-
! Done defining geographical window
outside_points &
- = convert_to_reg_coord(nearby_extr(s%extr_map(llc(1):urc(1), &
- llc(2):urc(2)), s%list, i), corner, step)
+ = nearby_extr(s%extr_map(llc(1):urc(1), llc(2):urc(2)), s%list, i)
- ! Shift the longitude of each point to a value close to the
- ! longitude of the target extremum:
+ ! Shift the x coordinate of each point to a value close to the
+ ! coordinate of the target extremum:
if (periodic) outside_points(1, :) = outside_points(1, :) &
- + ceiling((corner_window(1) - outside_points(1, :)) / twopi) &
- * twopi
+ + ceiling((llc(1) - outside_points(1, :)) / real(nlon)) * nlon
call set_contours(s%list(i)%out_cont, s%list(i)%speed_cont, &
s%list(i)%max_speed, s%list(i)%cyclone, s%list(i)%extr, &
s%list(i)%innermost_level, step, ssh(llc(1):urc(1), llc(2):urc(2)), &
u(llc(1):urc(1), llc(2):urc(2)), v(llc(1):urc(1), llc(2):urc(2)), &
- llc, convert_to_ind(outside_points, corner, step))
+ llc, outside_points)
end do
call cpu_time(t1)
diff --git a/Inst_eddies/set_contours.f90 b/Inst_eddies/set_contours.f90
index 8b18750c..50b50746 100644
--- a/Inst_eddies/set_contours.f90
+++ b/Inst_eddies/set_contours.f90
@@ -46,7 +46,7 @@ contains
! corresponding to ssh(1, 1). A priori, this is not the corner of
! the grid of the NetCDF file.
- real, intent(in):: outside_points_proj(:, :) ! (2, :)
+ integer, intent(in):: outside_points_proj(:, :) ! (2, :)
! coordinates in projection space of all the significant extrema,
! except the target extremum
@@ -94,8 +94,8 @@ contains
innermost_level, abs(extr%ssh - innermost_level) < min_amp)
call get_1_outerm(out_cont, cont_list, cont_list_proj, n_cont, cyclone, &
- real(extr%coord_proj), innermost_level_2, outside_points_proj, ssh, &
- real(corner), step)
+ real(extr%coord_proj), innermost_level_2, real(outside_points_proj), &
+ ssh, real(corner), step)
! Done with outermost contour, now let us take care of the
! max-speed contour.
--
GitLab