From 7d99eff6bb3948039d8afa8b0f4e062f7ec2d5ce Mon Sep 17 00:00:00 2001
From: Lionel GUEZ <guez@lmd.ens.fr>
Date: Fri, 15 Dec 2017 13:44:15 +0100
Subject: [PATCH] In procedure get_snapshot, encountered rounding error that
led to out-of-range index so add min when calculating urc.
Bug fix in procedure good_contour. polygon_contains_point should not
be called if polyline is not closed.
Procedure polygon_contains_point became unpure so remove pure for
procedure inside_4 too.
In program test_good_contour, read corner and step from namelist
rather than from coordinates in NetCDF file. Makes it possible to use
the same NetCDF file with different coordinates. Also read
outside_points from a file rather than standard input so there is no
need to specify in advance the number of outside points.
Take into account the possibility that no outermost contour was found
in program test_set_outermost_contour. Also, use fixed names for input
files because it is not convenient to require a particular order of
these two files on the command line.
---
Sources/get_snapshot.f | 7 ++-
Sources/good_contour.f | 29 +++++------
Sources/inside_4.f | 2 +-
Tests/Stdin/good_contour_nml_1.txt | 1 +
...contour_in2.txt => good_contour_nml_2.txt} | 5 +-
...contour_in.txt => no_good_contour_nml.txt} | 3 --
...od_contour_in.txt => outside_points_1.csv} | 2 -
Tests/Stdin/outside_points_2.csv | 2 +
Tests/test_good_contour.f | 38 ++++++++-------
Tests/test_set_outermost_contour.f | 48 +++++++++----------
Tests/tests.json | 17 ++++---
11 files changed, 78 insertions(+), 76 deletions(-)
create mode 100644 Tests/Stdin/good_contour_nml_1.txt
rename Tests/Stdin/{good_contour_in2.txt => good_contour_nml_2.txt} (51%)
rename Tests/Stdin/{no_good_contour_in.txt => no_good_contour_nml.txt} (77%)
rename Tests/Stdin/{good_contour_in.txt => outside_points_1.csv} (56%)
create mode 100644 Tests/Stdin/outside_points_2.csv
diff --git a/Sources/get_snapshot.f b/Sources/get_snapshot.f
index 58af9c39..75e3b453 100644
--- a/Sources/get_snapshot.f
+++ b/Sources/get_snapshot.f
@@ -179,9 +179,12 @@ contains
llc(:, i) = floor(convert_to_ind(minval( &
s%list_vis(i)%outermost_contour%points, dim = 2), corner, &
step))
- urc(:, i) = ceiling(convert_to_ind(maxval( &
+
+ urc(:, i) = min(ceiling(convert_to_ind(maxval( &
s%list_vis(i)%outermost_contour%points, dim = 2), corner, &
- step))
+ step)), [nlon, nlat])
+ ! (min should have no effect except because of roundup error)
+
corner_window(:, i) = corner + (llc(:, i) - 1) * step
ind_targ_extr(:, i) = s%ind_extr(:, i) - llc(:, i) + 1
diff --git a/Sources/good_contour.f b/Sources/good_contour.f
index 13962e80..00bfc038 100644
--- a/Sources/good_contour.f
+++ b/Sources/good_contour.f
@@ -29,7 +29,7 @@ contains
integer j ! index in outside_points
integer n_cont ! number of contours
integer n_out ! number of outside points
-
+
!--------------------------------------------------------------
n_out = size(outside_points, 2)
@@ -37,20 +37,21 @@ contains
n_cont = size(contours)
found_good_contour = .false.
i = 1
-
+
do while (.not. found_good_contour .and. i <= n_cont)
- if (contours(i)%closed .and. polygon_contains_point(contours(i)%points, &
- inside_point)) then
- found_good_contour = .true. ! maybe
- ! Does contours(i) contain one of outside_points?
- j = 1
+ if (contours(i)%closed) then
+ if (polygon_contains_point(contours(i)%points, inside_point)) then
+ found_good_contour = .true. ! maybe
+ ! Does contours(i) contain one of outside_points?
+ j = 1
- do while (found_good_contour .and. j <= n_out)
- found_good_contour &
- = .not. polygon_contains_point(contours(i)%points, &
- outside_points(:, j))
- j = j + 1
- end do
+ do while (found_good_contour .and. j <= n_out)
+ found_good_contour &
+ = .not. polygon_contains_point(contours(i)%points, &
+ outside_points(:, j))
+ j = j + 1
+ end do
+ end if
end if
i = i + 1
@@ -61,7 +62,7 @@ contains
else
good_contour = null_polyline()
end if
-
+
end function good_contour
end module good_contour_m
diff --git a/Sources/inside_4.f b/Sources/inside_4.f
index a8a96a74..edf19835 100644
--- a/Sources/inside_4.f
+++ b/Sources/inside_4.f
@@ -4,7 +4,7 @@ module inside_4_m
contains
- pure logical function inside_4(distance, center, v)
+ logical function inside_4(distance, center, v)
! This procedure returns true if the four points center \pm
! (distance(1), 0) and center \pm (0, distance(2)) are inside
diff --git a/Tests/Stdin/good_contour_nml_1.txt b/Tests/Stdin/good_contour_nml_1.txt
new file mode 100644
index 00000000..0ebbd508
--- /dev/null
+++ b/Tests/Stdin/good_contour_nml_1.txt
@@ -0,0 +1 @@
+&MAIN_NML corner = -2., -2., step = 0.08, 0.08 /
diff --git a/Tests/Stdin/good_contour_in2.txt b/Tests/Stdin/good_contour_nml_2.txt
similarity index 51%
rename from Tests/Stdin/good_contour_in2.txt
rename to Tests/Stdin/good_contour_nml_2.txt
index dafb4655..8e682e45 100644
--- a/Tests/Stdin/good_contour_in2.txt
+++ b/Tests/Stdin/good_contour_nml_2.txt
@@ -1,5 +1,2 @@
! Case where one of the contours tested does not contain inside point.
-&MAIN_NML INSIDE_POINT=0.3,0 /
-2
--0.5,-1
--0.7,0
+&MAIN_NML INSIDE_POINT=0.3,0, corner = -2., -2., step = 0.08, 0.08/
diff --git a/Tests/Stdin/no_good_contour_in.txt b/Tests/Stdin/no_good_contour_nml.txt
similarity index 77%
rename from Tests/Stdin/no_good_contour_in.txt
rename to Tests/Stdin/no_good_contour_nml.txt
index b6e951b2..61f8cde0 100644
--- a/Tests/Stdin/no_good_contour_in.txt
+++ b/Tests/Stdin/no_good_contour_nml.txt
@@ -1,5 +1,2 @@
! Case where there is no good contour
&MAIN_NML INSIDE_POINT=0.,0. /
-2
--0.5, -1.
-0.3, 0.
diff --git a/Tests/Stdin/good_contour_in.txt b/Tests/Stdin/outside_points_1.csv
similarity index 56%
rename from Tests/Stdin/good_contour_in.txt
rename to Tests/Stdin/outside_points_1.csv
index 017d4aa4..b6fb92e0 100644
--- a/Tests/Stdin/good_contour_in.txt
+++ b/Tests/Stdin/outside_points_1.csv
@@ -1,4 +1,2 @@
-&MAIN_NML /
-2
-0.5, -1.
0.3, 0.
diff --git a/Tests/Stdin/outside_points_2.csv b/Tests/Stdin/outside_points_2.csv
new file mode 100644
index 00000000..3260f5c4
--- /dev/null
+++ b/Tests/Stdin/outside_points_2.csv
@@ -0,0 +1,2 @@
+-0.5, -1.
+-0.7, 0.
diff --git a/Tests/test_good_contour.f b/Tests/test_good_contour.f
index edf48f7c..926e74f7 100644
--- a/Tests/test_good_contour.f
+++ b/Tests/test_good_contour.f
@@ -2,6 +2,7 @@ program test_good_contour
use contour_531, only: polyline
use good_contour_m, only: good_contour
+ use jumble, only: new_unit, count_lines
use netcdf, only: nf90_nowrite
use netcdf95, only: nf95_open, nf95_inq_dimid, nf95_inquire_dimension, &
nf95_get_var, nf95_inq_varid, nf95_close
@@ -13,18 +14,19 @@ program test_good_contour
implicit none
character(len = :), allocatable:: filename
- integer nx, ny, n_out
- real x(2), y(2) ! just the first two values to get the corner and step
+ integer nx, ny, n_out, i
type(polyline) c
real:: level = 5.
real, allocatable:: Z(:, :)
real:: inside_point(2) = [- 0.7, 0.]
real, allocatable:: outside_points(:, :) ! (2, n_out)
+ real:: step(2) = [0.04363323129985824, 0.04363323129985824] ! in rad
+ real:: corner(2) = [0.002181662, - 0.72213] ! in rad
TYPE(shpfileobject) shphandle
integer field_number, shape_number
- integer length, status, ncid, varid, dimid
+ integer length, status, ncid, varid, dimid, unit
- namelist /main_nml/ inside_point, level
+ namelist /main_nml/ inside_point, level, step, corner
!-------------------------------------------------------------------
@@ -38,13 +40,9 @@ program test_good_contour
call nf95_inq_dimid(ncid, "x", dimid)
call nf95_inquire_dimension(ncid, dimid, nclen = nx)
- call nf95_inq_varid(ncid, "x", varid)
- call nf95_get_var(ncid, varid, x)
call nf95_inq_dimid(ncid, "y", dimid)
call nf95_inquire_dimension(ncid, dimid, nclen = ny)
- call nf95_inq_varid(ncid, "y", varid)
- call nf95_get_var(ncid, varid, y)
allocate(z(nx, ny))
@@ -58,18 +56,22 @@ program test_good_contour
read(unit = *, nml = main_nml)
write(unit = *, nml = main_nml)
- print *, "Number of outside points? "
- read *, n_out
+ print *, 'Reading from "outside_points.csv"...'
+ call new_unit(unit)
+ open(unit, file = "outside_points.csv", status = "old", action = "read", &
+ position = "rewind")
+ call count_lines(unit, n_out)
+ rewind(unit)
allocate(outside_points(2, n_out))
- if (n_out /= 0) then
- print *, "Enter outside points."
- read *, outside_points
- end if
- print *, "outside_points = ", outside_points
+ if (n_out /= 0) read (unit, fmt = *) outside_points
+ close(unit)
+ print *, "outside_points:"
+
+ do i = 1, n_out
+ print *, outside_points(:, i)
+ end do
- c = good_contour(corner = [x(1), y(1)], step = [x(2) - x(1), y(2) - y(1)], &
- z = z, level = level, inside_point = inside_point, &
- outside_points = outside_points)
+ c = good_contour(corner, step, z, level, inside_point, outside_points)
if (c%n_points /= 0) then
call shp_create_03("test_good_contour", shpt_polygon, shphandle)
diff --git a/Tests/test_set_outermost_contour.f b/Tests/test_set_outermost_contour.f
index 415e1e7d..9995a584 100644
--- a/Tests/test_set_outermost_contour.f
+++ b/Tests/test_set_outermost_contour.f
@@ -3,11 +3,10 @@ program test_set_outermost_contour
use, intrinsic:: ISO_FORTRAN_ENV
use derived_types, only: eddy
- use jumble, only: get_command_arg_dyn
use netcdf, only: nf90_nowrite
use netcdf95, only: nf95_open, nf95_close, nf95_inq_varid, nf95_get_var, &
nf95_gw_var
- use nr_util, only: pi, assert
+ use nr_util, only: pi
use set_outermost_contour_m, only: set_outermost_contour
use shapelib, only: shpt_polygon, shpfileobject, ftdouble, shpclose
use shapelib_03, only: shp_create_03, dbf_add_field_03, &
@@ -15,7 +14,6 @@ program test_set_outermost_contour
implicit none
- character(len = :), allocatable:: filename
integer ncid, varid
real longitude(2) ! in rad
real latitude(2) ! in rad
@@ -37,28 +35,23 @@ program test_set_outermost_contour
!----------------------------------------------------------------
- call assert(COMMAND_ARGUMENT_COUNT() == 2, &
- "Required arguments: ADT-file extr_map-file")
-
write(unit = error_unit, nml = main_nml)
write(unit = error_unit, fmt = *) "Enter namelist main_nml."
read(unit = *, nml = main_nml)
write(unit = *, nml = main_nml)
- ! No problem of degenerate time coordinate with extr_map-file so
+ ! No problem of degenerate time coordinate with "extr_map.nc" so
! read it first:
- call get_command_arg_dyn(2, filename)
- print *, "Reading from ", filename, "..."
- call nf95_open(filename, nf90_nowrite, ncid)
+ print *, 'Reading from "extr_map.nc"...'
+ call nf95_open("extr_map.nc", nf90_nowrite, ncid)
call nf95_inq_varid(ncid, "extr_map", varid)
call nf95_gw_var(ncid, varid, extr_map)
call nf95_close(ncid)
- call get_command_arg_dyn(1, filename)
- print *, "Reading from ", filename, "..."
- call nf95_open(filename, nf90_nowrite, ncid)
+ print *, 'Reading from "h.nc"...'
+ call nf95_open("h.nc", nf90_nowrite, ncid)
call nf95_inq_varid(ncid, "lon", varid)
call nf95_get_var(ncid, varid, longitude)
@@ -83,18 +76,21 @@ program test_set_outermost_contour
call set_outermost_contour(e, ind_targ_extr, innermost_level, extr_map, &
ssh, corner = [longitude(1), latitude(1)], step = step)
- print *, "e%outermost_contour%closed = ", e%outermost_contour%closed
- print *, "Radius of disk of equal area: ", &
- sqrt(e%outermost_contour%area / pi) / 1e3, "km"
-
- call shp_create_03("test_set_outermost_contour", 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.)
- call dbf_write_attribute_03(shphandle, shape_number, ifield = 0, &
- fieldvalue = e%outermost_contour%ssh)
- CALL shpclose(shphandle)
- print *, 'Created shapefile "test_set_outermost_contour".'
+ if (e%outermost_contour%closed) then
+ print *, "Radius of disk of equal area: ", &
+ sqrt(e%outermost_contour%area / pi) / 1e3, "km"
+
+ call shp_create_03("test_set_outermost_contour", 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.)
+ call dbf_write_attribute_03(shphandle, shape_number, ifield = 0, &
+ fieldvalue = e%outermost_contour%ssh)
+ CALL shpclose(shphandle)
+ print *, 'Created shapefile "test_set_outermost_contour".'
+ else
+ print *, "Could not find an outermost contour."
+ end if
end program test_set_outermost_contour
diff --git a/Tests/tests.json b/Tests/tests.json
index 98bf5e66..7fd07ffd 100644
--- a/Tests/tests.json
+++ b/Tests/tests.json
@@ -1,21 +1,26 @@
[
{
- "stdin" : "$stdin_dir/good_contour_in.txt",
+ "stdin" : "$stdin_dir/good_contour_nml_1.txt",
"args" : ["$compil_prod_dir/test_good_contour",
"$input_dir/example.nc"],
- "title" : "Good_contour"
+ "title" : "Good_contour",
+ "required": [["$stdin_dir/outside_points_1.csv", "outside_points.csv"]],
+ "description": "3 contours at that level. 2 contain the inside point, one of which contains the 2 outside points."
},
{
- "stdin" : "$stdin_dir/good_contour_in2.txt",
+ "stdin" : "$stdin_dir/good_contour_nml_2.txt",
"args" : ["$compil_prod_dir/test_good_contour",
"$input_dir/example.nc"],
- "title" : "Good_contour_2"
+ "title" : "Good_contour_2",
+ "required": [["$stdin_dir/outside_points_2.csv", "outside_points.csv"]],
+ "description": "Select another good contour."
},
{
- "stdin" : "$stdin_dir/no_good_contour_in.txt",
+ "stdin" : "$stdin_dir/no_good_contour_nml.txt",
"args" : ["$compil_prod_dir/test_good_contour",
"$input_dir/example.nc"],
- "title" : "No_good_contour"
+ "title" : "No_good_contour",
+ "required": [["$stdin_dir/outside_points_1.csv", "outside_points.csv"]]
},
{
"args" : ["$compil_prod_dir/test_local_extrema",
--
GitLab