In main units test_get_snapshot and test_set_all_outerm, read last

longitude and latitude instead of second longitude and latitude to
compute step. This divides the rounding error on step by a factor nlon
or nlat.

Tests/test_get_snapshot.f

@@ -30,10 +30,7 @@ program test_get_snapshot
   ! in m
 
   integer ncid, dimid, varid, nlon, nlat
-
-  ! Just the first two values to get the corner and step:
-  real longitude(2), latitude(2) ! in degrees
-
+  real lon_min, lon_max, lat_min, lat_max ! longitudes and latitudes, in degrees
   real step(2) ! longitude and latitude steps, in rad
   logical periodic ! grid is periodic in longitude
 
@@ -54,15 +51,20 @@ program test_get_snapshot
 
   call find_coord(ncid, dimid = dimid, varid = varid, std_name = "longitude")
   call nf95_inquire_dimension(ncid, dimid, nclen = nlon)
-  call nf95_get_var(ncid, varid, longitude)
+  call nf95_get_var(ncid, varid, lon_min)
+  call nf95_get_var(ncid, varid, lon_max, start = [nlon])
 
   call find_coord(ncid, dimid = dimid, varid = varid, std_name = "latitude")
   call nf95_inquire_dimension(ncid, dimid, nclen = nlat)
-  call nf95_get_var(ncid, varid, latitude)
+  call nf95_get_var(ncid, varid, lat_min)
+  call nf95_get_var(ncid, varid, lat_max, start = [nlat])
 
   call nf95_close(ncid)
   
-  step = [longitude(2) - longitude(1), latitude(2) - latitude(1)]  * deg_to_rad
+  call assert(lon_max > lon_min .and. lat_max > lat_min, &
+       "test_get_snapshot coordinate order")
+  step = [(lon_max - lon_min) / (nlon - 1), (lat_max - lat_min) / (nlat - 1)]  &
+       * deg_to_rad
 
   ! As we are requiring the grid spacing to be uniform, the value of
   ! "periodic" must be consistent with the values of step(1) and nlon:
@@ -72,7 +74,7 @@ program test_get_snapshot
 
   call get_snapshot(s, min_amp, m = 1, n_proc = 1, k_end = 1, max_delta = 4, &
        nlon = nlon, nlat = nlat, k = 1, max_radius = max_radius, &
-       corner = [longitude(1), latitude(1)] * deg_to_rad, step = step, &
+       corner = [lon_min, lat_min] * deg_to_rad, step = step, &
        periodic = periodic)
 
   call init_shapefiles(hshp_extremum, hshp_outermost, hshp_max_speed)

Tests/test_set_all_outerm.f

@@ -24,11 +24,8 @@ program test_set_all_outerm
   ! in m
 
   character(len = :), allocatable:: filename
-
   integer nlon, nlat
-
-  ! Just the first two values to get the corner and step:
-  real longitude(2), latitude(2) ! in degrees
+  real lon_min, lon_max, lat_min, lat_max ! longitudes and latitudes, in degrees
 
   integer, parameter:: max_radius(2) = [20, 20]
   ! maximum radius of an eddy in longitude and latitude, in number of
@@ -65,11 +62,13 @@ program test_set_all_outerm
 
   call find_coord(ncid, dimid = dimid, varid = varid, std_name = "longitude")
   call nf95_inquire_dimension(ncid, dimid, nclen = nlon)
-  call nf95_get_var(ncid, varid, longitude)
+  call nf95_get_var(ncid, varid, lon_min)
+  call nf95_get_var(ncid, varid, lon_max, start = [nlon])
 
   call find_coord(ncid, dimid = dimid, varid = varid, std_name = "latitude")
   call nf95_inquire_dimension(ncid, dimid, nclen = nlat)
-  call nf95_get_var(ncid, varid, latitude)
+  call nf95_get_var(ncid, varid, lat_min)
+  call nf95_get_var(ncid, varid, lat_max, start = [nlat])
 
   allocate(ssh(nlon, nlat))
   call nf95_inq_varid(ncid, "adt", varid)
@@ -79,11 +78,14 @@ program test_set_all_outerm
 
   call nf95_close(ncid)
 
-  step = [longitude(2) - longitude(1), latitude(2) - latitude(1)]  * deg_to_rad
+  call assert(lon_max > lon_min .and. lat_max > lat_min, &
+       "test_set_all_outerm coordinate order")
+  step = [(lon_max - lon_min) / (nlon - 1), (lat_max - lat_min) / (nlat - 1)] &
+       * deg_to_rad
   call assert(2 * max_radius(1) * step(1) < pi, &
        "test_set_all_outerm max_radius")
   call set_all_outerm(s, min_amp, max_radius, &
-       corner = [longitude(1), latitude(1)] * deg_to_rad, step = step, &
+       corner = [lon_min, lat_min] * deg_to_rad, step = step, &
        periodic = .false., ssh = ssh)
 
   call shp_create_03("extremum_1", shpt_point, hshp_extremum)

Tests/test_spherical_polygon_area.f

@@ -2,12 +2,14 @@ program test_spherical_polygon_area
 
   use, intrinsic:: ISO_C_BINDING
 
+  ! Libraries:
   use gpc_f, only: polygon, shpobj2pol
   use jumble, only: get_command_arg_dyn
-  use nr_util, only: pi
+  use nr_util, only: deg_to_rad
   use shapelib, only: shpfileobject, shpobject, shpclose, shpdestroyobject, &
        shpgetinfo
   use shapelib_03, only: shp_open_03, shp_read_object_03
+  
   use spherical_polygon_area_m, only: spherical_polygon_area
 
   implicit none
@@ -17,7 +19,6 @@ program test_spherical_polygon_area
   TYPE(shpobject) psobject
   type(POLYGON), allocatable:: p(:)
   type(POLYGON) merged_p
-  real, parameter:: deg2rad = pi / 180.
   integer ishape, j, nentities, shapetype, dbffieldcount, dbfrecordcount
   real(c_double) minbound(4), maxbound(4)
 
@@ -52,7 +53,7 @@ program test_spherical_polygon_area
   end do
 
   forall (j = 1:merged_p%nparts) &
-       merged_p%part(j)%points = merged_p%part(j)%points * deg2rad
+       merged_p%part(j)%points = merged_p%part(j)%points * deg_to_rad
   print *, "Area = ", spherical_polygon_area(merged_p) / 1e6, "km2"
 
 end program test_spherical_polygon_area

Tests/test_spherical_polyline_area.f

 program test_spherical_polyline_area
 
+  ! Libraries:
   use contour_531, only: polyline
-  use nr_util, only: pi
+  use nr_util, only: deg_to_rad
+  
   use spherical_polyline_area_m, only: spherical_polyline_area
 
   implicit none
 
   real points(2, 4)
-  real, parameter:: deg2rad = pi / 180.
 
   !-------------------------------------------------------------------------
 
   points(:, 1) = [0., 0.]
-  points(:, 2) = [30. * deg2rad, 0.]
-  points(:, 3) = [0., 30. * deg2rad]
+  points(:, 2) = [30. * deg_to_rad, 0.]
+  points(:, 3) = [0., 30. * deg_to_rad]
   points(:, 4) = points(:, 1)
   
   print *, "Area = ", spherical_polyline_area(polyline(n_points = 4, &

Tests/test_successive_overlap.f

@@ -2,12 +2,14 @@ program test_successive_overlap
 
   use, intrinsic:: ISO_FORTRAN_ENV
 
-  use derived_types, only: snapshot
+  ! Libraries:
   use jumble, only: new_unit
   use nr_util, only: deg_to_rad
-  use read_snapshot_m, only: read_snapshot
   use shapelib, only: shpfileobject, shpclose
   use shapelib_03, only: shp_open_03
+
+  use derived_types, only: snapshot
+  use read_snapshot_m, only: read_snapshot
   use successive_overlap_m, only: successive_overlap
 
   implicit none

Tests/test_weight.f

 program test_weight
 
+  ! Libraries:
+  use nr_util, only: pi, deg_to_rad
+  
   use derived_types, only: eddy
-  use nr_util, only: pi
   use weight_m, only: weight
 
   implicit none
 
-  real, parameter:: deg_over_rad = pi / 180.
   type(eddy) e1, e2
 
   !---------------------------------------------------------------------
@@ -14,12 +15,12 @@ program test_weight
   e1%max_speed = 0.15
   e1%speed_cont%n_points = 4
   e1%speed_cont%area = pi * 5e4**2
-  e1%coord_extr(2) = - 20. * deg_over_rad
+  e1%coord_extr(2) = - 20. * deg_to_rad
 
   e2%max_speed = 0.1
   e2%speed_cont%n_points = 4
   e2%speed_cont%area = pi * 6e4**2
-  e2%coord_extr(2) = - 25. * deg_over_rad
+  e2%coord_extr(2) = - 25. * deg_to_rad
   
   print *, "weight = ", weight(e1, e2)