From 2a6cdf7dae07a6ca131eb873b2470505b4772a83 Mon Sep 17 00:00:00 2001 From: Lionel GUEZ <guez@lmd.ens.fr> Date: Wed, 20 Dec 2017 18:24:19 +0100 Subject: [PATCH] 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. --- Sources/get_1_outerm.f | 58 ++++++++++++++++++--------------------- Sources/set_all_outerm.f | 12 ++++---- Tests/test_get_1_outerm.f | 22 +++++++-------- 3 files changed, 44 insertions(+), 48 deletions(-) diff --git a/Sources/get_1_outerm.f b/Sources/get_1_outerm.f index a4e163d4..662a30f5 100644 --- a/Sources/get_1_outerm.f +++ b/Sources/get_1_outerm.f @@ -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 diff --git a/Sources/set_all_outerm.f b/Sources/set_all_outerm.f index e9836604..94f030eb 100644 --- a/Sources/set_all_outerm.f +++ b/Sources/set_all_outerm.f @@ -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)), & diff --git a/Tests/test_get_1_outerm.f b/Tests/test_get_1_outerm.f index 08ba1113..ae3a1aaf 100644 --- a/Tests/test_get_1_outerm.f +++ b/Tests/test_get_1_outerm.f @@ -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 -- GitLab