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

    ! Libraries:
    use jumble, only: rad_to_deg
    use shapelib, only: shpt_polygon
    use shapelib_03, only: shp_append_object_03, dbf_write_attribute_03
    use contour_531, only: convert_to_reg_coord

    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 input_ssh_m, only: corner_whole => corner
    use mean_speed_m, only: mean_speed
    use spher_polyline_area_m, only: spher_polyline_area
    use cont_list_m, only: write_cont_list, hshp_cont_list, ifield_ssh, &
         ifield_speed

    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. cont_list%area has its default initialization value and the
    ! order of points in the polyline components (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,
    ! corresponding to ssh(1, 1), in rad. A priori, this is not the
    ! corner of the grid of the NetCDF file.

    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, ishape
    real corner_proj(2), extr_coord_proj(2)

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

    corner_proj = (corner - corner_whole) / step + 1.
    extr_coord_proj = (extr_coord - corner_whole) / step + 1.
    i_outer = n_cont
    if (n_cont >= 2) call complete_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 = convert_to_reg_coord(good_contour(corner_proj, &
            ssh, cont_list(i)%ssh, extr_coord_proj), corner_whole, step)
    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

    if (write_cont_list) then
       do i = 1, n_cont
          call shp_append_object_03(ishape, hshp_cont_list, shpt_polygon, &
               cont_list(i)%points * rad_to_deg)
          call dbf_write_attribute_03(hshp_cont_list, ishape, ifield_ssh, &
               fieldvalue = cont_list(i)%ssh)

          if (ieee_is_nan(speed(i))) then
             call dbf_write_attribute_03(hshp_cont_list, ishape, ifield_speed, &
                  missing_speed)
             ! (Cannot write NaN to dbf file.)
          else
             call dbf_write_attribute_03(hshp_cont_list, ishape, ifield_speed, &
                  speed(i))
          end if
       end do
    end if

    i = maxloc(abs(speed), dim = 1, mask = .not. IEEE_IS_NAN(speed))
    ! (The speed may also 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