diff --git a/CMakeLists.txt b/CMakeLists.txt
index f60a728bec5f475b90a42dbdb028ac20096f91ef..e02593a0c84328633656d19de506b4ea3df9e1b0 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -56,4 +56,12 @@ target_link_libraries(extraction_eddies ${contour_531_LIBRARY}
include(Tests/CMakeLists.txt)
include(ConfigureCompilerFlags)
+
+if(CMAKE_Fortran_COMPILER_ID MATCHES GNU)
+ # Gfortran bug with quiet NaNs
+ # https://gcc.gnu.org/bugzilla/show_bug.cgi?id=70817
+ string(REPLACE "-ffpe-trap=invalid" "" CMAKE_Fortran_FLAGS_DEBUG
+ ${CMAKE_Fortran_FLAGS_DEBUG})
+endif()
+
include(TAGS)
diff --git a/Documentation_texfol/documentation.tex b/Documentation_texfol/documentation.tex
index 7b44ad308e135f2061ba8e04b32e62439695e16c..af8c26f08195b4be8b9f33718868c4fe842d6361 100644
--- a/Documentation_texfol/documentation.tex
+++ b/Documentation_texfol/documentation.tex
@@ -1433,8 +1433,8 @@ un seul processus.
Tous les arguments sont des données. e1, e2 : scalaires de type
eddy. j1, j2, indices de position dans la fenêtre correspondant à e1,
e2. j : indice de position dans la fenêtre auquel on veut
-interpoler. k, i : indices de date et de tourbillon. Cf. algorithme
-(\ref{alg:interpolate_eddy}).
+interpoler. k, i : indices de date et de tourbillon pour le tourbillon
+interpolé. Cf. algorithme (\ref{alg:interpolate_eddy}).
\begin{algorithm}[htbp]
\begin{algorithmic}
\STATE e = interpolation entre e1 et e2
@@ -1497,30 +1497,27 @@ flow(j)\%list\_vis\%delta\_in, flow\%number\_eddies. Cf. algorithme
\STATE weight\_delta = weight(flow(j - delta)\%list\_vis(i1),
flow(j)\%list\_vis(i2))
- \STATE flow(j - delta + 1)\%number\_eddies += 1
+ \STATE i\_pred = i1
- \STATE appel de interpolate\_eddy(flow(j -
- delta)\%list\_vis(i1), flow(j)\%list\_vis(i2), j - delta, j, j - delta + 1,
- k - delta + 1, flow(j - delta + 1)\%number\_eddies)
-
- \STATE écrire dans graph.txt : k - delta, i1, k - delta + 1,
- flow(j - delta + 1)\%number\_eddies, weight\_delta
-
- \FOR{j\_interp = j - delta + 2 \TO j - 1}
+ \FOR{j\_interp = j - delta + 1 \TO j - 1}
+ \STATE \COMMENT{compte d'un successeur fictif en j\_interp}
+
\STATE flow(j\_interp)\%number\_eddies += 1
\STATE appel de interpolate\_eddy(flow(j -
delta)\%list\_vis(i1), flow(j)\%list\_vis(i2), j - delta, j, j\_interp, k -
j + j\_interp, flow(j\_interp)\%number\_eddies)
- \STATE écrire dans graph.txt : k - j + j\_interp - 1,
- flow(j\_interp - 1)\%number\_eddies, k - j + j\_interp,
- flow(j\_interp)\%number\_eddies, weight\_delta
+ \STATE écrire dans graph.txt : k - j + j\_interp - 1, i\_pred, k -
+ j + j\_interp, flow(j\_interp)\%number\_eddies, weight\_delta
+
+ \STATE i\_pred = flow(j\_interp)\%number\_eddies
+
\ENDFOR
- \STATE écrire dans graph.txt : k - 1, flow(j - 1)\%number\_eddies, k,
- i2, weight\_delta
+ \STATE écrire dans graph.txt : k - 1, i\_pred, k, i2,
+ weight\_delta
\STATE flow(j - delta)\%list\_vis(i1)\%delta\_out = delta
diff --git a/Tests/CMakeLists.txt b/Tests/CMakeLists.txt
index 8b360eae526f280575f97d2c795c188e046ae92e..d51232cd644feeeed8e173cc6067a08a26d95423 100644
--- a/Tests/CMakeLists.txt
+++ b/Tests/CMakeLists.txt
@@ -30,9 +30,9 @@ target_link_libraries(test_nearby_extr ${contour_531_LIBRARY}
# test_successive_overlap
add_executable(test_successive_overlap derived_types.f90
- successive_overlap.f90 read_snapshot.f90 spher_polygon_area.f90
+ overlap.f90 read_snapshot.f90 spher_polygon_area.f90
spher_polyline_area.f90 weight.f90 read_eddy.f90
- read_field_indices.f90 candidate_overlap.f90
+ read_field_indices.f90 candidate_overlap.f90 interpolate_eddy.f90
${CMAKE_CURRENT_LIST_DIR}/test_successive_overlap.f90)
target_include_directories(test_successive_overlap PRIVATE
diff --git a/Tests/test_successive_overlap.f90 b/Tests/test_successive_overlap.f90
index 68aeb4f728933525aa4aeaeba4429e61889731cb..a5e4cc6e661da6e91d1af1bc7b6fff2f80891d78 100644
--- a/Tests/test_successive_overlap.f90
+++ b/Tests/test_successive_overlap.f90
@@ -11,7 +11,7 @@ program test_successive_overlap
use derived_types, only: snapshot
use read_field_indices_m, only: read_field_indices
use read_snapshot_m, only: read_snapshot
- use successive_overlap_m, only: successive_overlap
+ use overlap_m, only: overlap
implicit none
@@ -98,8 +98,8 @@ program test_successive_overlap
// '"successor eddy subscript"'
write(unit_edgelist, fmt = *) "k1 i1 k2 i2 weight"
- call successive_overlap(flow, unit_edgelist, nlon, nlat, periodic, dist_lim, &
- j = 2, k = k)
+ call overlap(flow, unit_edgelist, nlon, nlat, periodic, dist_lim, j = 2, &
+ k = k, delta = 1)
close(unit_edgelist)
print *, 'Created file "edgelist.csv".'
diff --git a/candidate_overlap.f90 b/candidate_overlap.f90
index 950601a282bc6c62b93ceca21a27266b37b306a5..b676a862395ff9751ab12743e2837ee018bac30a 100644
--- a/candidate_overlap.f90
+++ b/candidate_overlap.f90
@@ -4,7 +4,7 @@ module candidate_overlap_m
contains
- function candidate_overlap(extr_map, list_vis, cyclone)
+ function candidate_overlap(extr_map, list_vis, cyclone, delta_out)
! Find the eddies in extr_map that are valid and have a given
! cyclonicity.
@@ -22,6 +22,7 @@ contains
! cyclone to be defined.
logical, intent(in):: cyclone ! cyclonicity of the target extremum
+ integer, intent(in):: delta_out
!---------------------------------------------------------------------
diff --git a/interpolate_eddy.f90 b/interpolate_eddy.f90
new file mode 100644
index 0000000000000000000000000000000000000000..c224af41ba5b19cb9c706350cb858b854f956467
--- /dev/null
+++ b/interpolate_eddy.f90
@@ -0,0 +1,11 @@
+module interpolate_eddy_m
+
+ implicit none
+
+contains
+
+ subroutine interpolate_eddy
+
+ end subroutine interpolate_eddy
+
+end module interpolate_eddy_m
diff --git a/successive_overlap.f90 b/overlap.f90
similarity index 56%
rename from successive_overlap.f90
rename to overlap.f90
index e755c0ab18abc61cd39712ba75f1e8097a53f503..278a182d9548c845cfe7e8878ea26093a3138d84 100644
--- a/successive_overlap.f90
+++ b/overlap.f90
@@ -1,15 +1,15 @@
-module successive_overlap_m
+module overlap_m
implicit none
contains
- subroutine successive_overlap(flow, unit_edgelist, nlon, nlat, periodic, &
- dist_lim, j, k)
+ subroutine overlap(flow, unit_edgelist, nlon, nlat, periodic, &
+ dist_lim, j, k, delta)
- ! Finds edges between flow(j - 1) and flow(j), corresponding to
- ! dates k - 1 and k. Writes these edges to unit_edgelist. Updates
- ! flow(j - 1)%list_vis%delta_out and flow(j)%list_vis%delta_in.
+ ! Finds edges between flow(j - delta) and flow(j), corresponding to
+ ! dates k - delta and k. Writes these edges to unit_edgelist. Updates
+ ! flow(j - delta)%list_vis%delta_out and flow(j)%list_vis%delta_in.
! Libraries:
use contour_531, only: polyline
@@ -18,6 +18,7 @@ contains
use candidate_overlap_m, only: candidate_overlap
use derived_types, only: snapshot
+ use interpolate_eddy_m, only: interpolate_eddy
use spher_polygon_area_m, only: spher_polygon_area
use spher_polyline_area_m, only: spher_polyline_area
use weight_m, only: weight
@@ -35,12 +36,14 @@ contains
! position in time window, between 2 and max_delta + 1
integer, intent(in):: k ! date index
+ integer, intent(in):: delta
! Local:
- integer i1, i2, l, n_select, m
+ integer i1, i2, l, n_select, m, i_pred, j_interp
type(polyline) polyline_1, polyline_2
type(polygon) res_pol
+ real w
integer, allocatable:: selection(:)
! identifying numbers of a selection of extrema
@@ -51,12 +54,12 @@ contains
!-----------------------------------------------------------------------
- do i1 = 1, flow(j - 1)%number_vis_extr
- if (flow(j - 1)%list_vis(i1)%valid) then
+ do i1 = 1, flow(j - delta)%number_vis_extr
+ if (flow(j - delta)%list_vis(i1)%valid) then
! Define the geographical window around each eddy extremum:
- llc = flow(j - 1)%ind_extr(:, i1) - dist_lim
- urc = flow(j - 1)%ind_extr(:, i1) + dist_lim
+ llc = flow(j - delta)%ind_extr(:, i1) - dist_lim
+ urc = flow(j - delta)%ind_extr(:, i1) + dist_lim
llc(2) = max(llc(2), 1)
urc(2) = min(urc(2), nlat)
@@ -68,12 +71,13 @@ contains
selection = candidate_overlap(flow(j)%extr_map(llc(1):urc(1), &
llc(2):urc(2)), flow(j)%list_vis, &
- flow(j - 1)%list_vis(i1)%cyclone)
+ flow(j - delta)%list_vis(i1)%cyclone, &
+ flow(j - delta)%list_vis(i1)%delta_out)
n_select = size(selection)
if (n_select /= 0) polyline_1 &
- = merge(flow(j - 1)%list_vis(i1)%speed_cont%polyline, &
- flow(j - 1)%list_vis(i1)%out_cont%polyline, &
- flow(j - 1)%list_vis(i1)%speed_cont%n_points /= 0)
+ = merge(flow(j - delta)%list_vis(i1)%speed_cont%polyline, &
+ flow(j - delta)%list_vis(i1)%out_cont%polyline, &
+ flow(j - delta)%list_vis(i1)%speed_cont%n_points /= 0)
DO l = 1, n_select
i2 = selection(l)
@@ -83,12 +87,12 @@ contains
! Shift the longitudes of polyline_2 to values close to the
! longitude of extremum i1:
- m = floor((flow(j - 1)%list_vis(i1)%coord_extr(1) &
+ m = floor((flow(j - delta)%list_vis(i1)%coord_extr(1) &
- flow(j)%list_vis(i2)%coord_extr(1)) / twopi + 0.5)
polyline_2%points(1, :) = polyline_2%points(1, :) + m * twopi
- call gpc_polygon_clip_f(GPC_INT, polygon(nparts = 1, &
- part = [polyline_1], hole = [.false.]), &
+ call gpc_polygon_clip_f(GPC_INT, &
+ polygon(nparts = 1, part = [polyline_1], hole = [.false.]), &
polygon(nparts = 1, part = [polyline_2], hole = [.false.]), &
res_pol)
@@ -97,16 +101,31 @@ contains
if (spher_polygon_area(res_pol) >= 0.25 &
* min(spher_polyline_area(polyline_1), &
spher_polyline_area(polyline_2))) then
- write(unit_edgelist, fmt = *) k - 1, i1, k, i2, &
- weight(flow(j - 1)%list_vis(i1), flow(j)%list_vis(i2))
- flow(j - 1)%list_vis(i1)%delta_out = 1
- flow(j)%list_vis(i2)%delta_in = 1
+ w = weight(flow(j - delta)%list_vis(i1), &
+ flow(j)%list_vis(i2))
+ i_pred = i1
+
+ do j_interp = j - delta + 1, j - 1
+ ! Count one more fictitious eddy at j_interp!
+ flow(j_interp)%number_eddies &
+ = flow(j_interp)%number_eddies + 1
+
+ call interpolate_eddy
+ write(unit_edgelist, fmt = *) k - j + j_interp - 1, &
+ i_pred, k - j + j_interp, &
+ flow(j_interp)%number_eddies, w
+ i_pred = flow(j_interp)%number_eddies
+ end do
+
+ write(unit_edgelist, fmt = *) k - 1, i_pred, k, i2, w
+ flow(j - delta)%list_vis(i1)%delta_out = delta
+ flow(j)%list_vis(i2)%delta_in = delta
end if
end if
end DO
end if
end do
- end subroutine successive_overlap
+ end subroutine overlap
-end module successive_overlap_m
+end module overlap_m
diff --git a/read_snapshot.f90 b/read_snapshot.f90
index 4438f1bc706b24c0b8b4c4579e3e3410eed434d9..2b24996b6ff570dc8b172bea7782eac337067da6 100644
--- a/read_snapshot.f90
+++ b/read_snapshot.f90
@@ -11,7 +11,7 @@ contains
! Libraries:
use contour_531, only: convert_to_ind
- use nr_util, only: assert, deg_to_rad
+ use nr_util, only: assert
use shapelib, only: shpfileobject
use shapelib_03, only: shp_get_info_03