Open several SHPC in program `eddy_graph`

But reading from them is not here yet.

Common/read_snapshot.f90

@@ -14,7 +14,7 @@ contains
     use read_eddy_m, only: read_eddy
 
     type(snapshot), intent(out):: s ! completely defined
-    TYPE(shpc), intent(in):: hshp
+    TYPE(shpc), intent(in):: hshp(:)
 
     integer, intent(in):: nlon, nlat
     ! size of ssh array in input NetCDF, assuming no repeated point if
@@ -37,18 +37,18 @@ contains
 
     !---------------------------------------------------------------------
 
-    if (k == hshp%d0) then
+    if (k == hshp(1)%d0) then
        ishape_first = 0
     else
-       ! {k > hshp%d0}
-       ishape_first = hshp%ishape_last(k - 1) + 1
+       ! {k > hshp(1)%d0}
+       ishape_first = hshp(1)%ishape_last(k - 1) + 1
     end if
 
-    s%number_extr = hshp%ishape_last(k) - ishape_first + 1
+    s%number_extr = hshp(1)%ishape_last(k) - ishape_first + 1
     allocate(s%list(s%number_extr))
 
-    do ishape = ishape_first, hshp%ishape_last(k)
-       call read_eddy(e, date_read, eddy_i, hshp, ishape)
+    do ishape = ishape_first, hshp(1)%ishape_last(k)
+       call read_eddy(e, date_read, eddy_i, hshp(1), ishape)
 
        ! Check that all the eddies have the same date index:
        if (date_read /= k) then

Inst_eddies/Tests/test_nearby_extr.f90

@@ -44,7 +44,7 @@ program test_nearby_extr
   close(unit)
 
   call shpc_open(hshp, trim(shpc_dir), pszaccess = "rb")
-  call read_snapshot(s, hshp, nlon, nlat, k = hshp%d0, &
+  call read_snapshot(s, [hshp], nlon, nlat, k = hshp%d0, &
        corner = corner_deg * deg_to_rad, step = step_deg * deg_to_rad, &
        copy = 0)
   CALL shpc_close(hshp)

Overlap/Tests/test_get_dispatch_snap.f90

@@ -97,7 +97,7 @@ program test_get_dispatch_snap
   end if
   
   call get_snapshot(s, nlon, nlat, corner * deg_to_rad, step * deg_to_rad, &
-       copy, hshp, k, k_end, rank, n_proc, max_delta = 1)
+       copy, [hshp], k, k_end, rank, n_proc, max_delta = 1)
   CALL shpc_close(hshp)
   call dispatch_snapshot(s, unit_isolated, rank, k_begin, max_delta = 1, k = k)
 

Overlap/Tests/test_overlap.f90

@@ -74,9 +74,9 @@ program test_overlap
        action = "write")
   write(unit, fmt = *) '{"e_overestim": ', e_overestim, '}'
   close(unit)
-  call read_snapshot(flow(1), hshp, nlon, nlat, k_test_1, corner, step, copy)
-  call read_snapshot(flow(max_delta + 1), hshp, nlon, nlat, k_test_2, corner, &
-       step, copy)
+  call read_snapshot(flow(1), [hshp], nlon, nlat, k_test_1, corner, step, copy)
+  call read_snapshot(flow(max_delta + 1), [hshp], nlon, nlat, k_test_2, &
+       corner, step, copy)
   CALL shpc_close(hshp)
   print *, "Enter flow(1)%list%delta_out (array with ", &
        flow(1)%number_extr, " values):"

Overlap/Tests/test_read_snapshot.f90

@@ -62,7 +62,7 @@ program test_read_snapshot
   copy = merge(dist_lim, 0, periodic)
 
   call shpc_open(hshp, trim(shpc_dir), pszaccess = "rb")
-  call read_snapshot(s, hshp, nlon, nlat, k = hshp%d0, &
+  call read_snapshot(s, [hshp], nlon, nlat, k = hshp%d0, &
        corner = corner_deg * deg_to_rad, step = step_deg * deg_to_rad, &
        copy = copy)
   CALL shpc_close(hshp)

Overlap/Tests/test_send_recv.f90

@@ -95,7 +95,7 @@ program test_send_recv
   call ezmpi_bcast(copy, root = 0)
   call ezmpi_bcast(n_dates, root = 0)
 
-  if (rank == 1) call read_snapshot(s, hshp, nlon, nlat, k = hshp%d0, &
+  if (rank == 1) call read_snapshot(s, [hshp], nlon, nlat, k = hshp%d0, &
        corner = corner_deg * deg_to_rad, step = step_deg * deg_to_rad, &
        copy = copy)
   CALL shpc_close(hshp)

Overlap/eddy_graph.f90

@@ -19,7 +19,7 @@ program eddy_graph
   
   implicit none
 
-  integer rank, n_proc, copy, delta, j
+  integer rank, n_proc, copy, delta, i, j
   integer:: n_dates = huge(0) ! number of dates to read
   integer k ! date index
   integer k_begin, k_end, k_end_main_loop
@@ -49,8 +49,9 @@ program eddy_graph
 
   logical periodic ! grid is periodic in longitude
   integer unit_isolated, unit
+  integer n_shpc ! number of input SHPC directories
   integer e_overestim ! over-estimation of the number of eddies at each date
-  TYPE(shpc) hshp
+  TYPE(shpc), allocatable:: hshp(:) ! (n_shpc)
   type(snapshot), allocatable:: flow(:)
   character(len = 30) file
 
@@ -64,20 +65,36 @@ program eddy_graph
   call mpi_init
   call MPI_Comm_rank(MPI_Comm_world, rank)
   call MPI_Comm_size(MPI_Comm_world, n_proc)
-  call get_command_arg_dyn(1, shpc_dir, "Required argument: SHPC-directory")
-  call shpc_open(hshp, trim(shpc_dir), pszaccess = "rb")
 
   if (rank == 0) then
+     n_shpc = COMMAND_ARGUMENT_COUNT()
+     call assert(n_shpc /= 0, &
+          "Required arguments: SHPC-directory [SHPC-directory] ...")
+  end if
+
+  call ezmpi_bcast(n_shpc, root = 0)
+  allocate(hshp(n_shpc))
+
+  do i = 1, n_shpc
+     call get_command_arg_dyn(i, shpc_dir)
+     call shpc_open(hshp(i), trim(shpc_dir), pszaccess = "rb")
+  end do
+
+  if (rank == 0) then
+     ! Assuming grid_nml.txt is the same in all input SHPC:
      call new_unit(unit)
-     open(unit, file = shpc_dir // "/grid_nml.txt", status = "old", &
+     open(unit, file = hshp(1)%dir // "/grid_nml.txt", status = "old", &
           action = "read", position = "rewind")
      read(unit, nml = grid_nml)
      close(unit)
+
      write(unit = error_unit, nml = main_nml)
      write(unit = error_unit, fmt = *) "Enter namelist main_nml."
      read(unit = *, nml = main_nml)
      write(unit = *, nml = main_nml)
      call assert(max_delta >= 1, "eddy_graph max_delta")
+     call assert(n_dates == huge(0) .or. n_shpc == 1, "If we have several " &
+          // "SHPC directories, then we read all the dates from them")
 
      ! As we are requiring the grid spacing to be uniform, the value of
      ! "periodic" may be deduced from the values of step_deg(1) and nlon:
@@ -87,13 +104,14 @@ program eddy_graph
           "eddy_graph dist_lim")
      copy = merge(dist_lim, 0, periodic)
 
-     if (n_dates == huge(0)) n_dates = size(hshp%ishape_last)
+     if (n_dates == huge(0)) &
+          n_dates = sum([(size(hshp(i)%ishape_last), i = 1, n_shpc)])
      call assert(n_dates >= max_delta + 1, &
           "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([hshp%ishape_last(hshp%d0) + 1, &
-          ediff1d(hshp%ishape_last)])
+     e_overestim = maxval([(hshp(i)%ishape_last(hshp(i)%d0) + 1, &
+          ediff1d(hshp(i)%ishape_last), i = 1, n_shpc)])
      corner = corner_deg * deg_to_rad
      step = step_deg * deg_to_rad
   end if
@@ -125,13 +143,13 @@ program eddy_graph
   call open_edge_file(rank)
   ! We do not write the title line. That will be handled by eddy_graph.sh.
 
-  k_begin = hshp%d0 + (rank * n_dates) / n_proc
+  k_begin = hshp(1)%d0 + (rank * n_dates) / n_proc
 
   if (rank < n_proc - 1) then
-     k_end = hshp%d0 + ((rank + 1) * n_dates) / n_proc + max_delta - 1
+     k_end = hshp(1)%d0 + ((rank + 1) * n_dates) / n_proc + max_delta - 1
      k_end_main_loop = k_end - max_delta + 1
   else
-     k_end = hshp%d0 + n_dates - 1
+     k_end = hshp(1)%d0 + n_dates - 1
      k_end_main_loop = k_end
   end if
 
@@ -191,7 +209,10 @@ program eddy_graph
 
   ! 5. Back matter
 
-  CALL shpc_close(hshp)
+  do i = 1, n_shpc
+     CALL shpc_close(hshp(i))
+  end do
+
   close(unit_isolated)
   close(unit_edge)
   call mpi_finalize 

Overlap/get_snapshot.f90

@@ -19,7 +19,7 @@ contains
 
     real, intent(in):: step(:) ! (2) longitude and latitude steps, in rad
     integer, intent(in):: copy
-    TYPE(shpc), intent(in):: hshp
+    TYPE(shpc), intent(in):: hshp(:)
     integer, intent(in):: k ! date index
     integer, intent(in):: k_end ! last date index analyzed by this MPI process
     integer, intent(in):: rank ! of MPI process