get_1_outerm.f

module get_1_outerm_m

  implicit none

contains

  type(ssh_contour) function get_1_outerm(ssh_extremum, cyclone, coord_extr, &
       i, innermost_level, extr_map, ssh, corner, step)

    ! Gets one outermost contour, if it can find one. Method:
    ! dichotomy on level of ssh. If the procedure cannot find an
    ! outermost contour, it return a null polyline, zero area and ssh
    ! equal to ssh extremum. An outermost contour is not found if, and
    ! only if, there is no good contour at innermost level.

    use contour_531, only: polyline, convert_to_reg_coord
    use derived_types, only: ssh_contour
    use good_contour_m, only: good_contour
    use jumble, only: argwhere, pack_indices
    use nr_util, only: assert
    use outermost_possible_level_m, only: outermost_possible_level
    use spherical_polygon_area_m, only: spherical_polygon_area

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

    integer, intent(in):: i
    ! identifying number of the target extremum, that is, the extremum
    ! at coord_extr

    real, intent(in):: innermost_level
    ! ssh level of the innermost contour around the target extremum, in m

    integer, intent(in):: extr_map(:, :)
    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, allocatable:: nearby_extr(:, :) ! (2, :) longitude and
    ! latitude, in rad, of all the significant extrema, except the
    ! target extremum

    real level_try, level_good, level_bad ! in m
    type(polyline) tentative_contour
    real, parameter:: precision = 5e-4 ! in m
    logical mask(size(ssh, 1), size(ssh, 2))

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

    nearby_extr = convert_to_reg_coord(argwhere(extr_map > 0 &
         .and. extr_map /= i), corner, step)
    get_1_outerm%polyline = good_contour(corner, step, ssh, innermost_level, &
         coord_extr, nearby_extr)

    if (get_1_outerm%n_points == 0) then
       get_1_outerm%ssh = ssh_extremum
       get_1_outerm%area = 0.
    else
       level_good = innermost_level
       mask = ssh /= huge(0.)
       level_try = merge(maxval(ssh, mask), minval(ssh, mask), cyclone)
       call assert((level_try - ssh_extremum) &
            / (level_good - ssh_extremum) > 1., "get_1_outerm level_try")
       tentative_contour = good_contour(corner, step, ssh, level_try, &
            coord_extr, nearby_extr)

       if (tentative_contour%n_points /= 0) then
          get_1_outerm%polyline = tentative_contour
          get_1_outerm%ssh = level_try
       else
          ! {no good contour at level_try}
          level_bad = level_try

          do while (abs(level_bad - level_good) > precision)
             level_try = (level_good + level_bad) / 2.
             tentative_contour = good_contour(corner, step, ssh, level_try, &
                  coord_extr, nearby_extr)
             if (tentative_contour%n_points /= 0) then
                level_good = level_try
                get_1_outerm%polyline = tentative_contour
             else
                level_bad = level_try
             end if

             ! {get_1_outerm%polyline == good_contour(. . .,
             ! level_good) and get_1_outerm%%n_points /= 0 and
             ! no good contour at level_bad}
          end do
          
          get_1_outerm%ssh = level_good
       end if

       get_1_outerm%area = spherical_polygon_area(get_1_outerm%polyline)
    end if

  end function get_1_outerm

end module get_1_outerm_m