Skip to content
Snippets Groups Projects
Commit 0640e81d authored by Lionel GUEZ's avatar Lionel GUEZ
Browse files

Use associate constructs to shorten expressions

parent 0e1f9399
No related branches found
No related tags found
No related merge requests found
......@@ -109,39 +109,41 @@ program inst_eddies
! max-speed contours.
do i = 1, s%number_extr
if (s%list(i)%valid) then
! Restrict the field to the outermost contour:
llc = floor(convert_to_ind(minval(s%list(i)%out_cont%points, &
dim = 2), corner, step))
urc = ceiling(convert_to_ind(maxval( &
s%list(i)%out_cont%points, dim = 2), corner, step))
! Should have no effect except because of roundup error:
urc(2) = min(urc(2), nlat)
if (.not. periodic) urc(1) = min(urc(1), nlon)
corner_window = corner + (llc - 1) * step
outside_points = nearby_extr(s%extr_map(llc(1):urc(1), &
llc(2):urc(2)), s%list, i)
if (periodic) outside_points(1, :) = outside_points(1, :) &
+ ceiling((corner_window(1) - outside_points(1, :)) / twopi) &
* twopi
! (Shift the longitude of each point to a value close to
! the longitude of the target extremum.)
call set_max_speed(s%list(i), s%list(i)%extr%coord_proj - llc + 1, &
outside_points, ssh(llc(1):urc(1), llc(2):urc(2)), &
u(llc(1):urc(1), llc(2):urc(2)), &
v(llc(1):urc(1), llc(2):urc(2)), corner_window, step)
else
s%list(i)%speed_cont = null_ssh_contour()
s%list(i)%max_speed = missing_speed
s%list(i)%radius4 = 0
end if
associate (e => s%list(i))
if (e%valid) then
! Restrict the field to the outermost contour:
llc = floor(convert_to_ind(minval(e%out_cont%points, dim = 2), &
corner, step))
urc = ceiling(convert_to_ind(maxval(e%out_cont%points, dim = 2), &
corner, step))
! Should have no effect except because of roundup error:
urc(2) = min(urc(2), nlat)
if (.not. periodic) urc(1) = min(urc(1), nlon)
corner_window = corner + (llc - 1) * step
outside_points = nearby_extr(s%extr_map(llc(1):urc(1), &
llc(2):urc(2)), s%list, i)
if (periodic) outside_points(1, :) = outside_points(1, :) &
+ ceiling((corner_window(1) - outside_points(1, :)) / twopi) &
* twopi
! (Shift the longitude of each point to a value close to
! the longitude of the target extremum.)
call set_max_speed(e, e%extr%coord_proj - llc + 1, outside_points, &
ssh(llc(1):urc(1), llc(2):urc(2)), &
u(llc(1):urc(1), llc(2):urc(2)), &
v(llc(1):urc(1), llc(2):urc(2)), corner_window, step)
else
e%speed_cont = null_ssh_contour()
e%max_speed = missing_speed
e%radius4 = 0
end if
end associate
end do
call cpu_time(t1)
......
......@@ -52,72 +52,73 @@ contains
!-----------------------------------------------------------------------
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)%list(i1)%extr%coord_proj - dist_lim
urc = flow(j - delta)%list(i1)%extr%coord_proj + 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)%extr%coord(1) &
- flow(j)%list(i2)%extr%coord(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) >= min_inters &
* min(abs(spher_polyline_area(polyline_1)), &
abs(spher_polyline_area(polyline_2)))) then
write(unit_edge, fmt = *) (k - delta) * e_overestim + i1, &
k * e_overestim + 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
associate (e1 => flow(j - delta)%list(i1))
test_valid: if (e1%valid) then
! Define the geographical window around each eddy extremum:
llc = e1%extr%coord_proj - dist_lim
urc = e1%extr%coord_proj + 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, e1%delta_out, delta)
n_select = size(selection)
if (n_select /= 0) then
if (e1%speed_cont%n_points /= 0) then
polyline_1 = e1%speed_cont%polyline
else
polyline_1 = e1%out_cont%polyline
end if
end if
DO l = 1, n_select
i2 = selection(l)
associate (e2 => flow(j)%list(i2))
! Assertion: {e1%delta_out >= delta .or. e2%delta_in
! >= delta}
if (e2%speed_cont%n_points /= 0) then
polyline_2 = e2%speed_cont%polyline
else
polyline_2 = e2%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((e1%extr%coord(1) - e2%extr%coord(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) >= min_inters &
* min(abs(spher_polyline_area(polyline_1)), &
abs(spher_polyline_area(polyline_2)))) then
write(unit_edge, fmt = *) &
(k - delta) * e_overestim + i1, k * e_overestim + i2
e1%delta_out = min(e1%delta_out, delta)
e2%delta_in = min(e2%delta_in, delta)
end if
end if
end associate
end DO
end if test_valid
end associate
end do loop_i1
end subroutine overlap
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment