-
Lionel GUEZ authored
Preparing for using hshpc as argument of `write_aux`.
Lionel GUEZ authoredPreparing for using hshpc as argument of `write_aux`.
inst_eddies.f90 7.73 KiB
program inst_eddies
use, intrinsic:: ieee_arithmetic, only: IEEE_SUPPORT_DATATYPE, &
ieee_support_nan, ieee_value, IEEE_QUIET_NAN
! Libraries:
use contour_531, only: convert_to_ind
use netcdf, only: nf90_nowrite
use netcdf95, only: nf95_open, find_coord, nf95_inquire_dimension, &
nf95_get_var, nf95_close, nf95_get_att
use nr_util, only: assert, deg_to_rad, twopi, pi
use config_m, only: config, max_radius_deg, min_radius
use derived_types, only: snapshot, null_ssh_contour, missing_speed, shpc
use get_var_m, only: get_var
use nearby_extr_m, only: nearby_extr
use set_all_outerm_m, only: set_all_outerm
use set_max_speed_m, only: set_max_speed
use shpc_close_m, only: shpc_close
use shpc_create_m, only: shpc_create
use shpc_open_m, only: shpc_open
use write_aux_m, only: write_aux
use write_eddy_m, only: write_eddy
implicit none
type(snapshot) s
TYPE(shpc) hshpc_cyclo, hshpc_anti
integer i, days_1950, n_cyclo, n_anti
integer ncid, dimid, varid, copy
integer nlon, nlat
! size of ssh array in input NetCDF, assuming no repeated point if
! the grid is global
real lon_max, lat_max ! longitude and latitude, in degrees
real step_deg(2), step(2) ! longitude and latitude steps, in degrees and rad
logical periodic ! grid is periodic in longitude
real time
integer max_radius(2)
! maximum radius of an eddy in longitude and latitude, in number of
! grid points
real, allocatable:: ssh(:, :)
! (1 - max_radius(1):nlon + max_radius(1), nlat) if periodic, 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
real corner_deg(2), corner(2)
! longitude and latitude of the corner of the whole grid, in degrees
! and in rad
! Window around each extremum:
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
character(len = 30) time_unit logical exist
!--------------------------------------------------------------
call assert(IEEE_SUPPORT_DATATYPE(0.), ieee_support_nan(0.), &
"inst_eddies: not enough IEEE support")
call config ! We need main_nml before looking at the input files.
call nf95_open("h.nc", nf90_nowrite, ncid)
call find_coord(ncid, varid = varid, std_name = "time")
call nf95_get_att(ncid, varid, name = "units", values = time_unit)
call assert(time_unit == "days since 1950-01-01 00:00:00" &
.or. time_unit == "days since 1950-1-1 00:00:00", "main: bad time unit")
! We are assuming there is a single date in the input file:
call nf95_get_var(ncid, varid, time)
days_1950 = nint(time)
call assert(abs(time - days_1950) < 0.1, "main: bad time value")
call find_coord(ncid, dimid = dimid, varid = varid, std_name = "longitude")
call nf95_inquire_dimension(ncid, dimid, nclen = nlon)
call nf95_get_var(ncid, varid, corner_deg(1))
call nf95_get_var(ncid, varid, lon_max, start = [nlon])
call find_coord(ncid, dimid = dimid, varid = varid, std_name = "latitude")
call nf95_inquire_dimension(ncid, dimid, nclen = nlat)
call nf95_get_var(ncid, varid, corner_deg(2))
call nf95_get_var(ncid, varid, lat_max, start = [nlat])
call assert(lon_max > corner_deg(1) .and. lat_max > corner_deg(2), &
"inst_eddies coordinate order")
step_deg = [(lon_max - corner_deg(1)) / (nlon - 1), &
(lat_max - corner_deg(2)) / (nlat - 1)]
step = step_deg * deg_to_rad
! As we are requiring the grid spacing to be uniform, the value of
! "periodic" may be deduced from the values of step(1) and nlon:
periodic = nint(twopi / step(1)) == nlon
print *, "inst_eddies: periodic = ", periodic
call assert(2 * max_radius_deg(1) < 180., "inst_eddies max_radius_deg")
! (Let us require this even if not periodic. This is clearer.)
max_radius = nint(max_radius_deg / step_deg)
copy = merge(max_radius(1), 0, periodic)
allocate(ssh(1 - copy:nlon + copy, nlat), u(1 - copy:nlon + copy, nlat), &
v(1 - copy:nlon + copy, nlat))
corner = corner_deg * deg_to_rad
! Read ssh, u and v:
call get_var(periodic, max_radius(1), ssh, ncid, nlon, name = "adt", &
new_fill_value = huge(0.))
! (We cannot keep the original fill value because Contour_531
! works with an upper limit for valid values.)
call nf95_close(ncid)
call nf95_open("uv.nc", nf90_nowrite, ncid)
call get_var(periodic, max_radius(1), u, ncid, nlon, name = "ugos", &
new_fill_value = ieee_value(0., IEEE_QUIET_NAN))
call get_var(periodic, max_radius(1), v, ncid, nlon, name = "vgos", &
new_fill_value = ieee_value(0., IEEE_QUIET_NAN))
! (We will need quiet NaNs rather the original fill values for u and
! v when we compute the max-speed contours and when we search the
! ssh of max-speed contours.)
call nf95_close(ncid)
call set_all_outerm(s, max_radius, step, periodic, ssh, corner, &
min_area = pi * (min_radius * 1e3)**2)
! Done with outermost contours, now let us take care of the
! max-speed contours.
do i = 1, s%number_vis_extr
if (s%list_vis(i)%valid) then
! Restrict the field to the outermost contour:
llc = floor(convert_to_ind(minval(s%list_vis(i)%out_cont%points, &
dim = 2), corner, step))
urc = ceiling(convert_to_ind(maxval( &
s%list_vis(i)%out_cont%points, dim = 2), corner, step))
! Should have no effect except because of roundup error:
urc(2) = min(urc(2), nlat)
if (.not. periodic) urc(1) = min(urc(1), nlon)
corner_window = corner + (llc - 1) * step
outside_points = nearby_extr(s%extr_map(llc(1):urc(1), &
llc(2):urc(2)), s%list_vis, i)
if (periodic) outside_points(1, :) = outside_points(1, :) &
+ ceiling((corner_window(1) - outside_points(1, :)) / twopi) &
* twopi
! (Shift the longitude of each point to a value close to
! the longitude of the target extremum.)
call set_max_speed(s%list_vis(i), s%ind_extr(:, i) - llc + 1, &
outside_points, 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, step)
else
s%list_vis(i)%speed_cont = null_ssh_contour()
s%list_vis(i)%max_speed = missing_speed
s%list_vis(i)%radius4 = 0
end if
end do
! Write snapshot:
inquire(file = "SHPC_cyclo/extremum.shp", exist = exist)
if (exist) then
call shpc_open(hshpc_cyclo, shpc_dir = "SHPC_cyclo", rank = 0, &
pszaccess = "rb+")
call shpc_open(hshpc_anti, shpc_dir = "SHPC_anti", rank = 0, &
pszaccess = "rb+")
else
call shpc_create(hshpc_cyclo, shpc_dir = "SHPC_cyclo", cyclone = .true.)
call shpc_create(hshpc_anti, shpc_dir = "SHPC_anti", cyclone = .false.)
end if
n_cyclo = 0
n_anti = 0
do i = 1, s%number_vis_extr
if (s%list_vis(i)%cyclone) then
n_cyclo = n_cyclo + 1
call write_eddy(s%list_vis(i), hshpc_cyclo, days_1950, n_cyclo)
else
n_anti = n_anti + 1
call write_eddy(s%list_vis(i), hshpc_anti, days_1950, n_anti)
end if
end do
call write_aux(corner_deg, step_deg, nlon, nlat, n_cyclo, &
shpc_dir = "SHPC_cyclo")
call write_aux(corner_deg, step_deg, nlon, nlat, n_anti, &
shpc_dir = "SHPC_anti") CALL shpc_close(hshpc_cyclo)
CALL shpc_close(hshpc_anti)
print *, 'inst_eddies: Created shapefile collections in SHPC_cyclo and ', &
'SHPC_anti.'
end program inst_eddies