From c0f6b78f7c0bc915c3390c6499e5dea7fa6eded4 Mon Sep 17 00:00:00 2001
From: Lionel GUEZ <guez@lmd.ens.fr>
Date: Thu, 23 Dec 2021 19:49:50 +0100
Subject: [PATCH] Open several SHPC in program `eddy_graph`
But reading from them is not here yet.
---
Common/read_snapshot.f90 | 14 ++++----
Inst_eddies/Tests/test_nearby_extr.f90 | 2 +-
Overlap/Tests/test_get_dispatch_snap.f90 | 2 +-
Overlap/Tests/test_overlap.f90 | 6 ++--
Overlap/Tests/test_read_snapshot.f90 | 2 +-
Overlap/Tests/test_send_recv.f90 | 2 +-
Overlap/eddy_graph.f90 | 45 +++++++++++++++++-------
Overlap/get_snapshot.f90 | 2 +-
8 files changed, 48 insertions(+), 27 deletions(-)
diff --git a/Common/read_snapshot.f90 b/Common/read_snapshot.f90
index 06a1a2a1..2bfb033a 100644
--- a/Common/read_snapshot.f90
+++ b/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
diff --git a/Inst_eddies/Tests/test_nearby_extr.f90 b/Inst_eddies/Tests/test_nearby_extr.f90
index d48d5dfb..bcafa9a7 100644
--- a/Inst_eddies/Tests/test_nearby_extr.f90
+++ b/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)
diff --git a/Overlap/Tests/test_get_dispatch_snap.f90 b/Overlap/Tests/test_get_dispatch_snap.f90
index 2d018663..6cba2a6c 100644
--- a/Overlap/Tests/test_get_dispatch_snap.f90
+++ b/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)
diff --git a/Overlap/Tests/test_overlap.f90 b/Overlap/Tests/test_overlap.f90
index a120d407..30dc2624 100644
--- a/Overlap/Tests/test_overlap.f90
+++ b/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):"
diff --git a/Overlap/Tests/test_read_snapshot.f90 b/Overlap/Tests/test_read_snapshot.f90
index f3959b16..aba590bb 100644
--- a/Overlap/Tests/test_read_snapshot.f90
+++ b/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)
diff --git a/Overlap/Tests/test_send_recv.f90 b/Overlap/Tests/test_send_recv.f90
index d7550728..463e5e1a 100644
--- a/Overlap/Tests/test_send_recv.f90
+++ b/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)
diff --git a/Overlap/eddy_graph.f90 b/Overlap/eddy_graph.f90
index cc5b78c9..c03f5cdb 100644
--- a/Overlap/eddy_graph.f90
+++ b/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
diff --git a/Overlap/get_snapshot.f90 b/Overlap/get_snapshot.f90
index ce957639..68c37206 100644
--- a/Overlap/get_snapshot.f90
+++ b/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
--
GitLab