input_ssh.f90

module input_ssh_m

  implicit none

  real, protected:: corner(2) = huge(0.)
  ! longitude and latitude of the corner of the whole grid, 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

  private get_var

contains

  subroutine input_ssh(step, 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
    use netcdf95, only: nf95_open, nf95_find_coord, nf95_inquire_dimension, &
         nf95_get_var, nf95_close, nf95_nowrite

    use config_m, only: max_radius

    real, intent(out):: step(:) ! (2) longitude and latitude steps, in rad

    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

    !----------------------------------------------------------------------

    call nf95_open(ssh_fname, nf95_nowrite, ncid)
    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)]
    call new_unit(unit)
    open(unit, file = "grid_nml.txt", status = "replace", action = "write")
    write(unit, nml = grid_nml)
    close(unit)
    step = step_deg * deg_to_rad
    print *, "input_ssh: periodic = ?"
    read *, periodic
    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 = corner_deg * deg_to_rad

    ! 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