diff --git a/Inst_eddies/CMakeLists.txt b/Inst_eddies/CMakeLists.txt
index d3808689e87d0b36a45f3e7a01c6a67be99e2d1d..f3cd88697b56f951650a4b9015023e2b8567d551 100644
--- a/Inst_eddies/CMakeLists.txt
+++ b/Inst_eddies/CMakeLists.txt
@@ -3,7 +3,7 @@ add_executable(inst_eddies inst_eddies.f90 local_extrema.f90 set_max_speed.f90
get_1_outerm.f90 max_speed_contour_ssh.f90 good_contour.f90 mean_speed.f90
inside_4.f90 set_all_outerm.f90 nearby_extr.f90 get_var.f90 config.f90
input_ssh.f90 shpc_create.f90 write_snapshot.f90 ccw_orient.f90
- write_eddy.f90)
+ write_eddy.f90 set_all_extr.f90)
target_link_libraries(inst_eddies PRIVATE Contour_531::contour_531
Geometry::geometry NetCDF95::netcdf95 Shapelib_03::shapelib_03
Numer_Rec_95::numer_rec_95 Jumble::jumble NetCDF_Fortran::netcdff)
@@ -24,7 +24,7 @@ target_sources(test_get_1_outerm PRIVATE get_1_outerm.f90 good_contour.f90
target_sources(test_set_all_outerm PRIVATE set_all_outerm.f90 local_extrema.f90
get_1_outerm.f90 good_contour.f90 nearby_extr.f90 get_var.f90 config.f90
input_ssh.f90 shpc_create.f90 write_snapshot.f90 ccw_orient.f90
- write_eddy.f90)
+ write_eddy.f90 set_all_extr.f90)
target_sources(test_good_contour PRIVATE good_contour.f90)
target_sources(test_inside_4 PRIVATE inside_4.f90)
target_sources(test_mean_speed PRIVATE mean_speed.f90)
diff --git a/Inst_eddies/Tests/test_set_all_outerm.f90 b/Inst_eddies/Tests/test_set_all_outerm.f90
index ee48f0a229069b4f88f14cffac82b1a8471f708e..3336eb01124a7db2b2f7432c26d20f4dbf0c542b 100644
--- a/Inst_eddies/Tests/test_set_all_outerm.f90
+++ b/Inst_eddies/Tests/test_set_all_outerm.f90
@@ -7,6 +7,7 @@ program test_set_all_outerm
use derived_types, only: snapshot, shpc_slice_handler, null_ssh_contour, &
missing_speed
use input_ssh_m, only: input_ssh
+ use set_all_extr_m, only: set_all_extr
use set_all_outerm_m, only: set_all_outerm
use shpc_close_m, only: shpc_close
use shpc_create_m, only: shpc_create
@@ -33,11 +34,17 @@ program test_set_all_outerm
logical periodic ! grid is periodic in longitude
TYPE(shpc_slice_handler) hshp_cyclo, hshp_anti
+ real, allocatable:: innermost_level(:) ! (s%number_extr)
+ ! SSH level of the innermost contour around each extremum, in
+ ! m. By construction, innermost_level < extremum for a maximum, >
+ ! extremum for a minimum.
+
!--------------------------------------------------------------
call config
call input_ssh(corner, step, nlon, nlat, periodic, ssh, u, v)
- call set_all_outerm(s, step, periodic, ssh, corner, &
+ call set_all_extr(s, innermost_level, step, periodic, ssh, corner)
+ call set_all_outerm(s, step, periodic, ssh, corner, innermost_level, &
min_area = pi * (min_radius * 1e3)**2)
do i = 1, s%number_extr
diff --git a/Inst_eddies/inst_eddies.f90 b/Inst_eddies/inst_eddies.f90
index f0d51db36e861eff06e493d76e9636fdecc020b3..865eab67b7c597a3cd9e9d2e37f8e8ebdc803cb7 100644
--- a/Inst_eddies/inst_eddies.f90
+++ b/Inst_eddies/inst_eddies.f90
@@ -12,6 +12,7 @@ program inst_eddies
shpc_slice_handler
use input_ssh_m, only: input_ssh
use nearby_extr_m, only: nearby_extr
+ use set_all_extr_m, only: set_all_extr
use set_all_outerm_m, only: set_all_outerm
use set_max_speed_m, only: set_max_speed
use shpc_close_m, only: shpc_close
@@ -56,6 +57,11 @@ program inst_eddies
! latitude, in rad, of all the significant extrema, except the
! target extremum
+ real, allocatable:: innermost_level(:) ! (s%number_extr)
+ ! SSH level of the innermost contour around each extremum, in
+ ! m. By construction, innermost_level < extremum for a maximum, >
+ ! extremum for a minimum.
+
logical exist
real t0, t1 ! CPU times, in s
integer:: date = 0, slice = 0
@@ -100,7 +106,8 @@ program inst_eddies
write(unit, fmt = *) "CPU time for preamble, before computation:", t1 - t0, &
"s"
t0 = t1
- call set_all_outerm(s, step, periodic, ssh, corner, &
+ call set_all_extr(s, innermost_level, step, periodic, ssh, corner)
+ call set_all_outerm(s, step, periodic, ssh, corner, innermost_level, &
min_area = pi * (min_radius * 1e3)**2)
! Done with outermost contours, now let us take care of the
diff --git a/Inst_eddies/set_all_extr.f90 b/Inst_eddies/set_all_extr.f90
new file mode 100644
index 0000000000000000000000000000000000000000..6eba52244b20f194330f49f40587168e365b55f1
--- /dev/null
+++ b/Inst_eddies/set_all_extr.f90
@@ -0,0 +1,80 @@
+module set_all_extr_m
+
+ implicit none
+
+contains
+
+ subroutine set_all_extr(s, innermost_level, step, periodic, ssh, corner)
+
+ ! This procedures finds all the extrema of ssh. Not a function
+ ! because "s" is not completely defined.
+
+ use derived_types, only: snapshot
+ use input_ssh_m, only: max_radius
+ use local_extrema_m, only: local_extrema
+
+ type(snapshot), intent(out):: s
+ ! Define only s%extr_map, s%ind_extr, s%number_extr,
+ ! s%list%coord_extr, s%list%ssh_extr, s%list%cyclone
+
+ real, allocatable, intent(out):: innermost_level(:) ! (s%number_extr)
+ ! SSH level of the innermost contour around each extremum, in
+ ! m. By construction, innermost_level < extremum for a maximum, >
+ ! extremum for a minimum.
+
+ real, intent(in):: step(:) ! (2) longitude and latitude steps, in rad
+ logical, intent(in):: periodic ! grid is periodic in longitude
+
+ real, intent(in):: ssh(1 - merge(max_radius(1), 0, periodic):, :)
+ ! (1 - max_radius(1):nlon + max_radius(1), nlat) if the grid is periodic
+ ! in longitude, else (nlon, nlat). Sea-surface height, in m.
+
+ real, intent(in):: corner(:) ! (2) longitude and latitude of the
+ ! corner of the global grid, in rad
+
+ ! Local:
+
+ integer nlon, nlat
+
+ logical, allocatable:: cyclone(:) ! (s%number_extr)
+ integer i, copy
+
+ !--------------------------------------------------------------
+
+ copy = merge(max_radius(1), 0, periodic)
+ nlon = size(ssh, 1) - 2 * copy
+ nlat = size(ssh, 2)
+ allocate(s%extr_map(1 - copy:nlon + copy, nlat))
+ copy = merge(1, 0, periodic)
+ call local_extrema(s%extr_map(1 - copy:nlon + copy, :), s%ind_extr, &
+ innermost_level, cyclone, ssh(1 - copy:nlon + copy, :), periodic, &
+ my_lbound = [1 - copy, 1])
+ ! Assertion: \forall k, s%ind_extr(1, k) \in \{2 - copy, \dots,
+ ! nlon + copy - 1\} and s%ind_extr(2, k) \in \{2, \dots, nlat - 1\}
+
+ ! Assertion: \forall k, s%ind_extr(1, k) \in \{1, \dots, nlon\} and
+ ! s%ind_extr(2, k) \in \{2, \dots, nlat - 1\}
+
+ if (periodic) then
+ ! Extend in longitude:
+ s%extr_map(1 - max_radius(1):- 1, :) &
+ = s%extr_map(nlon + 1 - max_radius(1):nlon - 1, :)
+ s%extr_map(nlon + 2:nlon + max_radius(1), :) &
+ = s%extr_map(2:max_radius(1), :)
+ end if
+
+ s%number_extr = size(s%ind_extr, 2)
+ allocate(s%list(s%number_extr))
+
+ forall (i = 1:s%number_extr)
+ s%list(i)%coord_extr = corner + (s%ind_extr(:, i) - 1) * step
+ ! (Even when periodic, this is within the original NetCDF grid,
+ ! that is the grid without duplicated longitudes.)
+
+ s%list(i)%ssh_extr = ssh(s%ind_extr(1, i), s%ind_extr(2, i))
+ s%list(i)%cyclone = cyclone(i)
+ end forall
+
+ end subroutine set_all_extr
+
+end module set_all_extr_m
diff --git a/Inst_eddies/set_all_outerm.f90 b/Inst_eddies/set_all_outerm.f90
index 99f13c4b2b0fdb88dd3f2630e26f9aa7534d1480..eaa543b97ef3fe00a35549acbb9df9ff78ba481c 100644
--- a/Inst_eddies/set_all_outerm.f90
+++ b/Inst_eddies/set_all_outerm.f90
@@ -4,11 +4,11 @@ module set_all_outerm_m
contains
- subroutine set_all_outerm(s, step, periodic, ssh, corner, min_area)
+ subroutine set_all_outerm(s, step, periodic, ssh, corner, innermost_level, &
+ min_area)
- ! Extraction of eddies: find extrema and set all outermost
- ! contours in the snapshot "s". Not a function because "s" is not
- ! completely defined.
+ ! This procedure sets all outermost contours in the snapshot
+ ! "s".
! Libraries:
use jumble, only: argwhere
@@ -20,11 +20,9 @@ contains
use get_1_outerm_m, only: get_1_outerm
use input_ssh_m, only: max_radius
use nearby_extr_m, only: nearby_extr
- use local_extrema_m, only: local_extrema
- type(snapshot), intent(out):: s
- ! Specifically: define everything in s except
- ! s%list%speed_cont, s%list%max_speed and s%list%radius4.
+ type(snapshot), intent(inout):: s
+ ! Define s%list%valid, s%list%out_cont
real, intent(in):: step(:) ! (2) longitude and latitude steps, in rad
logical, intent(in):: periodic ! grid is periodic in longitude
@@ -36,21 +34,19 @@ contains
real, intent(in):: corner(:) ! (2) longitude and latitude of the
! corner of the global grid, in rad
+ real, intent(in):: innermost_level(:) ! (s%number_extr)
+ ! SSH level of the innermost contour around each extremum, in
+ ! m. By construction, innermost_level < extremum for a maximum, >
+ ! extremum for a minimum.
+
real, intent(in):: min_area
! minimum area of an outermost contour, in m2
-
+
! Local:
integer nlon, nlat
-
- real, allocatable:: innermost_level(:) ! (s%number_extr)
- ! SSH level of the innermost contour around each extremum, in
- ! m. By construction, innermost_level < extremum for a maximum, >
- ! extremum for a minimum.
-
real innermost_level_2
- logical, allocatable:: cyclone(:) ! (s%number_extr)
- integer i, l, copy
+ integer i, l
integer n_cycl ! number of cyclones
integer, allocatable:: sorted_extr(:) ! (s%number_extr)
@@ -75,51 +71,24 @@ contains
!--------------------------------------------------------------
- copy = merge(max_radius(1), 0, periodic)
- nlon = size(ssh, 1) - 2 * copy
+ if (.not. periodic) nlon = size(ssh, 1)
nlat = size(ssh, 2)
- allocate(s%extr_map(1 - copy:nlon + copy, nlat))
- copy = merge(1, 0, periodic)
- call local_extrema(s%extr_map(1 - copy:nlon + copy, :), s%ind_extr, &
- innermost_level, cyclone, ssh(1 - copy:nlon + copy, :), periodic, &
- my_lbound = [1 - copy, 1])
- ! Assertion: \forall k, s%ind_extr(1, k) \in \{2 - copy, \dots,
- ! nlon + copy - 1\} and s%ind_extr(2, k) \in \{2, \dots, nlat - 1\}
-
- ! Assertion: \forall k, s%ind_extr(1, k) \in \{1, \dots, nlon\} and
- ! s%ind_extr(2, k) \in \{2, \dots, nlat - 1\}
-
- if (periodic) then
- ! Extend in longitude:
- s%extr_map(1 - max_radius(1):- 1, :) &
- = s%extr_map(nlon + 1 - max_radius(1):nlon - 1, :)
- s%extr_map(nlon + 2:nlon + max_radius(1), :) &
- = s%extr_map(2:max_radius(1), :)
- end if
-
- s%number_extr = size(s%ind_extr, 2)
- allocate(s%list(s%number_extr), sorted_extr(s%number_extr))
- forall (i = 1:s%number_extr)
- s%list(i)%coord_extr = corner + (s%ind_extr(:, i) - 1) * step
- ! (Even when periodic, this is within the original NetCDF grid,
- ! that is the grid without duplicated longitudes.)
-
- s%list(i)%ssh_extr = ssh(s%ind_extr(1, i), s%ind_extr(2, i))
- s%list(i)%cyclone = cyclone(i)
+ allocate(sorted_extr(s%number_extr))
+ forall (i = 1:s%number_extr)
s%list(i)%valid = .true.
! must be intialized to true because it is used in nearby_extr
end forall
! Double sort of extrema on cyclonicity and SSH value:
- selection = argwhere(cyclone)
+ selection = argwhere(s%list%cyclone)
n_cycl = size(selection)
selection = selection(indexx(s%list(selection)%ssh_extr))
sorted_extr(:n_cycl) = selection(n_cycl:1:- 1) ! descending order of ssh
- selection = argwhere(.not. cyclone)
+ selection = argwhere(.not. s%list%cyclone)
selection = selection(indexx(s%list(selection)%ssh_extr))
sorted_extr(n_cycl + 1:) = selection