diff --git a/Analysis/plot_snapshot.py b/Analysis/plot_snapshot.py
index 6dcf8ca3795f74970340f8116206873384f0d49b..b3fe08d3eecf8c6a04a034666398b5789e36b55d 100755
--- a/Analysis/plot_snapshot.py
+++ b/Analysis/plot_snapshot.py
@@ -53,8 +53,8 @@ else:
m_s_iterShapes = reader_m_s.iterShapes()
with netCDF4.Dataset("h.nc") as f:
- longitude = f.variables["lon"][:]
- latitude = f.variables["lat"][:]
+ longitude = f.variables["longitude"][:]
+ latitude = f.variables["latitude"][:]
if args.window:
l = input("llcrnrlon, llcrnrlat, urcrnrlon, urcrnrlat = ? ")\
diff --git a/Documentation_texfol/Graphiques/call_graph.gv b/Documentation_texfol/Graphiques/call_graph.gv
index 0689b3c6bd6e621e6ba313a3aabb2d1e7b500f71..1729d3fa5a733569893792038fce50cc6d3c0a1d 100644
--- a/Documentation_texfol/Graphiques/call_graph.gv
+++ b/Documentation_texfol/Graphiques/call_graph.gv
@@ -3,7 +3,6 @@ digraph call_graph
main -> {get_snapshot successive_overlap non_successive_overlap};
main -> dispatch_snapshot;
non_successive_overlap -> {interpolate_eddy weight};
- dispatch_snapshot -> write_eddy;
interpolate_eddy -> write_eddy;
successive_overlap -> weight;
}
\ No newline at end of file
diff --git a/Documentation_texfol/Graphiques/input_output.odg b/Documentation_texfol/Graphiques/input_output.odg
index 001bc6e7995dea4372b2484bcac344ed828e94ed..cde843e63cfdced69511234f1567462a14a1731d 100644
Binary files a/Documentation_texfol/Graphiques/input_output.odg and b/Documentation_texfol/Graphiques/input_output.odg differ
diff --git a/Documentation_texfol/documentation.tex b/Documentation_texfol/documentation.tex
index eaf07eaba0736a2deffe42b888daf727aa1c3958..62d9abb243bb7bb93664e66f16f2d7141da601d4 100644
--- a/Documentation_texfol/documentation.tex
+++ b/Documentation_texfol/documentation.tex
@@ -992,8 +992,8 @@ fichiers de sortie suivants :
\end{description}
Les entrées-sorties sont dans : l'algorithme principal directement
-(lecture de longitude et latitude), \verb+get_snapshot+ (lecture de
-SSH, u, v), \verb+successive_overlap+ (écriture dans
+(lecture de longitude et latitude), \verb+get_snapshot+ (lecture des
+shapefiles), \verb+successive_overlap+ (écriture dans
\verb+edgelist_$m.txt+, arcs entre \verb+k_begin+(m) et
\verb+k_end_main_loop+(m)), \verb+non_successive_overlap+ directement
(écriture dans \verb+edgelist_$m.txt+), \verb+write_eddy+ (écriture
diff --git a/GNUmakefile b/GNUmakefile
index 1281706f821c7ad1b11cc4e39035a4e23dd2a1df..f195c3f922803432edc8533c84219f52bae3e6d1 100644
--- a/GNUmakefile
+++ b/GNUmakefile
@@ -13,7 +13,7 @@ src_test_get_1_outerm = good_contour.f test_get_1_outerm.f derived_types.f get_1
src_test_set_max_speed = test_set_max_speed.f derived_types.f set_max_speed.f good_contour.f max_speed_contour_ssh.f mean_speed.f spherical_polyline_area.f inside_4.f
-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_polyline_area.f mean_speed.f inside_4.f set_all_outerm.f derived_types.f init_shapefiles.f nearby_extr.f get_var.f
+src_test_get_snapshot = test_get_snapshot.f write_eddy.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_polyline_area.f mean_speed.f inside_4.f set_all_outerm.f derived_types.f init_shapefiles.f nearby_extr.f get_var.f
src_test_set_all_outerm = test_set_all_outerm.f derived_types.f set_all_outerm.f local_extrema.f get_1_outerm.f good_contour.f spherical_polyline_area.f nearby_extr.f get_var.f
@@ -23,7 +23,7 @@ src_test_spherical_polygon_area = spherical_polygon_area.f test_spherical_polygo
src_test_read_eddy = test_read_eddy.f derived_types.f init_shapefiles.f read_eddy.f write_eddy.f read_field_indices.f
-src_test_read_snapshot = test_read_snapshot.f derived_types.f dispatch_snapshot.f init_shapefiles.f read_snapshot.f write_eddy.f send_snapshot.f read_eddy.f read_field_indices.f
+src_test_read_snapshot = test_read_snapshot.f derived_types.f init_shapefiles.f read_snapshot.f write_eddy.f read_eddy.f read_field_indices.f
src_test_successive_overlap = test_successive_overlap.f derived_types.f successive_overlap.f read_snapshot.f spherical_polygon_area.f spherical_polyline_area.f weight.f read_eddy.f read_field_indices.f
diff --git a/Tests/test_get_snapshot.f b/Tests/test_get_snapshot.f
index a5443bf4f9255b00ec707bcb5ad78d6009685c9b..958917c49aea9c899ca7d84a2145329fa8e758d0 100644
--- a/Tests/test_get_snapshot.f
+++ b/Tests/test_get_snapshot.f
@@ -1,19 +1,23 @@
program test_get_snapshot
use, intrinsic:: ieee_arithmetic, only: IEEE_SUPPORT_DATATYPE, &
- ieee_support_nan
+ ieee_support_nan, ieee_value, IEEE_QUIET_NAN
use, intrinsic:: ISO_FORTRAN_ENV
! Libraries:
+ use contour_531, only: convert_to_ind
use netcdf, only: nf90_nowrite
use netcdf95, only: nf95_open, find_coord, nf95_inquire_dimension, &
nf95_get_var, nf95_close
use nr_util, only: assert, deg_to_rad, twopi, pi
use shapelib, only: shpfileobject, shpclose
- use derived_types, only: snapshot
- use get_snapshot_m, only: get_snapshot
+ use derived_types, only: snapshot, null_ssh_contour, missing_speed
+ use get_var_m, only: get_var
use init_shapefiles_m, only: init_shapefiles
+ use nearby_extr_m, only: nearby_extr
+ use set_all_outerm_m, only: set_all_outerm
+ use set_max_speed_m, only: set_max_speed
use write_eddy_m, only: write_eddy
implicit none
@@ -22,7 +26,7 @@ program test_get_snapshot
TYPE(shpfileobject) hshp_extremum ! shapefile extremum_$m
TYPE(shpfileobject) hshp_outermost ! shapefile outermost_contour_$m
TYPE(shpfileobject) hshp_max_speed ! shapefile max_speed_contour_$m
- integer i
+ integer i, k
real:: min_amp = 0.
! minimum amplitude of ssh, between extremum and outermost contour,
@@ -31,7 +35,7 @@ program test_get_snapshot
real:: min_radius = 20.
! minimum radius of the equal-area disk for an outermost contour, in km
- integer ncid, dimid, varid, nlon, nlat
+ integer ncid, dimid, varid, nlon, nlat, copy
real lon_min, lon_max, lat_min, lat_max ! longitudes and latitudes, in degrees
real step(2) ! longitude and latitude steps, in rad
logical periodic ! grid is periodic in longitude
@@ -40,6 +44,29 @@ program test_get_snapshot
! maximum radius of an eddy in longitude and latitude, in number of
! grid points
+ real, allocatable:: ssh(:, :)
+ ! (1 - max_radius(1):nlon + max_radius(1), nlat) if periodic, else
+ ! (nlon, nlat) sea-surface height, in m
+
+ real, allocatable:: u(:, :), v(:, :)
+ ! (1 - max_radius(1):nlon + max_radius(1), nlat) if periodic, else
+ ! (nlon, nlat) wind, in m s-1
+
+ real corner(2)
+ ! longitude and latitude of the corner of the whole grid, in rad
+
+ ! Window around each extremum:
+
+ 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
+
namelist /main_nml/ min_amp, max_radius, min_radius
!--------------------------------------------------------------
@@ -77,10 +104,71 @@ program test_get_snapshot
print *, "periodic = ", periodic
call assert(2 * max_radius(1) * step(1) < pi, "test_get_snapshot max_radius")
- call get_snapshot(s, min_amp, min_area = pi * (min_radius * 1e3)**2, m = 1, &
- n_proc = 1, k_end = 1, max_delta = 4, nlon = nlon, nlat = nlat, k = 1, &
- max_radius = max_radius, corner = [lon_min, lat_min] * deg_to_rad, &
- step = step, periodic = periodic)
+ copy = merge(max_radius(1), 0, periodic)
+ allocate(ssh(1 - copy:nlon + copy, nlat), u(1 - copy:nlon + copy, nlat), &
+ v(1 - copy:nlon + copy, nlat))
+ k = 1
+ corner = [lon_min, lat_min] * deg_to_rad
+
+ ! Read ssh, u and v at date k:
+
+ call nf95_open("h.nc", nf90_nowrite, ncid)
+ call get_var(periodic, max_radius(1), ssh, ncid, nlon, k, name = "adt", &
+ new_fill_value = huge(0.))
+ ! (We cannot keep the original fill value because Contour_531
+ ! works with an upper limit for valid values.)
+ call nf95_close(ncid)
+
+ call nf95_open("uv.nc", nf90_nowrite, ncid)
+ call get_var(periodic, max_radius(1), u, ncid, nlon, k, name = "u", &
+ new_fill_value = ieee_value(0., IEEE_QUIET_NAN))
+ call get_var(periodic, max_radius(1), v, ncid, nlon, k, name = "v", &
+ new_fill_value = ieee_value(0., IEEE_QUIET_NAN))
+ ! (We will need quiet NaNs rather the original fill values for u
+ ! and v when we compute the max-speed contours.)
+ call nf95_close(ncid)
+
+ call set_all_outerm(s, min_amp, max_radius, step, periodic, ssh, &
+ corner, min_area = pi * (min_radius * 1e3)**2)
+
+ ! Done with outermost contours, now let us take care of the
+ ! max-speed contours.
+
+ do i = 1, s%number_vis_extr
+ if (s%list_vis(i)%valid) then
+ ! Restrict the field to the outermost contour:
+
+ llc = floor(convert_to_ind(minval(s%list_vis(i)%out_cont%points, &
+ dim = 2), corner, step))
+
+ urc = ceiling(convert_to_ind(maxval( &
+ s%list_vis(i)%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)
+
+ corner_window = corner + (llc - 1) * step
+
+ outside_points = nearby_extr(s%extr_map(llc(1):urc(1), &
+ llc(2):urc(2)), s%list_vis, i)
+
+ if (periodic) outside_points(1, :) = outside_points(1, :) &
+ + ceiling((corner_window(1) - outside_points(1, :)) / twopi) &
+ * twopi
+ ! (Shift the longitude of each point to a value close to
+ ! the longitude of the target extremum.)
+
+ call set_max_speed(s%list_vis(i), s%ind_extr(:, i) - llc + 1, &
+ outside_points, 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)
+ else
+ s%list_vis(i)%speed_cont = null_ssh_contour()
+ s%list_vis(i)%max_speed = missing_speed
+ s%list_vis(i)%radius4 = 0
+ end if
+ end do
call init_shapefiles(hshp_extremum, hshp_outermost, hshp_max_speed)
diff --git a/depend.mk b/depend.mk
index 67e5337c7080c3bd3709772f0247c8691e36b78e..2d4a9fe194725831133c7062eaaf13f5f63510ec 100644
--- a/depend.mk
+++ b/depend.mk
@@ -1,17 +1,13 @@
-dispatch_snapshot.o : send_snapshot.o derived_types.o
get_1_outerm.o : spherical_polyline_area.o good_contour.o derived_types.o
-get_snapshot.o : set_all_outerm.o set_max_speed.o receive_snapshot.o nearby_extr.o get_var.o derived_types.o
nearby_extr.o : derived_types.o
read_eddy.o : read_field_indices.o derived_types.o
read_snapshot.o : read_field_indices.o read_eddy.o derived_types.o
-receive_snapshot.o : derived_types.o
-send_snapshot.o : derived_types.o
set_all_outerm.o : local_extrema.o nearby_extr.o get_1_outerm.o derived_types.o
set_max_speed.o : spherical_polyline_area.o mean_speed.o max_speed_contour_ssh.o inside_4.o good_contour.o derived_types.o
spherical_polygon_area.o : spherical_polyline_area.o
successive_overlap.o : weight.o spherical_polyline_area.o spherical_polygon_area.o derived_types.o
test_get_1_outerm.o : get_1_outerm.o derived_types.o
-test_get_snapshot.o : write_eddy.o init_shapefiles.o get_snapshot.o derived_types.o
+test_get_snapshot.o : write_eddy.o set_max_speed.o set_all_outerm.o nearby_extr.o init_shapefiles.o get_var.o derived_types.o
test_local_extrema.o : local_extrema.o
test_nearby_extr.o : nearby_extr.o derived_types.o
test_read_eddy.o : write_eddy.o read_field_indices.o read_eddy.o init_shapefiles.o derived_types.o
diff --git a/get_snapshot.f b/get_snapshot.f
index 48515b1a5e19bae2b4e6ec1f44771c67c4ce10b6..d8b6fc91a694551926f613459e9656fe6a5f1f19 100644
--- a/get_snapshot.f
+++ b/get_snapshot.f
@@ -4,33 +4,20 @@ module get_snapshot_m
contains
- subroutine get_snapshot(s, min_amp, min_area, m, n_proc, k_end, max_delta, &
- nlon, nlat, k, max_radius, corner, step, periodic)
+ subroutine get_snapshot(s, m, n_proc, k_end, max_delta, nlon, nlat, k, corner)
use, intrinsic:: ieee_arithmetic, only: ieee_value, IEEE_QUIET_NAN
! Libraries:
- use contour_531, only: convert_to_ind, null_polyline
- use netcdf, only: nf90_nowrite
- use netcdf95, only: nf95_open, nf95_close
- use nr_util, only: twopi
-
- use derived_types, only: snapshot, null_ssh_contour, missing_speed
- use get_var_m, only: get_var
- use nearby_extr_m, only: nearby_extr
+ use shapelib, only: shpfileobject, shpclose
+ use shapelib_03, only: shp_open_03
+
+ use derived_types, only: snapshot
+ use read_snapshot_m, only: read_snapshot
use receive_snapshot_m, only: receive_snapshot
- use set_max_speed_m, only: set_max_speed
- use set_all_outerm_m, only: set_all_outerm
type(snapshot), intent(out):: s
- real, intent(in):: min_amp
- ! minimum amplitude of ssh, between extremum and outermost
- ! contour, in m
-
- real, intent(in):: min_area
- ! minimum area of an outermost contour, in m2
-
integer, intent(in):: m ! rank of MPI process
integer, intent(in):: n_proc ! number of MPI processes
integer, intent(in):: k_end ! last date index analyzed by this MPI process
@@ -42,98 +29,28 @@ contains
integer, intent(in):: nlon, nlat
integer, intent(in):: k ! date index
- integer, intent(in):: max_radius(:) ! (2) maximum radius of an
- ! eddy in longitude and latitude, in number of grid points
-
real, intent(in):: corner(:) ! (2) longitude and latitude of the
! corner of the whole grid, in rad
-
- real, intent(in):: step(:) ! (2) longitude and latitude steps, in rad
- logical, intent(in):: periodic ! grid is periodic in longitude
! Local:
-
- integer ncid
- real ssh(1 - merge(max_radius(1), 0, periodic):nlon + merge(max_radius(1), &
- 0, periodic), nlat) ! sea-surface height, in m
- real, dimension(1 - merge(max_radius(1), 0, periodic):nlon &
- + merge(max_radius(1), 0, periodic), nlat):: u, v ! wind, in m s-1
- integer i
-
- ! Window around each extremum:
-
- 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
+ TYPE(shpfileobject) hshp_extremum ! shapefile extremum_1
+ TYPE(shpfileobject) hshp_outermost ! shapefile outermost_contour_1
+ TYPE(shpfileobject) hshp_max_speed ! shapefile max_speed_contour_1
!--------------------------------------------------------------
if (m == n_proc .or. k <= k_end - max_delta) then
- ! Read ssh, u and v at date k:
-
- call nf95_open("h.nc", nf90_nowrite, ncid)
- call get_var(periodic, max_radius(1), ssh, ncid, nlon, k, name = "adt", &
- new_fill_value = huge(0.))
- ! (We cannot keep the original fill value because Contour_531
- ! works with an upper limit for valid values.)
- call nf95_close(ncid)
-
- call nf95_open("uv.nc", nf90_nowrite, ncid)
- call get_var(periodic, max_radius(1), u, ncid, nlon, k, name = "u", &
- new_fill_value = ieee_value(0., IEEE_QUIET_NAN))
- call get_var(periodic, max_radius(1), v, ncid, nlon, k, name = "v", &
- new_fill_value = ieee_value(0., IEEE_QUIET_NAN))
- ! (We will need quiet NaNs rather the original fill values for u
- ! and v when we compute the max-speed contours.)
- call nf95_close(ncid)
-
- call set_all_outerm(s, min_amp, max_radius, step, periodic, ssh, &
- corner, min_area)
-
- ! Done with outermost contours, now let us take care of the
- ! max-speed contours.
-
- do i = 1, s%number_vis_extr
- if (s%list_vis(i)%valid) then
- ! Restrict the field to the outermost contour:
-
- llc = floor(convert_to_ind(minval(s%list_vis(i)%out_cont%points, &
- dim = 2), corner, step))
-
- urc = ceiling(convert_to_ind(maxval( &
- s%list_vis(i)%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)
-
- corner_window = corner + (llc - 1) * step
-
- outside_points = nearby_extr(s%extr_map(llc(1):urc(1), &
- llc(2):urc(2)), s%list_vis, i)
-
- if (periodic) outside_points(1, :) = outside_points(1, :) &
- + ceiling((corner_window(1) - outside_points(1, :)) / twopi) &
- * twopi
- ! (Shift the longitude of each point to a value close to
- ! the longitude of the target extremum.)
-
- call set_max_speed(s%list_vis(i), s%ind_extr(:, i) - llc + 1, &
- outside_points, 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)
- else
- s%list_vis(i)%speed_cont = null_ssh_contour()
- s%list_vis(i)%max_speed = missing_speed
- s%list_vis(i)%radius4 = 0
- end if
- end do
+ ! Read at date k:
+ call shp_open_03("extremum", pszaccess = "rb", hshp = hshp_extremum)
+ call shp_open_03("outermost_contour", pszaccess = "rb", &
+ hshp = hshp_outermost)
+ call shp_open_03("max_speed_contour", pszaccess = "rb", &
+ hshp = hshp_max_speed)
+ call read_snapshot(s, hshp_extremum, hshp_outermost, hshp_max_speed, &
+ corner, nlon, nlat)
+ CALL shpclose(hshp_extremum)
+ CALL shpclose(hshp_outermost)
+ CALL shpclose(hshp_max_speed)
s%number_eddies = s%number_vis_extr
else