Move essentially all the content of procedure get_snapshot to main

unit test_get_snapshot. get_snapshot reads from shapefiles instead of
netCDF files.

Bug fix: programs test_get_snapshot and test_read_snapshot do not need
dispatch_snapshot.f and send_snapshot.f.

Procedures dispatch_snapshot, get_snapshot, send_snapshot and
receive_snapshot are no longer referenced, for now. They will be
called by the overlap program.

Analysis/plot_snapshot.py

@@ -53,8 +53,8 @@ else:
     m_s_iterShapes = reader_m_s.iterShapes()
 
 with netCDF4.Dataset("h.nc") as f:
-    longitude = f.variables["lon"][:]
-    latitude = f.variables["lat"][:]
+    longitude = f.variables["longitude"][:]
+    latitude = f.variables["latitude"][:]
     
 if args.window:
     l = input("llcrnrlon, llcrnrlat, urcrnrlon, urcrnrlat = ? ")\

Documentation_texfol/Graphiques/call_graph.gv

@@ -3,7 +3,6 @@ digraph call_graph
 	main -> {get_snapshot successive_overlap non_successive_overlap};
 	main -> dispatch_snapshot;
 	non_successive_overlap -> {interpolate_eddy weight};
-	dispatch_snapshot -> write_eddy;
 	interpolate_eddy -> write_eddy;
 	successive_overlap -> weight;
 }
\ No newline at end of file

Documentation_texfol/documentation.tex

@@ -992,8 +992,8 @@ fichiers de sortie suivants :
 \end{description}
 
 Les entrées-sorties sont dans : l'algorithme principal directement
-(lecture de longitude et latitude), \verb+get_snapshot+ (lecture de
-SSH, u, v), \verb+successive_overlap+ (écriture dans
+(lecture de longitude et latitude), \verb+get_snapshot+ (lecture des
+shapefiles), \verb+successive_overlap+ (écriture dans
 \verb+edgelist_$m.txt+, arcs entre \verb+k_begin+(m) et
 \verb+k_end_main_loop+(m)), \verb+non_successive_overlap+ directement
 (écriture dans \verb+edgelist_$m.txt+), \verb+write_eddy+ (écriture

GNUmakefile

@@ -13,7 +13,7 @@ src_test_get_1_outerm = good_contour.f test_get_1_outerm.f derived_types.f get_1
 
 src_test_set_max_speed = test_set_max_speed.f derived_types.f set_max_speed.f good_contour.f max_speed_contour_ssh.f mean_speed.f spherical_polyline_area.f inside_4.f
 
-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_polyline_area.f mean_speed.f inside_4.f set_all_outerm.f derived_types.f init_shapefiles.f nearby_extr.f get_var.f
+src_test_get_snapshot = test_get_snapshot.f write_eddy.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_polyline_area.f mean_speed.f inside_4.f set_all_outerm.f derived_types.f init_shapefiles.f nearby_extr.f get_var.f
 
 src_test_set_all_outerm = test_set_all_outerm.f derived_types.f set_all_outerm.f local_extrema.f get_1_outerm.f good_contour.f spherical_polyline_area.f nearby_extr.f get_var.f
 
@@ -23,7 +23,7 @@ src_test_spherical_polygon_area = spherical_polygon_area.f test_spherical_polygo
 
 src_test_read_eddy = test_read_eddy.f derived_types.f init_shapefiles.f read_eddy.f write_eddy.f read_field_indices.f
 
-src_test_read_snapshot = test_read_snapshot.f derived_types.f dispatch_snapshot.f init_shapefiles.f read_snapshot.f write_eddy.f send_snapshot.f read_eddy.f read_field_indices.f
+src_test_read_snapshot = test_read_snapshot.f derived_types.f init_shapefiles.f read_snapshot.f write_eddy.f read_eddy.f read_field_indices.f
 
 src_test_successive_overlap = test_successive_overlap.f derived_types.f successive_overlap.f read_snapshot.f spherical_polygon_area.f spherical_polyline_area.f weight.f read_eddy.f read_field_indices.f
 

Tests/test_get_snapshot.f

 program test_get_snapshot
 
   use, intrinsic:: ieee_arithmetic, only: IEEE_SUPPORT_DATATYPE, &
-       ieee_support_nan
+       ieee_support_nan, ieee_value, IEEE_QUIET_NAN
   use, intrinsic:: ISO_FORTRAN_ENV
 
   ! Libraries:
+  use contour_531, only: convert_to_ind
   use netcdf, only: nf90_nowrite
   use netcdf95, only: nf95_open, find_coord, nf95_inquire_dimension, &
        nf95_get_var, nf95_close
   use nr_util, only: assert, deg_to_rad, twopi, pi
   use shapelib, only: shpfileobject, shpclose
 
-  use derived_types, only: snapshot
-  use get_snapshot_m, only: get_snapshot
+  use derived_types, only: snapshot, null_ssh_contour, missing_speed
+  use get_var_m, only: get_var
   use init_shapefiles_m, only: init_shapefiles
+  use nearby_extr_m, only: nearby_extr
+  use set_all_outerm_m, only: set_all_outerm
+  use set_max_speed_m, only: set_max_speed
   use write_eddy_m, only: write_eddy
 
   implicit none
@@ -22,7 +26,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 max_speed_contour_$m
-  integer i
+  integer i, k
 
   real:: min_amp = 0.
   ! minimum amplitude of ssh, between extremum and outermost contour,
@@ -31,7 +35,7 @@ program test_get_snapshot
   real:: min_radius = 20.
   ! minimum radius of the equal-area disk for an outermost contour, in km
 
-  integer ncid, dimid, varid, nlon, nlat
+  integer ncid, dimid, varid, nlon, nlat, copy
   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
@@ -40,6 +44,29 @@ program test_get_snapshot
   ! maximum radius of an eddy in longitude and latitude, in number of
   ! grid points
 
+  real, allocatable:: ssh(:, :)
+  ! (1 - max_radius(1):nlon + max_radius(1), nlat) if periodic, else
+  ! (nlon, nlat) sea-surface height, in m
+
+  real, allocatable:: u(:, :), v(:, :)
+  ! (1 - max_radius(1):nlon + max_radius(1), nlat) if periodic, else
+  ! (nlon, nlat) wind, in m s-1
+
+  real corner(2)
+  ! longitude and latitude of the corner of the whole grid, in rad
+
+  ! Window around each extremum:
+
+  integer llc(2) ! indices in global grid of lower left corner
+  integer urc(2) ! indices in global grid of upper right corner
+
+  real corner_window(2) ! longitude and latitude of the window around each
+  ! extremum, in rad
+
+  real, allocatable:: outside_points(:, :) ! (2, :) longitude and
+  ! latitude, in rad, of all the significant extrema, except the
+  ! target extremum
+
   namelist /main_nml/ min_amp, max_radius, min_radius
 
   !--------------------------------------------------------------
@@ -77,10 +104,71 @@ program test_get_snapshot
   print *, "periodic = ", periodic
   call assert(2 * max_radius(1) * step(1) < pi, "test_get_snapshot max_radius")
 
-  call get_snapshot(s, min_amp, min_area = pi * (min_radius * 1e3)**2, m = 1, &
-       n_proc = 1, k_end = 1, max_delta = 4, nlon = nlon, nlat = nlat, k = 1, &
-       max_radius = max_radius, corner = [lon_min, lat_min] * deg_to_rad, &
-       step = step, periodic = periodic)
+  copy = merge(max_radius(1), 0, periodic)
+  allocate(ssh(1 - copy:nlon + copy, nlat), u(1 - copy:nlon + copy, nlat), &
+       v(1 - copy:nlon + copy, nlat))
+  k = 1
+  corner = [lon_min, lat_min] * deg_to_rad
+
+  ! Read ssh, u and v at date k:
+
+  call nf95_open("h.nc", nf90_nowrite, ncid)
+  call get_var(periodic, max_radius(1), ssh, ncid, nlon, k, name = "adt", &
+       new_fill_value = huge(0.))
+  ! (We cannot keep the original fill value because Contour_531
+  ! works with an upper limit for valid values.)
+  call nf95_close(ncid)
+
+  call nf95_open("uv.nc", nf90_nowrite, ncid)
+  call get_var(periodic, max_radius(1), u, ncid, nlon, k, name = "u", &
+       new_fill_value = ieee_value(0., IEEE_QUIET_NAN))
+  call get_var(periodic, max_radius(1), v, ncid, nlon, k, name = "v", &
+       new_fill_value = ieee_value(0., IEEE_QUIET_NAN))
+  ! (We will need quiet NaNs rather the original fill values for u
+  ! and v when we compute the max-speed contours.)
+  call nf95_close(ncid)
+
+  call set_all_outerm(s, min_amp, max_radius, step, periodic, ssh, &
+       corner, min_area = pi * (min_radius * 1e3)**2)
+
+  ! Done with outermost contours, now let us take care of the
+  ! max-speed contours.
+
+  do i = 1, s%number_vis_extr
+     if (s%list_vis(i)%valid) then
+        ! Restrict the field to the outermost contour:
+
+        llc = floor(convert_to_ind(minval(s%list_vis(i)%out_cont%points, &
+             dim = 2), corner, step))
+
+        urc = ceiling(convert_to_ind(maxval( &
+             s%list_vis(i)%out_cont%points, dim = 2), corner, step))
+
+        ! Should have no effect except because of roundup error:
+        urc(2) = min(urc(2), nlat)
+        if (.not. periodic) urc(1) = min(urc(1), nlon)
+
+        corner_window = corner + (llc - 1) * step
+
+        outside_points = nearby_extr(s%extr_map(llc(1):urc(1), &
+             llc(2):urc(2)), s%list_vis, i)
+
+        if (periodic) outside_points(1, :) = outside_points(1, :) &
+             + ceiling((corner_window(1) - outside_points(1, :)) / twopi) &
+             * twopi
+        ! (Shift the longitude of each point to a value close to
+        ! the longitude of the target extremum.)
+
+        call set_max_speed(s%list_vis(i), s%ind_extr(:, i) - llc + 1, &
+             outside_points, ssh(llc(1):urc(1), llc(2):urc(2)), &
+             u(llc(1):urc(1), llc(2):urc(2)), &
+             v(llc(1):urc(1), llc(2):urc(2)), corner_window, step)
+     else
+        s%list_vis(i)%speed_cont = null_ssh_contour()
+        s%list_vis(i)%max_speed = missing_speed
+        s%list_vis(i)%radius4 = 0
+     end if
+  end do
 
   call init_shapefiles(hshp_extremum, hshp_outermost, hshp_max_speed)
 

depend.mk

-dispatch_snapshot.o : send_snapshot.o derived_types.o 
 get_1_outerm.o : spherical_polyline_area.o good_contour.o derived_types.o 
-get_snapshot.o : set_all_outerm.o set_max_speed.o receive_snapshot.o nearby_extr.o get_var.o derived_types.o 
 nearby_extr.o : derived_types.o 
 read_eddy.o : read_field_indices.o derived_types.o 
 read_snapshot.o : read_field_indices.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 nearby_extr.o get_1_outerm.o derived_types.o 
 set_max_speed.o : spherical_polyline_area.o mean_speed.o max_speed_contour_ssh.o inside_4.o good_contour.o derived_types.o 
 spherical_polygon_area.o : spherical_polyline_area.o 
 successive_overlap.o : weight.o spherical_polyline_area.o spherical_polygon_area.o derived_types.o 
 test_get_1_outerm.o : get_1_outerm.o derived_types.o 
-test_get_snapshot.o : write_eddy.o init_shapefiles.o get_snapshot.o derived_types.o 
+test_get_snapshot.o : write_eddy.o set_max_speed.o set_all_outerm.o nearby_extr.o init_shapefiles.o get_var.o derived_types.o 
 test_local_extrema.o : local_extrema.o 
 test_nearby_extr.o : nearby_extr.o derived_types.o 
 test_read_eddy.o : write_eddy.o read_field_indices.o read_eddy.o init_shapefiles.o derived_types.o 

get_snapshot.f

@@ -4,33 +4,20 @@ module get_snapshot_m
 
 contains
 
-  subroutine get_snapshot(s, min_amp, min_area, m, n_proc, k_end, max_delta, &
-       nlon, nlat, k, max_radius, corner, step, periodic)
+  subroutine get_snapshot(s, m, n_proc, k_end, max_delta, nlon, nlat, k, corner)
 
     use, intrinsic:: ieee_arithmetic, only: ieee_value, IEEE_QUIET_NAN
 
     ! Libraries:
-    use contour_531, only: convert_to_ind, null_polyline
-    use netcdf, only: nf90_nowrite
-    use netcdf95, only: nf95_open, nf95_close
-    use nr_util, only: twopi
-
-    use derived_types, only: snapshot, null_ssh_contour, missing_speed
-    use get_var_m, only: get_var
-    use nearby_extr_m, only: nearby_extr
+    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 receive_snapshot_m, only: receive_snapshot
-    use set_max_speed_m, only: set_max_speed
-    use set_all_outerm_m, only: set_all_outerm
 
     type(snapshot), intent(out):: s
 
-    real, intent(in):: min_amp
-    ! minimum amplitude of ssh, between extremum and outermost
-    ! contour, in m
-
-    real, intent(in):: min_area
-    ! minimum area of an outermost contour, in m2
-    
     integer, intent(in):: m ! rank of MPI process
     integer, intent(in):: n_proc ! number of MPI processes
     integer, intent(in):: k_end ! last date index analyzed by this MPI process
@@ -42,98 +29,28 @@ contains
     integer, intent(in):: nlon, nlat
     integer, intent(in):: k ! date index
 
-    integer, intent(in):: max_radius(:) ! (2) maximum radius of an
-    ! eddy in longitude and latitude, in number of grid points
-
     real, intent(in):: corner(:) ! (2) longitude and latitude of the
     ! corner of the whole grid, in rad
-    
-    real, intent(in):: step(:) ! (2) longitude and latitude steps, in rad
-    logical, intent(in):: periodic ! grid is periodic in longitude
 
     ! Local:
-
-    integer ncid
-    real ssh(1 - merge(max_radius(1), 0, periodic):nlon + merge(max_radius(1), &
-         0, periodic), nlat) ! sea-surface height, in m
-    real, dimension(1 - merge(max_radius(1), 0, periodic):nlon &
-         + merge(max_radius(1), 0, periodic), nlat):: u, v ! wind, in m s-1
-    integer i
-
-    ! Window around each extremum:
-    
-    integer llc(2) ! indices in global grid of lower left corner
-    integer urc(2) ! indices in global grid of upper right corner
-
-    real corner_window(2) ! longitude and latitude of the window around each
-    ! extremum, in rad
-
-    real, allocatable:: outside_points(:, :) ! (2, :) longitude and
-    ! latitude, in rad, of all the significant extrema, except the
-    ! target extremum
+    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
 
     !--------------------------------------------------------------
 
     if (m == n_proc .or. k <= k_end - max_delta) then
-       ! Read ssh, u and v at date k:
-
-       call nf95_open("h.nc", nf90_nowrite, ncid)
-       call get_var(periodic, max_radius(1), ssh, ncid, nlon, k, name = "adt", &
-            new_fill_value = huge(0.))
-       ! (We cannot keep the original fill value because Contour_531
-       ! works with an upper limit for valid values.)
-       call nf95_close(ncid)
-
-       call nf95_open("uv.nc", nf90_nowrite, ncid)
-       call get_var(periodic, max_radius(1), u, ncid, nlon, k, name = "u", &
-            new_fill_value = ieee_value(0., IEEE_QUIET_NAN))
-       call get_var(periodic, max_radius(1), v, ncid, nlon, k, name = "v", &
-            new_fill_value = ieee_value(0., IEEE_QUIET_NAN))
-       ! (We will need quiet NaNs rather the original fill values for u
-       ! and v when we compute the max-speed contours.)
-       call nf95_close(ncid)
-
-       call set_all_outerm(s, min_amp, max_radius, step, periodic, ssh, &
-            corner, min_area)
-
-       ! Done with outermost contours, now let us take care of the
-       ! max-speed contours.
-
-       do i = 1, s%number_vis_extr
-          if (s%list_vis(i)%valid) then
-             ! Restrict the field to the outermost contour:
-
-             llc = floor(convert_to_ind(minval(s%list_vis(i)%out_cont%points, &
-                  dim = 2), corner, step))
-
-             urc = ceiling(convert_to_ind(maxval( &
-                  s%list_vis(i)%out_cont%points, dim = 2), corner, step))
-
-             ! Should have no effect except because of roundup error:
-             urc(2) = min(urc(2), nlat)
-             if (.not. periodic) urc(1) = min(urc(1), nlon)
-
-             corner_window = corner + (llc - 1) * step
-             
-             outside_points = nearby_extr(s%extr_map(llc(1):urc(1), &
-                  llc(2):urc(2)), s%list_vis, i)
-
-             if (periodic) outside_points(1, :) = outside_points(1, :) &
-                  + ceiling((corner_window(1) - outside_points(1, :)) / twopi) &
-                  * twopi
-             ! (Shift the longitude of each point to a value close to
-             ! the longitude of the target extremum.)
-             
-             call set_max_speed(s%list_vis(i), s%ind_extr(:, i) - llc + 1, &
-                  outside_points, ssh(llc(1):urc(1), llc(2):urc(2)), &
-                  u(llc(1):urc(1), llc(2):urc(2)), &
-                  v(llc(1):urc(1), llc(2):urc(2)), corner_window, step)
-          else
-             s%list_vis(i)%speed_cont = null_ssh_contour()
-             s%list_vis(i)%max_speed = missing_speed
-             s%list_vis(i)%radius4 = 0
-          end if
-       end do
+       ! Read at date k:
+       call shp_open_03("extremum", pszaccess = "rb", hshp = hshp_extremum)
+       call shp_open_03("outermost_contour", pszaccess = "rb", &
+            hshp = hshp_outermost)
+       call shp_open_03("max_speed_contour", pszaccess = "rb", &
+            hshp = hshp_max_speed)
+       call read_snapshot(s, hshp_extremum, hshp_outermost, hshp_max_speed, &
+            corner, nlon, nlat)
+       CALL shpclose(hshp_extremum)
+       CALL shpclose(hshp_outermost)
+       CALL shpclose(hshp_max_speed)
 
        s%number_eddies = s%number_vis_extr
     else