-
GUEZ Lionel authored
In the case of non-uniform longitude-latitude grid: we read the whole coordinate arrays from the NetCDF file; we have a different way of inverting the projection, in `invert_proj`; the program `eddy_graph` does not use `corner_whole` and step, which are not defined, it uses instead the shapefile `extr_proj`.
GUEZ Lionel authoredIn the case of non-uniform longitude-latitude grid: we read the whole coordinate arrays from the NetCDF file; we have a different way of inverting the projection, in `invert_proj`; the program `eddy_graph` does not use `corner_whole` and step, which are not defined, it uses instead the shapefile `extr_proj`.
input_ssh.f90 7.01 KiB
module input_ssh_m
implicit none
real, protected:: corner_whole(2) = huge(0.)
! longitude and latitude of the corner of the whole grid, in rad
real, protected:: step(2) = huge(0.) ! longitude and latitude steps, in rad
integer, protected:: nlon = - huge(0), nlat = - huge(0)
! size of ssh array in input NetCDF, assuming no repeated point if
! the grid is global
logical, protected:: uniform_lon_lat = .true.
! grid is cartesian and uniform in longitude and latitude
real, allocatable, save, protected:: longitude(:) ! (nlon)
! allocated only for non-uniform grid, coordinate for the whole grid, in rad
real, allocatable, save, protected:: latitude(:) ! (nlat)
! allocated only for non-uniform grid, coordinate for the whole grid, in rad
real, allocatable, save, protected:: lon_steps(:) ! (nlon - 1)
! allocated only for non-uniform grid, longitude steps, in rad
real, allocatable, save, protected:: lat_steps(:) ! (nlat - 1)
! allocated only for non-uniform grid, longitude steps, in rad
private get_var
contains
subroutine input_ssh(periodic, ssh, u, v, ssh_fname, u_fname, v_fname)
use, intrinsic:: ieee_arithmetic, only: ieee_value, IEEE_QUIET_NAN
! Libraries:
use jumble, only: new_unit, assert, deg_to_rad, ediff1d
use netcdf95, only: nf95_open, nf95_find_coord, nf95_inquire_dimension, &
nf95_get_var, nf95_close, nf95_nowrite, nf95_gw_var
use config_m, only: max_radius
logical, intent(out):: periodic ! grid is periodic in longitude
! If periodic, we assume that in the NetCDF file, the last
! longitude does not repeat the first longitude modulo 360
! degrees.
real, allocatable, intent(out):: ssh(:, :)
! (1 - max_radius(1):nlon + max_radius(1), nlat) if periodic, else
! (nlon, nlat) sea-surface height, in m
real, allocatable, intent(out):: u(:, :), v(:, :)
! (1 - max_radius(1):nlon + max_radius(1), nlat) if periodic, else
! (nlon, nlat) wind, in m s-1
character(len = *), intent(in):: ssh_fname, u_fname, v_fname ! file names
! Local:
integer unit
integer ncid, dimid, varid
real:: corner_deg(2) = huge(0.)
! longitude and latitude of the corner of the whole grid, in degrees
real:: step_deg(2) = huge(0.) ! longitude and latitude steps, in degrees
real lon_max, lat_max ! longitude and latitude, in degrees
integer copy
namelist /grid_nml/ corner_deg, step_deg, nlon, nlat, periodic, & uniform_lon_lat
namelist /input_ssh_nml/ periodic, uniform_lon_lat
!----------------------------------------------------------------------
periodic = .false. ! default value
print *, "input_ssh: Enter input_ssh_nml:"
read(unit = *, nml = input_ssh_nml)
call nf95_open(ssh_fname, nf95_nowrite, ncid)
if (uniform_lon_lat) then
call nf95_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 nf95_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), &
"input_ssh: coordinate order")
step_deg = [(lon_max - corner_deg(1)) / (nlon - 1), &
(lat_max - corner_deg(2)) / (nlat - 1)]
step = step_deg * deg_to_rad
corner_whole = corner_deg * deg_to_rad
else
call nf95_find_coord(ncid, varid = varid, std_name = "longitude")
call nf95_gw_var(ncid, varid, longitude)
longitude = longitude * deg_to_rad
call nf95_find_coord(ncid, varid = varid, std_name = "latitude")
call nf95_gw_var(ncid, varid, latitude)
latitude = latitude * deg_to_rad
call assert(longitude(2) > longitude(1) &
.and. latitude(2) > latitude(1), "input_ssh: coordinate order")
nlon = size(longitude)
nlat = size(latitude)
lon_steps = ediff1d(longitude)
lat_steps = ediff1d(latitude)
end if
call new_unit(unit)
open(unit, file = "grid_nml.txt", status = "replace", action = "write")
write(unit, nml = grid_nml)
close(unit)
if (periodic) call assert(4 * max_radius(1) < nlon, "input_ssh: max_radius")
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))
! Read ssh, u and v:
call get_var(periodic, ssh, ncid, 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.)
if (u_fname /= ssh_fname) then
call nf95_close(ncid)
call nf95_open(u_fname, nf95_nowrite, ncid)
end if
call get_var(periodic, u, ncid, name = "ugos", &
new_fill_value = ieee_value(0., IEEE_QUIET_NAN))
if (v_fname /= u_fname) then
call nf95_close(ncid)
call nf95_open(v_fname, nf95_nowrite, ncid)
end if
call get_var(periodic, v, ncid, 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)
end subroutine input_ssh
!*********************************************************************
subroutine get_var(periodic, values, ncid, name, new_fill_value)
! Read a NetCDF variable, change the missing value and extend it
! in longitude if periodic.
! Libraries:
use netcdf95, only: nf95_inq_varid, nf95_get_var, nf95_get_missing
use config_m, only: max_radius
logical, intent(in):: periodic ! grid is periodic in longitude
real, intent(out):: values(1 - merge(max_radius(1), 0, periodic):, :)
! (1 - merge(max_radius(1), 0, periodic):nlon + merge(max_radius(1),
! 0, periodic), nlat)
! ssh, u or v. We cannot place this argument first because the
! declaration references periodic.
integer, intent(in):: ncid
character(len = *), intent(in):: name ! of NetCDF variable
real, intent(in):: new_fill_value
! Local:
integer varid
real Fill_Value
!-------------------------------------------------------------------------
call nf95_inq_varid(ncid, name, varid)
call nf95_get_var(ncid, varid, values(1:nlon, :))
call nf95_get_missing(ncid, varid, Fill_Value)
! Change the missing value:
where (values(1:nlon, :) == Fill_Value) values(1:nlon, :) = new_fill_value
if (periodic) then
! Extend in longitude:
values(1 - max_radius(1):0, :) &
= values(nlon + 1 - max_radius(1):nlon, :)
values(nlon + 1:nlon + max_radius(1), :) = values(1:max_radius(1), :)
end if
end subroutine get_var
end module input_ssh_m