get_1_outerm.f90

module get_1_outerm_m

  implicit none

contains

  type(ssh_contour) function get_1_outerm(cyclone, coord_extr, &
       innermost_level_2, outside_points, ssh, corner, step)

    ! This procedure gets one outermost good contour, if it can find
    ! one. If the procedure cannot find an outermost contour, it
    ! returns a null ssh_contour. The procedure also finds some inner
    ! good contours. The number of good contours stored, including the
    ! outermost contour, is lower than or equal to n_max_cont. If we
    ! find more than n_max_cont - 1 good contours excluding the
    ! outermost contour, then we store only the first n_max_cont - 1
    ! found, because they provide a better sample of the ssh range. If
    ! the procedure cannot find an outermost contour, it does not
    ! store inner contours. An outermost contour is not found if, and
    ! only if, there is no good contour at innermost_level_2. If an
    ! outermost contour is found then its level is farther from the
    ! extremum than innermost_level_2. In the array of stored
    ! contours, the contours are in monotonic order of ssh: ascending
    ! if the extremum is a minimum, descending if the extremum is a
    ! maximum. The area is only computed for the outermost contour.

    ! Method: dichotomy on level of ssh.

    ! 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. The
    ! longitudes of outside points and the target extremum should be
    ! in the interval [corner(1), corner(1) + step(1) * (size(ssh, 1)
    ! - 1)] as there will be no automatic shifting by a mulitple of
    ! 360 degrees.

    ! Libraries:
    use contour_531, only: polyline
    use jumble, only: assert

    use ccw_orient_m, only: ccw_orient
    use derived_types, only: ssh_contour, null_ssh_contour
    use good_contour_m, only: good_contour
    use spher_polyline_area_m, only: spher_polyline_area

    logical, intent(in):: cyclone
    real, intent(in):: coord_extr(:) ! (2)

    real, intent(in):: innermost_level_2
    ! ssh level of the innermost contour to consider around the target
    ! extremum, in m. Assume: innermost_level_2 < extremum for a
    ! maximum, > extremum for a minimum.

    real, intent(in):: outside_points(:, :) ! (2, :) longitude and
    ! latitude, in rad, of all the significant extrema, except the
    ! target extremum

    real, intent(in):: ssh(:, :) ! in m

    real, intent(in):: corner(:) ! (2)
    ! longitude and latitude corresponding to ssh(1,1), in rad

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

    ! Local:
    real level_try, level_good, level_bad ! in m
    real, parameter:: precision = 5e-4 ! in m
    logical mask(size(ssh, 1), size(ssh, 2))
    integer, parameter:: n_max_cont = 31

    type(ssh_contour) cont_list(n_max_cont + 1)
    ! Contour list. The outermost contour is element n_cont. The
    ! contours are in monotonic order of ssh.

    integer n_cont ! number of good contours found and stored

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

    n_cont = 0
    cont_list(1)%polyline = good_contour(corner, step, ssh, innermost_level_2, &
         coord_extr, outside_points)

    if (cont_list(1)%n_points == 0) then
       get_1_outerm = null_ssh_contour()
    else
       n_cont = 1
       cont_list(1)%ssh = innermost_level_2
       level_good = innermost_level_2
       mask = ssh /= huge(0.)
       level_try = merge(maxval(ssh, mask), minval(ssh, mask), cyclone)
       call assert(merge(level_try > level_good, level_try < level_good, &
            cyclone), "get_1_outerm level_try")
       cont_list(2)%polyline = good_contour(corner, step, ssh, level_try, &
            coord_extr, outside_points)

       if (cont_list(2)%n_points /= 0) then
          ! This can only happen if there are degenerate extrema
          ! everywhere around the path of contour. So it should
          ! probably never happen.
          n_cont = 2
          cont_list(2)%ssh = level_try
       else
          ! Assertion: {no good contour at maxval or minval}
          level_bad = level_try

          do while (abs(level_bad - level_good) > precision)
             level_try = (level_good + level_bad) / 2.
             cont_list(n_cont + 1)%polyline = good_contour(corner, step, ssh, &
                  level_try, coord_extr, outside_points)

             if (cont_list(n_cont + 1)%n_points /= 0) then
                cont_list(n_cont + 1)%ssh = level_try

                if (n_cont == n_max_cont) then
                   ! Replace the previous good contour found
                   cont_list(n_cont) = cont_list(n_cont + 1)
                else
                   n_cont = n_cont + 1
                end if

                level_good = level_try
             else
                level_bad = level_try
             end if

             ! Assertion: cont_list(n_cont)%polyline == good_contour(. . .,
             ! level_good) and cont_list(n_cont)%%n_points /= 0 and
             ! no good contour at level_bad
          end do
       end if

       cont_list(n_cont)%area = spher_polyline_area(cont_list(n_cont)%polyline)
       call ccw_orient(cont_list(n_cont))
       get_1_outerm = cont_list(n_cont)
    end if

  end function get_1_outerm

end module get_1_outerm_m