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

@@ -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
 

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

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

Tests/Stdin/good_contour_nml_1.txt

+&MAIN_NML corner = -2., -2., step = 0.08, 0.08 /

Tests/Stdin/good_contour_in2.txt → Tests/Stdin/good_contour_nml_2.txt

 ! 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/

Tests/Stdin/no_good_contour_in.txt → Tests/Stdin/no_good_contour_nml.txt

 ! Case where there is no good contour
 &MAIN_NML INSIDE_POINT=0.,0. /
-2
--0.5, -1.
-0.3, 0.

Tests/Stdin/good_contour_in.txt → Tests/Stdin/outside_points_1.csv

-&MAIN_NML /
-2
 -0.5, -1.
 0.3, 0.

Tests/Stdin/outside_points_2.csv

+-0.5, -1.
+-0.7, 0.

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)

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

Tests/tests.json

 [
     {
-	"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",