set_contours.f90

module set_contours_m

  implicit none

contains

  subroutine set_contours(out_cont, speed_cont, max_speed, cyclone, extr, &
       innermost_level, step, ssh, u, v, corner, outside_points)

    ! This procedure sets contours in an eddy.

    ! Libraries:
    use contour_531, only: polyline

    use config_m, only: min_amp
    use derived_types, only: extremum, ssh_contour
    use get_1_outerm_m, only: get_1_outerm
    use set_max_speed_m, only: set_max_speed

    type(ssh_contour), intent(out):: out_cont ! outermost contour

    type(ssh_contour), intent(out):: speed_cont
    ! contour with maximum average azimuthal speed

    real, intent(out):: max_speed
    ! Average of azimuthal speed on out_cont 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. Not defined if out_cont is not closed.

    logical, intent(in):: cyclone
    type(extremum), intent(in):: extr

    real, intent(in):: innermost_level
    ! SSH level of the innermost contour around each extremum, in
    ! m. By construction, innermost_level < extremum for a maximum, >
    ! extremum for a minimum.

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

    real, intent(in), dimension(:, :):: ssh, u, v
    ! sea-surface height, in m, and wind, in m s-1

    integer, intent(in):: corner(:) ! (2)
    ! Coordinates in projection space of the corner of the grid,
    ! corresponding to ssh(1, 1). A priori, this is not the corner of
    ! the grid of the NetCDF file.

    integer, intent(in):: outside_points(:, :) ! (2, :)
    ! coordinates in projection space of all the significant extrema,
    ! except the target extremum

    ! Local:

    real innermost_level_2

    ! Window around each extremum:

    integer llc(2) ! indices in grid of lower left corner
    integer urc(2) ! indices in grid of upper right corner

    integer corner_window(2)
    ! coordinates in projection space of corner of the window around
    ! each extremum

    integer, parameter:: n_max_cont = 31 ! must be >= 3

    type(polyline) cont_list(n_max_cont)
    ! Contour list, with points given by longitude and latitude, in
    ! radians. Defined only for subscripts 1:n_cont. The order of
    ! points in each contour (clockwise or counter-clockwise) is not
    ! specified.

    type(polyline) cont_list_proj(n_max_cont)
    ! Contour list, with points given by projection
    ! coordinates. Defined only for subscripts 1:n_cont. The contour
    ! with a given subscript in cont_list_proj corresponds to the
    ! contour with the same subscript in cont_list. More precisely,
    ! the point with subscript j in contour i of cont_list_proj
    ! corresponds to the point with subscript j in contour i of
    ! cont_list (this is important for mean_speed). For i <= n_cont,
    ! cont_list_proj(i)%area has its default initialization value and
    ! the order of points in cont_list_proj(i) (clockwise or
    ! counter-clockwise) is not specified.

    real cont_list_ssh(n_max_cont)

    integer n_cont
    ! number of good contours found and stored in cont_list, 0 <=
    ! n_cont <= n_max_cont

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

    ! Only look for good contours with amplitudes >= min_amp:
    innermost_level_2 = merge(extr%ssh + merge(min_amp, - min_amp, cyclone), &
         innermost_level, abs(extr%ssh - innermost_level) < min_amp)

    call get_1_outerm(out_cont, cont_list, cont_list_proj, cont_list_ssh, &
         n_cont, cyclone, real(extr%coord_proj), innermost_level_2, &
         real(outside_points), ssh, real(corner), step)

    ! Done with outermost contour, now let us take care of the
    ! max-speed contour.
    if (out_cont%closed) then
       ! assert n_cont >= 1
       ! {begin} Restrict the field to the outermost contour:

       corner_window = floor(minval(cont_list_proj(n_cont)%points, dim = 2))
       llc = corner_window - corner + 1
       urc = ceiling(maxval(cont_list_proj(n_cont)%points, dim = 2)) - corner &
            + 1

       ! Should have no effect except because of roundup error:
       urc(2) = min(urc(2), size(ssh, 2))
       urc(1) = min(urc(1), size(ssh, 1))

       ! {end} Done restricting field

       call set_max_speed(speed_cont, max_speed, cont_list, cont_list_proj, &
            cont_list_ssh, n_cont, extr%coord, real(extr%coord_proj), &
            ssh(llc(1):urc(1), llc(2):urc(2)), &
            u(llc(1):urc(1), llc(2):urc(2)), v(llc(1):urc(1), llc(2):urc(2)), &
            real(corner_window), step)
    end if

  end subroutine set_contours

end module set_contours_m