module get_1_outerm_m
implicit none
contains
subroutine get_1_outerm(out_cont, cont_list, n_cont, cyclone, coord_extr, &
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, 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 only computed for the
! outermost contour. On return, the array
! abs(ediff1d(cont_list(: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.
! 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
type(ssh_contour), intent(out):: cont_list(:) ! (n_max_cont)
! Contour list. 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. For i <= n_cont - 1,
! cont_list(i)%area has its default initialization value and the
! order of points in cont_list(i)%polyline (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):: coord_extr(:) ! (2) longitude and latitude, in rad
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), coord_extr_proj(2)
real outside_points_proj(2, size(outside_points, 2))
!-----------------------------------------------------------------
corner_proj = (corner - corner_whole) / step + 1.
coord_extr_proj = (coord_extr - corner_whole) / step + 1.
outside_points_proj = convert_to_ind(outside_points, corner_whole, step)
n_max_cont = size(cont_list)
n_cont = 0
cont_list(1)%polyline = convert_to_reg_coord(good_contour(corner_proj, &
ssh, innermost_level_2, coord_extr_proj, outside_points_proj), &
corner_whole, step)
test_n_points: if (cont_list(1)%n_points == 0) then
out_cont = null_ssh_contour()
else test_n_points
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 = convert_to_reg_coord(good_contour(corner_proj, &
ssh, level_try, coord_extr_proj, outside_points_proj), &
corner_whole, step)
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 &
= convert_to_reg_coord(good_contour(corner_proj, ssh, &
level_try, coord_extr_proj, outside_points_proj), &
corner_whole, step)
if (cont_list(n_cont + 1)%n_points /= 0) then
cont_list(n_cont + 1)%ssh = level_try
if (n_cont == n_max_cont - 1) then
! Replace the previous good contour found
cont_list(n_cont) = cont_list(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(n_cont)%polyline == good_contour(. . .,
! level_good) and cont_list(n_cont)%%n_points /= 0 and
! no good contour at level_bad
end do
! Assertion: n_cont <= n_max_cont - 1
end if
cont_list(n_cont)%area = spher_polyline_area(cont_list(n_cont)%polyline)
call ccw_orient(cont_list(n_cont))
if (cont_list(n_cont)%area < min_area) then
out_cont = null_ssh_contour()
n_cont = 0
else
out_cont = cont_list(n_cont)
end if
end if test_n_points
end subroutine get_1_outerm
end module get_1_outerm_m