diff --git a/Inst_eddies/inst_eddies.f90 b/Inst_eddies/inst_eddies.f90
index a01f0c55bf33578148a08463d2a7abb54192a897..0ead454cc1f4277cf23b0cceddd731e46a23e2fb 100644
--- a/Inst_eddies/inst_eddies.f90
+++ b/Inst_eddies/inst_eddies.f90
@@ -7,9 +7,10 @@ program inst_eddies
use, intrinsic:: ISO_FORTRAN_ENV
! Libraries:
- use jumble, only: assert, new_unit
+ use jumble, only: assert, new_unit, pi, argwhere
+ use numer_rec_95, only: indexx
- use config_m, only: config
+ use config_m, only: config, min_radius
use derived_types, only: snapshot, shpc_slice_handler
use input_ssh_m, only: input_ssh
use set_all_extr_m, only: set_all_extr
@@ -23,7 +24,7 @@ program inst_eddies
type(snapshot) s
TYPE(shpc_slice_handler) hshpc_cyclo, hshpc_anti
- integer iostat
+ integer iostat, i, l
character(len = 200) iomsg
integer unit
@@ -47,6 +48,18 @@ program inst_eddies
integer:: date = 0, slice = 0
namelist /dates_nml/ date, slice
+ integer, allocatable:: sorted_extr(:) ! (s%number_extr)
+ ! Sorted identifying number of extrema: first, minima sorted in
+ ! order of decreasing SSH, and second, maxima sorted in order of
+ ! increasing SSH.
+
+ real min_area ! minimum area of an outermost contour, in m2
+
+ integer, allocatable:: selection(:)
+ ! identifying numbers of a selection of extrema
+
+ integer n_cycl ! number of cyclones
+
!--------------------------------------------------------------
call assert(IEEE_SUPPORT_DATATYPE(0.), ieee_support_nan(0.), &
@@ -96,7 +109,33 @@ program inst_eddies
"s"
t0 = t1
call set_all_extr(s, step, periodic, ssh, corner)
- call set_all_contours(s, step, periodic, ssh, u, v, corner)
+
+ allocate(sorted_extr(s%number_extr))
+ min_area = pi * (min_radius * 1e3)**2
+
+ forall (i = 1:s%number_extr)
+ s%list(i)%out_cont%closed = .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(s%list%cyclone)
+ n_cycl = size(selection)
+ selection = selection(indexx(s%list(selection)%extr%ssh))
+ sorted_extr(:n_cycl) = selection(n_cycl:1:- 1) ! descending order of ssh
+
+ selection = argwhere(.not. s%list%cyclone)
+ selection = selection(indexx(s%list(selection)%extr%ssh))
+ sorted_extr(n_cycl + 1:) = selection
+
+ ! Done sorting
+
+ loop_extr: do l = 1, s%number_extr
+ i = sorted_extr(l)
+ call set_all_contours(s, step, periodic, ssh, u, v, corner, min_area, i)
+ end do loop_extr
+
call cpu_time(t1)
write(unit, fmt = *) "CPU time for computation, before output:", t1 - t0, "s"
t0 = t1
diff --git a/Inst_eddies/set_all_contours.f90 b/Inst_eddies/set_all_contours.f90
index b6cbc9eeac46eac869d4ba549d40169f2b77dca9..c1be9601afdb855acfada911bdd626506f0dffce 100644
--- a/Inst_eddies/set_all_contours.f90
+++ b/Inst_eddies/set_all_contours.f90
@@ -4,16 +4,15 @@ module set_all_contours_m
contains
- subroutine set_all_contours(s, step, periodic, ssh, u, v, corner)
+ subroutine set_all_contours(s, step, periodic, ssh, u, v, corner, min_area, i)
! This procedure sets all contours in the snapshot "s".
! Libraries:
use contour_531, only: convert_to_ind
- use jumble, only: argwhere, twopi, pi
- use numer_rec_95, only: indexx
+ use jumble, only: twopi
- use config_m, only: min_amp, min_radius
+ 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
@@ -37,21 +36,13 @@ contains
real, intent(in):: corner(:) ! (2) longitude and latitude of the
! corner of the global grid, in rad
+ real, intent(in):: min_area ! minimum area of an outermost contour, in m2
+ integer, intent(in):: i
+
! Local:
- real min_area ! minimum area of an outermost contour, in m2
integer nlon, nlat
real innermost_level_2
- integer i, l
- integer n_cycl ! number of cyclones
-
- integer, allocatable:: sorted_extr(:) ! (s%number_extr)
- ! Sorted identifying number of extrema: first, minima sorted in
- ! order of decreasing SSH, and second, maxima sorted in order of
- ! increasing SSH.
-
- integer, allocatable:: selection(:)
- ! identifying numbers of a selection of extrema
! Window around each extremum:
@@ -83,88 +74,64 @@ contains
if (.not. periodic) nlon = size(ssh, 1)
nlat = size(ssh, 2)
- allocate(sorted_extr(s%number_extr))
- min_area = pi * (min_radius * 1e3)**2
- forall (i = 1:s%number_extr)
- s%list(i)%out_cont%closed = .true.
- ! must be intialized to true because it is used in nearby_extr
- end forall
+ associate (e => s%list(i))
+ ! Define the geographical window around each eddy extremum:
- ! Double sort of extrema on cyclonicity and SSH value:
+ llc = e%extr%coord_proj - max_radius
+ urc = e%extr%coord_proj + max_radius
- selection = argwhere(s%list%cyclone)
- n_cycl = size(selection)
- selection = selection(indexx(s%list(selection)%extr%ssh))
- sorted_extr(:n_cycl) = selection(n_cycl:1:- 1) ! descending order of ssh
+ llc(2) = max(llc(2), 1)
+ urc(2) = min(urc(2), nlat)
- selection = argwhere(.not. s%list%cyclone)
- selection = selection(indexx(s%list(selection)%extr%ssh))
- sorted_extr(n_cycl + 1:) = selection
+ if (.not. periodic) then
+ llc(1) = max(llc(1), 1)
+ urc(1) = min(urc(1), nlon)
+ end if
- ! Done sorting
+ corner_window = corner + (llc - 1) * step
- loop_extr: do l = 1, s%number_extr
- i = sorted_extr(l)
+ ! 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)
- associate (e => s%list(i))
- ! Define the geographical window around each eddy extremum:
+ outside_points = nearby_extr(s%extr_map(llc(1):urc(1), &
+ llc(2):urc(2)), s%list, i)
- llc = e%extr%coord_proj - max_radius
- urc = e%extr%coord_proj + max_radius
+ ! 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
- 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
+ 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)
- ! 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)
+ ! Done with outermost contour, now let us take care of the
+ ! max-speed contour.
+ if (e%out_cont%closed) then
+ ! Restrict the field to the outermost contour:
- outside_points = nearby_extr(s%extr_map(llc(1):urc(1), &
- llc(2):urc(2)), s%list, i)
+ llc = floor(convert_to_ind(minval(e%out_cont%points, dim = 2), &
+ corner, step))
- ! 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
+ urc = ceiling(convert_to_ind(maxval(e%out_cont%points, dim = 2), &
+ corner, step))
- 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)
-
- ! Done with outermost contour, now let us take care of the
- ! max-speed contour.
- if (e%out_cont%closed) then
- ! Restrict the field to the outermost contour:
-
- llc = floor(convert_to_ind(minval(e%out_cont%points, dim = 2), &
- corner, step))
-
- urc = ceiling(convert_to_ind(maxval(e%out_cont%points, dim = 2), &
- corner, step))
-
- ! Should have no effect except because of roundup error:
- urc(2) = min(urc(2), nlat)
- if (.not. periodic) urc(1) = min(urc(1), nlon)
+ ! Should have no effect except because of roundup error:
+ urc(2) = min(urc(2), nlat)
+ if (.not. periodic) urc(1) = min(urc(1), nlon)
- corner_window = corner + (llc - 1) * step
+ corner_window = corner + (llc - 1) * step
- call set_max_speed(e, cont_list, n_cont, &
- 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, step)
- end if
- end associate
- end do loop_extr
+ call set_max_speed(e, cont_list, n_cont, &
+ 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, step)
+ end if
+ end associate
end subroutine set_all_contours