From 92a6f46683a37f2156e052aae69f459aa22f867d Mon Sep 17 00:00:00 2001 From: Lionel GUEZ <guez@lmd.ens.fr> Date: Fri, 28 Sep 2018 17:10:26 +0200 Subject: [PATCH] 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 | 18 ++++++++++-------- Tests/test_set_all_outerm.f | 18 ++++++++++-------- Tests/test_spherical_polygon_area.f | 7 ++++--- Tests/test_spherical_polyline_area.f | 9 +++++---- Tests/test_successive_overlap.f | 6 ++++-- Tests/test_weight.f | 9 +++++---- 6 files changed, 38 insertions(+), 29 deletions(-) diff --git a/Tests/test_get_snapshot.f b/Tests/test_get_snapshot.f index ec6f7554..28c0200c 100644 --- a/Tests/test_get_snapshot.f +++ b/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) diff --git a/Tests/test_set_all_outerm.f b/Tests/test_set_all_outerm.f index 0e8bec64..4789fc25 100644 --- a/Tests/test_set_all_outerm.f +++ b/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) diff --git a/Tests/test_spherical_polygon_area.f b/Tests/test_spherical_polygon_area.f index c093831b..3c338753 100644 --- a/Tests/test_spherical_polygon_area.f +++ b/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 diff --git a/Tests/test_spherical_polyline_area.f b/Tests/test_spherical_polyline_area.f index b77e28eb..e53649a8 100644 --- a/Tests/test_spherical_polyline_area.f +++ b/Tests/test_spherical_polyline_area.f @@ -1,19 +1,20 @@ 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, & diff --git a/Tests/test_successive_overlap.f b/Tests/test_successive_overlap.f index 53bd292b..add638bf 100644 --- a/Tests/test_successive_overlap.f +++ b/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 diff --git a/Tests/test_weight.f b/Tests/test_weight.f index c2ef9914..1e8123ce 100644 --- a/Tests/test_weight.f +++ b/Tests/test_weight.f @@ -1,12 +1,13 @@ 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) -- GitLab