-
Lionel GUEZ authoredLionel GUEZ authored
eddy_graph.f90 7.09 KiB
program eddy_graph
use, intrinsic:: ISO_FORTRAN_ENV
! Libraries:
use ezmpi, only: ezmpi_bcast
use jumble, only: get_command_arg_dyn, read_column, new_unit
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
use shapelib_03, only: dbf_read_attribute_03
use derived_types, only: snapshot, shpc
use dispatch_snapshot_m, only: dispatch_snapshot
use get_snapshot_m, only: get_snapshot
use overlap_m, only: overlap
use unit_edge_m, only: open_edge_file, unit_edge
use shpc_open_m, only: shpc_open
use shpc_close_m, only: shpc_close
use shpc_create_m, only: shpc_create
use write_overlap_m, only: init_interpolated_eddy
implicit none
integer rank, n_proc, copy, delta, j
integer:: n_dates = huge(0) ! number of dates to read
integer d_init, k_begin, k_end, k, k_end_main_loop
character(len = :), allocatable:: shpc_dir
integer:: max_delta = 1
real:: corner_deg(2) = [huge(0.), huge(0.)]
! longitude and latitude of the corner of the whole grid in input
! NetCDF, in degrees
real corner(2)
! longitude and latitude of the corner of the whole grid in input
! NetCDF, in rad
real:: step_deg(2) = [huge(0.), huge(0.)]
! longitude and latitude steps, in degrees
real step(2) ! longitude and latitude steps, in rad
integer:: nlon = - 1, nlat = - 1
! size of ssh array in input NetCDF, assuming no repeated point if
! the grid is global
integer:: dist_lim = 12
! We look for an overlapping eddy at dist_lim (in grid points) of
! the first extremum.
logical periodic ! grid is periodic in longitude
integer, allocatable:: ishape_last(:)
! shape index (0-based) in the collection of shapefiles of the last
! shape at a given date index
integer unit_isolated, unit_number_eddies, unit
TYPE(shpc) hshp
TYPE(shpc) hshp_interp ! for interpolated eddies
type(snapshot), allocatable:: flow(:)
character(len = 30) file
namelist /grid_nml/ corner_deg, step_deg, nlon, nlat
namelist /main_nml/ max_delta, dist_lim, n_dates
!-------------------------------------------------------------------------
! 1. Front matter
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")
if (rank == 0) then
call new_unit(unit)
open(unit, file = shpc_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")
! 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:
periodic = nint(360. / step_deg(1)) == nlon
print *, "periodic = ", periodic
if (periodic) call assert(2 * dist_lim * step_deg(1) < 180., &
"eddy_graph dist_lim")
copy = merge(dist_lim, 0, periodic)
call read_column(ishape_last, &
file = trim(shpc_dir) // "/ishape_last.txt", last = n_dates)
n_dates = size(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)")
corner = corner_deg * deg_to_rad
step = step_deg * deg_to_rad
end if
call ezmpi_bcast(max_delta, root = 0)
call ezmpi_bcast(corner, root = 0)
call ezmpi_bcast(step, root = 0)
call ezmpi_bcast(nlon, root = 0)
call ezmpi_bcast(nlat, root = 0)
call ezmpi_bcast(dist_lim, root = 0)
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(n_dates))
call ezmpi_bcast(ishape_last, root = 0)
call shpc_open(hshp, trim(shpc_dir), rank)
if (rank == 0) call dbf_read_attribute_03(d_init, hshp%extremum, &
hshp%extr_date, ishape = 0)
call ezmpi_bcast(d_init, root = 0)
! Open output files:
call new_unit(unit_isolated)
write(unit = file, fmt = '("isolated_nodes_", i0, ".txt")') rank
open(unit_isolated, file = trim(file), status = "replace", action = "write")
call new_unit(unit_number_eddies)
write(unit = file, fmt = '("number_eddies_", i0, ".csv")') rank
open(unit_number_eddies, file = trim(file), status = "replace", &
action = "write")
call open_edge_file(rank)
! We do not write the title line. That will be handled by eddy_graph.sh.
write(unit = file, fmt = '("SHPC_", i0)') rank
call shpc_create(hshp_interp, shpc_dir = trim(file), cyclone = hshp%cyclone)
call init_interpolated_eddy
k_begin = d_init + (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_main_loop = k_end - max_delta + 1
else
k_end = d_init + n_dates - 1
k_end_main_loop = k_end
end if
allocate(flow(max_delta + 1))
! 2. Prologue:
do k = k_begin, k_begin + max_delta
call get_snapshot(flow(k - k_begin + 1), nlon, nlat, ishape_last, corner, &
step, copy, hshp, d_init, k, k_end, rank, n_proc, max_delta) ! (read)
end do
do delta = 1, max_delta
do k = k_begin + delta, k_begin + max_delta
call overlap(flow, nlon, nlat, periodic, dist_lim, hshp_interp, k, &
delta, j = k - k_begin + 1)
end do
end do
! 3. Main loop
do k = k_begin + max_delta + 1, k_end_main_loop
call dispatch_snapshot(flow(1), unit_isolated, unit_number_eddies, rank, &
k_begin, max_delta, k = k - max_delta - 1)
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)
do delta = 1, max_delta
call overlap(flow, nlon, nlat, periodic, dist_lim, hshp_interp, k, &
delta, j = max_delta + 1)
end do
end do
! 4. Epilogue
do k = k_end_main_loop + 1, k_end
! {rank < n_proc - 1 and k >= k_end - max_delta + 2}
call dispatch_snapshot(flow(1), unit_isolated, unit_number_eddies, rank, &
k_begin, max_delta, k = k - max_delta - 1)
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)
! (reception)
! Stitching:
do delta = k - k_end + max_delta, max_delta
call overlap(flow, nlon, nlat, periodic, dist_lim, hshp_interp, k, &
delta, j = max_delta + 1)
end do
end do
do j = 1, max_delta + 1
call dispatch_snapshot(flow(j), unit_isolated, unit_number_eddies, rank, &
k_begin, max_delta, k = k_end - max_delta - 1 + j)
end do
! 5. Back matter
CALL shpc_close(hshp_interp)
CALL shpc_close(hshp)
close(unit_isolated)
close(unit_number_eddies) close(unit_edge)
call mpi_finalize
end program eddy_graph