overlap.f90

module overlap_m

  implicit none

contains

  subroutine overlap(flow, nlon, nlat, periodic, dist_lim, e_overestim, k, &
       delta, j)

    ! This procedure finds edges between flow(j - delta) and flow(j),
    ! corresponding to date indices k - delta and k. It writes these
    ! edges. It updates flow(j - delta)%list%delta_out and
    ! flow(j)%list%delta_in.

    ! Libraries:
    use contour_531, only: polyline
    use gpc_f, only: gpc_polygon_clip_f, GPC_INT, polygon
    use nr_util, only: twopi

    use candidate_overlap_m, only: candidate_overlap
    use derived_types, only: snapshot
    use spher_polygon_area_m, only: spher_polygon_area
    use spher_polyline_area_m, only: spher_polyline_area
    use unit_edge_m, only: unit_edge
    use weight_m, only: weight

    type(snapshot), intent(inout):: flow(:) ! (max_delta + 1)
    integer, intent(in):: nlon, nlat
    logical, intent(in):: periodic ! grid is periodic in longitude

    integer, intent(in):: dist_lim
    ! We look for an overlapping eddy at dist_lim (in grid points) of
    ! the extremum of a given eddy.

    integer, intent(in):: e_overestim
    ! over-estimation of the number of eddies at each date
    
    integer, intent(in):: k ! date index (0-based)
    integer, intent(in):: delta ! between 1 and max_delta

    integer, intent(in):: j
    ! position in time window, between 2 and max_delta + 1

    ! Local:

    integer i1 ! eddy index of predecessor
    integer i2 ! eddy index of successor
    integer l, n_select
    type(polyline) polyline_1, polyline_2
    type(polygon) res_pol

    integer, allocatable:: selection(:)
    ! identifying numbers of a selection of extrema

    ! Window around each extremum:
    integer llc(2) ! indices in global grid of lower left corner
    integer urc(2) ! indices in global grid of upper right corner

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

    loop_i1: do i1 = 1, flow(j - delta)%number_extr
       test_valid: if (flow(j - delta)%list(i1)%valid) then
          ! Define the geographical window around each eddy extremum:

          llc = flow(j - delta)%ind_extr(:, i1) - dist_lim
          urc = flow(j - delta)%ind_extr(:, i1) + dist_lim

          llc(2) = max(llc(2), 1)
          urc(2) = min(urc(2), nlat)

          if (.not. periodic) then
             llc(1) = max(llc(1), 1)
             urc(1) = min(urc(1), nlon)
          end if

          ! Pre-select potential successors:
          selection = candidate_overlap(flow(j)%extr_map(llc(1):urc(1), &
               llc(2):urc(2)), flow(j)%list, &
               flow(j - delta)%list(i1)%delta_out, delta)
          
          n_select = size(selection)

          if (n_select /= 0) then
             if (flow(j - delta)%list(i1)%speed_cont%n_points /= 0) then
                polyline_1 = flow(j - delta)%list(i1)%speed_cont%polyline
             else
                polyline_1 = flow(j - delta)%list(i1)%out_cont%polyline
             end if
          end if

          DO l = 1, n_select
             i2 = selection(l)
             ! Assertion: {flow(j - delta)%list(i1)%delta_out >=
             ! delta .or. flow(j)%list(i2)%delta_in >= delta}

             if (flow(j)%list(i2)%speed_cont%n_points /= 0) then
                polyline_2 = flow(j)%list(i2)%speed_cont%polyline
             else
                polyline_2 = flow(j)%list(i2)%out_cont%polyline
             end if

             ! Shift the longitudes of polyline_2 to values close to the
             ! longitude of extremum i1:
             polyline_2%points(1, :) = polyline_2%points(1, :) &
                  + floor((flow(j - delta)%list(i1)%coord_extr(1) &
                  - flow(j)%list(i2)%coord_extr(1)) / twopi + 0.5) * twopi

             call gpc_polygon_clip_f(GPC_INT, &
                  polygon(nparts = 1, part = [polyline_1], hole = [.false.]), &
                  polygon(nparts = 1, part = [polyline_2], hole = [.false.]), &
                  res_pol)

             if (res_pol%nparts /= 0) then
                ! polyline_1 and polyline_2 overlap
                if (spher_polygon_area(res_pol) >= 0.25 &
                     * min(spher_polyline_area(polyline_1), &
                     spher_polyline_area(polyline_2))) then
                   write(unit_edge, fmt = *) (k - delta) * e_overestim + i1, &
                        k * e_overestim + i2, &
                        weight(flow(j - delta)%list(i1), &
                        flow(j)%list(i2))
                   flow(j - delta)%list(i1)%delta_out &
                        = min(flow(j - delta)%list(i1)%delta_out, delta)
                   flow(j)%list(i2)%delta_in &
                        = min(flow(j)%list(i2)%delta_in, delta)
                end if
             end if
          end DO
       end if test_valid
    end do loop_i1

  end subroutine overlap
  
end module overlap_m