From ad047254f7842119ba6973114795ace0b44effe5 Mon Sep 17 00:00:00 2001
From: Lionel GUEZ <guez@lmd.ens.fr>
Date: Wed, 14 Aug 2019 16:00:11 +0200
Subject: [PATCH] Circumvent gfortran bug in `CMakeLists.txt`.
Prepare for overlap of eddies at non-successive dates:
- add argument `delta_out` to procedure `candidate_overlap` (not yet
used);
- add procedure `interpolate_eddy` (empty for now);
- rename `successive_overlap` to overlap because the procedure will handle both
cases, successive dates and non-successive dates;
- generalize procedure overlap to handle delta between dates.
---
CMakeLists.txt | 8 +++
Documentation_texfol/documentation.tex | 29 +++++------
Tests/CMakeLists.txt | 4 +-
Tests/test_successive_overlap.f90 | 6 +--
candidate_overlap.f90 | 3 +-
interpolate_eddy.f90 | 11 +++++
successive_overlap.f90 => overlap.f90 | 67 +++++++++++++++++---------
read_snapshot.f90 | 2 +-
8 files changed, 83 insertions(+), 47 deletions(-)
create mode 100644 interpolate_eddy.f90
rename successive_overlap.f90 => overlap.f90 (56%)
diff --git a/CMakeLists.txt b/CMakeLists.txt
index f60a728b..e02593a0 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 7b44ad30..af8c26f0 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 8b360eae..d51232cd 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 68aeb4f7..a5e4cc6e 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 950601a2..b676a862 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 00000000..c224af41
--- /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 e755c0ab..278a182d 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 4438f1bc..2b24996b 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
--
GitLab