-
Lionel GUEZ authoredLionel GUEZ authored
inst_eddies.f90 6.40 KiB
program inst_eddies
! The current directory at run-time should contain directory SHPC.
use, intrinsic:: ieee_arithmetic, only: IEEE_SUPPORT_DATATYPE, &
ieee_support_nan
use, intrinsic:: ISO_FORTRAN_ENV
! Libraries:
use jumble, only: assert, new_unit, pi, argwhere, twopi
use numer_rec_95, only: indexx
use config_m, only: config, min_radius
use derived_types, only: snapshot, shpc_slice_handler
use input_ssh_m, only: input_ssh, max_radius
use nearby_extr_m, only: nearby_extr
use set_all_extr_m, only: set_all_extr
use set_contours_m, only: set_contours
use shpc_close_m, only: shpc_close
use shpc_create_m, only: shpc_create
use shpc_open_m, only: shpc_open
use write_snapshot_m, only: write_snapshot
implicit none
type(snapshot) s
TYPE(shpc_slice_handler) hshpc_cyclo, hshpc_anti
integer iostat, i, l, nlon, nlat
character(len = 200) iomsg
integer unit
real corner(2)
! longitude and latitude of the corner of the whole grid, in rad
real step(2) ! longitude and latitude steps, in rad
logical periodic ! grid is periodic in longitude
real, allocatable:: ssh(:, :)
! (1 - max_radius(1):nlon + max_radius(1), nlat) if the grid is periodic
! in longitude, else (nlon, nlat). Sea-surface height, in m.
real, allocatable:: u(:, :), v(:, :)
! (1 - max_radius(1):nlon + max_radius(1), nlat) if periodic, else
! (nlon, nlat) wind, in m s-1
logical exist
real t0, t1 ! CPU times, in s
character(len = 30) timings_fname ! file name
integer:: date = 0, slice = 0
namelist /dates_nml/ date, slice
integer, allocatable:: sorted_extr(:) ! (s%number_extr)
! Sorted identifying number of extrema: first, minima sorted in
! order of decreasing SSH, and second, maxima sorted in order of
! increasing SSH.
real min_area ! minimum area of an outermost contour, in m2
integer, allocatable:: selection(:)
! identifying numbers of a selection of extrema
integer n_cycl ! number of cyclones
integer llc(2) ! indices in global grid of lower left corner
integer urc(2) ! indices in global grid of upper right corner
real corner_window(2) ! longitude and latitude of the window around each
! extremum, in rad
real, allocatable:: outside_points(:, :) ! (2, :) longitude and
! latitude, in rad, of all the significant extrema, except the ! target extremum
!--------------------------------------------------------------
call assert(IEEE_SUPPORT_DATATYPE(0.), ieee_support_nan(0.), &
"inst_eddies: not enough IEEE support")
call cpu_time(t0)
print *, "inst_eddies: Enter dates_nml:"
read(unit = *, nml = dates_nml)
write(timings_fname, fmt = '("SHPC/Slice_", i0, "/timings.txt")') slice
inquire(file = timings_fname, exist = exist)
call new_unit(unit)
if (exist) then
open(unit, file = timings_fname, STATUS = "old", ACTION = "WRITE", &
position = "append")
else
OPEN(unit, file = timings_fname, STATUS = "new", ACTION = "WRITE", &
iostat = iostat, iomsg = iomsg)
if (iostat /= 0) then
write(ERROR_UNIT, fmt = "(a, i0, a)") "Is the directory SHPC/Slice_", &
slice, " created?"
write(ERROR_UNIT, fmt = *) trim(iomsg)
stop 1
end if
end if
call config
call input_ssh(corner, step, periodic, ssh, u, v)
! Open or create now so we test the existence of directories:
call shpc_open(hshpc_cyclo, shpc_dir = "SHPC", cyclone = .true., &
slice = slice, pszaccess = "rb+", iostat = iostat)
if (iostat == 0) then
call shpc_open(hshpc_anti, shpc_dir = "SHPC", cyclone = .false., &
slice = slice, pszaccess = "rb+")
else
print *, "inst_eddies: Opening failed so we will create a slice..."
call shpc_create(hshpc_cyclo, shpc_dir = "SHPC", cyclone = .true., &
slice = slice, grid_lon_lat = .true.)
call shpc_create(hshpc_anti, shpc_dir = "SHPC", cyclone = .false., &
slice = slice, grid_lon_lat = .true.)
end if
call cpu_time(t1)
write(unit, fmt = *) "CPU time for preamble, before computation:", t1 - t0, &
"s"
t0 = t1
call set_all_extr(s, step, periodic, ssh, corner)
allocate(sorted_extr(s%number_extr))
min_area = pi * (min_radius * 1e3)**2
forall (i = 1:s%number_extr)
s%list(i)%out_cont%closed = .true.
! must be intialized to true because it is used in nearby_extr
end forall
! Double sort of extrema on cyclonicity and SSH value:
selection = argwhere(s%list%cyclone)
n_cycl = size(selection)
selection = selection(indexx(s%list(selection)%extr%ssh))
sorted_extr(:n_cycl) = selection(n_cycl:1:- 1) ! descending order of ssh
selection = argwhere(.not. s%list%cyclone)
selection = selection(indexx(s%list(selection)%extr%ssh)) sorted_extr(n_cycl + 1:) = selection
! Done sorting
if (.not. periodic) nlon = size(ssh, 1)
nlat = size(ssh, 2)
loop_extr: do l = 1, s%number_extr
i = sorted_extr(l)
associate (e => s%list(i))
! Define the geographical window around each eddy extremum:
llc = e%extr%coord_proj - max_radius
urc = e%extr%coord_proj + max_radius
llc(2) = max(llc(2), 1)
urc(2) = min(urc(2), nlat)
if (.not. periodic) then
llc(1) = max(llc(1), 1)
urc(1) = min(urc(1), nlon)
end if
corner_window = corner + (llc - 1) * step
end associate
outside_points = nearby_extr(s%extr_map(llc(1):urc(1), llc(2):urc(2)), &
s%list, i)
! Shift the longitude of each point to a value close to the
! longitude of the target extremum:
if (periodic) outside_points(1, :) = outside_points(1, :) &
+ ceiling((corner_window(1) - outside_points(1, :)) / twopi) &
* twopi
call set_contours(s, step, ssh(llc(1):urc(1), llc(2):urc(2)), &
u(llc(1):urc(1), llc(2):urc(2)), v(llc(1):urc(1), llc(2):urc(2)), &
corner_window, min_area, i, outside_points)
end do loop_extr
call cpu_time(t1)
write(unit, fmt = *) "CPU time for computation, before output:", t1 - t0, "s"
t0 = t1
print *, "inst_eddies: s%number_extr = ", s%number_extr
call write_snapshot(s, hshpc_cyclo, hshpc_anti, date)
CALL shpc_close(hshpc_cyclo)
CALL shpc_close(hshpc_anti)
if (iostat == 0) then
print *, 'inst_eddies: Appended to shapefiles in one slice of SHPC.'
else
print *, 'inst_eddies: Created shapefiles in one slice of SHPC.'
end if
call cpu_time(t1)
write(unit, fmt = *) "CPU time for output:", t1 - t0, "s"
close(unit)
end program inst_eddies