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

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

Tests/read_eddy.f

+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

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
 

Tests/test_read_eddy.f

+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

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

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
   

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"
     }
 ]

depend.mk

 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 

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