Write node indices

Write node indices instead of dates and eddy indices in
edgelist. Compute `e_overestim` from `ishape_last`. Remove header
lines of edgelist. See commit 3370733a. Replace argument d of procedure
overlap by arguments `e_overestim` and k.

Overlap/Tests/test_overlap.f90

@@ -3,7 +3,7 @@ program test_overlap
   use, intrinsic:: ISO_FORTRAN_ENV
 
   ! Libraries:
-  use jumble, only: get_command_arg_dyn, read_column, new_unit
+  use jumble, only: get_command_arg_dyn, read_column, new_unit, ediff1d
   use mpi_f08, only: mpi_init, mpi_finalize, MPI_Comm_rank, MPI_Comm_world, &
        MPI_Comm_size, mpi_abort
   use nr_util, only: deg_to_rad, assert
@@ -36,6 +36,7 @@ program test_overlap
 
   logical periodic ! grid is periodic in longitude
   integer:: nlon = - 1, nlat = - 1, max_delta = 1
+  integer e_overestim ! over-estimation of the number of eddies at each date
 
   integer:: dist_lim = 12
   ! We look for an overlapping eddy at dist_lim (in grid points) of
@@ -78,6 +79,7 @@ program test_overlap
   step = step_deg * deg_to_rad
   allocate(flow(max_delta + 1))
   call read_column(ishape_last, file = trim(shpc_dir) // "/ishape_last.txt")
+  e_overestim = maxval([ishape_last(1) + 1, ediff1d(ishape_last)])
   call shpc_open(hshp, trim(shpc_dir), rank = 0)
   call dbf_read_attribute_03(d_init, hshp%extremum, hshp%extr_date, ishape = 0)
   call read_snapshot(flow(1), hshp, nlon, nlat, d_init, k_test_1, corner, &
@@ -92,15 +94,8 @@ program test_overlap
        flow(max_delta + 1)%number_vis_extr, " values)):"
   read *, flow(max_delta + 1)%list_vis%delta_in
   call open_edge_file(rank = 0)
-
-  ! Title lines:
-  write(unit_edge, fmt = "(1x, a)") '"predecessor date subscript" ' &
-       // '"predecessor eddy subscript" "successor date subscript" ' &
-       // '"successor eddy subscript"'
-  write(unit_edge, fmt = *) "k1 i1 k2 i2 weight"
-
-  call overlap(flow, nlon, nlat, periodic, dist_lim, d = d_init + k_test_2, &
-       delta = max_delta, j = max_delta + 1)
+  call overlap(flow, nlon, nlat, periodic, dist_lim, e_overestim, &
+       k = k_test_2, delta = max_delta, j = max_delta + 1)
   close(unit_edge)
   print *, 'Created file "edgelist.csv".'
 

Overlap/eddy_graph.f90

@@ -4,7 +4,7 @@ program eddy_graph
 
   ! Libraries:
   use ezmpi, only: ezmpi_bcast
-  use jumble, only: get_command_arg_dyn, read_column, new_unit
+  use jumble, only: get_command_arg_dyn, read_column, new_unit, ediff1d
   use mpi_f08, only: mpi_init, mpi_finalize, MPI_Comm_rank, MPI_Comm_world, &
        MPI_Comm_size
   use nr_util, only: assert, deg_to_rad
@@ -20,7 +20,7 @@ program eddy_graph
   
   implicit none
 
-  integer rank, n_proc, copy, delta, j, d
+  integer rank, n_proc, copy, delta, j
   integer:: n_dates = huge(0) ! number of dates to read
   integer k ! date index ((0-based)
   integer d_init, k_begin, k_end, k_end_main_loop
@@ -55,6 +55,7 @@ program eddy_graph
   ! shape at a given date index
 
   integer unit_isolated, unit
+  integer e_overestim ! over-estimation of the number of eddies at each date
   TYPE(shpc) hshp
   type(snapshot), allocatable:: flow(:)
   character(len = 30) file
@@ -97,6 +98,7 @@ program eddy_graph
           "eddy_graph: n_dates should be >= max_delta + 1")
      call assert(n_proc <= n_dates / (max_delta + 1), &
           "eddy_graph: n_proc should be <= n_dates / (max_delta + 1)")
+     e_overestim = maxval([ishape_last(1) + 1, ediff1d(ishape_last)])
      corner = corner_deg * deg_to_rad
      step = step_deg * deg_to_rad
   end if
@@ -112,6 +114,7 @@ program eddy_graph
   call ezmpi_bcast(n_dates, root = 0)
   if (rank /= 0) allocate(ishape_last(n_dates))
   call ezmpi_bcast(ishape_last, root = 0)
+  call ezmpi_bcast(e_overestim, root = 0)
 
   call shpc_open(hshp, trim(shpc_dir), rank)
   if (rank == 0) call dbf_read_attribute_03(d_init, hshp%extremum, &
@@ -148,8 +151,8 @@ program eddy_graph
 
   do delta = 1, max_delta
      do k = k_begin + delta, k_begin + max_delta
-        call overlap(flow, nlon, nlat, periodic, dist_lim, d = d_init + k, &
-             delta = delta, j = k - k_begin + 1)
+        call overlap(flow, nlon, nlat, periodic, dist_lim, e_overestim, k, &
+             delta, j = k - k_begin + 1)
      end do
   end do
 
@@ -161,11 +164,10 @@ program eddy_graph
      flow(:max_delta) = flow(2:)
      call get_snapshot(flow(max_delta + 1), nlon, nlat, ishape_last, corner, &
           step, copy, hshp, d_init, k, k_end, rank, n_proc, max_delta)
-     d = d_init + k
      
      do delta = 1, max_delta
-        call overlap(flow, nlon, nlat, periodic, dist_lim, d, delta, &
-             j = max_delta + 1)
+        call overlap(flow, nlon, nlat, periodic, dist_lim, e_overestim, k, &
+             delta, j = max_delta + 1)
      end do
   end do
 
@@ -181,12 +183,9 @@ program eddy_graph
      ! (reception)
 
      ! Stitching:
-
-     d = d_init + k
-
      do delta = k - k_end + max_delta, max_delta
-        call overlap(flow, nlon, nlat, periodic, dist_lim, d, delta, &
-             j = max_delta + 1)
+        call overlap(flow, nlon, nlat, periodic, dist_lim, e_overestim, k, &
+             delta, j = max_delta + 1)
      end do
   end do
 

Overlap/eddy_graph_in.sh

@@ -29,9 +29,7 @@ done
 export LC_NUMERIC=C # for the sort command
 
 cat >edgelist.csv <<EOF
- "predecessor date" "predecessor eddy subscript" "successor date" "successor eddy subscript"
- d1 i1 d2 i2 weight
-$(sort --key=1 --key=2 --numeric-sort edgelist_no_header.csv)
+$(sort --key=1 --numeric-sort edgelist_no_header.csv)
 EOF
 rm edgelist_no_header.csv 
 

Overlap/overlap.f90

@@ -4,10 +4,11 @@ module overlap_m
 
 contains
 
-  subroutine overlap(flow, nlon, nlat, periodic, dist_lim, d, delta, j)
+  subroutine overlap(flow, nlon, nlat, periodic, dist_lim, e_overestim, k, &
+       delta, j)
 
     ! This procedure finds edges between flow(j - delta) and flow(j),
-    ! corresponding to dates d - delta and d. It writes these
+    ! corresponding to date indices k - delta and k. It writes these
     ! edges. It updates flow(j - delta)%list_vis%delta_out and
     ! flow(j)%list_vis%delta_in.
 
@@ -32,7 +33,10 @@ contains
     ! We look for an overlapping eddy at dist_lim (in grid points) of
     ! the extremum of a given eddy.
 
-    integer, intent(in):: d ! date, in days since 1950-1-1
+    integer, intent(in):: e_overestim
+    ! over-estimation of the number of eddies at each date
+    
+    integer, intent(in):: k ! date index (0-based)
     integer, intent(in):: delta ! between 1 and max_delta
 
     integer, intent(in):: j
@@ -113,7 +117,8 @@ contains
                 if (spher_polygon_area(res_pol) >= 0.25 &
                      * min(spher_polyline_area(polyline_1), &
                      spher_polyline_area(polyline_2))) then
-                   write(unit_edge, fmt = *) d - delta, i1, d, i2, &
+                   write(unit_edge, fmt = *) (k - delta) * e_overestim + i1, &
+                        k * e_overestim + i2, &
                         weight(flow(j - delta)%list_vis(i1), &
                         flow(j)%list_vis(i2))
                    flow(j - delta)%list_vis(i1)%delta_out &