max_speed_contour_ssh.f90

module max_speed_contour_ssh_m

  implicit none

contains

  pure real function max_speed_contour_ssh(ssh, u, v, ind_extr, radius)

    ! Choose an SSH level by examining the values of velocities inside
    ! radius.

    use, intrinsic:: IEEE_ARITHMETIC, only: IEEE_IS_NAN

    use derived_types, only: missing_ssh

    real, intent(in):: ssh(:, :), u(:, :), v(:, :)
    integer, intent(in):: ind_extr(:) ! (2) indices of the extremum in ssh

    integer, intent(in):: radius
    ! ind_extr + radius - 1 in all four directions is inside outermost
    ! contour, ind_extr + radius may be outside. radius should be >=
    ! 2.

    ! Local:

    real v_azim(4, radius - 1)
    ! azimuthal velocity
    ! First dimension: direction, index 1 for east, 2 for north, etc.
    ! Second dimension: grid index, starting one point away from the extremum

    integer l, direction
    integer, parameter:: east = 1, north = 2, west = 3, south = 4

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

    ! Assuming that ascending indices in u and v correspond to
    ! ascending longitude and latitude:
    v_azim(east, :) = v(ind_extr(1) + 1:ind_extr(1) + radius - 1, &
         ind_extr(2))
    v_azim(north, :) = - u(ind_extr(1), ind_extr(2) + 1:ind_extr(2) &
         + radius - 1)
    v_azim(west, :) = - v(ind_extr(1) - 1:ind_extr(1) - radius + 1:- 1, &
         ind_extr(2))
    v_azim(south, :) = u(ind_extr(1), ind_extr(2) - 1:ind_extr(2) - radius &
         + 1:- 1)

    if (any(ieee_is_nan(v_azim))) then
       ! This is strange, we are inside an outermost contour, and yet
       ! there is a missing value. Well, we will go on with the
       ! outermost contour as max-speed contour.
       max_speed_contour_ssh = missing_ssh
    else
       if (radius >= 3) then
          l = maxloc(abs(sum(v_azim, dim = 1) / 4.), dim = 1)
       else
          l = 1
       end if
       
       direction = maxloc(abs(v_azim(:, l)), dim = 1)

       select case (direction)
       case(east)
          max_speed_contour_ssh = ssh(ind_extr(1) + l, ind_extr(2))
       case(north)
          max_speed_contour_ssh = ssh(ind_extr(1), ind_extr(2) + l)
       case(west)
          max_speed_contour_ssh = ssh(ind_extr(1) - l, ind_extr(2))
       case(south)
          max_speed_contour_ssh = ssh(ind_extr(1), ind_extr(2) - l)
       end select
    end if

  end function max_speed_contour_ssh
  
end module max_speed_contour_ssh_m