module set_all_outerm_m

  implicit none

contains

  subroutine set_all_outerm(s, min_amp, max_radius, corner, step, ssh)

    ! Extraction of eddies: find extrema and set all outermost
    ! contours in snapshot. Not a function because snapshot is not
    ! completely defined.

    use derived_types, only: snapshot
    use get_1_outerm_m, only: get_1_outerm
    use local_extrema_m, only: local_extrema

    type(snapshot), intent(out):: s
    ! Specifically: define everything in s except
    ! s%list_vis%max_speed_contour, s%list_vis%max_speed and
    ! s%number_eddies.

    real, intent(in):: min_amp ! minimum amplitude of ssh, between
    ! extremum and outermost contour, in m

    integer, intent(in):: max_radius(:) ! (2) maximum radius of an
    ! eddy in longitude and latitude, in number of grid points

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

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

    real, intent(in):: ssh(:, :) ! sea-surface height, in m

    ! Local:

    real, allocatable:: innermost_level(:) ! (s%number_vis_eddies)
    ! level of innermost contour, for each extremum

    logical, allocatable:: cyclone(:) ! (s%number_vis_eddies)
    integer i
    logical, allocatable:: flat_extr(:) ! (s%number_vis_eddies) flat extremum
    logical, allocatable:: noise_around(:) ! (s%number_vis_eddies)

    ! Window around each extremum:

    integer, allocatable:: llc(:, :) ! (2, s%number_vis_eddies)
    ! indices in global grid of lower left corner

    integer, allocatable:: urc(:, :) ! (2, s%number_vis_eddies)
    ! indices in global grid of upper right corner

    real, allocatable:: corner_window(:, :) ! (2, s%number_vis_eddies)
    ! longitude and latitude, in rad

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

    allocate(s%extr_map(size(ssh, 1), size(ssh, 2)))
    call local_extrema(ssh, s%extr_map, s%ind_extr, innermost_level, cyclone)
    s%number_vis_eddies = size(s%ind_extr, 2)
    allocate(s%list_vis(s%number_vis_eddies), flat_extr(s%number_vis_eddies), &
         noise_around(s%number_vis_eddies), llc(2, s%number_vis_eddies), &
         urc(2, s%number_vis_eddies), corner_window(2, s%number_vis_eddies))

    forall (i = 1:s%number_vis_eddies)
       s%list_vis(i)%coord_extr = corner + (s%ind_extr(:, i) - 1) * step
       s%list_vis(i)%ssh_extremum = ssh(s%ind_extr(1, i), s%ind_extr(2, i))
       s%list_vis(i)%cyclone = cyclone(i)
       s%list_vis(i)%interpolated = .false.
    end forall

    ! Define the geographical window around each eddy extremum:
    forall (i = 1:s%number_vis_eddies)
       llc(:, i) = max(s%ind_extr(:, i) - max_radius, 1)
       urc(:, i) = min(s%ind_extr(:, i) + max_radius, &
            [size(ssh, 1), size(ssh, 2)])
       corner_window(:, i) = corner + (llc(:, i) - 1) * step
    end forall

    if (min_amp == 0.) then
       flat_extr = .false.
    else
       flat_extr = abs(innermost_level - s%list_vis%ssh_extremum) < min_amp
    end if

    do i = 1, s%number_vis_eddies
       if (flat_extr(i)) then
          s%list_vis(i)%outermost_contour &
               = get_1_outerm(s%list_vis(i)%ssh_extremum, &
               s%list_vis(i)%cyclone, s%list_vis(i)%coord_extr, i, &
               innermost_level(i), &
               s%extr_map(llc(1, i):urc(1, i), llc(2, i):urc(2, i)), &
               ssh(llc(1, i):urc(1, i), llc(2, i):urc(2, i)), &
               corner_window(:, i), step)
          if (s%list_vis(i)%outermost_contour%n_points == 0) then
             s%list_vis(i)%suff_amp = .false.
          else
             s%list_vis(i)%suff_amp &
                  = abs(s%list_vis(i)%outermost_contour%ssh &
                  - s%list_vis(i)%ssh_extremum) >= min_amp
          end if
       end if
    end do

    ! We must modify s%extr_map in this loop, separate from the
    ! previous loop, so we do not influence the batch of
    ! get_1_outerm.
    do i = 1, s%number_vis_eddies
       if (flat_extr(i) .and. .not. s%list_vis(i)%suff_amp) &
            s%extr_map(s%ind_extr(1, i), s%ind_extr(2, i)) &
            = - s%extr_map(s%ind_extr(1, i), s%ind_extr(2, i))
    end do
    ! (no forall because of gfortran internal compiler error)

    ! Define noise_around:
    if (min_amp /= 0.) &
         forall &
         (i = 1:s%number_vis_eddies, &
         flat_extr(i) .and. s%list_vis(i)%suff_amp) &
         noise_around(i) &
         = any(s%extr_map(llc(1, i):urc(1, i), llc(2, i):urc(2, i)) < 0)

    s%list_vis%twice = flat_extr .and. s%list_vis%suff_amp .and. noise_around

    do i = 1, s%number_vis_eddies
       if (s%list_vis(i)%suff_amp .and. noise_around(i) &
            .or. .not. flat_extr(i)) then
          s%list_vis(i)%outermost_contour &
               = get_1_outerm(s%list_vis(i)%ssh_extremum, &
               s%list_vis(i)%cyclone, s%list_vis(i)%coord_extr, i, &
               innermost_level(i), &
               s%extr_map(llc(1, i):urc(1, i), llc(2, i):urc(2, i)), &
               ssh(llc(1, i):urc(1, i), llc(2, i):urc(2, i)), &
               corner_window(:, i), step)
          s%list_vis(i)%suff_amp &
               = s%list_vis(i)%outermost_contour%n_points /= 0
       end if
    end do

    ! We must modify s%extr_map in this loop, separate from the
    ! previous loop, so we do not influence the batch of
    ! get_1_outerm.
    do i = 1, s%number_vis_eddies
       if (.not. flat_extr(i) .and. .not. s%list_vis(i)%suff_amp) &
            s%extr_map(s%ind_extr(1, i), s%ind_extr(2, i)) &
            = - s%extr_map(s%ind_extr(1, i), s%ind_extr(2, i))
    end do
    ! (no forall because of gfortran internal compiler error)

  end subroutine set_all_outerm

end module set_all_outerm_m