From bd45e4cf48ad86552b690ad82a8d11db66ca86e8 Mon Sep 17 00:00:00 2001
From: Lionel GUEZ <guez@lmd.ens.fr>
Date: Thu, 19 Apr 2018 16:45:14 +0200
Subject: [PATCH] Add procedure read_eddy and corresponding test. This is in
preparation for upcoming procedure read_snapshot.
Bug fix in procedure rossby_number: e%speed_cont could be null.
---
GNUmakefile | 22 ++----
Tests/read_eddy.f | 111 ++++++++++++++++++++++++++++
Tests/test_get_snapshot.f | 2 +-
Tests/test_read_eddy.f | 40 ++++++++++
Tests/test_spherical_polygon_area.f | 22 +++---
Tests/test_weight.f | 2 +
Tests/tests.json | 9 ++-
depend.mk | 2 +
weight.f | 16 +++-
9 files changed, 193 insertions(+), 33 deletions(-)
create mode 100644 Tests/read_eddy.f
create mode 100644 Tests/test_read_eddy.f
diff --git a/GNUmakefile b/GNUmakefile
index 75ad5b7f..492d9ee0 100644
--- a/GNUmakefile
+++ b/GNUmakefile
@@ -19,7 +19,9 @@ src_test_weight = test_weight.f weight.f derived_types.f
src_test_spherical_polygon_area = spherical_polygon_area.f test_spherical_polygon_area.f spherical_polyline_area.f
-sources := $(sort ${src_test_local_extrema} ${src_test_get_1_outerm} ${src_test_set_max_speed} ${src_test_get_snapshot} ${src_test_set_all_outerm} ${src_test_weight} ${src_test_spherical_polygon_area}) test_good_contour.f test_inside_4.f test_max_speed_contour_ssh.f test_mean_speed.f test_spherical_polyline_area.f
+src_test_read_eddy = test_read_eddy.f derived_types.f init_shapefiles.f read_eddy.f write_eddy.f
+
+sources := $(sort ${src_test_local_extrema} ${src_test_get_1_outerm} ${src_test_set_max_speed} ${src_test_get_snapshot} ${src_test_set_all_outerm} ${src_test_weight} ${src_test_spherical_polygon_area} ${src_test_read_eddy}) test_good_contour.f test_inside_4.f test_max_speed_contour_ssh.f test_mean_speed.f test_spherical_polyline_area.f
lib_list = GPC_F contour_531 numer_rec_95 shapelib_03 netcdf95 geometry jumble netcdff fortrangis shp fortranc nr_util
@@ -32,15 +34,15 @@ obj_test_get_snapshot := $(src_test_get_snapshot:.f=.o)
obj_test_set_all_outerm := $(src_test_set_all_outerm:.f=.o)
obj_test_weight := $(src_test_weight:.f=.o)
obj_test_spherical_polygon_area := $(src_test_spherical_polygon_area:.f=.o)
+obj_test_read_eddy := $(src_test_read_eddy:.f=.o)
objects := $(sources:.f=.o)
-execut = test_good_contour test_inside_4 test_get_1_outerm test_local_extrema test_max_speed_contour_ssh test_mean_speed test_set_max_speed test_get_snapshot test_set_all_outerm test_weight test_spherical_polyline_area test_spherical_polygon_area
+execut = test_good_contour test_inside_4 test_get_1_outerm test_local_extrema test_max_speed_contour_ssh test_mean_speed test_set_max_speed test_get_snapshot test_set_all_outerm test_weight test_spherical_polyline_area test_spherical_polygon_area test_read_eddy
# 3. Compiler-dependent part
mode = debug
include ${general_compiler_options_dir}/${FC}_${mode}.mk
-comma=,
ifeq (${FC},gfortran)
# gfortran bug:
@@ -49,9 +51,6 @@ endif
# 4. Rules
-SHELL = bash
-LINK.o = $(FC) $(LDFLAGS) $(TARGET_ARCH)
-
%.o: %.f
@echo "Building $@..."
@$(COMPILE.f) $(OUTPUT_OPTION) $<
@@ -60,9 +59,9 @@ LINK.o = $(FC) $(LDFLAGS) $(TARGET_ARCH)
@echo "Linking $@..."
@$(LINK.o) $^ $(LOADLIBES) $(LDLIBS) -o $@
-.DELETE_ON_ERROR:
.PHONY: all clean clobber depend
all: ${execut} log
+include ${general_compiler_options_dir}/settings.mk
test_get_1_outerm: ${obj_test_get_1_outerm}
test_set_max_speed: ${obj_test_set_max_speed}
test_max_speed_contour_ssh: max_speed_contour_ssh.o
@@ -75,24 +74,17 @@ test_set_all_outerm: ${obj_test_set_all_outerm}
test_weight: ${obj_test_weight}
test_spherical_polyline_area: spherical_polyline_area.o
test_spherical_polygon_area: ${obj_test_spherical_polygon_area}
+test_read_eddy: ${obj_test_read_eddy}
depend ${makefile_dir}/depend.mk:
makedepf90 -free -Wmissing -Wconfused $(addprefix -I, ${VPATH}) -nosrc $(addprefix -u , ${lib_list} shapelib netcdf intrinsic) ${sources} >${makefile_dir}/depend.mk
-TAGS: ${sources}
- ctags -e --language-force=fortran $^
-
clean:
rm -f ${execut} ${objects} log
clobber: clean
rm -f *.mod ${makefile_dir}/depend.mk TAGS
-log:
- hostname >$@
- ${FC} ${version_flag} >>$@ 2>&1
- @echo -e "\nFC = ${FC}\n\nCPPFLAGS = ${CPPFLAGS}\n\nFFLAGS = ${FFLAGS}\n\nLDLIBS = ${LDLIBS}\n\nLDFLAGS = ${LDFLAGS}" >>$@
-
ifneq ($(MAKECMDGOALS), clobber)
include ${makefile_dir}/depend.mk
endif
diff --git a/Tests/read_eddy.f b/Tests/read_eddy.f
new file mode 100644
index 00000000..f15327c8
--- /dev/null
+++ b/Tests/read_eddy.f
@@ -0,0 +1,111 @@
+module read_eddy_m
+
+ implicit none
+
+contains
+
+ subroutine read_eddy(e, k, i, hshp_extremum, hshp_outermost, hshp_max_speed, &
+ ishape)
+
+ use, intrinsic:: ISO_C_BINDING
+
+ use contour_531, only: null_polyline
+ use derived_types, only: eddy
+ use gpc_f, only: shpobj2pol, polygon
+ use nr_util, only: pi
+ use shapelib, only: shpfileobject, shpobject, dbfreadattribute, &
+ shpdestroyobject
+ use shapelib_03, only: shp_read_object_03
+
+ type(eddy), intent(out):: e
+ integer, intent(out):: k ! date index
+ integer, intent(out):: i ! eddy index
+
+ TYPE(shpfileobject), intent(inout):: hshp_extremum ! shapefile extremum_$m
+
+ TYPE(shpfileobject), intent(inout):: hshp_outermost
+ ! shapefile outermost_contour_$m
+
+ TYPE(shpfileobject), intent(inout):: hshp_max_speed
+ ! shapefile x_speed_contour_$m
+
+ integer, intent(in):: ishape
+
+ ! Local:
+ integer interpolated, cyclone, suff_amp, twice
+ TYPE(shpobject) psobject
+ type(polygon) p
+
+ real(c_double) attr
+ ! We need this intermediary variable because the kind of the
+ ! argument to dbfreadattribute cannot be the default kind
+
+ real, parameter:: deg_to_rad = pi / 180.
+
+ !---------------------------------------------------------------------
+
+ call dbfreadattribute(hshp_extremum, ishape, ifield = 0, attr = attr)
+ e%ssh_extr = attr
+
+ call dbfreadattribute(hshp_extremum, ishape, ifield = 1, attr = k)
+ call dbfreadattribute(hshp_extremum, ishape, ifield = 2, attr = i)
+
+ call dbfreadattribute(hshp_extremum, ishape, ifield = 3, &
+ attr = interpolated)
+ e%interpolated = .false.
+
+ call dbfreadattribute(hshp_extremum, ishape, ifield = 4, attr = cyclone)
+ e%cyclone = cyclone == 1
+
+ call dbfreadattribute(hshp_extremum, ishape, ifield = 5, attr = suff_amp)
+ e%suff_amp = suff_amp == 1
+
+ call dbfreadattribute(hshp_extremum, ishape, ifield = 6, attr = attr)
+ e%max_speed = attr
+
+ call shp_read_object_03(hshp_extremum, ishape, psobject)
+ e%coord_extr = [psobject%padfx(1), psobject%padfy(1)] * deg_to_rad
+ call shpdestroyobject(psobject)
+
+ call dbfreadattribute(hshp_outermost, ishape, ifield = 0, attr = attr)
+ e%out_cont%area = attr
+
+ call dbfreadattribute(hshp_outermost, ishape, ifield = 1, attr = attr)
+ e%out_cont%ssh = attr
+
+ call dbfreadattribute(hshp_outermost, ishape, ifield = 4, attr = twice)
+ e%twice = twice == 1
+
+ call dbfreadattribute(hshp_outermost, ishape, ifield = 5, attr = e%radius4)
+
+ call shp_read_object_03(hshp_outermost, ishape, psobject)
+ p = shpobj2pol(psobject)
+ call shpdestroyobject(psobject)
+
+ if (p%nparts == 0) then
+ e%out_cont%polyline = null_polyline()
+ else
+ e%out_cont%polyline = p%part(1)
+ e%out_cont%polyline%points = e%out_cont%polyline%points * deg_to_rad
+ end if
+
+ call dbfreadattribute(hshp_max_speed, ishape, ifield = 0, attr = attr)
+ e%speed_cont%area = attr
+
+ call dbfreadattribute(hshp_max_speed, ishape, ifield = 1, attr = attr)
+ e%speed_cont%ssh = attr
+
+ call shp_read_object_03(hshp_max_speed, ishape, psobject)
+ p = shpobj2pol(psobject)
+ call shpdestroyobject(psobject)
+
+ if (p%nparts == 0) then
+ e%speed_cont%polyline = null_polyline()
+ else
+ e%speed_cont%polyline = p%part(1)
+ e%speed_cont%polyline%points = e%speed_cont%polyline%points * deg_to_rad
+ end if
+
+ end subroutine read_eddy
+
+end module read_eddy_m
diff --git a/Tests/test_get_snapshot.f b/Tests/test_get_snapshot.f
index 36253bcb..5ad99eec 100644
--- a/Tests/test_get_snapshot.f
+++ b/Tests/test_get_snapshot.f
@@ -21,7 +21,7 @@ program test_get_snapshot
TYPE(shpfileobject) hshp_extremum ! shapefile extremum_$m
TYPE(shpfileobject) hshp_outermost ! shapefile outermost_contour_$m
- TYPE(shpfileobject) hshp_max_speed ! shapefile x_speed_contour_$m
+ TYPE(shpfileobject) hshp_max_speed ! shapefile max_speed_contour_$m
integer unit_isolated, unit_number_eddies
diff --git a/Tests/test_read_eddy.f b/Tests/test_read_eddy.f
new file mode 100644
index 00000000..73c48b42
--- /dev/null
+++ b/Tests/test_read_eddy.f
@@ -0,0 +1,40 @@
+program test_read_eddy
+
+ use derived_types, only: eddy
+ use init_shapefiles_m, only: init_shapefiles
+ use read_eddy_m, only: read_eddy
+ use shapelib, only: shpfileobject, shpclose
+ use shapelib_03, only: shp_open_03
+ use write_eddy_m, only: write_eddy
+
+ implicit none
+
+ type(eddy) e
+ integer k, i
+ TYPE(shpfileobject) hshp_extremum ! shapefile extremum_1
+ TYPE(shpfileobject) hshp_outermost ! shapefile outermost_contour_1
+ TYPE(shpfileobject) hshp_max_speed ! shapefile max_speed_contour_1
+
+ !-------------------------------------------------------------------------
+
+ call shp_open_03("extremum_1_old", pszaccess = "rb", hshp = hshp_extremum)
+ call shp_open_03("outermost_contour_1_old", pszaccess = "rb", &
+ hshp = hshp_outermost)
+ call shp_open_03("max_speed_contour_1_old", pszaccess = "rb", &
+ hshp = hshp_max_speed)
+ call read_eddy(e, k, i, hshp_extremum, hshp_outermost, hshp_max_speed, &
+ ishape = 0)
+ CALL shpclose(hshp_extremum)
+ CALL shpclose(hshp_outermost)
+ CALL shpclose(hshp_max_speed)
+
+ call init_shapefiles(hshp_extremum, hshp_outermost, hshp_max_speed)
+ call write_eddy(e, k, i, hshp_extremum, hshp_outermost, hshp_max_speed)
+ CALL shpclose(hshp_extremum)
+ print *, 'Created shapefile "extremum_1".'
+ CALL shpclose(hshp_outermost)
+ print *, 'Created shapefile "outermost_contour_1".'
+ CALL shpclose(hshp_max_speed)
+ print *, 'Created shapefile "max_speed_contour_1".'
+
+end program test_read_eddy
diff --git a/Tests/test_spherical_polygon_area.f b/Tests/test_spherical_polygon_area.f
index 07dc1272..c093831b 100644
--- a/Tests/test_spherical_polygon_area.f
+++ b/Tests/test_spherical_polygon_area.f
@@ -18,7 +18,7 @@ program test_spherical_polygon_area
type(POLYGON), allocatable:: p(:)
type(POLYGON) merged_p
real, parameter:: deg2rad = pi / 180.
- integer i, j, nentities, shapetype, dbffieldcount, dbfrecordcount
+ integer ishape, j, nentities, shapetype, dbffieldcount, dbfrecordcount
real(c_double) minbound(4), maxbound(4)
!---------------------------------------------------------------------
@@ -27,11 +27,11 @@ program test_spherical_polygon_area
call shp_open_03(filename, "rb", hshp)
call shpgetinfo(hshp, nentities, shapetype, minbound, maxbound, &
dbffieldcount, dbfrecordcount)
- allocate(p(nentities))
+ allocate(p(0:nentities - 1))
- do i = 1, nentities
- call shp_read_object_03(hshp, i - 1, psobject)
- p(i) = shpobj2pol(psobject)
+ do ishape = 0, nentities - 1
+ call shp_read_object_03(hshp, ishape, psobject)
+ p(ishape) = shpobj2pol(psobject)
call shpdestroyobject(psobject)
end do
@@ -40,19 +40,19 @@ program test_spherical_polygon_area
allocate(merged_p%part(merged_p%nparts), merged_p%hole(merged_p%nparts))
j = 1
- do i = 1, nentities
- merged_p%part(j:j + p(i)%nparts - 1) = p(i)%part
+ do ishape = 0, nentities - 1
+ merged_p%part(j:j + p(ishape)%nparts - 1) = p(ishape)%part
! Assume that the first part of each shape is the exterior ring,
! and the following parts are holes (GeoJSON convention):
merged_p%hole(j) = .false.
- merged_p%hole(j + 1:j + p(i)%nparts - 1) = .true.
+ merged_p%hole(j + 1:j + p(ishape)%nparts - 1) = .true.
- j = j + p(i)%nparts
+ j = j + p(ishape)%nparts
end do
- forall (i = 1:merged_p%nparts) &
- merged_p%part(i)%points = merged_p%part(i)%points * deg2rad
+ forall (j = 1:merged_p%nparts) &
+ merged_p%part(j)%points = merged_p%part(j)%points * deg2rad
print *, "Area = ", spherical_polygon_area(merged_p) / 1e6, "km2"
end program test_spherical_polygon_area
diff --git a/Tests/test_weight.f b/Tests/test_weight.f
index d3ed9005..c2ef9914 100644
--- a/Tests/test_weight.f
+++ b/Tests/test_weight.f
@@ -12,10 +12,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
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
diff --git a/Tests/tests.json b/Tests/tests.json
index 76907d48..1b90bbe8 100644
--- a/Tests/tests.json
+++ b/Tests/tests.json
@@ -105,9 +105,7 @@
"required": [["$input_dir/2006_01/h_region_1.nc", "h.nc"],
["$input_dir/2006_01/uv_region_1.nc", "uv.nc"],
["$input_dir/extr_map_negative_2_8.nc", "extr_map.nc"],
- "$tests_old_dir/Get_1_outerm_noise_2_8/test_get_1_outerm.shp",
- "$tests_old_dir/Get_1_outerm_noise_2_8/test_get_1_outerm.shx",
- "$tests_old_dir/Get_1_outerm_noise_2_8/test_get_1_outerm.dbf"],
+ "$tests_old_dir/Get_1_outerm_noise_2_8/test_get_1_outerm.*"],
"input": "&MAIN_NML IND_TARG_EXTR=19,11 /\n"
},
{
@@ -204,5 +202,10 @@
"$input_dir/tri_square_hole"],
"title": "Area_multi_polygon",
"description": "Area of a multipolygon, with a hole in one polygon."
+ },
+ {
+ "args": "$compil_prod_dir/test_read_eddy",
+ "required": ["$input_dir/Test_read_eddy/*"],
+ "title": "Read_eddy"
}
]
diff --git a/depend.mk b/depend.mk
index b711aad2..ff7be99d 100644
--- a/depend.mk
+++ b/depend.mk
@@ -1,6 +1,7 @@
dispatch_snapshot.o : write_eddy.o send_snapshot.o derived_types.o
get_1_outerm.o : spherical_polyline_area.o outermost_possible_level.o good_contour.o derived_types.o
get_snapshot.o : set_all_outerm.o set_max_speed.o receive_snapshot.o derived_types.o
+read_eddy.o : derived_types.o
receive_snapshot.o : derived_types.o
send_snapshot.o : derived_types.o
set_all_outerm.o : local_extrema.o get_1_outerm.o derived_types.o
@@ -9,6 +10,7 @@ spherical_polygon_area.o : spherical_polyline_area.o
test_get_1_outerm.o : get_1_outerm.o derived_types.o
test_get_snapshot.o : init_shapefiles.o get_snapshot.o dispatch_snapshot.o derived_types.o
test_local_extrema.o : local_extrema.o
+test_read_eddy.o : write_eddy.o read_eddy.o init_shapefiles.o derived_types.o
test_set_all_outerm.o : set_all_outerm.o derived_types.o
test_set_max_speed.o : set_max_speed.o derived_types.o
test_spherical_polygon_area.o : spherical_polygon_area.o
diff --git a/weight.f b/weight.f
index 88dc0b2e..a4b20f16 100644
--- a/weight.f
+++ b/weight.f
@@ -24,7 +24,9 @@ contains
pure real function rossby_number(e)
- ! Computes the Rossby number for the maximum-speed contour.
+ ! Computes the Rossby number for the maximum-speed
+ ! contour. Assumes that e%out_cont is not null. e%max_speed can be
+ ! a quiet NaN, and then the result is a quiet NaN.
use derived_types, only: eddy
use nr_util, only: pi
@@ -35,11 +37,19 @@ contains
real, parameter:: omega = 2 * pi / 86164.
! angular speed of rotation of the Earth, in rad s-1
+
+ real radius ! of the disk of same area
!----------------------------------------------------------
- rossby_number = abs(e%max_speed / (sqrt(e%speed_cont%area / pi) &
- * 2. * omega * sin(e%coord_extr(2))))
+ if (e%speed_cont%n_points /= 0) then
+ radius = sqrt(e%speed_cont%area / pi)
+ else
+ radius = sqrt(e%out_cont%area / pi)
+ end if
+
+ rossby_number = abs(e%max_speed &
+ / (radius * 2. * omega * sin(e%coord_extr(2))))
end function rossby_number
--
GitLab