mean_speed.f90

module mean_speed_m

  implicit none

contains

  real function mean_speed(p, points_xy, center)

    ! Interpolates the wind at each point of the input polyline and
    ! computes the mean azimuthal speed at interpolation points. A
    ! positive value of `mean_speed` means counterclockwise
    ! circulation. Interpolation is in projection space.

    ! We assume (and do not check) that the longitude of the center
    ! has been brought modulo 2 pi in the longitude range of the
    ! polyline `p` and that the polyline is in the longitude range of
    ! the global grid corresponding to u and v.

    use, intrinsic:: iso_fortran_env, only: error_unit

    ! Libraries:
    use contour_531, only: polyline
    use numer_rec_95, only: bilinear_interp2_reg, bilinear_interp_reg
    use jumble, only: rad_to_deg

    use input_ssh_m, only: corner_u, corner_v, single_grid_uv, u, v

    type(polyline), intent(in):: p
    ! Contour, with points given by longitude and latitude, in
    ! radians. Should be closed.

    real, intent(in):: points_xy(:, :) ! (2, p%n_points)
    ! Points of the contour `p` given by projection coordinates.

    real, intent(in):: center(:) ! (2)
    ! Longitude and latitude, in rad. Azimuthal speed is computed with
    ! respect to this point.

    ! Local:

    real ui(p%n_points - 1), vi(p%n_points - 1) ! (ni)
    ! speed interpolated at polygon points

    real x, y ! coordinates in tangent plane, divided by Earth radius
    real v_azim(p%n_points - 1) ! (ni) azimuthal speed at polygon points
    integer j
    integer ni ! number of interpolation points

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

    if (p%n_points <= 1) then
       write(error_unit, fmt = *) &
            "mean_speed: requires non-null, closed contour"
       write(error_unit, fmt = *) "center (in °) = ", center * rad_to_deg
       write(error_unit, fmt = *) "p%n_points = ", p%n_points
       write(error_unit, fmt = *) "Aborting."
       stop 1
    end if

    ni = p%n_points - 1

    if (single_grid_uv) then
       call bilinear_interp2_reg(corner_u, [1., 1.], u, v, points_xy(:, :ni), &
            ui, vi)
    else
       ui = bilinear_interp_reg(corner_u, [1., 1.], u, points_xy(:, :ni))
       vi = bilinear_interp_reg(corner_v, [1., 1.], v, points_xy(:, :ni))
    end if

    do j = 1, ni
       x = cos(center(2)) * (p%points(1, j) - center(1))
       y = p%points(2, j) - center(2)
       v_azim(j) = (x * vi(j) - y * ui(j)) / sqrt(x**2 + y**2)
    end do

    mean_speed = sum(v_azim) / ni

  end function mean_speed

end module mean_speed_m