-
Lionel GUEZ authoredLionel GUEZ authored
inst_eddies.f90 7.51 KiB
program inst_eddies
use, intrinsic:: ieee_arithmetic, only: IEEE_SUPPORT_DATATYPE, &
ieee_support_nan, ieee_value, IEEE_QUIET_NAN
use, intrinsic:: ISO_FORTRAN_ENV
! Libraries:
use contour_531, only: convert_to_ind
use jumble, only: new_unit
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 derived_types, only: snapshot, null_ssh_contour, missing_speed, shp_tr
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 shp_tr_close_m, only: shp_tr_close
use shp_tr_create_m, only: shp_tr_create
use write_eddy_m, only: write_eddy
implicit none
type(snapshot) s
TYPE(shp_tr) hshp
integer i, k, unit
real:: min_amp = 0.001
! minimum amplitude of ssh, between extremum and outermost contour,
! in m
real:: min_radius = 25.
! minimum radius of the equal-area disk for an outermost contour, in km
integer ncid, dimid, varid, copy, iostat
integer nlon, nlat
! size of ssh array in input NetCDF, assuming no repeated point if
! the grid is global
real lon_min, lon_max, lat_min, lat_max ! longitudes and latitudes, in degrees
real step(2) ! longitude and latitude steps, in rad
logical periodic ! grid is periodic in longitude
real time
integer:: max_radius(2) = [20, 20]
! 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(2)
! longitude and latitude of the corner of the whole grid, 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
character(len = 200) iomsg
namelist /main_nml/ min_amp, max_radius, min_radius
!--------------------------------------------------------------
call assert(IEEE_SUPPORT_DATATYPE(0.), ieee_support_nan(0.), &
"inst_eddies: not enough IEEE support")
! Trick to verify that the directory "SHP_triplet" exists:
call new_unit(unit)
open(unit, file = "SHP_triplet/extremum.dbf", status = "replace", &
action = "write", iostat = iostat, iomsg = iomsg)
if (iostat /= 0) then
write(unit = error_unit, fmt = *) trim(iomsg)
write(unit = error_unit, fmt = *) &
'Have you created the directory "SHP_triplet"?'
stop 1
end if
close(unit)
! We need this before looking at the input files:
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 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)
k = nint(time)
call assert(abs(time - k) < 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, lon_min)
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, lat_min)
call nf95_get_var(ncid, varid, lat_max, start = [nlat])
call assert(lon_max > lon_min .and. lat_max > lat_min, &
"inst_eddies coordinate order")
step = [(lon_max - lon_min) / (nlon - 1), (lat_max - lat_min) / (nlat - 1)] &
* 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 *, "periodic = ", periodic
call assert(2 * max_radius(1) * step(1) < pi, "inst_eddies max_radius")
! (Let us require this even if not periodic. This is clearer.)
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 = [lon_min, lat_min] * 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, min_amp, 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
call shp_tr_create(hshp, shp_tr_dir = "SHP_triplet")
! Write snapshot:
do i = 1, s%number_vis_extr
call write_eddy(s%list_vis(i), hshp, k, i)
end do
print *, "Number of extrema:", s%number_vis_extr
CALL shp_tr_close(hshp)
print *, 'Created shapefiles in SHP_triplet.'
end program inst_eddies