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

Make get_1_outerm a function instead of a subroutine. The argument e 

with some components intent(in) and one component intent(out) was less
clear. Also, replace argument ind_targ_extr by argument i. This is
simpler both in get_1_outerm and in set_all_outerm.

There is now no need to use an eddy in test_get_1_outerm.
parent bce9ff0a
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Sources/get_1_outerm.f
+ 27
− 31
View file @ 2a6cdf7d
......@@ -4,25 +4,26 @@ module get_1_outerm_m
contains
subroutine get_1_outerm(e, ind_targ_extr, innermost_level, &
extr_map, ssh, corner, step)
type(ssh_contour) function get_1_outerm(ssh_extremum, cyclone, coord_extr, &
i, innermost_level, extr_map, ssh, corner, step)
! Defines e%outermost_contour. Method: dichotomy on level of ssh.
! Get one outermost contour. Method: dichotomy on level of ssh.
use contour_531, only: polyline, convert_to_reg_coord
use derived_types, only: eddy
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
type(eddy), intent(inout):: e
! Intent(in): e%ssh_extremum, e%cyclone, e%coord_extr
! Intetn(out): e%outermost_contour
real, intent(in):: ssh_extremum
logical, intent(in):: cyclone
real, intent(in):: coord_extr(:) ! (2)
integer, intent(in):: ind_targ_extr(:)
! (2) indices of the target extremum, that is, the extremum of eddy e
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
......@@ -38,8 +39,6 @@ contains
! Local:
integer i ! identifying number of the eddy
real, allocatable:: nearby_extr(:, :) ! (2, :) longitude and
! latitude, in rad, of all the significant extrema, except the
! target extremum
......@@ -50,26 +49,24 @@ contains
!-----------------------------------------------------------------
i = extr_map(ind_targ_extr(1), ind_targ_extr(2))
nearby_extr = convert_to_reg_coord(argwhere(extr_map > 0 &
.and. extr_map /= i), corner, step)
e%outermost_contour%polyline = good_contour(corner, step, ssh, &
innermost_level, e%coord_extr, nearby_extr)
get_1_outerm%polyline = good_contour(corner, step, ssh, innermost_level, &
coord_extr, nearby_extr)
if (e%outermost_contour%n_points == 0) then
e%outermost_contour%ssh = e%ssh_extremum
e%outermost_contour%area = 0.
if (get_1_outerm%n_points == 0) then
get_1_outerm%ssh = ssh_extremum
get_1_outerm%area = 0.
else
level_good = innermost_level
level_try = merge(maxval(ssh), minval(ssh), e%cyclone)
call assert((level_try - e%ssh_extremum) &
/ (level_good - e%ssh_extremum) > 1., &
"get_1_outerm level_try")
level_try = merge(maxval(ssh), minval(ssh), 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, &
e%coord_extr, nearby_extr)
coord_extr, nearby_extr)
if (tentative_contour%n_points /= 0) then
e%outermost_contour%polyline = tentative_contour
get_1_outerm%polyline = tentative_contour
level_good = level_try
else
! {no good contour at level_try}
......@@ -78,25 +75,24 @@ contains
do while (abs(level_bad - level_good) > precision)
level_try = (level_good + level_bad) / 2.
tentative_contour = good_contour(corner, step, ssh, level_try, &
e%coord_extr, nearby_extr)
coord_extr, nearby_extr)
if (tentative_contour%n_points /= 0) then
level_good = level_try
e%outermost_contour%polyline = tentative_contour
get_1_outerm%polyline = tentative_contour
else
level_bad = level_try
end if
! {e%outermost_contour%polyline == good_contour(. . .,
! level_good) and e%outermost_contour%%n_points /= 0 and
! {get_1_outerm%polyline == good_contour(. . .,
! level_good) and get_1_outerm%%n_points /= 0 and
! no good contour at level_bad}
end do
end if
e%outermost_contour%ssh = level_good
e%outermost_contour%area &
= spherical_polygon_area(e%outermost_contour%polyline)
get_1_outerm%ssh = level_good
get_1_outerm%area = spherical_polygon_area(get_1_outerm%polyline)
end if
end subroutine get_1_outerm
end function get_1_outerm
end module get_1_outerm_m
Sources/set_all_outerm.f
+ 6
− 6
View file @ 2a6cdf7d
......@@ -55,9 +55,6 @@ contains
real corner_window(2, s%number_vis_eddies) ! longitude and latitude, in rad
integer ind_targ_extr(2, s%number_vis_eddies)
! indices in the window of the target extremum
!--------------------------------------------------------------
! Define the geographical window around each eddy extremum:
......@@ -66,7 +63,6 @@ contains
urc(:, i) = min(s%ind_extr(:, i) + max_radius, &
[size(ssh, 1), size(ssh, 2)])
corner_window(:, i) = corner + (llc(:, i) - 1) * step
ind_targ_extr(:, i) = s%ind_extr(:, i) - llc(:, i) + 1
end forall
if (min_amp == 0.) then
......@@ -77,7 +73,9 @@ contains
do i = 1, s%number_vis_eddies
if (flat_extr(i)) then
call get_1_outerm(s%list_vis(i), ind_targ_extr(:, i), &
s%list_vis(i)%outermost_contour &
= get_1_outerm(s%list_vis(i)%ssh_extremum, &
s%list_vis(i)%cyclone, s%list_vis(i)%coord_extr, i, &
innermost_level(i), &
s%extr_map(llc(1, i):urc(1, i), llc(2, i):urc(2, i)), &
ssh(llc(1, i):urc(1, i), llc(2, i):urc(2, i)), &
......@@ -115,7 +113,9 @@ contains
do i = 1, s%number_vis_eddies
if (s%list_vis(i)%suff_amp .and. noise_around(i) &
.or. .not. flat_extr(i)) then
call get_1_outerm(s%list_vis(i), ind_targ_extr(:, i), &
s%list_vis(i)%outermost_contour &
= get_1_outerm(s%list_vis(i)%ssh_extremum, &
s%list_vis(i)%cyclone, s%list_vis(i)%coord_extr, i, &
innermost_level(i), &
s%extr_map(llc(1, i):urc(1, i), llc(2, i):urc(2, i)), &
ssh(llc(1, i):urc(1, i), llc(2, i):urc(2, i)), &
......
Tests/test_get_1_outerm.f
+ 11
− 11
View file @ 2a6cdf7d
......@@ -2,7 +2,7 @@ program test_get_1_outerm
use, intrinsic:: ISO_FORTRAN_ENV
use derived_types, only: eddy
use derived_types, only: ssh_contour
use netcdf, only: nf90_nowrite
use netcdf95, only: nf95_open, nf95_close, nf95_inq_varid, nf95_get_var, &
nf95_gw_var
......@@ -26,7 +26,7 @@ program test_get_1_outerm
real:: innermost_level = 0.2933! level of innermost contour, for each extremum
logical:: cyclone = .true.
type(eddy) e
type(ssh_contour) outermost_contour
TYPE(shpfileobject) shphandle
integer field_number, shape_number
real:: ssh_extremum = 0.2759
......@@ -70,23 +70,23 @@ program test_get_1_outerm
call nf95_close(ncid)
step = [longitude(2) - longitude(1), latitude(2) - latitude(1)]
e%coord_extr = [longitude(1), latitude(1)] + (ind_targ_extr - 1) * step
e%cyclone = cyclone
e%ssh_extremum = ssh_extremum
call get_1_outerm(e, ind_targ_extr, innermost_level, extr_map, &
ssh, corner = [longitude(1), latitude(1)], step = step)
outermost_contour = get_1_outerm(ssh_extremum, cyclone, &
coord_extr = [longitude(1), latitude(1)] + (ind_targ_extr - 1) * step, &
i = extr_map(ind_targ_extr(1), ind_targ_extr(2)), &
innermost_level = innermost_level, extr_map = extr_map, ssh = ssh, &
corner = [longitude(1), latitude(1)], step = step)
if (e%outermost_contour%closed) then
if (outermost_contour%closed) then
print *, "Radius of disk of equal area: ", &
sqrt(e%outermost_contour%area / pi) / 1e3, "km"
sqrt(outermost_contour%area / pi) / 1e3, "km"
call shp_create_03("test_get_1_outerm", shpt_polygon, shphandle)
call dbf_add_field_03(field_number, shphandle, 'level', ftdouble, &
nwidth = 13, ndecimals = 6)
call shp_append_simple_object_03(shape_number, shphandle, shpt_polygon, &
e%outermost_contour%points / pi * 180.)
outermost_contour%points / pi * 180.)
call dbf_write_attribute_03(shphandle, shape_number, ifield = 0, &
fieldvalue = e%outermost_contour%ssh)
fieldvalue = outermost_contour%ssh)
CALL shpclose(shphandle)
print *, 'Created shapefile "test_get_1_outerm".'
else
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment