-
Lionel GUEZ authored
Write longitude and latitude if `leave_trace` and not `uniform_lon_lat`.
Lionel GUEZ authoredWrite longitude and latitude if `leave_trace` and not `uniform_lon_lat`.
input_ssh.f90 9.04 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
real, protected:: corner_u(2) = huge(0.), corner_v(2) = huge(0.)
! Positions in projection space x, y of points corresponding to u(1
! - copy, 1) and v(1 - copy, 1).
logical, protected:: single_grid_uv = .true. ! same grid for u and v
real, allocatable, save, protected:: u(:, :), v(:, :)
! (1 - max_radius(1):nlon + max_radius(1), nlat) if periodic, else
! (nlon, nlat). Wind, in m s-1.
private get_var
contains
subroutine input_ssh(periodic, ssh, ssh_fname, u_fname, v_fname, leave_trace)
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
character(len = *), intent(in):: ssh_fname, u_fname, v_fname ! file names
logical, intent(in):: leave_trace
! Local:
integer unit
integer ncid, dimid, varid
integer:: nc_time = 1
! Index, 1-based, in the NetCDF time coordinate, if any, of the
! time slice to read. If there is no time coordinate, then leave
! nc_time to its default 1 value. If there is a time coordinate
! then we assume that the dimensions of the ssh, u and v NetCDF
! variables are longitude, latitude and time, in Fortran order.
real:: corner_deg(2) = huge(0.)
! longitude and latitude of the corner of the whole grid, in degrees
real:: stag_u(2) = 0., stag_v(2) = 0.
! Staggering for u and v. That is, position in projection space x, y of
! point corresponding to u(1, 1) or v(1, 1) respective to point
! corresponding to ssh(1, 1).
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, stag_u, stag_v
namelist /input_ssh_nml/ periodic, uniform_lon_lat, nc_time, stag_u, stag_v
!----------------------------------------------------------------------
periodic = .false. ! default value
print *, "input_ssh: Enter input_ssh_nml:"
read(unit = *, nml = input_ssh_nml)
single_grid_uv = all(stag_u == stag_v)
call nf95_open(ssh_fname, nf95_nowrite, ncid)
if (leave_trace) call new_unit(unit)
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)
call nf95_find_coord(ncid, varid = varid, std_name = "latitude")
call nf95_gw_var(ncid, varid, latitude)
if (leave_trace) then
open(unit, file = "longitude.txt", status = "replace", &
action = "write")
write(unit, fmt = *) longitude
close(unit)
open(unit, file = "latitude.txt", status = "replace", &
action = "write")
write(unit, fmt = *) latitude
close(unit)
end if
longitude = longitude * deg_to_rad
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
if (leave_trace) then
open(unit, file = "grid_nml.txt", status = "replace", action = "write")
write(unit, nml = grid_nml)
close(unit)
end if
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))
corner_u = [1. - copy + stag_u(1), 1. + stag_u(2)]
if (.not. single_grid_uv) corner_v = [1. - copy + stag_v(1), 1. + stag_v(2)]
! Read ssh, u and v:
call get_var(periodic, ssh, ncid, nc_time, 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, nc_time, 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, nc_time, 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, nc_time, 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
integer, intent(in):: nc_time ! Index, 1-based, in the NetCDF time coordinate, if any, of the
! time slice to read. If there is no time coordinate, then nc_time
! should be equal to 1. If there is a time coordinate then we
! assume that the dimension of the NetCDF variable to read are
! longitude, latitude and time, in Fortran order.
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)
if (nc_time == 1) then
call nf95_get_var(ncid, varid, values(1:nlon, :))
else
call nf95_get_var(ncid, varid, values(1:nlon, :), &
start = [1, 1, nc_time])
end if
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