From 10f339d21d75ec993948225a00a1c9798f0d587c Mon Sep 17 00:00:00 2001
From: Lionel GUEZ <guez@lmd.ens.fr>
Date: Mon, 2 Sep 2019 16:11:18 +0200
Subject: [PATCH] Move the inner part of overlap to `write_overlap`

Add `.cmake` files. The idea is to have a build system independent of
the environment. This allows us to have a customized version of
`ConfigureCompilerFlags.cmake` without the `-ffpe-trap=invalid`
option.

Move the inner part of procedure overlap to a separate procedure
`write_overlap`. overlap was too complex with 6 nested constructs. The
downside is that the initialization of `e` is repeated for every
non-successive overlap. If necessary, we can later add an
initialization on first entry (but an initialization at declaration is
not possible). `write_overlap` has some consistency: all the arguments
are `intent(in)`, scalars or small vectors of intrinsic type, and
`write_overlap` only writes to files. All modifications to flow are
done outside of `write_overlap`.
---
 CMakeLists.txt               | 11 +----
 ConfigureCompilerFlags.cmake | 47 +++++++++++++++++++
 FindNetCDF.cmake             | 41 +++++++++++++++++
 TAGS.cmake                   | 61 ++++++++++++++++++++++++
 Tests/CMakeLists.txt         |  1 +
 overlap.f90                  | 84 ++++++++++------------------------
 write_overlap.f90            | 89 ++++++++++++++++++++++++++++++++++++
 7 files changed, 265 insertions(+), 69 deletions(-)
 create mode 100644 ConfigureCompilerFlags.cmake
 create mode 100644 FindNetCDF.cmake
 create mode 100644 TAGS.cmake
 create mode 100644 write_overlap.f90

diff --git a/CMakeLists.txt b/CMakeLists.txt
index e02593a0..c1df1125 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,5 +1,6 @@
 cmake_minimum_required(VERSION 3.14)
 project(Detection_eddies Fortran)
+set(CMAKE_MODULE_PATH ${CMAKE_CURRENT_LIST_DIR})
 
 find_path(jumble_INCLUDE_DIR jumble.mod)
 find_library(jumble_LIBRARY jumble)
@@ -55,13 +56,5 @@ target_link_libraries(extraction_eddies ${contour_531_LIBRARY}
   ${shapelib_LIBRARY} ${geometry_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(ConfigureCompilerFlags)  
 include(TAGS)
diff --git a/ConfigureCompilerFlags.cmake b/ConfigureCompilerFlags.cmake
new file mode 100644
index 00000000..139e87d5
--- /dev/null
+++ b/ConfigureCompilerFlags.cmake
@@ -0,0 +1,47 @@
+if(CMAKE_Fortran_COMPILER_ID MATCHES GNU)
+  # Fortran language options:
+  string(APPEND CMAKE_Fortran_FLAGS " -std=f2003")
+
+  # Error and warning options:
+  string(APPEND CMAKE_Fortran_FLAGS
+    " -fmax-errors=1 -pedantic -Wall -Wcharacter-truncation -Wunused-parameter -Wno-conversion -Wimplicit-interface -Wimplicit-procedure -Wno-integer-division")
+
+  # Debugging options:
+  set(CMAKE_Fortran_FLAGS_DEBUG "-fbacktrace -g -ffpe-trap=zero,overflow")
+  # (We cannot use -ffpe-trap=invalid because of Gfortran bug with
+  # quiet NaNs: https://gcc.gnu.org/bugzilla/show_bug.cgi?id=70817.)
+
+  # Code generation options:
+  string(APPEND CMAKE_Fortran_FLAGS_DEBUG
+    " -fcheck=bounds,do,mem,pointer,recursion -finit-real=nan -O0")
+
+  # Optimization options:
+  set(CMAKE_Fortran_FLAGS_RELEASE -O3)
+
+  # Hardware model options:
+  string(APPEND CMAKE_Fortran_FLAGS_RELEASE " -mcmodel=medium")
+elseif(CMAKE_Fortran_COMPILER_ID MATCHES Intel)
+  # Language:
+  string(APPEND CMAKE_Fortran_FLAGS
+    " -noaltparam -stand f03 -standard-semantics -assume nostd_mod_proc_name")
+  string(APPEND CMAKE_Fortran_FLAGS_DEBUG
+    " -check bounds,format,output_conversion,pointers,stack,uninit")
+
+  # Data:
+  string(APPEND CMAKE_Fortran_FLAGS " -auto -mcmodel=medium")
+  string(APPEND CMAKE_Fortran_FLAGS_DEBUG " -init=arrays,minus_huge,snan")
+  
+  # Compiler diagnostics:
+  string(APPEND CMAKE_Fortran_FLAGS
+    " -warn declarations,general,stderrors,truncated_source,uncalled,unused,usage -traceback -diag-error-limit 1")
+  
+  # Optimization:
+  string(APPEND CMAKE_Fortran_FLAGS_DEBUG " -O0")
+  
+  # Floating point:
+  string(APPEND CMAKE_Fortran_FLAGS_DEBUG " -fp-stack-check -fpe-all=0")
+  
+  # Debug:
+  string(APPEND CMAKE_Fortran_FLAGS_DEBUG
+    " -debug full -debug-parameters all -ftrapuv")
+endif()
diff --git a/FindNetCDF.cmake b/FindNetCDF.cmake
new file mode 100644
index 00000000..05b4e84a
--- /dev/null
+++ b/FindNetCDF.cmake
@@ -0,0 +1,41 @@
+find_path(NetCDF_INCLUDE_DIR
+  NAMES netcdf.h
+  DOC "netcdf include directories")
+##mark_as_advanced(NetCDF_INCLUDE_DIR)
+
+find_library(NetCDF_LIBRARY
+  NAMES netcdf
+  DOC "netcdf library")
+##mark_as_advanced(NetCDF_LIBRARY)
+
+if (NetCDF_INCLUDE_DIR)
+  file(STRINGS "${NetCDF_INCLUDE_DIR}/netcdf_meta.h" _netcdf_version_lines
+    REGEX "#define[ \t]+NC_VERSION_(MAJOR|MINOR|PATCH|NOTE)")
+  string(REGEX REPLACE ".*NC_VERSION_MAJOR *\([0-9]*\).*" "\\1" _netcdf_version_major "${_netcdf_version_lines}")
+  string(REGEX REPLACE ".*NC_VERSION_MINOR *\([0-9]*\).*" "\\1" _netcdf_version_minor "${_netcdf_version_lines}")
+  string(REGEX REPLACE ".*NC_VERSION_PATCH *\([0-9]*\).*" "\\1" _netcdf_version_patch "${_netcdf_version_lines}")
+  string(REGEX REPLACE ".*NC_VERSION_NOTE *\"\([^\"]*\)\".*" "\\1" _netcdf_version_note "${_netcdf_version_lines}")
+  set(NetCDF_VERSION "${_netcdf_version_major}.${_netcdf_version_minor}.${_netcdf_version_patch}${_netcdf_version_note}")
+  unset(_netcdf_version_major)
+  unset(_netcdf_version_minor)
+  unset(_netcdf_version_patch)
+  unset(_netcdf_version_note)
+  unset(_netcdf_version_lines)
+endif ()
+
+include(FindPackageHandleStandardArgs)
+find_package_handle_standard_args(NetCDF
+  REQUIRED_VARS NetCDF_LIBRARY NetCDF_INCLUDE_DIR
+  VERSION_VAR NetCDF_VERSION)
+
+if (NetCDF_FOUND)
+  set(NetCDF_INCLUDE_DIRS "${NetCDF_INCLUDE_DIR}")
+  set(NetCDF_LIBRARIES "${NetCDF_LIBRARY}")
+
+  if (NOT TARGET NetCDF::NetCDF)
+    add_library(NetCDF::NetCDF UNKNOWN IMPORTED)
+    set_target_properties(NetCDF::NetCDF PROPERTIES
+      IMPORTED_LOCATION "${NetCDF_LIBRARY}"
+      INTERFACE_INCLUDE_DIRECTORIES "${NetCDF_INCLUDE_DIR}")
+  endif ()
+endif ()
diff --git a/TAGS.cmake b/TAGS.cmake
new file mode 100644
index 00000000..4cdf2953
--- /dev/null
+++ b/TAGS.cmake
@@ -0,0 +1,61 @@
+#
+# (C) Copyright 2018 Tillmann Heidsieck
+#
+# SPDX-License-Identifier: MIT
+#
+
+find_program(CTAGS ctags)
+
+function(tags_add_dir _dir)
+	get_property(_subdirlist DIRECTORY "${_dir}" PROPERTY SUBDIRECTORIES)
+	foreach(_sd ${_subdirlist})
+		tags_add_dir(${_sd})
+	endforeach()
+
+	get_property(_targets_list DIRECTORY "${_dir}" PROPERTY
+	  BUILDSYSTEM_TARGETS)
+	foreach(_tgt ${_targets_list})
+	  get_property(_t_files  TARGET ${_tgt} PROPERTY SOURCES)
+	  list(FILTER _t_files EXCLUDE REGEX TARGET_OBJECTS)
+	  get_property(_t_defs   TARGET ${_tgt} PROPERTY
+	    COMPILE_DEFINITIONS)
+	  
+	  list(APPEND _all_source_files ${_t_files})
+	  string(REGEX REPLACE "${CMAKE_SOURCE_DIR}/" "" _t_files
+	    "${_t_files}")
+	  string(REGEX REPLACE ";" "\n" _t_files "${_t_files}")
+	  
+	  string(REGEX REPLACE ";" "\n" _t_defs "${_t_defs}")
+	  string(REGEX REPLACE "^([^=]+)$" "\\1=1" _t_defs "${_t_defs}")
+	endforeach()
+	set(_all_source_files ${_all_source_files} PARENT_SCOPE)
+	set(_all_definitions ${_all_definitions} PARENT_SCOPE)
+endfunction()
+
+if(CTAGS)
+	set(SOURCES_LIST "${CMAKE_BINARY_DIR}/sources_list")
+	set(_all_source_files "")
+	set(_all_definitions "")
+
+	file(REMOVE "${SOURCES_LIST}")
+
+	tags_add_dir(${CMAKE_SOURCE_DIR})
+
+	list(SORT _all_source_files)
+	list(REMOVE_DUPLICATES _all_source_files)
+	string(REGEX REPLACE "${CMAKE_SOURCE_DIR}/" "" _asf
+	  "${_all_source_files}")
+	string(REGEX REPLACE ";" "\n" _asf "${_asf}")
+	file(APPEND ${SOURCES_LIST} "\n${_asf}")
+
+	string(REGEX REPLACE ";" "\n" _all_definitions "${_all_definitions}")
+	string(REGEX REPLACE "^([^=]+)$" "\\1=1" _all_definitions
+	  "${_all_definitions}")
+
+	add_custom_target(TAGS DEPENDS ${CMAKE_SOURCE_DIR}/TAGS)
+	add_custom_command(OUTPUT ${CMAKE_SOURCE_DIR}/TAGS
+		COMMAND ctags -e --language-force=fortran -L ${SOURCES_LIST}
+		DEPENDS ${_all_source_files}
+		WORKING_DIRECTORY ${CMAKE_SOURCE_DIR}
+		VERBATIM)
+endif()
diff --git a/Tests/CMakeLists.txt b/Tests/CMakeLists.txt
index 6a8a3c1d..386ae43e 100644
--- a/Tests/CMakeLists.txt
+++ b/Tests/CMakeLists.txt
@@ -33,6 +33,7 @@ add_executable(test_successive_overlap derived_types.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 write_eddy.f90 init_shapefiles.f90
+  write_overlap.f90
   ${CMAKE_CURRENT_LIST_DIR}/test_successive_overlap.f90)
 
 target_include_directories(test_successive_overlap PRIVATE
diff --git a/overlap.f90 b/overlap.f90
index bbcae6a4..900d53d6 100644
--- a/overlap.f90
+++ b/overlap.f90
@@ -7,23 +7,25 @@ contains
   subroutine overlap(flow, unit_edgelist, nlon, nlat, periodic, dist_lim, &
        hshp_extremum, hshp_outermost, hshp_max_speed, j, k, delta)
 
-    ! 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.
+    ! Finds edges between flow(j - delta) and flow(j), corresponding
+    ! to dates k - delta and k. Calls write_overlap for these
+    ! edges. Updates flow(j - delta)%list_vis%delta_out,
+    ! flow(j)%list_vis%delta_in and flow(j - delta + 1:j -
+    ! 1)%number_eddies.
 
     ! Libraries:
     use contour_531, only: polyline
     use gpc_f, only: gpc_polygon_clip_f, GPC_INT, polygon
     use nr_util, only: twopi
+    use shapelib, only: shpfileobject
 
     use candidate_overlap_m, only: candidate_overlap
-    use derived_types, only: snapshot, eddy, null_ssh_contour, missing_speed
+    use derived_types, only: snapshot
     use interpolate_eddy_m, only: interpolate_eddy
-    use shapelib, only: shpfileobject
     use spher_polygon_area_m, only: spher_polygon_area
     use spher_polyline_area_m, only: spher_polyline_area
     use weight_m, only: weight
-    use write_eddy_m, only: write_eddy
+    use write_overlap_m, only: write_overlap
 
     type(snapshot), intent(inout):: flow(:) ! (max_delta + 1)
     integer, intent(in):: unit_edgelist ! logical unit for edgelist file
@@ -49,11 +51,9 @@ contains
 
     ! Local:
 
-    integer i1, i2, l, n_select, m, i_pred, j_interp, k_interp
+    integer i1, i2, l, n_select, m
     type(polyline) polyline_1, polyline_2
     type(polygon) res_pol
-    real w, w_interp, coord_extr_2(2)
-    type(eddy) e
 
     integer, allocatable:: selection(:)
     ! identifying numbers of a selection of extrema
@@ -64,16 +64,6 @@ contains
 
     !-----------------------------------------------------------------------
 
-    if (delta >= 2) then
-       ! Set some null components of the interpolated eddy:
-       e%out_cont = null_ssh_contour()
-       e%speed_cont = null_ssh_contour()
-       e%max_speed = missing_speed
-       e%valid = .false.
-       e%interpolated = .true.
-       e%radius4 = 0
-    end if
-
     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:
@@ -121,47 +111,21 @@ contains
                 if (spher_polygon_area(res_pol) >= 0.25 &
                      * min(spher_polyline_area(polyline_1), &
                      spher_polyline_area(polyline_2))) then
-                   w = weight(flow(j - delta)%list_vis(i1), &
-                        flow(j)%list_vis(i2))
-                   i_pred = i1
-                   
-                   if (delta >= 2) then
-                      ! Prepare interpolation:
-                      
-                      e%cyclone = flow(j - delta)%list_vis(i1)%cyclone
-
-                      ! Define coord_extr_2 with the longitude of
-                      ! extremum i2 shifted to a value close to the
-                      ! longitude of extremum i1:
-                      coord_extr_2(1) = flow(j)%list_vis(i2)%coord_extr(1) &
-                           + m * twopi
-                      coord_extr_2(2) = flow(j)%list_vis(i2)%coord_extr(2)
-                   end if
-
-                   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
-
-                      ! Interpolate extremum:
-                      w_interp = (j_interp - j + delta) / real(delta)
-                      e%coord_extr = (1 - w_interp) &
-                           * flow(j - delta)%list_vis(i1)%coord_extr &
-                           + w_interp * coord_extr_2
-                      e%ssh_extr = (1 - w_interp) &
-                           * flow(j - delta)%list_vis(i1)%ssh_extr + w_interp &
-                           * flow(j)%list_vis(i2)%ssh_extr
-
-                      k_interp = k - j + j_interp
-                      call write_eddy(e, hshp_extremum, hshp_outermost, &
-                           hshp_max_speed, k_interp, &
-                           flow(j_interp)%number_eddies)
-                      write(unit_edgelist, fmt = *) k_interp - 1, i_pred, &
-                           k_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
+                   ! Count one more fictitious eddy between j - delta
+                   ! + 1 and j - 1:
+                   flow(j - delta + 1:j - 1)%number_eddies &
+                        = flow(j - delta + 1:j - 1)%number_eddies + 1
+
+                   call write_overlap(k, delta, i1, i2, &
+                        flow(j - delta + 1:j - 1)%number_eddies, &
+                        flow(j - delta)%list_vis(i1)%cyclone, &
+                        flow(j - delta)%list_vis(i1)%coord_extr, &
+                        flow(j - delta)%list_vis(i1)%ssh_extr, &
+                        flow(j)%list_vis(i2)%coord_extr, &
+                        flow(j)%list_vis(i2)%ssh_extr, m, hshp_extremum, &
+                        hshp_outermost, hshp_max_speed, unit_edgelist, &
+                        weight(flow(j - delta)%list_vis(i1), &
+                        flow(j)%list_vis(i2)))
                    flow(j - delta)%list_vis(i1)%delta_out = delta
                    flow(j)%list_vis(i2)%delta_in = delta
                 end if
diff --git a/write_overlap.f90 b/write_overlap.f90
new file mode 100644
index 00000000..69d570dc
--- /dev/null
+++ b/write_overlap.f90
@@ -0,0 +1,89 @@
+module write_overlap_m
+
+  implicit none
+
+contains
+
+  subroutine write_overlap(k, delta, i1, i2, i_interp, cyclone, coord_extr_1, &
+       ssh_extr_1, coord_extr_2, ssh_extr_2, m, hshp_extremum, hshp_outermost, &
+       hshp_max_speed, unit_edgelist, w)
+
+    ! Libraries:
+    use nr_util, only: twopi
+    use shapelib, only: shpfileobject
+
+    use derived_types, only: eddy, null_ssh_contour, missing_speed
+    use write_eddy_m, only: write_eddy
+
+    integer, intent(in):: k ! date index
+    integer, intent(in):: delta
+    integer, intent(in):: i1
+    integer, intent(in):: i2
+    integer, intent(in):: i_interp(:) ! (delta - 1)
+    logical, intent(in):: cyclone
+    real, intent(in):: coord_extr_1(:) ! (2)
+    real, intent(in):: ssh_extr_1
+
+    real, intent(in):: coord_extr_2(:) ! (2)
+    ! (Note: value attribute not allowed for an array.)
+
+    real, intent(in):: ssh_extr_2
+    integer, intent(in):: m
+    TYPE(shpfileobject), intent(in):: hshp_extremum ! shapefile extremum
+
+    TYPE(shpfileobject), intent(in):: hshp_outermost
+    ! shapefile outermost_contour
+
+    TYPE(shpfileobject), intent(in):: hshp_max_speed ! shapefile x_speed_contour
+    integer, intent(in):: unit_edgelist ! logical unit for edgelist file
+    real, intent(in):: w
+
+    ! Local:
+    type(eddy) e
+    integer i_pred, j_interp, k_interp
+    real w_interp
+    real coord_extr_2_local(2)
+
+    !-------------------------------------------------------------------
+
+    if (delta >= 2) then
+       ! Set some null components of the interpolated eddy:
+       e%out_cont = null_ssh_contour()
+       e%speed_cont = null_ssh_contour()
+       e%max_speed = missing_speed
+       e%valid = .false.
+       e%interpolated = .true.
+       e%radius4 = 0
+
+       ! Prepare interpolation:
+
+       e%cyclone = cyclone
+
+       ! Define coord_extr_2_local with the longitude of extremum i2 shifted
+       ! to a value close to the longitude of extremum i1:
+       coord_extr_2_local(1) = coord_extr_2(1) + m * twopi
+       coord_extr_2_local(2) = coord_extr_2(2)
+    end if
+
+    i_pred = i1
+
+    do j_interp = 1, delta - 1
+       ! Interpolate extremum:
+       w_interp = j_interp / real(delta)
+       e%coord_extr = (1 - w_interp) * coord_extr_1 &
+            + w_interp * coord_extr_2_local
+       e%ssh_extr = (1 - w_interp) * ssh_extr_1 + w_interp * ssh_extr_2
+
+       k_interp = k - delta + j_interp
+       call write_eddy(e, hshp_extremum, hshp_outermost, hshp_max_speed, &
+            k_interp, i_interp(j_interp))
+       write(unit_edgelist, fmt = *) k_interp - 1, i_pred, k_interp, &
+            i_interp(j_interp), w
+       i_pred = i_interp(j_interp)
+    end do
+
+    write(unit_edgelist, fmt = *) k - 1, i_pred, k, i2, w
+
+  end subroutine write_overlap
+
+end module write_overlap_m
-- 
GitLab