input_ssh.f90

module input_ssh_m

  implicit none

  integer, protected:: max_radius(2)
  ! maximum radius of an eddy in longitude and latitude, in number of
  ! grid points

contains

  subroutine input_ssh(corner, step, periodic, ssh, u, v)

    use, intrinsic:: ieee_arithmetic, only: ieee_value, IEEE_QUIET_NAN

    ! Libraries:
    use jumble, only: new_unit, assert, deg_to_rad, twopi, get_command_arg_dyn
    use netcdf95, only: nf95_open, nf95_find_coord, nf95_inquire_dimension, &
         nf95_get_var, nf95_close, nf95_nowrite

    use config_m, only: max_radius_deg
    use get_var_m, only: get_var

    real, intent(out):: corner(:) ! (2)
    ! longitude and latitude of the corner of the whole grid, in rad

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

    logical, intent(out):: periodic ! grid is periodic in longitude

    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

    ! Local:

    integer nlon, nlat
    ! size of ssh array in input NetCDF, assuming no repeated point if
    ! the grid is global

    logical exist
    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
    integer nlon_new, nlat_new
    real corner_deg_new(2), step_deg_new(2)
    real lon_max, lat_max ! longitude and latitude, in degrees
    integer copy
    character(len = :), allocatable:: fname ! file name
    namelist /grid_nml/ corner_deg, step_deg, nlon, nlat

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

    inquire(file = "SHPC/grid_nml.txt", exist = exist)
    call get_command_arg_dyn(1, fname, "Required arguments: h-file uv-file")
    call nf95_open(fname, nf95_nowrite, ncid)
    call nf95_find_coord(ncid, dimid = dimid, varid = varid, &
         std_name = "longitude")
    call new_unit(unit)

    if (exist) then
       print *, "input_ssh: Found SHPC/grid_nml.txt"
       nlon = - huge(0)
       nlat = - huge(0)
       open(unit, file = "SHPC/grid_nml.txt", status = "old", action = "read", &
            position = "rewind")
       read(unit, nml = grid_nml)
       close(unit)
       call nf95_inquire_dimension(ncid, dimid, nclen = nlon_new)
       call assert(nlon_new == nlon, "input_ssh: nlon")
       call nf95_get_var(ncid, varid, corner_deg_new(1))
       call assert(corner_deg_new(1) == corner_deg(1), &
            "input_ssh: corner_deg(1)")
    else
       print *, "input_ssh: SHPC/grid_nml.txt not found. Will be created..."
       call nf95_inquire_dimension(ncid, dimid, nclen = nlon)
       call nf95_get_var(ncid, varid, corner_deg(1))
    end if

    call nf95_get_var(ncid, varid, lon_max, start = [nlon])
    call nf95_find_coord(ncid, dimid = dimid, varid = varid, &
         std_name = "latitude")

    if (exist) then
       call nf95_inquire_dimension(ncid, dimid, nclen = nlat_new)
       call assert(nlat_new == nlat, "input_ssh: nlat")
       call nf95_get_var(ncid, varid, corner_deg_new(2))
       call assert(corner_deg_new(2) == corner_deg(2), &
            "input_ssh: corner_deg(2)")
    else
       call nf95_inquire_dimension(ncid, dimid, nclen = nlat)
       call nf95_get_var(ncid, varid, corner_deg(2))
    end if

    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")

    if (exist) then
       step_deg_new = [(lon_max - corner_deg(1)) / (nlon - 1), &
            (lat_max - corner_deg(2)) / (nlat - 1)]
       call assert(all(step_deg_new == step_deg), "input_ssh: step_deg")
    else
       step_deg = [(lon_max - corner_deg(1)) / (nlon - 1), &
            (lat_max - corner_deg(2)) / (nlat - 1)]
       open(unit, file = "SHPC/grid_nml.txt", status = "replace", &
            action = "write")
       write(unit, nml = grid_nml)
       close(unit)
    end if

    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 *, "input_ssh: periodic = ", periodic

    call assert(2 * max_radius_deg(1) < 180., "input_ssh: 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 get_command_arg_dyn(2, fname, "Required arguments: h-file uv-file")
    call nf95_open(fname, nf95_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)

  end subroutine input_ssh

end module input_ssh_m