Skip to content
Snippets Groups Projects
Commit 6c937242 authored by Lionel GUEZ's avatar Lionel GUEZ
Browse files

Create subroutines `config_graph` and `read_grid`

Which read namelists. Motivations: diminish the number of arguments of
procedure overlap; lighten the main program unit. We conditionally
compile MPI calls in `read_grid` because we want to use variables of
module`read_grid_m` in procedure `read_snapshot`, and `read_snapshot`
is called from sequential program `test_nearby_extr`.
parent a56f5269
No related branches found
No related tags found
No related merge requests found
module read_grid_m
implicit none
integer, protected:: nlon = - 1, nlat = - 1
! size of ssh array in input NetCDF, assuming no repeated point if
! the grid is global
logical, protected:: periodic ! grid is periodic in longitude
real, protected:: corner(2)
! longitude and latitude of the corner of the whole grid in input
! NetCDF, in rad
real, protected:: step(2) ! longitude and latitude steps, in rad
contains
subroutine read_grid(rank, shpc_dir)
! Libraries:
#ifndef SEQUENTIAL
use ezmpi, only: ezmpi_bcast
#endif
use jumble, only: new_unit, deg_to_rad
integer, intent(in):: rank
character(len = *), intent(in):: shpc_dir
! Local:
integer unit
real:: corner_deg(2) = [huge(0.), huge(0.)]
! longitude and latitude of the corner of the whole grid in input
! NetCDF, in degrees
real:: step_deg(2) = [huge(0.), huge(0.)]
! longitude and latitude steps, in degrees
namelist /grid_nml/ corner_deg, step_deg, nlon, nlat
!-----------------------------------------------------------------
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)
! 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
corner = corner_deg * deg_to_rad
step = step_deg * deg_to_rad
end if
#ifndef SEQUENTIAL
call ezmpi_bcast(nlon, root = 0)
call ezmpi_bcast(nlat, root = 0)
call ezmpi_bcast(periodic, root = 0)
call ezmpi_bcast(corner, root = 0)
call ezmpi_bcast(step, root = 0)
#endif
end subroutine read_grid
end module read_grid_m
module config_graph_m
implicit none
integer, protected:: copy
integer, protected:: max_delta = 1
integer, protected:: dist_lim = 12
! We look for an overlapping eddy at dist_lim (in grid points) of
! the first extremum.
logical, protected:: cyclone = .true.
contains
subroutine config_graph(rank)
use, intrinsic:: ISO_FORTRAN_ENV
! Libraries:
use ezmpi, only: ezmpi_bcast
use jumble, only: assert, pi
use read_grid_m, only: periodic, step
integer, intent(in):: rank
! Local:
namelist /config_graph_nml/ max_delta, dist_lim, cyclone
!-----------------------------------------------------------------
if (rank == 0) then
! config_graph_nml:
write(unit = error_unit, nml = config_graph_nml)
write(unit = error_unit, fmt = *) "Enter namelist config_graph_nml."
read(unit = *, nml = config_graph_nml)
write(unit = *, nml = config_graph_nml)
call assert(max_delta >= 1, "config_graph max_delta")
if (periodic) call assert(2 * dist_lim * step(1) < PI, &
"config_graph dist_lim")
copy = merge(dist_lim, 0, periodic)
end if
call ezmpi_bcast(max_delta, root = 0)
call ezmpi_bcast(dist_lim, root = 0)
call ezmpi_bcast(cyclone, root = 0)
call ezmpi_bcast(copy, root = 0)
end subroutine config_graph
end module config_graph_m
...@@ -5,50 +5,28 @@ program eddy_graph ...@@ -5,50 +5,28 @@ program eddy_graph
! Libraries: ! Libraries:
use ezmpi, only: ezmpi_bcast use ezmpi, only: ezmpi_bcast
use jumble, only: get_command_arg_dyn, new_unit, ediff1d, assert, & use jumble, only: get_command_arg_dyn, new_unit, ediff1d, assert, &
deg_to_rad, read_opcol read_opcol
use mpi_f08, only: mpi_init, mpi_finalize, MPI_Comm_rank, MPI_Comm_world, & use mpi_f08, only: mpi_init, mpi_finalize, MPI_Comm_rank, MPI_Comm_world, &
MPI_Comm_size MPI_Comm_size
use shapelib_03, only: dbf_read_attribute_03 use shapelib_03, only: dbf_read_attribute_03
use config_graph_m, only: config_graph
use derived_types, only: snapshot, shpc_slice_handler, shpc_slice_meta use derived_types, only: snapshot, shpc_slice_handler, shpc_slice_meta
use dispatch_snapshot_m, only: dispatch_snapshot use dispatch_snapshot_m, only: dispatch_snapshot
use get_snapshot_m, only: get_snapshot use get_snapshot_m, only: get_snapshot
use overlap_m, only: overlap use overlap_m, only: overlap
use read_grid_m, only: read_grid
use unit_edge_m, only: open_edge_file, unit_edge use unit_edge_m, only: open_edge_file, unit_edge
use shpc_open_m, only: shpc_open use shpc_open_m, only: shpc_open
use shpc_close_m, only: shpc_close use shpc_close_m, only: shpc_close
implicit none implicit none
integer rank, n_proc, copy, delta, i, j, iostat integer rank, n_proc, delta, i, j, iostat
integer:: n_dates = huge(0) ! number of dates to read integer:: n_dates = huge(0) ! number of dates to read
integer k ! date index integer k ! date index
integer k_begin, k_end, k_end_main_loop integer k_begin, k_end, k_end_main_loop
character(len = :), allocatable:: shpc_dir 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 unit_isolated, unit integer unit_isolated, unit
integer n_slices ! number of slices in the input SHPC integer n_slices ! number of slices in the input SHPC
integer e_overestim ! over-estimation of the number of eddies at each date integer e_overestim ! over-estimation of the number of eddies at each date
...@@ -56,10 +34,8 @@ program eddy_graph ...@@ -56,10 +34,8 @@ program eddy_graph
type(shpc_slice_meta), allocatable:: ssm(:) ! (n_slices) type(shpc_slice_meta), allocatable:: ssm(:) ! (n_slices)
type(snapshot), allocatable:: flow(:) ! (max_delta + 1) type(snapshot), allocatable:: flow(:) ! (max_delta + 1)
character(len = 30) file character(len = 30) file
logical:: cyclone = .true.
namelist /grid_nml/ corner_deg, step_deg, nlon, nlat namelist /main_nml/ n_dates
namelist /main_nml/ max_delta, dist_lim, n_dates, cyclone
!------------------------------------------------------------------------- !-------------------------------------------------------------------------
...@@ -69,16 +45,13 @@ program eddy_graph ...@@ -69,16 +45,13 @@ program eddy_graph
call MPI_Comm_rank(MPI_Comm_world, rank) call MPI_Comm_rank(MPI_Comm_world, rank)
call MPI_Comm_size(MPI_Comm_world, n_proc) call MPI_Comm_size(MPI_Comm_world, n_proc)
call get_command_arg_dyn(1, shpc_dir, "Required argument: SHPC-directory") call get_command_arg_dyn(1, shpc_dir, "Required argument: SHPC-directory")
call read_grid(rank, shpc_dir)
call config_graph(rank)
if (rank == 0) then 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)
! n_slices: ! n_slices:
call new_unit(unit)
open(unit, file = shpc_dir // "/n_slices.txt", status = "old", & open(unit, file = shpc_dir // "/n_slices.txt", status = "old", &
action = "read", position = "rewind", iostat = iostat) action = "read", position = "rewind", iostat = iostat)
...@@ -96,24 +69,11 @@ program eddy_graph ...@@ -96,24 +69,11 @@ program eddy_graph
write(unit = error_unit, fmt = *) "Enter namelist main_nml." write(unit = error_unit, fmt = *) "Enter namelist main_nml."
read(unit = *, nml = main_nml) read(unit = *, nml = main_nml)
write(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_slices == 1, "If we have several " & call assert(n_dates == huge(0) .or. n_slices == 1, "If we have several " &
// "slices then we read all the dates from them") // "slices 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:
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)
corner = corner_deg * deg_to_rad
step = step_deg * deg_to_rad
end if end if
call ezmpi_bcast(n_slices, root = 0) call ezmpi_bcast(n_slices, root = 0)
call ezmpi_bcast(cyclone, root = 0)
allocate(hshpc(n_slices), ssm(n_slices)) allocate(hshpc(n_slices), ssm(n_slices))
do i = 1, n_slices do i = 1, n_slices
...@@ -139,14 +99,6 @@ program eddy_graph ...@@ -139,14 +99,6 @@ program eddy_graph
close(unit) close(unit)
end if 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) call ezmpi_bcast(n_dates, root = 0)
call ezmpi_bcast(e_overestim, root = 0) call ezmpi_bcast(e_overestim, root = 0)
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment