set_max_speed.F90

module set_max_speed_m

  implicit none

contains

  subroutine set_max_speed(speed_cont, max_speed, cont_list, n_cont, &
       extr_coord, ssh, u, v, corner, step)

    ! This procedure defines speed_cont and max_speed. On return,
    ! speed_cont may be a null ssh contour and max_speed may be set to
    ! missing_speed. max_speed is never NaN.

    ! The length of the interval of longitudes of the domain, step(1)
    ! * (size(ssh, 1) - 1), should be < 180 degrees, so that the
    ! geometrical processing done here is non-ambiguous.

    use, intrinsic:: IEEE_ARITHMETIC, only: IEEE_IS_NAN, ieee_value, &
         IEEE_QUIET_NAN

    use ccw_orient_m, only: ccw_orient
    use complete_ssh_m, only: complete_ssh
    use derived_types, only: eddy, null_ssh_contour, ssh_contour, missing_speed
    use good_contour_m, only: good_contour
    use mean_speed_m, only: mean_speed
    use spher_polyline_area_m, only: spher_polyline_area

    type(ssh_contour), intent(out):: speed_cont

    real, intent(out):: max_speed
     ! Average of azimuthal speed on cont_list(i_outer) or speed_cont,
     ! in m s-1. Positive means counterclockwise circulation. If
     ! speed_cont is not a null ssh contour then max_speed is the
     ! speed on speed_cont.

    type(ssh_contour), intent(inout):: cont_list(:) ! (n_max_cont)
    ! Contour list. Defined only for subscripts 1:n_cont. On entry,
    ! the outermost contour is element with subscript n_cont. The
    ! contours between 1 and i_outer are in monotonic order of
    ! ssh. For i /= i_outer, cont_list(i)%area has its default
    ! initialization value and the order of points in
    ! cont_list(i)%polyline (clockwise or counter-clockwise) is not
    ! specified.

    integer, intent(inout):: n_cont
    ! Number of good contours found and stored in cont_list. On entry,
    ! 1 <= n_cont <= n_max_cont - 1 (not 0 because set_max_speed is
    ! only called when there is an outermost contour). On return, 1 <=
    ! n_cont <= n_max_cont.

    real, intent(in):: extr_coord(:)
    ! (2) longitude and latitude of extremum , in rad

    real, intent(in):: ssh(:, :), u(:, :), v(:, :)
    ! The domain should be the bounding box of out_cont, because we do
    ! not exclude other extrema.

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

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

    ! Local:

    real, allocatable:: speed(:) ! (n_cont) speed on the contour
    integer i, i_outer

    integer sort_ssh(size(cont_list)) ! (n_max_cont)
    ! Sorts cont_list%ssh in ascending order. Defined only for
    ! subscripts 1:n_cont.

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

    i_outer = n_cont
    if (n_cont >= 2) call complete_ssh(sort_ssh, cont_list%ssh, n_cont)

    ! Now find the contours associated to the new values of SSH:
    do i = i_outer + 1, n_cont
       cont_list(i)%polyline = good_contour(corner, step, ssh, &
            cont_list(i)%ssh, extr_coord)
    end do

    allocate(speed(n_cont))

    do i = 1, n_cont
       if (cont_list(i)%closed) then
          speed(i) = mean_speed(u, v, cont_list(i)%polyline, extr_coord, &
               corner, step)
       else
          ! It is possible that good_contour returned a null polyline
          ! if there is a missing value of SSH inside the outermost
          ! contour.
          speed(i) = ieee_value(0., IEEE_QUIET_NAN)
       end if
    end do

#ifdef DEBUG
    print *, "set_max_speed: speed(sort_ssh) = ", speed(sort_ssh(:n_cont))
    print *, "set_max_speed: cont_list%ssh(sort_ssh) = ", &
         cont_list(sort_ssh(:n_cont))%ssh
#endif

    i = maxloc(abs(speed), dim = 1, mask = .not. IEEE_IS_NAN(speed))
    ! (The speed may be NaN when we are near land.)

    if (i == 0) then
       ! All speed values are NaN
       speed_cont = null_ssh_contour()
       max_speed = missing_speed
    else
       max_speed = speed(i)

       if (i == i_outer) then
          ! Maximum speed on the outermost contour
          speed_cont = null_ssh_contour()
       else
          speed_cont = cont_list(i)
          speed_cont%area = spher_polyline_area(speed_cont%polyline)
          call ccw_orient(speed_cont)
       end if
    end if

  end subroutine set_max_speed

end module set_max_speed_m