diff --git a/GNUmakefile b/GNUmakefile
index 492d9ee072adfb8b45d5c1433eb149188deda9ff..2c198fefcf2088720a30b54ec02f31c852a38439 100644
--- a/GNUmakefile
+++ b/GNUmakefile
@@ -21,7 +21,9 @@ 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
 
-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
+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
+
+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} ${src_test_read_snapshot}) 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
 
@@ -35,9 +37,11 @@ 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)
+obj_test_read_snapshot := $(src_test_read_snapshot:.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 test_read_eddy
+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 test_read_snapshot
 
 # 3. Compiler-dependent part
 
@@ -75,6 +79,7 @@ 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}
+test_read_snapshot: ${obj_test_read_snapshot}
 
 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
diff --git a/Tests/read_eddy.f b/Tests/read_eddy.f
index f15327c8c3b40c7622d3659197f9928ce85c1a60..0aa2ff8c38e0710391fe31c32007118d83e86d6a 100644
--- a/Tests/read_eddy.f
+++ b/Tests/read_eddy.f
@@ -29,7 +29,7 @@ contains
     TYPE(shpfileobject), intent(inout):: hshp_max_speed
     ! shapefile x_speed_contour_$m
 
-    integer, intent(in):: ishape
+    integer, intent(in):: ishape ! 0-based
 
     ! Local:
     integer interpolated, cyclone, suff_amp, twice
diff --git a/Tests/read_snapshot.f b/Tests/read_snapshot.f
new file mode 100644
index 0000000000000000000000000000000000000000..3b0086239f42c09f847470deefe55ba1e85a0802
--- /dev/null
+++ b/Tests/read_snapshot.f
@@ -0,0 +1,62 @@
+module read_snapshot_m
+
+  implicit none
+
+contains
+
+  subroutine read_snapshot(s, hshp_extremum, hshp_outermost, hshp_max_speed)
+
+    use, intrinsic:: ISO_C_BINDING
+
+    use contour_531, only: convert_to_ind
+    use derived_types, only: snapshot, eddy
+    use nr_util, only: pi
+    use read_eddy_m, only: read_eddy
+    use shapelib, only: shpfileobject, shpgetinfo
+
+    type(snapshot), intent(out):: s
+    TYPE(shpfileobject), intent(inout):: hshp_extremum ! shapefile extremum_1
+
+    TYPE(shpfileobject), intent(inout):: hshp_outermost
+    ! shapefile outermost_contour_1
+    
+    TYPE(shpfileobject), intent(inout):: hshp_max_speed
+    ! shapefile max_speed_contour_1
+
+    ! Local:
+    integer ishape, k, i
+    integer shapetype, dbffieldcount, dbfrecordcount
+    real(c_double) minbound(4), maxbound(4)
+    type(eddy) e
+    real, parameter:: deg_to_rad = pi / 180.
+
+    !---------------------------------------------------------------------
+
+    call shpgetinfo(hshp_extremum, s%number_vis_eddies, shapetype, minbound, &
+         maxbound, dbffieldcount, dbfrecordcount)
+    allocate(s%list_vis(s%number_vis_eddies))
+
+    do ishape = 0, s%number_vis_eddies - 1
+       call read_eddy(e, k, i, hshp_extremum, hshp_outermost, hshp_max_speed, &
+            ishape)
+       s%list_vis(i) = e
+    end do
+
+    allocate(s%ind_extr(2, s%number_vis_eddies))
+
+    forall (i = 1:s%number_vis_eddies) s%ind_extr(:, i) &
+         = nint(convert_to_ind(s%list_vis(i)%coord_extr, &
+         corner = [0.125, -59.875] * deg_to_rad, &
+         step = [0.25, 0.25] * deg_to_rad))
+
+    allocate(s%extr_map(120, 180))
+    s%extr_map = 0
+    forall (i = 1:s%number_vis_eddies) &
+         s%extr_map(s%ind_extr(1, i), s%ind_extr(2, i)) &
+         = merge(i, -i, s%list_vis(i)%suff_amp)
+
+    s%number_eddies = s%number_vis_eddies
+
+  end subroutine read_snapshot
+  
+end module read_snapshot_m
diff --git a/Tests/test_read_snapshot.f b/Tests/test_read_snapshot.f
new file mode 100644
index 0000000000000000000000000000000000000000..1eeacf646d5d2a42ec4264573a89411b18b7ec19
--- /dev/null
+++ b/Tests/test_read_snapshot.f
@@ -0,0 +1,60 @@
+program test_read_snapshot
+
+  use derived_types, only: snapshot
+  use dispatch_snapshot_m, only: dispatch_snapshot
+  use init_shapefiles_m, only: init_shapefiles
+  use jumble, only: new_unit
+  use read_snapshot_m, only: read_snapshot
+  use shapelib, only: shpfileobject, shpclose
+  use shapelib_03, only: shp_open_03
+
+  implicit none
+
+  type(snapshot) s
+  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
+  integer unit_isolated, unit_number_eddies
+
+  !-------------------------------------------------------------------------
+
+  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_snapshot(s, hshp_extremum, hshp_outermost, hshp_max_speed)
+  CALL shpclose(hshp_extremum)
+  CALL shpclose(hshp_outermost)
+  CALL shpclose(hshp_max_speed)
+
+  call init_shapefiles(hshp_extremum, hshp_outermost, hshp_max_speed)
+  
+  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".'
+
+end program test_read_snapshot
diff --git a/Tests/tests.json b/Tests/tests.json
index 1b90bbe89cbdc23461d9ec42c4cc99ac148a860f..15e2a0f67c463fa315a0672bf9e9ad39463090d3 100644
--- a/Tests/tests.json
+++ b/Tests/tests.json
@@ -207,5 +207,10 @@
 	"args": "$compil_prod_dir/test_read_eddy",
 	"required": ["$input_dir/Test_read_eddy/*"],
 	"title": "Read_eddy"
+    },
+    {
+	"args": "$compil_prod_dir/test_read_snapshot",
+	"required": ["$input_dir/Test_read_eddy/*"],
+	"title": "Read_snapshot"
     }
 ]
diff --git a/depend.mk b/depend.mk
index ff7be99d0fbabd89bae7afde8cd3259caefd0ce9..93ccfd2b324f70e35001cbcce8c98fa39ef5ea23 100644
--- a/depend.mk
+++ b/depend.mk
@@ -2,6 +2,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 
+read_snapshot.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 
@@ -11,6 +12,7 @@ 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_read_snapshot.o : read_snapshot.o init_shapefiles.o dispatch_snapshot.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