From e1b9ac7a01559c3c311a922927baf1ef7d4f0532 Mon Sep 17 00:00:00 2001
From: Lionel GUEZ <guez@lmd.ens.fr>
Date: Thu, 25 May 2023 22:00:58 +0200
Subject: [PATCH] Move some computation to main program unit
Move the computation of the search window for the outermost contour
and `outside_points` from procedure `set_contours` to main program
unit. The motivation is to have a clearer interface of procedure
`set_contours`, without `s` and `i`.
---
Inst_eddies/inst_eddies.f90 | 50 +++++++++++++++++++++++++++++++++---
Inst_eddies/set_contours.f90 | 47 ++++++++-------------------------
2 files changed, 56 insertions(+), 41 deletions(-)
diff --git a/Inst_eddies/inst_eddies.f90 b/Inst_eddies/inst_eddies.f90
index 4a267803..a5768244 100644
--- a/Inst_eddies/inst_eddies.f90
+++ b/Inst_eddies/inst_eddies.f90
@@ -7,12 +7,13 @@ program inst_eddies
use, intrinsic:: ISO_FORTRAN_ENV
! Libraries:
- use jumble, only: assert, new_unit, pi, argwhere
+ use jumble, only: assert, new_unit, pi, argwhere, twopi
use numer_rec_95, only: indexx
use config_m, only: config, min_radius
use derived_types, only: snapshot, shpc_slice_handler
- use input_ssh_m, only: input_ssh
+ 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
use shpc_close_m, only: shpc_close
@@ -24,7 +25,7 @@ program inst_eddies
type(snapshot) s
TYPE(shpc_slice_handler) hshpc_cyclo, hshpc_anti
- integer iostat, i, l
+ integer iostat, i, l, nlon, nlat
character(len = 200) iomsg
integer unit
@@ -59,6 +60,16 @@ program inst_eddies
! identifying numbers of a selection of extrema
integer n_cycl ! number of cyclones
+ 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
+
!--------------------------------------------------------------
@@ -133,7 +144,38 @@ program inst_eddies
loop_extr: do l = 1, s%number_extr
i = sorted_extr(l)
- call set_contours(s, step, periodic, ssh, u, v, corner, min_area, i)
+ if (.not. periodic) nlon = size(ssh, 1)
+ nlat = size(ssh, 2)
+
+ associate (e => s%list(i))
+ ! Define the geographical window around each eddy extremum:
+
+ llc = e%extr%coord_proj - max_radius
+ urc = e%extr%coord_proj + max_radius
+
+ llc(2) = max(llc(2), 1)
+ urc(2) = min(urc(2), nlat)
+
+ if (.not. periodic) then
+ llc(1) = max(llc(1), 1)
+ urc(1) = min(urc(1), nlon)
+ end if
+
+ corner_window = corner + (llc - 1) * step
+ end associate
+
+ outside_points = 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:
+ if (periodic) outside_points(1, :) = outside_points(1, :) &
+ + ceiling((corner_window(1) - outside_points(1, :)) / twopi) &
+ * twopi
+
+ call set_contours(s, 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)), &
+ corner_window, min_area, i, outside_points)
end do loop_extr
call cpu_time(t1)
diff --git a/Inst_eddies/set_contours.f90 b/Inst_eddies/set_contours.f90
index 9f538292..7c809d51 100644
--- a/Inst_eddies/set_contours.f90
+++ b/Inst_eddies/set_contours.f90
@@ -4,19 +4,17 @@ module set_contours_m
contains
- subroutine set_contours(s, step, periodic, ssh, u, v, corner, min_area, i)
+ subroutine set_contours(s, step, ssh, u, v, corner, min_area, i, &
+ outside_points)
! This procedure sets all contours in the snapshot "s".
! Libraries:
use contour_531, only: convert_to_ind
- use jumble, only: twopi
use config_m, only: min_amp
use derived_types, only: snapshot, null_ssh_contour, ssh_contour
use get_1_outerm_m, only: get_1_outerm
- use input_ssh_m, only: max_radius
- use nearby_extr_m, only: nearby_extr
use set_max_speed_m, only: set_max_speed
type(snapshot), intent(inout):: s
@@ -25,9 +23,8 @@ contains
! s%list%innermost_level, s%extr_map should be defined on entry.
real, intent(in):: step(:) ! (2) longitude and latitude steps, in rad
- logical, intent(in):: periodic ! grid is periodic in longitude
- real, intent(in), dimension(merge(1 - max_radius(1), 1, periodic):, :):: &
+ real, intent(in), dimension(:, :):: &
ssh, u, v
! (1 - max_radius(1):nlon + max_radius(1), nlat) if the grid is
! periodic in longitude, else (nlon, nlat). Sea-surface height, in
@@ -39,6 +36,10 @@ contains
real, intent(in):: min_area ! minimum area of an outermost contour, in m2
integer, intent(in):: i
+ real, intent(in):: outside_points(:, :) ! (2, :) longitude and
+ ! latitude, in rad, of all the significant extrema, except the
+ ! target extremum
+
! Local:
integer nlon, nlat
@@ -52,10 +53,6 @@ contains
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, parameter:: n_max_cont = 31 ! must be >= 3
type(ssh_contour) cont_list(n_max_cont)
@@ -72,42 +69,18 @@ contains
!--------------------------------------------------------------
- if (.not. periodic) nlon = size(ssh, 1)
+ nlon = size(ssh, 1)
nlat = size(ssh, 2)
associate (e => s%list(i))
- ! Define the geographical window around each eddy extremum:
-
- llc = e%extr%coord_proj - max_radius
- urc = e%extr%coord_proj + max_radius
-
- llc(2) = max(llc(2), 1)
- urc(2) = min(urc(2), nlat)
-
- if (.not. periodic) then
- llc(1) = max(llc(1), 1)
- urc(1) = min(urc(1), nlon)
- end if
-
- corner_window = corner + (llc - 1) * step
-
! Only look for good contours with amplitudes >= min_amp:
innermost_level_2 = merge(e%extr%ssh + merge(min_amp, - min_amp, &
e%cyclone), e%innermost_level, &
abs(e%extr%ssh - e%innermost_level) < min_amp)
- outside_points = 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:
- if (periodic) outside_points(1, :) = outside_points(1, :) &
- + ceiling((corner_window(1) - outside_points(1, :)) / twopi) &
- * twopi
-
call get_1_outerm(e%out_cont, cont_list, n_cont, e%cyclone, &
e%extr%coord, innermost_level_2, outside_points, &
- ssh(llc(1):urc(1), llc(2):urc(2)), corner_window, step, min_area)
+ ssh, corner, step, min_area)
! Done with outermost contour, now let us take care of the
! max-speed contour.
@@ -122,7 +95,7 @@ contains
! Should have no effect except because of roundup error:
urc(2) = min(urc(2), nlat)
- if (.not. periodic) urc(1) = min(urc(1), nlon)
+ urc(1) = min(urc(1), nlon)
corner_window = corner + (llc - 1) * step
--
GitLab