set_all_extr.f90

module set_all_extr_m

  implicit none

contains

  subroutine set_all_extr(s, step, periodic, ssh, corner)

    ! This procedures finds all the extrema of ssh. Not a function
    ! because "s" is not completely defined.

    use derived_types, only: snapshot
    use input_ssh_m, only: max_radius
    use local_extrema_m, only: local_extrema

    type(snapshot), intent(out):: s
    ! Define only s%extr_map, s%list%extr, s%number_extr,
    ! s%list%cyclone, s%list%innermost_level.

    real, intent(in):: step(:) ! (2) longitude and latitude steps, in rad
    logical, intent(in):: periodic ! grid is periodic in longitude

    real, intent(in):: ssh(1 - merge(max_radius(1), 0, periodic):, :)
    ! (1 - max_radius(1):nlon + max_radius(1), nlat) if the grid is periodic
    ! in longitude, else (nlon, nlat). Sea-surface height, in m.

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

    ! Local:

    integer nlon, nlat
    logical, allocatable:: cyclone(:) ! (s%number_extr)
    integer i, copy

    integer, allocatable:: ind_extr(:, :) ! (2, n_extr)
    ! List of couple of subscripts of extrema of ssh, in the two
    ! dimensions.

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

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

    copy = merge(max_radius(1), 0, periodic)
    nlon = size(ssh, 1) - 2 * copy
    nlat = size(ssh, 2)
    allocate(s%extr_map(1 - copy:nlon + copy, nlat))
    copy = merge(1, 0, periodic)
    call local_extrema(s%extr_map(1 - copy:nlon + copy, :), ind_extr, &
         innermost_level, cyclone, ssh(1 - copy:nlon + copy, :), periodic, &
         my_lbound = [1 - copy, 1])
    ! Assertion: \forall k, ind_extr(1, k) \in \{2 - copy, \dots,
    ! nlon + copy - 1\} and ind_extr(2, k) \in \{2, \dots, nlat - 1\}

    ! Assertion: \forall k, ind_extr(1, k) \in \{1, \dots, nlon\} and
    ! ind_extr(2, k) \in \{2, \dots, nlat - 1\}

    if (periodic) then
       ! Extend in longitude:
       s%extr_map(1 - max_radius(1):- 1, :) &
            = s%extr_map(nlon + 1 - max_radius(1):nlon - 1, :)
       s%extr_map(nlon + 2:nlon + max_radius(1), :) &
            = s%extr_map(2:max_radius(1), :)
    end if

    s%number_extr = size(ind_extr, 2)
    allocate(s%list(s%number_extr))
    s%list%innermost_level = innermost_level

    forall (i = 1:s%number_extr)
       s%list(i)%extr%coord = corner + (ind_extr(:, i) - 1) * step
       ! (Even when periodic, this is within the original NetCDF grid,
       ! that is the grid without duplicated longitudes.)

       s%list(i)%extr%coord_proj = ind_extr(:, i)
       s%list(i)%extr%ssh = ssh(ind_extr(1, i), ind_extr(2, i))
       s%list(i)%cyclone = cyclone(i)
    end forall

  end subroutine set_all_extr

end module set_all_extr_m