diff --git a/Common/derived_types.f90 b/Common/derived_types.f90 index 0a065491c9a46f8bf2457b44b6d4dcd6b603c352..8867ebef7681255b6c00da43504e7e780ccf6184 100644 --- a/Common/derived_types.f90 +++ b/Common/derived_types.f90 @@ -75,6 +75,12 @@ module derived_types max_speed_eddy_index character(len = :), allocatable:: dir ! directory + + integer, allocatable:: ishape_last(:) ! (d0:) + ! shape index (0-based) in the collection of shapefiles of the last + ! shape at a given date index + + integer d0 ! first date end type Shpc private shpfileobject diff --git a/Common/read_snapshot.f90 b/Common/read_snapshot.f90 index 1ee8d43437c294d309f99dc5b73643b3a1e2b2ea..06a1a2a1f1de52999104952bcf485aaf0824e606 100644 --- a/Common/read_snapshot.f90 +++ b/Common/read_snapshot.f90 @@ -4,8 +4,7 @@ module read_snapshot_m contains - subroutine read_snapshot(s, hshp, nlon, nlat, d_init, k, corner, step, copy, & - ishape_last) + subroutine read_snapshot(s, hshp, nlon, nlat, k, corner, step, copy) ! Libraries: use contour_531, only: convert_to_ind @@ -21,7 +20,6 @@ contains ! size of ssh array in input NetCDF, assuming no repeated point if ! the grid is global - integer, intent(in):: d_init integer, intent(in):: k ! date index real, intent(in):: corner(:) ! (2) longitude and latitude of the @@ -29,7 +27,7 @@ contains real, intent(in):: step(:) ! (2) longitude and latitude steps, in rad - integer, intent(in):: copy, ishape_last(d_init:) + integer, intent(in):: copy ! Local: integer ishape, ishape_first @@ -39,17 +37,17 @@ contains !--------------------------------------------------------------------- - if (k == d_init) then + if (k == hshp%d0) then ishape_first = 0 else - ! {k > d_init} - ishape_first = ishape_last(k - 1) + 1 + ! {k > hshp%d0} + ishape_first = hshp%ishape_last(k - 1) + 1 end if - s%number_extr = ishape_last(k) - ishape_first + 1 + s%number_extr = hshp%ishape_last(k) - ishape_first + 1 allocate(s%list(s%number_extr)) - do ishape = ishape_first, ishape_last(k) + do ishape = ishape_first, hshp%ishape_last(k) call read_eddy(e, date_read, eddy_i, hshp, ishape) ! Check that all the eddies have the same date index: diff --git a/Inst_eddies/Tests/test_nearby_extr.f90 b/Inst_eddies/Tests/test_nearby_extr.f90 index 6fedeb86eeffd07a8e09939835ebd092d2d8d1f9..183e7405fd107bea9ad14f740a000db30498a8c3 100644 --- a/Inst_eddies/Tests/test_nearby_extr.f90 +++ b/Inst_eddies/Tests/test_nearby_extr.f90 @@ -16,7 +16,7 @@ program test_nearby_extr character(len = :), allocatable:: shpc_dir type(snapshot) s TYPE(shpc) hshp - integer d_init, unit, ishape_last, l + integer unit, l real corner_deg(2) ! longitude and latitude of the corner of the whole grid in input @@ -44,14 +44,15 @@ program test_nearby_extr close(unit) call shpc_open(hshp, trim(shpc_dir), pszaccess = "rb") - call dbf_read_attribute_03(d_init, hshp%extremum, hshp%extr_date, ishape = 0) + call dbf_read_attribute_03(hshp%d0, hshp%extremum, hshp%extr_date, ishape = 0) open(unit, file = shpc_dir // "/ishape_last.txt", status = "old", & action = "read", position = "rewind") - read(unit, fmt = *) ishape_last ! first date + allocate(hshp%ishape_last(hshp%d0:hshp%d0)) + read(unit, fmt = *) hshp%ishape_last ! first date close(unit) - call read_snapshot(s, hshp, nlon, nlat, d_init, k = d_init, & + call read_snapshot(s, hshp, nlon, nlat, k = hshp%d0, & corner = corner_deg * deg_to_rad, step = step_deg * deg_to_rad, & - copy = 0, ishape_last = [ishape_last]) + copy = 0) CALL shpc_close(hshp) nearby = nearby_extr(s%extr_map, s%list, i = 4) diff --git a/Overlap/Tests/test_get_dispatch_snap.f90 b/Overlap/Tests/test_get_dispatch_snap.f90 index 3e4d98182c76069f7dbc0d871f24c5581e554383..15da05bbf88c800cec4a8a525d719d93fe1397d2 100644 --- a/Overlap/Tests/test_get_dispatch_snap.f90 +++ b/Overlap/Tests/test_get_dispatch_snap.f90 @@ -23,7 +23,7 @@ program test_get_dispatch_snap character(len = :), allocatable:: shpc_dir type(snapshot) s TYPE(shpc) hshp - integer k_begin, d_init, copy, rank, n_proc, k_end, n_dates + integer k_begin, copy, rank, n_proc, k_end, n_dates integer unit_isolated real:: corner(2) = [0.125, - 59.875] @@ -42,11 +42,6 @@ program test_get_dispatch_snap ! the first extremum. integer:: k = 0 - - integer, allocatable:: ishape_last(:) ! (d_init:d_init + n_dates - 1) - ! shape index (0-based) in the collection of shapefiles of the last - ! shape at a given date index - namelist /main_nml/ corner, step, nlon, nlat, dist_lim, k !------------------------------------------------------------------------- @@ -80,11 +75,11 @@ program test_get_dispatch_snap call shpc_open(hshp, trim(shpc_dir), pszaccess = "rb") if (rank == 0) then - call dbf_read_attribute_03(d_init, hshp%extremum, hshp%extr_date, & + call dbf_read_attribute_03(hshp%d0, hshp%extremum, hshp%extr_date, & ishape = 0) - call read_column(ishape_last, & - file = trim(shpc_dir) // "/ishape_last.txt", my_lbound = d_init) - n_dates = size(ishape_last) + call read_column(hshp%ishape_last, & + file = trim(shpc_dir) // "/ishape_last.txt", my_lbound = hshp%d0) + n_dates = size(hshp%ishape_last) call new_unit(unit_isolated) open(unit_isolated, file = "isolated_nodes.txt", status = "replace", & action = "write") @@ -97,19 +92,19 @@ program test_get_dispatch_snap call ezmpi_bcast(k, root = 0) call ezmpi_bcast(copy, root = 0) call ezmpi_bcast(n_dates, root = 0) - if (rank /= 0) allocate(ishape_last(d_init:d_init + n_dates - 1)) - call ezmpi_bcast(ishape_last, root = 0) - call ezmpi_bcast(d_init, root = 0) - k_begin = d_init + (rank * n_dates) / n_proc + if (rank /= 0) allocate(hshp%ishape_last(hshp%d0:hshp%d0 + n_dates - 1)) + call ezmpi_bcast(hshp%ishape_last, root = 0) + call ezmpi_bcast(hshp%d0, root = 0) + k_begin = hshp%d0 + (rank * n_dates) / n_proc if (rank < n_proc - 1) then - k_end = d_init + ((rank + 1) * n_dates) / n_proc + k_end = hshp%d0 + ((rank + 1) * n_dates) / n_proc else - k_end = d_init + n_dates - 1 + k_end = hshp%d0 + n_dates - 1 end if - call get_snapshot(s, nlon, nlat, d_init, ishape_last, corner * deg_to_rad, & - step * deg_to_rad, copy, hshp, k, k_end, rank, n_proc, max_delta = 1) + 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) 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 f6f073f0ce0d9316de3ec93efdb6d76a198f6ba7..689b9b828ce6efb0244b9d60f413ba08629617d1 100644 --- a/Overlap/Tests/test_overlap.f90 +++ b/Overlap/Tests/test_overlap.f90 @@ -17,10 +17,8 @@ program test_overlap implicit none character(len = :), allocatable:: shpc_dir - integer d_init integer:: k_test_1 = 0, k_test_2 = 1 integer unit, i, copy, n_dates - integer, allocatable:: ishape_last(:) type(snapshot), allocatable:: flow(:) ! (max_delta + 1) TYPE(shpc) hshp @@ -67,21 +65,21 @@ program test_overlap step = step_deg * deg_to_rad allocate(flow(max_delta + 1)) call shpc_open(hshp, trim(shpc_dir), pszaccess = "rb") - call dbf_read_attribute_03(d_init, hshp%extremum, hshp%extr_date, ishape = 0) - call read_column(ishape_last, file = trim(shpc_dir) // "/ishape_last.txt", & - my_lbound = d_init) - n_dates = size(ishape_last) - call assert(d_init <= [k_test_1, k_test_2] .and. [k_test_1, k_test_2] & - < d_init + n_dates, "test_overlap k_test_1, k_test_2") - e_overestim = maxval([ishape_last(d_init) + 1, ediff1d(ishape_last)]) + call dbf_read_attribute_03(hshp%d0, hshp%extremum, hshp%extr_date, ishape = 0) + call read_column(hshp%ishape_last, & + file = trim(shpc_dir) // "/ishape_last.txt", my_lbound = hshp%d0) + n_dates = size(hshp%ishape_last) + call assert(hshp%d0 <= [k_test_1, k_test_2] .and. [k_test_1, k_test_2] & + < hshp%d0 + n_dates, "test_overlap k_test_1, k_test_2") + e_overestim = maxval([hshp%ishape_last(hshp%d0) + 1, & + ediff1d(hshp%ishape_last)]) open(unit, file = "node_id_param.json", status = "replace", & action = "write") write(unit, fmt = *) '{"e_overestim": ', e_overestim, '}' close(unit) - call read_snapshot(flow(1), hshp, nlon, nlat, d_init, k_test_1, corner, & - step, copy, ishape_last) - call read_snapshot(flow(max_delta + 1), hshp, nlon, nlat, d_init, k_test_2, & - corner, step, copy, ishape_last) + 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 a2b8ba54207d2189ae5320538890da78be03ed22..1eb2481cf1a7700c2898616f598f99e18753cc0e 100644 --- a/Overlap/Tests/test_read_snapshot.f90 +++ b/Overlap/Tests/test_read_snapshot.f90 @@ -18,7 +18,7 @@ program test_read_snapshot character(len = :), allocatable:: shpc_dir type(snapshot) s TYPE(shpc) hshp - integer d_init, copy, unit, ishape_last + integer copy, unit real:: corner_deg(2) = [0.125, - 59.875] ! longitude and latitude of the corner of the whole grid in input @@ -62,16 +62,17 @@ program test_read_snapshot copy = merge(dist_lim, 0, periodic) call shpc_open(hshp, trim(shpc_dir), pszaccess = "rb") - call dbf_read_attribute_03(d_init, hshp%extremum, hshp%extr_date, ishape = 0) + call dbf_read_attribute_03(hshp%d0, hshp%extremum, hshp%extr_date, ishape = 0) open(unit, file = shpc_dir // "/ishape_last.txt", status = "old", & action = "read", position = "rewind") - read(unit, fmt = *) ishape_last ! first date + allocate(hshp%ishape_last(hshp%d0:hshp%d0)) + read(unit, fmt = *) hshp%ishape_last ! first date close(unit) - call read_snapshot(s, hshp, nlon, nlat, d_init, k = d_init, & + call read_snapshot(s, hshp, nlon, nlat, k = hshp%d0, & corner = corner_deg * deg_to_rad, step = step_deg * deg_to_rad, & - copy = copy, ishape_last = [ishape_last]) + copy = copy) CALL shpc_close(hshp) - call write_snapshot(s, corner_deg, step_deg, nlon, nlat, copy, d_init) + call write_snapshot(s, corner_deg, step_deg, nlon, nlat, copy, hshp%d0) end program test_read_snapshot diff --git a/Overlap/Tests/test_send_recv.f90 b/Overlap/Tests/test_send_recv.f90 index 99a270a9d3f8299caf5f0bb7ee4099ae9bf701f0..bb952071272c7c317486e1c4707bc46701dad6c0 100644 --- a/Overlap/Tests/test_send_recv.f90 +++ b/Overlap/Tests/test_send_recv.f90 @@ -24,7 +24,7 @@ program test_send_recv character(len = :), allocatable:: shpc_dir type(snapshot) s - integer rank, n_proc, d_init, copy, n_dates, unit + integer rank, n_proc, copy, n_dates, unit logical flag INTEGER(KIND=MPI_ADDRESS_KIND) attribute_val @@ -43,10 +43,6 @@ program test_send_recv logical periodic ! grid is periodic in longitude TYPE(shpc) hshp - integer, allocatable:: ishape_last(:) ! (d_init:d_init + n_dates - 1) - ! shape index (0-based) in the collection of shapefiles of the last - ! shape at a given date index - !--------------------------------------------------------------------- call mpi_init @@ -89,11 +85,11 @@ program test_send_recv call shpc_open(hshp, trim(shpc_dir), pszaccess = "rb") if (rank == 0) then - call dbf_read_attribute_03(d_init, hshp%extremum, hshp%extr_date, & + call dbf_read_attribute_03(hshp%d0, hshp%extremum, hshp%extr_date, & ishape = 0) - call read_column(ishape_last, & - file = trim(shpc_dir) // "/ishape_last.txt", my_lbound = d_init) - n_dates = size(ishape_last) + call read_column(hshp%ishape_last, & + file = trim(shpc_dir) // "/ishape_last.txt", my_lbound = hshp%d0) + n_dates = size(hshp%ishape_last) end if call ezmpi_bcast(corner_deg, root = 0) @@ -102,24 +98,24 @@ program test_send_recv call ezmpi_bcast(nlat, root = 0) call ezmpi_bcast(copy, root = 0) call ezmpi_bcast(n_dates, root = 0) - if (rank == 1) allocate(ishape_last(d_init:d_init + n_dates - 1)) - call ezmpi_bcast(ishape_last, root = 0) + if (rank == 1) allocate(hshp%ishape_last(hshp%d0:hshp%d0 + n_dates - 1)) + call ezmpi_bcast(hshp%ishape_last, root = 0) - call ezmpi_bcast(d_init, root = 0) - if (rank == 1) call read_snapshot(s, hshp, nlon, nlat, d_init, k = d_init, & + call ezmpi_bcast(hshp%d0, root = 0) + 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, ishape_last = ishape_last) + copy = copy) CALL shpc_close(hshp) if (rank == 1) then - call send_snapshot(s, dest = 0, tag = d_init) + call send_snapshot(s, dest = 0, tag = hshp%d0) else ! rank == 0 call MPI_Comm_get_attr(MPI_Comm_world, MPI_TAG_UB, attribute_val, flag) call assert(flag, "test_send_recv MPI_Comm_get_attr MPI_TAG_UB") print *, "MPI_TAG_UB = ", attribute_val - call recv_snapshot(s, nlon, nlat, copy, source = 1, tag = d_init) - call write_snapshot(s, corner_deg, step_deg, nlon, nlat, copy, d_init) + call recv_snapshot(s, nlon, nlat, copy, source = 1, tag = hshp%d0) + call write_snapshot(s, corner_deg, step_deg, nlon, nlat, copy, hshp%d0) end if call mpi_finalize diff --git a/Overlap/eddy_graph.f90 b/Overlap/eddy_graph.f90 index c129e4576dd92f58e683ace756e37cf9e24d49a9..985fb873a61c82bc17babb7f5832a3db939519e0 100644 --- a/Overlap/eddy_graph.f90 +++ b/Overlap/eddy_graph.f90 @@ -23,7 +23,7 @@ program eddy_graph integer rank, n_proc, copy, delta, j integer:: n_dates = huge(0) ! number of dates to read integer k ! date index - integer d_init, k_begin, k_end, k_end_main_loop + integer k_begin, k_end, k_end_main_loop character(len = :), allocatable:: shpc_dir integer:: max_delta = 1 @@ -49,11 +49,6 @@ program eddy_graph ! the first extremum. logical periodic ! grid is periodic in longitude - - integer, allocatable:: ishape_last(:) ! (d_init:d_init + n_dates - 1) - ! shape index (0-based) in the collection of shapefiles of the last - ! 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 @@ -93,17 +88,18 @@ program eddy_graph "eddy_graph dist_lim") copy = merge(dist_lim, 0, periodic) - call dbf_read_attribute_03(d_init, hshp%extremum, hshp%extr_date, & + call dbf_read_attribute_03(hshp%d0, hshp%extremum, hshp%extr_date, & ishape = 0) - call read_column(ishape_last, & + call read_column(hshp%ishape_last, & file = trim(shpc_dir) // "/ishape_last.txt", last = n_dates, & - my_lbound = d_init) - if (n_dates == huge(0)) n_dates = size(ishape_last) + my_lbound = hshp%d0) + if (n_dates == huge(0)) n_dates = size(hshp%ishape_last) 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([ishape_last(d_init) + 1, ediff1d(ishape_last)]) + e_overestim = maxval([hshp%ishape_last(hshp%d0) + 1, & + ediff1d(hshp%ishape_last)]) corner = corner_deg * deg_to_rad step = step_deg * deg_to_rad end if @@ -117,10 +113,10 @@ program eddy_graph call ezmpi_bcast(periodic, root = 0) call ezmpi_bcast(copy, root = 0) call ezmpi_bcast(n_dates, root = 0) - if (rank /= 0) allocate(ishape_last(d_init:d_init + n_dates - 1)) - call ezmpi_bcast(ishape_last, root = 0) + if (rank /= 0) allocate(hshp%ishape_last(hshp%d0:hshp%d0 + n_dates - 1)) + call ezmpi_bcast(hshp%ishape_last, root = 0) call ezmpi_bcast(e_overestim, root = 0) - call ezmpi_bcast(d_init, root = 0) + call ezmpi_bcast(hshp%d0, root = 0) if (rank == 0) then open(unit, file = "node_id_param.json", status = "replace", & @@ -138,13 +134,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 = d_init + (rank * n_dates) / n_proc + k_begin = hshp%d0 + (rank * n_dates) / n_proc if (rank < n_proc - 1) then - k_end = d_init + ((rank + 1) * n_dates) / n_proc + max_delta - 1 + k_end = hshp%d0 + ((rank + 1) * n_dates) / n_proc + max_delta - 1 k_end_main_loop = k_end - max_delta + 1 else - k_end = d_init + n_dates - 1 + k_end = hshp%d0 + n_dates - 1 k_end_main_loop = k_end end if @@ -153,8 +149,8 @@ program eddy_graph ! 2. Prologue: do k = k_begin, k_begin + max_delta - call get_snapshot(flow(k - k_begin + 1), nlon, nlat, d_init, ishape_last, & - corner, step, copy, hshp, k, k_end, rank, n_proc, max_delta) ! (read) + call get_snapshot(flow(k - k_begin + 1), nlon, nlat, corner, step, copy, & + hshp, k, k_end, rank, n_proc, max_delta) ! (read) end do do delta = 1, max_delta @@ -170,8 +166,8 @@ program eddy_graph call dispatch_snapshot(flow(1), unit_isolated, rank, k_begin, max_delta, & k = k - max_delta - 1) flow(:max_delta) = flow(2:) - call get_snapshot(flow(max_delta + 1), nlon, nlat, d_init, ishape_last, & - corner, step, copy, hshp, k, k_end, rank, n_proc, max_delta) + call get_snapshot(flow(max_delta + 1), nlon, nlat, corner, step, copy, & + hshp, k, k_end, rank, n_proc, max_delta) do delta = 1, max_delta call overlap(flow, nlon, nlat, periodic, dist_lim, e_overestim, k, & @@ -186,8 +182,8 @@ program eddy_graph call dispatch_snapshot(flow(1), unit_isolated, rank, k_begin, max_delta, & k = k - max_delta - 1) flow(:max_delta) = flow(2:) - call get_snapshot(flow(max_delta + 1), nlon, nlat, d_init, ishape_last, & - corner, step, copy, hshp, k, k_end, rank, n_proc, max_delta) + call get_snapshot(flow(max_delta + 1), nlon, nlat, corner, step, copy, & + hshp, k, k_end, rank, n_proc, max_delta) ! (reception) ! Stitching: diff --git a/Overlap/get_snapshot.f90 b/Overlap/get_snapshot.f90 index 28b1d351b94d2a5ff314e12aafb4557188d98ca9..ce957639820ff528bcc3ddb77c44ef2cb2991fb8 100644 --- a/Overlap/get_snapshot.f90 +++ b/Overlap/get_snapshot.f90 @@ -4,8 +4,8 @@ module get_snapshot_m contains - subroutine get_snapshot(s, nlon, nlat, d_init, ishape_last, corner, step, & - copy, hshp, k, k_end, rank, n_proc, max_delta) + subroutine get_snapshot(s, nlon, nlat, corner, step, copy, hshp, k, k_end, & + rank, n_proc, max_delta) use derived_types, only: snapshot, shpc use read_snapshot_m, only: read_snapshot @@ -13,11 +13,6 @@ contains type(snapshot), intent(out):: s integer, intent(in):: nlon, nlat - integer, intent(in):: d_init ! first date in the collection of shapefiles - - integer, intent(in):: ishape_last(d_init:) - ! shape index (0-based) in the collection of shapefiles of the last - ! shape at a given date index real, intent(in):: corner(:) ! (2) longitude and latitude of the ! corner of the whole grid, in rad @@ -37,8 +32,7 @@ contains !-------------------------------------------------------------- if (rank == n_proc - 1 .or. k <= k_end - max_delta) then - call read_snapshot(s, hshp, nlon, nlat, d_init, k, corner, step, copy, & - ishape_last) + call read_snapshot(s, hshp, nlon, nlat, k, corner, step, copy) else call recv_snapshot(s, nlon, nlat, copy, source = rank + 1, tag = k) end if