get_1_outerm.f90

module get_1_outerm_m

  implicit none

contains

  subroutine get_1_outerm(out_cont, cont_list, cont_list_proj, n_cont, &
       cyclone, extr_coord_proj, innermost_level_2, outside_points, ssh, &
       corner, step)

    ! This procedure gets one outermost good contour with sufficient
    ! area, if it can find one. If the procedure cannot find an
    ! outermost contour, it returns a null ssh_contour. The procedure
    ! also returns some inner good contours. The maximum number of
    ! inner good contours returned is given by the size of cont_list -
    ! 1. If we find more inner good contours, then we return only the
    ! first ones found, because they provide a better sample of the
    ! ssh range. If the procedure cannot find an outermost contour, it
    ! does not return any inner contour. An outermost contour is not
    ! returned if there is no good contour at innermost_level_2 or if
    ! the outermost contour does not have sufficient area. If an
    ! outermost contour is found then its level is farther from the
    ! extremum than innermost_level_2. The outermost contour is
    ! duplicated in the array cont_list. In cont_list_proj, the
    ! contours, including the outermost contour, are in monotonic
    ! order of ssh: ascending if the extremum is a minimum, descending
    ! if the extremum is a maximum. The area is not computed for
    ! contours in cont_list_proj. On return, the array
    ! abs(ediff1d(cont_list_proj(:n_cont)%ssh)) is in descending
    ! order. (See notes for a proof.) On return, n_cont == 1 means we
    ! have only found one good contour at innermost_level_2 and it has
    ! sufficient area. On return, n_cont == 0 .eqv. out_cont ==
    ! null_polyline().

    ! 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, convert_to_reg_coord, convert_to_ind
    use jumble, only: assert

    use ccw_orient_m, only: ccw_orient
    use config_m, only: min_area
    use derived_types, only: ssh_contour, null_ssh_contour
    use good_contour_m, only: good_contour
    use input_ssh_m, only: corner_whole => corner
    use spher_polyline_area_m, only: spher_polyline_area

    type(ssh_contour), intent(out):: out_cont
    ! points given by longitude, latitude, in rad

    type(polyline), intent(out):: cont_list(:) ! (n_max_cont)
    ! Contour list, with points given by longitude and latitude, in
    ! radians. Defined only for subscripts 1:n_cont. If n_cont >= 1,
    ! the outermost contour is element with subscript n_cont. For i <=
    ! n_cont, the order of points in cont_list(i) (clockwise or
    ! counter-clockwise) is not specified.

    type(ssh_contour), intent(out):: cont_list_proj(:) ! (n_max_cont)
    ! Contour list, with points given by projection
    ! coordinates. Defined only for subscripts 1:n_cont. If n_cont >=
    ! 1, the outermost contour is element with subscript n_cont. The
    ! contours are in monotonic order of ssh. The contour with a given
    ! subscript in cont_list_proj corresponds to the contour with the
    ! same subscript in cont_list. More precisely, the point with
    ! subscript j in contour i of cont_list_proj corresponds to the
    ! point with subscript j in contour i of cont_list (this is
    ! important for mean_speed). For i <= n_cont,
    ! cont_list_proj(i)%area has its default initialization value and
    ! the order of points in cont_list_proj(i) (clockwise or
    ! counter-clockwise) is not specified.

    integer, intent(out):: n_cont
    ! number of good contours found and stored, 0 <= n_cont <= n_max_cont - 1

    logical, intent(in):: cyclone

    real, intent(in):: extr_coord_proj(:) ! (2)
    ! coordinates of extremum in projection space

    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. A
    ! priori, this is not the corner of the grid of the NetCDF file.

    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 n_max_cont ! >= 3
    integer i
    real corner_proj(2)
    real outside_points_proj(2, size(outside_points, 2))

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

    corner_proj = (corner - corner_whole) / step + 1.
    outside_points_proj = convert_to_ind(outside_points, corner_whole, step)
    n_max_cont = size(cont_list_proj)
    n_cont = 0
    cont_list_proj(1)%polyline = good_contour(corner_proj, ssh, &
         innermost_level_2, extr_coord_proj, outside_points_proj)

    test_n_points: if (.not. cont_list_proj(1)%closed) then
       ! cont_list_proj(1)%n_points == 0
       out_cont = null_ssh_contour()
    else test_n_points
       n_cont = 1
       cont_list_proj(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_proj(2)%polyline = good_contour(corner_proj, ssh, level_try, &
            extr_coord_proj, outside_points_proj)

       if (cont_list_proj(2)%closed) then
          ! cont_list_proj(2)%n_points /= 0
          ! 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_proj(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_proj(n_cont + 1)%polyline = good_contour(corner_proj, &
                  ssh, level_try, extr_coord_proj, outside_points_proj)

             if (cont_list_proj(n_cont + 1)%closed) then
                ! cont_list_proj(n_cont + 1)%n_points /= 0
                cont_list_proj(n_cont + 1)%ssh = level_try

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

                level_good = level_try
             else
                level_bad = level_try
             end if

             ! Assertion: cont_list_proj(n_cont)%polyline ==
             ! good_contour(. . ., level_good) and
             ! cont_list_proj(n_cont)%closed and no good contour at
             ! level_bad
          end do
          ! Assertion: n_cont <= n_max_cont - 1
       end if

       out_cont%polyline &
            = convert_to_reg_coord(cont_list_proj(n_cont)%polyline, &
            corner_whole, step)
       out_cont%area = spher_polyline_area(out_cont%polyline)

       if (abs(out_cont%area) < min_area) then
          out_cont = null_ssh_contour()
          n_cont = 0
       else
          ! n_cont >= 1
          forall (i = 1:n_cont - 1) cont_list(i) &
               = convert_to_reg_coord(cont_list_proj(i)%polyline, &
               corner_whole, step)
          cont_list(n_cont) = out_cont%polyline
          out_cont%ssh = cont_list_proj(n_cont)%ssh

          call ccw_orient(out_cont)
          ! (Note we have to orient out_cont after copying it to
          ! cont_list, so that cont_list corresponds to cont_list_proj
          ! point by point.)
       end if
    end if test_n_points

  end subroutine get_1_outerm

end module get_1_outerm_m