diff --git a/GNUmakefile b/GNUmakefile
index c7b8e2ce25df0eafd6edcdb30b0c230c2b1ee9cc..22d0f8d57204573adbcdb8585962d7de6e38d4bb 100644
--- a/GNUmakefile
+++ b/GNUmakefile
@@ -13,6 +13,8 @@ src_test_set_max_speed = test_set_max_speed.f derived_types.f set_max_speed.f go
src_test_get_snapshot = test_get_snapshot.f get_snapshot.f dispatch_snapshot.f write_eddy.f send_snapshot.f receive_snapshot.f local_extrema.f set_max_speed.f outermost_possible_level.f get_1_outerm.f max_speed_contour_ssh.f good_contour.f spherical_polygon_area.f mean_speed.f inside_4.f set_all_outerm.f
+src_test_set_all_outerm = test_set_all_outerm.f derived_types.f dispatch_snapshot.f set_all_outerm.f write_eddy.f send_snapshot.f local_extrema.f get_1_outerm.f good_contour.f spherical_polygon_area.f
+
sources := $(shell cat ${makefile_dir}/file_list)
lib_list = contour_531 numer_rec_95 shapelib_03 netcdf95 geometry jumble netcdff fortrangis shp fortranc nr_util
@@ -23,9 +25,10 @@ obj_test_local_extrema := $(src_test_local_extrema:.f=.o)
obj_test_get_1_outerm := $(src_test_get_1_outerm:.f=.o)
obj_test_set_max_speed := $(src_test_set_max_speed:.f=.o)
obj_test_get_snapshot := $(src_test_get_snapshot:.f=.o)
+obj_test_set_all_outerm := $(src_test_set_all_outerm:.f=.o)
objects := $(sources:.f=.o)
-execut = test_good_contour test_inside_4 test_get_1_outerm test_local_extrema test_max_speed_contour_ssh test_mean_speed test_set_max_speed test_get_snapshot
+execut = test_good_contour test_inside_4 test_get_1_outerm test_local_extrema test_max_speed_contour_ssh test_mean_speed test_set_max_speed test_get_snapshot test_set_all_outerm
# 3. Compiler-dependent part
@@ -47,6 +50,7 @@ test_inside_4: inside_4.o
test_local_extrema: ${obj_test_local_extrema}
test_mean_speed: mean_speed.o
test_get_snapshot: ${obj_test_get_snapshot}
+test_set_all_outerm: ${obj_test_set_all_outerm}
depend ${makefile_dir}/depend.mk:
makedepf90 -free -Wmissing -Wconfused $(addprefix -I, ${VPATH}) -nosrc $(addprefix -u , ${lib_list} shapelib netcdf) ${sources} >${makefile_dir}/depend.mk
diff --git a/Sources/get_snapshot.f b/Sources/get_snapshot.f
index 48c598b5be1c4722091b467ed96ad5c4f4427a77..f3d53289fc44e980fb8bfb0c51f08ac5d9d31cbd 100644
--- a/Sources/get_snapshot.f
+++ b/Sources/get_snapshot.f
@@ -9,7 +9,6 @@ contains
use contour_531, only: convert_to_ind, null_polyline
use derived_types, only: snapshot
- use local_extrema_m, only: local_extrema
use netcdf, only: nf90_nowrite
use netcdf95, only: nf95_open, nf95_close, nf95_inq_varid, nf95_get_var
use receive_snapshot_m, only: receive_snapshot
@@ -43,10 +42,6 @@ contains
real ssh(nlon, nlat) ! sea-surface height, in m
real u(nlon, nlat), v(nlon, nlat) ! wind, in m s-1
- real, allocatable:: innermost_level(:) ! (s%number_vis_eddies)
- ! level of innermost contour, for each extremum
-
- logical, allocatable:: cyclone(:) ! (s%number_vis_eddies)
integer i
! Window around each extremum:
@@ -70,22 +65,7 @@ contains
call nf95_get_var(ncid, varid, v, start = [1, 1, k])
call nf95_close(ncid)
- ! Extraction of eddies:
-
- allocate(s%extr_map(nlon, nlat))
- call local_extrema(ssh, s%extr_map, s%ind_extr, innermost_level, cyclone)
- s%number_vis_eddies = size(s%ind_extr, 2)
- allocate(s%list_vis(s%number_vis_eddies))
-
- forall (i = 1:s%number_vis_eddies)
- s%list_vis(i)%coord_extr = corner + (s%ind_extr(:, i) - 1) * step
- s%list_vis(i)%ssh_extremum = ssh(s%ind_extr(1, i), s%ind_extr(2, i))
- s%list_vis(i)%cyclone = cyclone(i)
- s%list_vis(i)%interpolated = .false.
- end forall
-
- call set_all_outerm(s, min_amp, max_radius, corner, step, &
- innermost_level, ssh)
+ call set_all_outerm(s, min_amp, max_radius, corner, step, ssh)
! Done with outermost contours, now let us take care of the
! max-speed contours.
diff --git a/Sources/set_all_outerm.f b/Sources/set_all_outerm.f
index 94f030eb02ea6adf3fadbef0b6fec79e99a4d672..051a15361a905f78e37efe9caf68c2ab925a04c8 100644
--- a/Sources/set_all_outerm.f
+++ b/Sources/set_all_outerm.f
@@ -4,24 +4,20 @@ module set_all_outerm_m
contains
- subroutine set_all_outerm(s, min_amp, max_radius, corner, step, &
- innermost_level, ssh)
+ subroutine set_all_outerm(s, min_amp, max_radius, corner, step, ssh)
- ! Set all outermost contours.
+ ! Extraction of eddies: find extrema and set all outermost
+ ! contours in snapshot. Not a function because snapshot is not
+ ! completely defined.
use derived_types, only: snapshot
use get_1_outerm_m, only: get_1_outerm
+ use local_extrema_m, only: local_extrema
- type(snapshot), intent(inout):: s
-
- ! Intent(in): s%number_vis_eddies, s%ind_extr,
- ! s%list_vis%ssh_extremum, s%list_vis%cyclone,
- ! s%list_vis%coord_extr
-
- ! Intent(out): s%list_vis%outermost_contour, s%list_vis%suff_amp,
- ! s%list_vis%twice
-
- ! Intent(inout): s%extr_map
+ type(snapshot), intent(out):: s
+ ! Specifically: define everything in s except
+ ! s%list_vis%max_speed_contour, s%list_vis%max_speed and
+ ! s%number_eddies.
real, intent(in):: min_amp ! minimum amplitude of ssh, between
! extremum and outermost contour, in m
@@ -34,29 +30,45 @@ contains
real, intent(in):: step(:) ! (2) longitude and latitude steps, in rad
- real, intent(in):: innermost_level(:) ! (s%number_vis_eddies)
- ! level of innermost contour, for each extremum
-
real, intent(in):: ssh(:, :) ! sea-surface height, in m
! Local:
+ real, allocatable:: innermost_level(:) ! (s%number_vis_eddies)
+ ! level of innermost contour, for each extremum
+
+ logical, allocatable:: cyclone(:) ! (s%number_vis_eddies)
integer i
- logical flat_extr(s%number_vis_eddies) ! flat extremum
- logical noise_around(s%number_vis_eddies)
+ logical, allocatable:: flat_extr(:) ! (s%number_vis_eddies) flat extremum
+ logical, allocatable:: noise_around(:) ! (s%number_vis_eddies)
! Window around each extremum:
- integer llc(2, s%number_vis_eddies) ! indices in global grid of
- ! lower left corner
+ integer, allocatable:: llc(:, :) ! (2, s%number_vis_eddies)
+ ! indices in global grid of lower left corner
- integer urc(2, s%number_vis_eddies) ! indices in global grid of
- ! upper right corner
+ integer, allocatable:: urc(:, :) ! (2, s%number_vis_eddies)
+ ! indices in global grid of upper right corner
- real corner_window(2, s%number_vis_eddies) ! longitude and latitude, in rad
+ real, allocatable:: corner_window(:, :) ! (2, s%number_vis_eddies)
+ ! longitude and latitude, in rad
!--------------------------------------------------------------
+ allocate(s%extr_map(size(ssh, 1), size(ssh, 2)))
+ call local_extrema(ssh, s%extr_map, s%ind_extr, innermost_level, cyclone)
+ s%number_vis_eddies = size(s%ind_extr, 2)
+ allocate(s%list_vis(s%number_vis_eddies), flat_extr(s%number_vis_eddies), &
+ noise_around(s%number_vis_eddies), llc(2, s%number_vis_eddies), &
+ urc(2, s%number_vis_eddies), corner_window(2, s%number_vis_eddies))
+
+ forall (i = 1:s%number_vis_eddies)
+ s%list_vis(i)%coord_extr = corner + (s%ind_extr(:, i) - 1) * step
+ s%list_vis(i)%ssh_extremum = ssh(s%ind_extr(1, i), s%ind_extr(2, i))
+ s%list_vis(i)%cyclone = cyclone(i)
+ s%list_vis(i)%interpolated = .false.
+ end forall
+
! Define the geographical window around each eddy extremum:
forall (i = 1:s%number_vis_eddies)
llc(:, i) = max(s%ind_extr(:, i) - max_radius, 1)
diff --git a/Tests/test_set_all_outerm.f b/Tests/test_set_all_outerm.f
new file mode 100644
index 0000000000000000000000000000000000000000..cb38946c0667a4c09dbd7b3ed98160ebc4a36bd7
--- /dev/null
+++ b/Tests/test_set_all_outerm.f
@@ -0,0 +1,162 @@
+program test_set_all_outerm
+
+ use contour_531, only: null_polyline
+ use derived_types, only: snapshot
+ use dispatch_snapshot_m, only: dispatch_snapshot
+ use jumble, only: new_unit, get_command_arg_dyn
+ use netcdf, only: nf90_nowrite
+ use netcdf95, only: nf95_open, nf95_close, nf95_inq_varid, nf95_get_var, &
+ find_coord, nf95_inquire_dimension
+ use nr_util, only: pi
+ use set_all_outerm_m, only: set_all_outerm
+ use shapelib, only: shpt_point, shpt_polygon, shpfileobject, ftdouble, &
+ shpclose, ftinteger
+ use shapelib_03, only: shp_create_03, dbf_add_field_03
+
+ implicit none
+
+ type(snapshot) s
+
+ real:: min_amp = 0.
+ ! minimum amplitude of ssh, between extremum and outermost contour,
+ ! in m
+
+ character(len = :), allocatable:: filename
+
+ integer nlon, nlat
+
+ ! Just the first two values to get the corner and step:
+ real longitude(2), latitude(2) ! in degrees
+
+ integer, parameter:: max_radius(2) = [20, 20]
+ ! maximum radius of an eddy in longitude and latitude, in number of
+ ! grid points
+
+ real, parameter:: deg_over_rad = pi / 180.
+
+ integer ncid, varid, dimid, i
+ real, allocatable:: ssh(:, :) ! (nlon, nlat) sea-surface height, in m
+
+ TYPE(shpfileobject) hshp_extremum
+ ! shapefile extremum_$m. The fields in the DBF file are, in that
+ ! order: ssh, date index, eddy index, interpolated, cyclone, suff
+ ! amplitude.
+
+ TYPE(shpfileobject) hshp_outermost
+ ! shapefile outermost_contour_$m. The fields in the DBF file are,
+ ! in that order: area, ssh, date index, eddy index, twice.
+
+ TYPE(shpfileobject) hshp_max_speed
+ ! shapefile x_speed_contour_$m. The fields in the DBF file are, in
+ ! that order: area, ssh, spped, date index, eddy index.
+
+ integer ifield, unit_isolated, unit_number_eddies
+
+ !--------------------------------------------------------------
+
+ call get_command_arg_dyn(1, filename)
+
+ print *, "min_amp = ? "
+ read *, min_amp
+ print *, "min_amp = ", min_amp
+
+ ! Read ssh:
+
+ call nf95_open(filename, nf90_nowrite, ncid)
+
+ call find_coord(ncid, dimid = dimid, varid = varid, std_name = "longitude")
+ call nf95_inquire_dimension(ncid, dimid, nclen = nlon)
+ call nf95_get_var(ncid, varid, longitude)
+
+ call find_coord(ncid, dimid = dimid, varid = varid, std_name = "latitude")
+ call nf95_inquire_dimension(ncid, dimid, nclen = nlat)
+ call nf95_get_var(ncid, varid, latitude)
+
+ allocate(ssh(nlon, nlat))
+ call nf95_inq_varid(ncid, "adt", varid)
+ call nf95_get_var(ncid, varid, ssh)
+
+ call nf95_close(ncid)
+
+ call set_all_outerm(s, min_amp, max_radius, &
+ corner = [longitude(1), latitude(1)] * deg_over_rad, step &
+ = [longitude(2) - longitude(1), latitude(2) - latitude(1)] &
+ * deg_over_rad, ssh = ssh)
+
+ do i = 1, s%number_vis_eddies
+ s%list_vis(i)%max_speed_contour%polyline = null_polyline()
+ s%list_vis(i)%max_speed_contour%ssh = s%list_vis(i)%ssh_extremum
+ s%list_vis(i)%max_speed_contour%area = 0.
+ s%list_vis(i)%max_speed = 0.
+ end do
+
+ s%number_eddies = s%number_vis_eddies
+
+ call shp_create_03("extremum_1", shpt_point, hshp_extremum)
+ call dbf_add_field_03(ifield, hshp_extremum, 'ssh', ftdouble, nwidth = 13, &
+ ndecimals = 6)
+ call dbf_add_field_03(ifield, hshp_extremum, 'date_index', ftinteger, &
+ nwidth = 4, ndecimals = 0)
+ call dbf_add_field_03(ifield, hshp_extremum, 'eddy_index', ftinteger, &
+ nwidth = 5, ndecimals = 0)
+ call dbf_add_field_03(ifield, hshp_extremum, 'interpolat', ftinteger, &
+ nwidth = 1, ndecimals = 0)
+ call dbf_add_field_03(ifield, hshp_extremum, 'cyclone', ftinteger, &
+ nwidth = 1, ndecimals = 0)
+ call dbf_add_field_03(ifield, hshp_extremum, 'suff_amp', ftinteger, &
+ nwidth = 1, ndecimals = 0)
+
+ call shp_create_03("outermost_contour_1", shpt_polygon, hshp_outermost)
+ call dbf_add_field_03(ifield, hshp_outermost, 'area', ftdouble, nwidth = 20, &
+ ndecimals = 6)
+ call dbf_add_field_03(ifield, hshp_outermost, 'ssh', ftdouble, nwidth = 13, &
+ ndecimals = 6)
+ call dbf_add_field_03(ifield, hshp_outermost, 'date_index', ftinteger, &
+ nwidth = 4, ndecimals = 0)
+ call dbf_add_field_03(ifield, hshp_outermost, 'eddy_index', ftinteger, &
+ nwidth = 5, ndecimals = 0)
+ call dbf_add_field_03(ifield, hshp_outermost, 'twice', ftinteger, &
+ nwidth = 1, ndecimals = 0)
+
+ call shp_create_03("max_speed_contour_1", shpt_polygon, hshp_max_speed)
+ call dbf_add_field_03(ifield, hshp_max_speed, 'area', ftdouble, nwidth = 20, &
+ ndecimals = 6)
+ call dbf_add_field_03(ifield, hshp_max_speed, 'ssh', ftdouble, nwidth = 13, &
+ ndecimals = 6)
+ call dbf_add_field_03(ifield, hshp_max_speed, 'speed', &
+ ftdouble, nwidth = 13, ndecimals = 6)
+ call dbf_add_field_03(ifield, hshp_max_speed, 'date_index', ftinteger, &
+ nwidth = 4, ndecimals = 0)
+ call dbf_add_field_03(ifield, hshp_max_speed, 'eddy_index', ftinteger, &
+ nwidth = 5, ndecimals = 0)
+
+ call new_unit(unit_isolated)
+ open(unit_isolated, file = "isolated_nodes_1.csv", status = "replace", &
+ action = "write")
+ write(unit_isolated, fmt = *) '"date index" "eddy index"' ! title line
+
+ call new_unit(unit_number_eddies)
+ open(unit_number_eddies, file = "number_eddies_1.csv", status = "replace", &
+ action = "write")
+ write(unit_number_eddies, fmt = *) '"date index" ' &
+ // '"number of visible eddies" "number of interpolated eddies"'
+ ! (title line)
+
+ call dispatch_snapshot(s, hshp_extremum, hshp_outermost, hshp_max_speed, &
+ unit_isolated, unit_number_eddies, k = 1, m = 1, k_begin = 1, &
+ max_delta = 4)
+
+ CALL shpclose(hshp_extremum)
+ print *, 'Created shapefile "extremum_1".'
+ CALL shpclose(hshp_outermost)
+ print *, 'Created shapefile "outermost_contour_1".'
+ CALL shpclose(hshp_max_speed)
+ print *, 'Created shapefile "max_speed_contour_1".'
+ close(unit_isolated)
+ print *, 'Created "isolated_nodes_1.csv".'
+ close(unit_number_eddies)
+ print *, 'Created "number_eddies_1.csv".'
+ print *, "Average number of points per outermost contour: ", &
+ sum(s%list_vis%outermost_contour%n_points) / real(s%number_vis_eddies)
+
+end program test_set_all_outerm
diff --git a/Tests/tests.json b/Tests/tests.json
index fdab93ef5b30dec62a7312e2c57ae297163a9ae0..b53e606ea9dec2f5a713926502a86a5320cc640d 100644
--- a/Tests/tests.json
+++ b/Tests/tests.json
@@ -101,6 +101,12 @@
"$tests_old_dir/Get_1_outerm_noise_2_8/test_get_1_outerm.dbf"],
"input": "&MAIN_NML IND_TARG_EXTR=19,11 /\n"
},
+ {
+ "args": ["$compil_prod_dir/test_set_all_outerm",
+ "$input_dir/h_comparison_region.nc"],
+ "title": "Set_all_outerm",
+ "input": "0.001\n"
+ },
{
"args": "$compil_prod_dir/test_get_snapshot",
"title": "Get_snapshot",
diff --git a/depend.mk b/depend.mk
index 36c32cb837fc177109da5cdea738b7cf7a52fba3..175beb861a5cc3a23b0bed6e72ec2c4256f89089 100644
--- a/depend.mk
+++ b/depend.mk
@@ -1,8 +1,8 @@
dispatch_snapshot.o : write_eddy.o send_snapshot.o derived_types.o
-get_snapshot.o : set_all_outerm.o set_max_speed.o receive_snapshot.o local_extrema.o derived_types.o
+get_snapshot.o : set_all_outerm.o set_max_speed.o receive_snapshot.o derived_types.o
receive_snapshot.o : derived_types.o
send_snapshot.o : derived_types.o
-set_all_outerm.o : get_1_outerm.o derived_types.o
+set_all_outerm.o : local_extrema.o get_1_outerm.o derived_types.o
set_max_speed.o : spherical_polygon_area.o mean_speed.o max_speed_contour_ssh.o inside_4.o good_contour.o derived_types.o
get_1_outerm.o : spherical_polygon_area.o outermost_possible_level.o good_contour.o derived_types.o
test_get_snapshot.o : get_snapshot.o dispatch_snapshot.o derived_types.o
@@ -11,6 +11,7 @@ test_inside_4.o : inside_4.o
test_local_extrema.o : local_extrema.o
test_max_speed_contour_ssh.o : max_speed_contour_ssh.o
test_mean_speed.o : mean_speed.o
+test_set_all_outerm.o : set_all_outerm.o dispatch_snapshot.o derived_types.o
test_set_max_speed.o : set_max_speed.o derived_types.o
test_get_1_outerm.o : get_1_outerm.o derived_types.o
write_eddy.o : derived_types.o
diff --git a/file_list b/file_list
index d7ad655a6b6a819bce443d1fde0aa3ca76028deb..73904baf99023dcbf10e1e62da5074354dbc8b65 100644
--- a/file_list
+++ b/file_list
@@ -19,6 +19,7 @@ test_inside_4.f
test_local_extrema.f
test_max_speed_contour_ssh.f
test_mean_speed.f
+test_set_all_outerm.f
test_set_max_speed.f
test_get_1_outerm.f
write_eddy.f