local_extrema.f90

module local_extrema_m

  implicit none

  private local_extremum

contains

  subroutine local_extrema(extr_map, ind_extr, innermost_level, local_min, &
       field, periodic, my_lbound)

    ! Finds local extrema in a field and closest value around each
    ! extremum. The extrema are output both in a map and as a
    ! list. The input field may be periodic in the first
    ! dimension. Then it is assumed that it has one duplicate column
    ! at each bound. If periodic then some extrema may be duplicated
    ! in extr_map but there is no duplication in ind_extr,
    ! innermost_level and local_min.

    ! Note that no coordinate grid is used here so there is no
    ! assumption on the grid underlying "field".

    use nr_util, only: assert_eq

    integer, intent(out):: extr_map(:, :) ! (m, n) Map of
    ! identification numbers of extrema. Identification numbers are
    ! subscripts in the second dimension of ind_extr. During a call to
    ! local_extrema, we use the array extr_map to temporarily record
    ! whether the extremum is a maximum or a minimum: extr_map is
    ! positive for a maximum, negative for a minimum, 0
    ! elsewhere. But, at the end of the call, all values of extr_map
    ! are set positive. If periodic then, on output, extr_map(1, :)
    ! will be equal to extr_map(m - 1, :) and extr_map(m, :) will be
    ! equal to extr_map(2, :).

    integer, intent(out), allocatable:: ind_extr(:, :) ! (2, n_extr)
    ! List of couple of subscripts of extrema of field, in the two
    ! dimensions of the field, considering that the lower bounds of
    ! field is given by my_lbound.

    real, intent(out), allocatable:: innermost_level(:) ! (n_extr)
    ! Level of innermost contour, for each extremum. By construction,
    ! innermost_level < extremum for a maximum, > extremum for a
    ! minimum.

    logical, intent(out), allocatable:: local_min(:) ! (n_extr)

    real, intent(in):: field(:, :) ! (m, n) Missing value is
    ! huge(0.). If periodic in the first dimension then the period
    ! should be m - 2, so field(1, :) should be equal to field(m - 1,
    ! :) and field(m, :) should be equal to field(2, :).

    logical, intent(in):: periodic ! field is periodic in the first dimension

    integer, intent(in):: my_lbound(:) ! (2)
    ! lower bounds of actual argument associated to field

    ! Local:

    integer n_extr, i, j, k, m, n

    real innermost_map(size(field, 1), size(field, 2)) ! (m, n) map of
    ! levels for innermost contours, defined only for extremum points

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

    m = assert_eq(size(field, 1), size(extr_map, 1), "local_extrema m")
    n = assert_eq(size(field, 2), size(extr_map, 2), "local_extrema n")

    extr_map = 0
    n_extr = 0

    ! No extremum on the borders of the second dimension of field so
    ! loop on columns 2 to n - 1.

    if (periodic) then
       do j = 2, n - 1
          ! Row 2:
          call local_extremum(innermost_map(2, j), extr_map(:3, j - 1:j + 1), &
               n_extr, field(:3, j - 1:j + 1))
          extr_map(m, j) = extr_map(2, j)

          do i = 3, m - 2
             call local_extremum(innermost_map(i, j), &
                  extr_map(i - 1:i + 1, j - 1:j + 1), n_extr, &
                  field(i - 1:i + 1, j - 1:j + 1))
          end do

          ! Row m - 1:
          call local_extremum(innermost_map(m - 1, j), &
               extr_map(m - 2:, j - 1:j + 1), n_extr, &
               field(m - 2:, j - 1:j + 1))
          extr_map(1, j) = extr_map(m - 1, j)
       end do
    else
       do j = 2, n - 1
          ! No extremum on the borders of the first dimension of field
          ! so loop on row 2 to n - 1:
          do i = 2, m - 1
             call local_extremum(innermost_map(i, j), &
                  extr_map(i - 1:i + 1, j - 1:j + 1), n_extr, &
                  field(i - 1:i + 1, j - 1:j + 1))
          end do
       end do
    end if

    ! Define ind_extr, considering that the lower bounds of field are
    ! 1 and 1. We will adjust it later for lower bounds my_lbound. We
    ! do not call pack_indices because we already know the number of
    ! extrema so the following algorithm is faster (stops at n_extr)
    ! and uses less memory.

    allocate(ind_extr(2, n_extr))

    i = 1
    j = 2

    do k = 1, n_extr
       do
          ! next (i, j)
          i = i + 1

          if (i == m) then
             ! Go to next column
             i = 2
             j = j + 1
          end if

          if (extr_map(i, j) /= 0) exit
       end do

       ind_extr(1, k) = i
       ind_extr(2, k) = j
       ! {abs(extr_map(ind_extr(1, k), ind_extr(2, k))) == k}
    end do

    ! Assertion: \forall k, ind_extr(1, k) \in \{2, \dots, m - 1\} et
    ! ind_extr(2, k) \in \{2, \dots, n - 1\}

    allocate(innermost_level(n_extr), local_min(n_extr))
    forall (k = 1:n_extr)
       innermost_level(k) = innermost_map(ind_extr(1, k), ind_extr(2, k))
       local_min(k) = extr_map(ind_extr(1, k), ind_extr(2, k)) < 0
    end forall

    ! Make all values of extr_map positive:
    forall (k = 1:n_extr, local_min(k)) &
         extr_map(ind_extr(1, k), ind_extr(2, k)) &
         = - extr_map(ind_extr(1, k), ind_extr(2, k))
    extr_map(1, :) = abs(extr_map(1, :))
    extr_map(m, :) = abs(extr_map(m, :))

    forall (k = 1:n_extr) ind_extr(:, k) = ind_extr(:, k) - 1 + my_lbound

  end subroutine local_extrema

  !*********************************************************************

  subroutine local_extremum(extr_around, local_map, n_extr, local_field)

    ! Look for an extremum at the center of a 3x3 domain.

    real, intent(out):: extr_around
    ! Value of extremum different from the center value. May not be
    ! defined, but guaranteed to be defined if the center is a local
    ! extremum.

    integer, intent(inout):: local_map(- 1:, - 1:) ! (3, 3) Map of
    ! identification numbers of extrema. "local_map" is positive for a
    ! maximum, negative for a minimum.  Only local_map(0, 0) is
    ! possibly updated.

    integer, intent(inout):: n_extr

    real, intent(in):: local_field(- 1:, - 1:) ! (3, 3) Missing value
    ! is huge(0.).

    ! Local:
    logical diff_center(3, 3) ! values different from the center value
    integer count_diff

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

    if (all(local_field /= huge(0.))) then
       diff_center = local_field /= local_field(0, 0)
       count_diff = count(diff_center) ! {between 0 and 8}

       if (count_diff /= 0) then
          ! {count_diff >= 1}
          extr_around = maxval(local_field, diff_center)

          if (local_field(0, 0) > extr_around) then
             ! {local maximum, possibly degenerate}
             if (count_diff == 8 .or. all(local_map <= 0)) then
                ! {not a degenerate maximum or none of the previously
                ! examined values is a maximum}
                n_extr = n_extr + 1
                local_map(0, 0) = n_extr
             end if
          else
             ! {local_field(0, 0) < maxval(local_field, diff_center)}
             extr_around = minval(local_field, diff_center)

             if (local_field(0, 0) < extr_around &
                  .and. (count_diff == 8 .or. all(local_map >= 0))) then
                ! {local minimum, not degenerate or none of previously
                ! examined values is a minimum}
                n_extr = n_extr + 1
                local_map(0, 0) = - n_extr
             end if
          end if
       end if
    end if

  end subroutine local_extremum

end module local_extrema_m