Move computation of extrema from procedure get_snapshot to procedure

set_all_outerm. Thus, we have the clearer intent(out) for argument s
of set_all_outerm. Also, this prepares a possible replacement of
set_all_outerm by a procedure from the Topology toolkit.

Create test program for procedure set_all_outerm.

GNUmakefile

@@ -13,6 +13,8 @@ src_test_set_max_speed = test_set_max_speed.f derived_types.f set_max_speed.f go
 
 src_test_get_snapshot = test_get_snapshot.f get_snapshot.f dispatch_snapshot.f write_eddy.f send_snapshot.f receive_snapshot.f local_extrema.f set_max_speed.f outermost_possible_level.f get_1_outerm.f max_speed_contour_ssh.f good_contour.f spherical_polygon_area.f mean_speed.f inside_4.f set_all_outerm.f
 
+src_test_set_all_outerm = test_set_all_outerm.f derived_types.f dispatch_snapshot.f set_all_outerm.f write_eddy.f send_snapshot.f local_extrema.f get_1_outerm.f good_contour.f spherical_polygon_area.f
+
 sources := $(shell cat ${makefile_dir}/file_list)
 
 lib_list = contour_531 numer_rec_95 shapelib_03 netcdf95 geometry jumble netcdff fortrangis shp fortranc nr_util
@@ -23,9 +25,10 @@ obj_test_local_extrema := $(src_test_local_extrema:.f=.o)
 obj_test_get_1_outerm := $(src_test_get_1_outerm:.f=.o)
 obj_test_set_max_speed := $(src_test_set_max_speed:.f=.o)
 obj_test_get_snapshot := $(src_test_get_snapshot:.f=.o)
+obj_test_set_all_outerm := $(src_test_set_all_outerm:.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
+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
 
 # 3. Compiler-dependent part
 
@@ -47,6 +50,7 @@ test_inside_4: inside_4.o
 test_local_extrema: ${obj_test_local_extrema}
 test_mean_speed: mean_speed.o
 test_get_snapshot: ${obj_test_get_snapshot}
+test_set_all_outerm: ${obj_test_set_all_outerm}
 
 depend ${makefile_dir}/depend.mk:
 	makedepf90 -free -Wmissing -Wconfused $(addprefix -I, ${VPATH}) -nosrc $(addprefix -u , ${lib_list} shapelib netcdf) ${sources} >${makefile_dir}/depend.mk

Sources/get_snapshot.f

@@ -9,7 +9,6 @@ contains
 
     use contour_531, only: convert_to_ind, null_polyline
     use derived_types, only: snapshot
-    use local_extrema_m, only: local_extrema
     use netcdf, only: nf90_nowrite
     use netcdf95, only: nf95_open, nf95_close, nf95_inq_varid, nf95_get_var
     use receive_snapshot_m, only: receive_snapshot
@@ -43,10 +42,6 @@ contains
     real ssh(nlon, nlat) ! sea-surface height, in m
     real u(nlon, nlat), v(nlon, nlat) ! wind, in m s-1
 
-    real, allocatable:: innermost_level(:) ! (s%number_vis_eddies)
-    ! level of innermost contour, for each extremum
-
-    logical, allocatable:: cyclone(:) ! (s%number_vis_eddies)
     integer i
 
     ! Window around each extremum:
@@ -70,22 +65,7 @@ contains
        call nf95_get_var(ncid, varid, v, start = [1, 1, k])
        call nf95_close(ncid)
 
-       ! Extraction of eddies:
-
-       allocate(s%extr_map(nlon, nlat))
-       call local_extrema(ssh, s%extr_map, s%ind_extr, innermost_level, cyclone)
-       s%number_vis_eddies = size(s%ind_extr, 2)
-       allocate(s%list_vis(s%number_vis_eddies))
-
-       forall (i = 1:s%number_vis_eddies)
-          s%list_vis(i)%coord_extr = corner + (s%ind_extr(:, i) - 1) * step
-          s%list_vis(i)%ssh_extremum = ssh(s%ind_extr(1, i), s%ind_extr(2, i))
-          s%list_vis(i)%cyclone = cyclone(i)
-          s%list_vis(i)%interpolated = .false.
-       end forall
-
-       call set_all_outerm(s, min_amp, max_radius, corner, step, &
-            innermost_level, ssh)
+       call set_all_outerm(s, min_amp, max_radius, corner, step, ssh)
 
        ! Done with outermost contours, now let us take care of the
        ! max-speed contours.

Sources/set_all_outerm.f

@@ -4,24 +4,20 @@ module set_all_outerm_m
 
 contains
 
-  subroutine set_all_outerm(s, min_amp, max_radius, corner, step, &
-       innermost_level, ssh)
+  subroutine set_all_outerm(s, min_amp, max_radius, corner, step, ssh)
 
-    ! Set all outermost contours.
+    ! Extraction of eddies: find extrema and set all outermost
+    ! contours in snapshot. Not a function because snapshot is not
+    ! completely defined.
 
     use derived_types, only: snapshot
     use get_1_outerm_m, only: get_1_outerm
+    use local_extrema_m, only: local_extrema
 
-    type(snapshot), intent(inout):: s
-    
-    ! Intent(in): s%number_vis_eddies, s%ind_extr,
-    ! s%list_vis%ssh_extremum, s%list_vis%cyclone,
-    ! s%list_vis%coord_extr
-
-    ! Intent(out): s%list_vis%outermost_contour, s%list_vis%suff_amp,
-    ! s%list_vis%twice
-
-    ! Intent(inout): s%extr_map
+    type(snapshot), intent(out):: s
+    ! Specifically: define everything in s except
+    ! s%list_vis%max_speed_contour, s%list_vis%max_speed and
+    ! s%number_eddies.
 
     real, intent(in):: min_amp ! minimum amplitude of ssh, between
     ! extremum and outermost contour, in m
@@ -34,29 +30,45 @@ contains
 
     real, intent(in):: step(:) ! (2) longitude and latitude steps, in rad
 
-    real, intent(in):: innermost_level(:) ! (s%number_vis_eddies)
-    ! level of innermost contour, for each extremum
-
     real, intent(in):: ssh(:, :) ! sea-surface height, in m
 
     ! Local:
 
+    real, allocatable:: innermost_level(:) ! (s%number_vis_eddies)
+    ! level of innermost contour, for each extremum
+
+    logical, allocatable:: cyclone(:) ! (s%number_vis_eddies)
     integer i
-    logical flat_extr(s%number_vis_eddies) ! flat extremum
-    logical noise_around(s%number_vis_eddies)
+    logical, allocatable:: flat_extr(:) ! (s%number_vis_eddies) flat extremum
+    logical, allocatable:: noise_around(:) ! (s%number_vis_eddies)
 
     ! Window around each extremum:
 
-    integer llc(2, s%number_vis_eddies) ! indices in global grid of
-    ! lower left corner
+    integer, allocatable:: llc(:, :) ! (2, s%number_vis_eddies)
+    ! indices in global grid of lower left corner
 
-    integer urc(2, s%number_vis_eddies) ! indices in global grid of
-    ! upper right corner
+    integer, allocatable:: urc(:, :) ! (2, s%number_vis_eddies)
+    ! indices in global grid of upper right corner
 
-    real corner_window(2, s%number_vis_eddies) ! longitude and latitude, in rad
+    real, allocatable:: corner_window(:, :) ! (2, s%number_vis_eddies)
+    ! longitude and latitude, in rad
 
     !--------------------------------------------------------------
 
+    allocate(s%extr_map(size(ssh, 1), size(ssh, 2)))
+    call local_extrema(ssh, s%extr_map, s%ind_extr, innermost_level, cyclone)
+    s%number_vis_eddies = size(s%ind_extr, 2)
+    allocate(s%list_vis(s%number_vis_eddies), flat_extr(s%number_vis_eddies), &
+         noise_around(s%number_vis_eddies), llc(2, s%number_vis_eddies), &
+         urc(2, s%number_vis_eddies), corner_window(2, s%number_vis_eddies))
+
+    forall (i = 1:s%number_vis_eddies)
+       s%list_vis(i)%coord_extr = corner + (s%ind_extr(:, i) - 1) * step
+       s%list_vis(i)%ssh_extremum = ssh(s%ind_extr(1, i), s%ind_extr(2, i))
+       s%list_vis(i)%cyclone = cyclone(i)
+       s%list_vis(i)%interpolated = .false.
+    end forall
+
     ! Define the geographical window around each eddy extremum:
     forall (i = 1:s%number_vis_eddies)
        llc(:, i) = max(s%ind_extr(:, i) - max_radius, 1)

Tests/test_set_all_outerm.f

+program test_set_all_outerm
+
+  use contour_531, only: null_polyline
+  use derived_types, only: snapshot
+  use dispatch_snapshot_m, only: dispatch_snapshot
+  use jumble, only: new_unit, get_command_arg_dyn
+  use netcdf, only: nf90_nowrite
+  use netcdf95, only: nf95_open, nf95_close, nf95_inq_varid, nf95_get_var, &
+       find_coord, nf95_inquire_dimension
+  use nr_util, only: pi
+  use set_all_outerm_m, only: set_all_outerm
+  use shapelib, only: shpt_point, shpt_polygon, shpfileobject, ftdouble, &
+       shpclose, ftinteger
+  use shapelib_03, only: shp_create_03, dbf_add_field_03
+
+  implicit none
+
+  type(snapshot) s
+
+  real:: min_amp = 0.
+  ! minimum amplitude of ssh, between extremum and outermost contour,
+  ! 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
+
+  integer, parameter:: max_radius(2) = [20, 20]
+  ! maximum radius of an eddy in longitude and latitude, in number of
+  ! grid points
+
+  real, parameter:: deg_over_rad = pi / 180.
+
+  integer ncid, varid, dimid, i
+  real, allocatable:: ssh(:, :) ! (nlon, nlat) sea-surface height, in m
+
+  TYPE(shpfileobject) hshp_extremum
+  ! shapefile extremum_$m. The fields in the DBF file are, in that
+  ! order: ssh, date index, eddy index, interpolated, cyclone, suff
+  ! amplitude.
+
+  TYPE(shpfileobject) hshp_outermost
+  ! shapefile outermost_contour_$m. The fields in the DBF file are,
+  ! in that order: area, ssh, date index, eddy index, twice.
+
+  TYPE(shpfileobject) hshp_max_speed
+  ! shapefile x_speed_contour_$m. The fields in the DBF file are, in
+  ! that order: area, ssh, spped, date index, eddy index.
+
+  integer ifield, unit_isolated, unit_number_eddies
+
+  !--------------------------------------------------------------
+
+  call get_command_arg_dyn(1, filename)
+
+  print *, "min_amp = ? "
+  read *, min_amp
+  print *, "min_amp = ", min_amp
+
+  ! Read ssh:
+
+  call nf95_open(filename, nf90_nowrite, ncid)
+
+  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 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)
+
+  allocate(ssh(nlon, nlat))
+  call nf95_inq_varid(ncid, "adt", varid)
+  call nf95_get_var(ncid, varid, ssh)
+
+  call nf95_close(ncid)
+
+  call set_all_outerm(s, min_amp, max_radius, &
+       corner = [longitude(1), latitude(1)] * deg_over_rad, step &
+       = [longitude(2) - longitude(1), latitude(2) - latitude(1)]  &
+       * deg_over_rad, ssh = ssh)
+
+  do i = 1, s%number_vis_eddies
+     s%list_vis(i)%max_speed_contour%polyline = null_polyline()
+     s%list_vis(i)%max_speed_contour%ssh = s%list_vis(i)%ssh_extremum
+     s%list_vis(i)%max_speed_contour%area = 0.
+     s%list_vis(i)%max_speed = 0.
+  end do
+
+  s%number_eddies = s%number_vis_eddies
+
+  call shp_create_03("extremum_1", shpt_point, hshp_extremum)
+  call dbf_add_field_03(ifield, hshp_extremum, 'ssh', ftdouble, nwidth = 13, &
+       ndecimals = 6)
+  call dbf_add_field_03(ifield, hshp_extremum, 'date_index', ftinteger, &
+       nwidth = 4, ndecimals = 0)
+  call dbf_add_field_03(ifield, hshp_extremum, 'eddy_index', ftinteger, &
+       nwidth = 5, ndecimals = 0)
+  call dbf_add_field_03(ifield, hshp_extremum, 'interpolat', ftinteger, &
+       nwidth = 1, ndecimals = 0)
+  call dbf_add_field_03(ifield, hshp_extremum, 'cyclone', ftinteger, &
+       nwidth = 1, ndecimals = 0)
+  call dbf_add_field_03(ifield, hshp_extremum, 'suff_amp', ftinteger, &
+       nwidth = 1, ndecimals = 0)
+
+  call shp_create_03("outermost_contour_1", shpt_polygon, hshp_outermost)
+  call dbf_add_field_03(ifield, hshp_outermost, 'area', ftdouble, nwidth = 20, &
+       ndecimals = 6)
+  call dbf_add_field_03(ifield, hshp_outermost, 'ssh', ftdouble, nwidth = 13, &
+       ndecimals = 6)
+  call dbf_add_field_03(ifield, hshp_outermost, 'date_index', ftinteger, &
+       nwidth = 4, ndecimals = 0)
+  call dbf_add_field_03(ifield, hshp_outermost, 'eddy_index', ftinteger, &
+       nwidth = 5, ndecimals = 0)
+  call dbf_add_field_03(ifield, hshp_outermost, 'twice', ftinteger, &
+       nwidth = 1, ndecimals = 0)
+
+  call shp_create_03("max_speed_contour_1", shpt_polygon, hshp_max_speed)
+  call dbf_add_field_03(ifield, hshp_max_speed, 'area', ftdouble, nwidth = 20, &
+       ndecimals = 6)
+  call dbf_add_field_03(ifield, hshp_max_speed, 'ssh', ftdouble, nwidth = 13, &
+       ndecimals = 6)
+  call dbf_add_field_03(ifield, hshp_max_speed, 'speed', &
+       ftdouble, nwidth = 13, ndecimals = 6)
+  call dbf_add_field_03(ifield, hshp_max_speed, 'date_index', ftinteger, &
+       nwidth = 4, ndecimals = 0)
+  call dbf_add_field_03(ifield, hshp_max_speed, 'eddy_index', ftinteger, &
+       nwidth = 5, ndecimals = 0)
+
+  call new_unit(unit_isolated)
+  open(unit_isolated, file = "isolated_nodes_1.csv", status = "replace", &
+       action = "write")
+  write(unit_isolated, fmt = *) '"date index" "eddy index"' ! title line
+
+  call new_unit(unit_number_eddies)
+  open(unit_number_eddies, file = "number_eddies_1.csv", status = "replace", &
+       action = "write")
+  write(unit_number_eddies, fmt = *) '"date index" ' &
+       // '"number of visible eddies" "number of interpolated eddies"'
+  ! (title line)
+
+  call dispatch_snapshot(s, hshp_extremum, hshp_outermost, hshp_max_speed, &
+       unit_isolated, unit_number_eddies, k = 1, m = 1, k_begin = 1, &
+       max_delta = 4)
+
+  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".'
+  close(unit_isolated)
+  print *, 'Created "isolated_nodes_1.csv".'
+  close(unit_number_eddies)
+  print *, 'Created "number_eddies_1.csv".'
+  print *, "Average number of points per outermost contour: ", &
+       sum(s%list_vis%outermost_contour%n_points) / real(s%number_vis_eddies)
+
+end program test_set_all_outerm

Tests/tests.json

@@ -101,6 +101,12 @@
 		     "$tests_old_dir/Get_1_outerm_noise_2_8/test_get_1_outerm.dbf"],
 	"input": "&MAIN_NML IND_TARG_EXTR=19,11 /\n"
     },
+    {
+	"args": ["$compil_prod_dir/test_set_all_outerm",
+		 "$input_dir/h_comparison_region.nc"],
+	"title": "Set_all_outerm",
+	"input": "0.001\n"
+    },
     {
 	"args": "$compil_prod_dir/test_get_snapshot",
 	"title": "Get_snapshot",

depend.mk

 dispatch_snapshot.o : write_eddy.o send_snapshot.o derived_types.o 
-get_snapshot.o : set_all_outerm.o set_max_speed.o receive_snapshot.o local_extrema.o derived_types.o 
+get_snapshot.o : set_all_outerm.o set_max_speed.o receive_snapshot.o derived_types.o 
 receive_snapshot.o : derived_types.o 
 send_snapshot.o : derived_types.o 
-set_all_outerm.o : get_1_outerm.o derived_types.o 
+set_all_outerm.o : local_extrema.o get_1_outerm.o derived_types.o 
 set_max_speed.o : spherical_polygon_area.o mean_speed.o max_speed_contour_ssh.o inside_4.o good_contour.o derived_types.o 
 get_1_outerm.o : spherical_polygon_area.o outermost_possible_level.o good_contour.o derived_types.o 
 test_get_snapshot.o : get_snapshot.o dispatch_snapshot.o derived_types.o 
@@ -11,6 +11,7 @@ test_inside_4.o : inside_4.o
 test_local_extrema.o : local_extrema.o 
 test_max_speed_contour_ssh.o : max_speed_contour_ssh.o 
 test_mean_speed.o : mean_speed.o 
+test_set_all_outerm.o : set_all_outerm.o dispatch_snapshot.o derived_types.o 
 test_set_max_speed.o : set_max_speed.o derived_types.o 
 test_get_1_outerm.o : get_1_outerm.o derived_types.o 
 write_eddy.o : derived_types.o 

file_list

@@ -19,6 +19,7 @@ test_inside_4.f
 test_local_extrema.f
 test_max_speed_contour_ssh.f
 test_mean_speed.f
+test_set_all_outerm.f
 test_set_max_speed.f
 test_get_1_outerm.f
 write_eddy.f