CarpetIOF5: Make generic reading work

20 files changed, 1103 insertions, 247 deletions

diff --git a/CarpetDev/CarpetIOF5/configuration.ccl b/CarpetDev/CarpetIOF5/configuration.ccl
index 09fb32efd..a08dea2c0 100644
--- a/CarpetDev/CarpetIOF5/configuration.ccl
+++ b/CarpetDev/CarpetIOF5/configuration.ccl
@@ -1,6 +1,7 @@
 # Configuration definition for thorn CarpetIOF5
 
 REQUIRES Carpet CarpetLib F5 HDF5
-OPTIONAL LoopControl
+
+OPTIONAL LoopControl MPI
 {
 }
diff --git a/CarpetDev/CarpetIOF5/par/iof5-checkpoint.par b/CarpetDev/CarpetIOF5/par/iof5-checkpoint.par
index 3306b73e6..f3bacfa89 100644
--- a/CarpetDev/CarpetIOF5/par/iof5-checkpoint.par
+++ b/CarpetDev/CarpetIOF5/par/iof5-checkpoint.par
@@ -49,3 +49,4 @@ InitBase::initial_data_setup_method = "init_all_levels"
 
 CarpetIOF5::checkpoint = yes
 IO::checkpoint_every   = 1024
+IO::checkpoint_dir     = $parfile
diff --git a/CarpetDev/CarpetIOF5/par/iof5-multipatch-input.par b/CarpetDev/CarpetIOF5/par/iof5-multipatch-input.par
new file mode 100644
index 000000000..202cb179b
--- /dev/null
+++ b/CarpetDev/CarpetIOF5/par/iof5-multipatch-input.par
@@ -0,0 +1,52 @@
+ActiveThorns = "
+        Boundary
+        Carpet
+        CarpetIOF5
+        CarpetLib
+        CarpetRegrid2
+        CartGrid3D
+        CoordBase
+        Coordinates
+        F5
+        Formaline
+        GSL
+        HDF5
+        InitBase
+        IOUtil
+        LoopControl
+        SymBase
+"
+
+
+
+Cactus::cctk_run_title             = "IOF5"
+Cactus::cctk_full_warnings         = yes
+Cactus::highlight_warning_messages = no
+
+Cactus::cctk_itlast = 0
+
+IO::out_dir = $parfile
+
+Carpet::domain_from_multipatch   = yes
+CartGrid3D::type                 = "multipatch"
+Coordinates::coordinate_system   = "Thornburg04"
+Coordinates::h_cartesian         = 0.25
+Coordinates::h_radial            = 0.25
+Coordinates::sphere_inner_radius = 1.0
+Coordinates::sphere_outer_radius = 3.0
+Coordinates::n_angular           = 8
+
+Carpet::max_refinement_levels = 10
+CarpetRegrid2::num_centres    =  1
+CarpetRegrid2::num_levels_1   =  3
+CarpetRegrid2::radius_1[1]    =  0.5
+CarpetRegrid2::radius_1[2]    =  0.2
+
+InitBase::initial_data_setup_method = "init_all_levels"
+
+IO::filereader_ID_dir   = "iof5-multipatch"
+IO::filereader_ID_files = "iof5-multipatch"
+IO::filereader_ID_vars  = "grid::coordinates"
+
+CarpetIOF5::out_every = 1
+CarpetIOF5::out_vars  = "grid::coordinates"
diff --git a/CarpetDev/CarpetIOF5/par/iof5-refined-cell-input.par b/CarpetDev/CarpetIOF5/par/iof5-refined-cell-input.par
new file mode 100644
index 000000000..322067721
--- /dev/null
+++ b/CarpetDev/CarpetIOF5/par/iof5-refined-cell-input.par
@@ -0,0 +1,62 @@
+ActiveThorns = "
+        Boundary
+        Carpet
+        CarpetIOF5
+        CarpetLib
+        CarpetRegrid2
+        CartGrid3D
+        CoordBase
+        F5
+        Formaline
+        GSL
+        HDF5
+        InitBase
+        IOUtil
+        LoopControl
+        SymBase
+"
+
+
+
+Cactus::cctk_run_title             = "IOF5"
+Cactus::cctk_full_warnings         = yes
+Cactus::highlight_warning_messages = no
+
+Cactus::cctk_itlast = 1024
+
+IO::out_dir = $parfile
+
+Carpet::domain_from_coordbase         = yes
+CartGrid3D::type                      = "coordbase"
+CoordBase::domainsize                 = "minmax"
+CoordBase::xmin                       = -1.0
+CoordBase::ymin                       = -1.0
+CoordBase::zmin                       = -1.0
+CoordBase::xmax                       = +1.0
+CoordBase::ymax                       = +1.0
+CoordBase::zmax                       = +1.0
+CoordBase::dx                         =  0.25
+CoordBase::dy                         =  0.25
+CoordBase::dz                         =  0.25
+CoordBase::boundary_staggered_x_lower = yes
+CoordBase::boundary_staggered_y_lower = yes
+CoordBase::boundary_staggered_z_lower = yes
+CoordBase::boundary_staggered_x_upper = yes
+CoordBase::boundary_staggered_y_upper = yes
+CoordBase::boundary_staggered_z_upper = yes
+
+Carpet::refinement_centering  = "cell"
+Carpet::max_refinement_levels = 10
+CarpetRegrid2::num_centres    =  1
+CarpetRegrid2::num_levels_1   =  3
+CarpetRegrid2::radius_1[1]    =  0.5
+CarpetRegrid2::radius_1[2]    =  0.2
+
+InitBase::initial_data_setup_method = "init_all_levels"
+
+IO::filereader_ID_dir   = "iof5-refined-cell"
+IO::filereader_ID_files = "iof5-refined-cell"
+IO::filereader_ID_vars  = "grid::coordinates"
+
+CarpetIOF5::out_every = 256
+CarpetIOF5::out_vars  = "grid::coordinates"
diff --git a/CarpetDev/CarpetIOF5/par/iof5-refined-input.par b/CarpetDev/CarpetIOF5/par/iof5-refined-input.par
new file mode 100644
index 000000000..30917f401
--- /dev/null
+++ b/CarpetDev/CarpetIOF5/par/iof5-refined-input.par
@@ -0,0 +1,55 @@
+ActiveThorns = "
+        Boundary
+        Carpet
+        CarpetIOF5
+        CarpetLib
+        CarpetRegrid2
+        CartGrid3D
+        CoordBase
+        F5
+        Formaline
+        GSL
+        HDF5
+        InitBase
+        IOUtil
+        LoopControl
+        SymBase
+"
+
+
+
+Cactus::cctk_run_title             = "IOF5"
+Cactus::cctk_full_warnings         = yes
+Cactus::highlight_warning_messages = no
+
+Cactus::cctk_itlast = 1024
+
+IO::out_dir = $parfile
+
+Carpet::domain_from_coordbase = yes
+CartGrid3D::type              = "coordbase"
+CoordBase::domainsize         = "minmax"
+CoordBase::xmin               = -1.0
+CoordBase::ymin               = -1.0
+CoordBase::zmin               = -1.0
+CoordBase::xmax               = +1.0
+CoordBase::ymax               = +1.0
+CoordBase::zmax               = +1.0
+CoordBase::dx                 =  0.25
+CoordBase::dy                 =  0.25
+CoordBase::dz                 =  0.25
+
+Carpet::max_refinement_levels = 10
+CarpetRegrid2::num_centres    =  1
+CarpetRegrid2::num_levels_1   =  3
+CarpetRegrid2::radius_1[1]    =  0.5
+CarpetRegrid2::radius_1[2]    =  0.2
+
+InitBase::initial_data_setup_method = "init_all_levels"
+
+IO::filereader_ID_dir   = "iof5-refined"
+IO::filereader_ID_files = "iof5-refined"
+IO::filereader_ID_vars  = "grid::coordinates"
+
+CarpetIOF5::out_every = 256
+CarpetIOF5::out_vars  = "grid::coordinates"
diff --git a/CarpetDev/CarpetIOF5/par/iof5-uniform-fragment-input.par b/CarpetDev/CarpetIOF5/par/iof5-uniform-fragment-input.par
new file mode 100644
index 000000000..3593ea50f
--- /dev/null
+++ b/CarpetDev/CarpetIOF5/par/iof5-uniform-fragment-input.par
@@ -0,0 +1,46 @@
+ActiveThorns = "
+        Boundary
+        Carpet
+        CarpetIOF5
+        CarpetLib
+        CartGrid3D
+        CoordBase
+        F5
+        Formaline
+        GSL
+        HDF5
+        InitBase
+        IOUtil
+        LoopControl
+        SymBase
+"
+
+
+
+Cactus::cctk_run_title             = "IOF5"
+Cactus::cctk_full_warnings         = yes
+Cactus::highlight_warning_messages = no
+
+Cactus::cctk_itlast = 2
+
+IO::out_dir = $parfile
+
+Carpet::domain_from_coordbase = yes
+CartGrid3D::type              = "coordbase"
+CoordBase::domainsize         = "minmax"
+CoordBase::xmin               = +1.0
+CoordBase::ymin               = +1.0
+CoordBase::zmin               = +1.0
+CoordBase::xmax               = +3.0
+CoordBase::ymax               = +3.0
+CoordBase::zmax               = +3.0
+CoordBase::dx                 =  0.25
+CoordBase::dy                 =  0.25
+CoordBase::dz                 =  0.25
+
+IO::filereader_ID_dir   = "iof5-uniform-fragment"
+IO::filereader_ID_files = "iof5-uniform-fragment"
+IO::filereader_ID_vars  = "grid::coordinates"
+
+CarpetIOF5::out_every = 2
+CarpetIOF5::out_vars  = "grid::coordinates"
diff --git a/CarpetDev/CarpetIOF5/par/iof5-uniform-fragment.par b/CarpetDev/CarpetIOF5/par/iof5-uniform-fragment.par
new file mode 100644
index 000000000..6c363bb52
--- /dev/null
+++ b/CarpetDev/CarpetIOF5/par/iof5-uniform-fragment.par
@@ -0,0 +1,44 @@
+ActiveThorns = "
+        Boundary
+        Carpet
+        CarpetIOF5
+        CarpetLib
+        CartGrid3D
+        CoordBase
+        F5
+        Formaline
+        GSL
+        HDF5
+        InitBase
+        IOUtil
+        LoopControl
+        SymBase
+"
+
+
+
+Cactus::cctk_run_title             = "IOF5"
+Cactus::cctk_full_warnings         = yes
+Cactus::highlight_warning_messages = no
+
+Cactus::cctk_itlast = 2
+
+IO::out_dir = $parfile
+
+Carpet::domain_from_coordbase = yes
+CartGrid3D::type              = "coordbase"
+CoordBase::domainsize         = "minmax"
+CoordBase::xmin               = -1.0
+CoordBase::ymin               = -1.0
+CoordBase::zmin               = -1.0
+CoordBase::xmax               = +1.0
+CoordBase::ymax               = +1.0
+CoordBase::zmax               = +1.0
+CoordBase::dx                 =  0.25
+CoordBase::dy                 =  0.25
+CoordBase::dz                 =  0.25
+
+CarpetIOF5::fragment_contiguous_components = yes
+
+CarpetIOF5::out_every = 2
+CarpetIOF5::out_vars  = "grid::coordinates"
diff --git a/CarpetDev/CarpetIOF5/par/iof5-uniform-input.par b/CarpetDev/CarpetIOF5/par/iof5-uniform-input.par
new file mode 100644
index 000000000..7dfbefd80
--- /dev/null
+++ b/CarpetDev/CarpetIOF5/par/iof5-uniform-input.par
@@ -0,0 +1,46 @@
+ActiveThorns = "
+        Boundary
+        Carpet
+        CarpetIOF5
+        CarpetLib
+        CartGrid3D
+        CoordBase
+        F5
+        Formaline
+        GSL
+        HDF5
+        InitBase
+        IOUtil
+        LoopControl
+        SymBase
+"
+
+
+
+Cactus::cctk_run_title             = "IOF5"
+Cactus::cctk_full_warnings         = yes
+Cactus::highlight_warning_messages = no
+
+Cactus::cctk_itlast = 2
+
+IO::out_dir = $parfile
+
+Carpet::domain_from_coordbase = yes
+CartGrid3D::type              = "coordbase"
+CoordBase::domainsize         = "minmax"
+CoordBase::xmin               = +1.0
+CoordBase::ymin               = +1.0
+CoordBase::zmin               = +1.0
+CoordBase::xmax               = +3.0
+CoordBase::ymax               = +3.0
+CoordBase::zmax               = +3.0
+CoordBase::dx                 =  0.25
+CoordBase::dy                 =  0.25
+CoordBase::dz                 =  0.25
+
+IO::filereader_ID_dir   = "iof5-uniform"
+IO::filereader_ID_files = "iof5-uniform"
+IO::filereader_ID_vars  = "grid::coordinates"
+
+CarpetIOF5::out_every = 2
+CarpetIOF5::out_vars  = "grid::coordinates"
diff --git a/CarpetDev/CarpetIOF5/par/iof5-uniform-separate-fragment-input.par b/CarpetDev/CarpetIOF5/par/iof5-uniform-separate-fragment-input.par
new file mode 100644
index 000000000..55f69a326
--- /dev/null
+++ b/CarpetDev/CarpetIOF5/par/iof5-uniform-separate-fragment-input.par
@@ -0,0 +1,46 @@
+ActiveThorns = "
+        Boundary
+        Carpet
+        CarpetIOF5
+        CarpetLib
+        CartGrid3D
+        CoordBase
+        F5
+        Formaline
+        GSL
+        HDF5
+        InitBase
+        IOUtil
+        LoopControl
+        SymBase
+"
+
+
+
+Cactus::cctk_run_title             = "IOF5"
+Cactus::cctk_full_warnings         = yes
+Cactus::highlight_warning_messages = no
+
+Cactus::cctk_itlast = 2
+
+IO::out_dir = $parfile
+
+Carpet::domain_from_coordbase = yes
+CartGrid3D::type              = "coordbase"
+CoordBase::domainsize         = "minmax"
+CoordBase::xmin               = +1.0
+CoordBase::ymin               = +1.0
+CoordBase::zmin               = +1.0
+CoordBase::xmax               = +3.0
+CoordBase::ymax               = +3.0
+CoordBase::zmax               = +3.0
+CoordBase::dx                 =  0.25
+CoordBase::dy                 =  0.25
+CoordBase::dz                 =  0.25
+
+IO::filereader_ID_dir   = "iof5-uniform-separate-fragment"
+IO::filereader_ID_files = "iof5-uniform-separate-fragment"
+IO::filereader_ID_vars  = "grid::coordinates"
+
+CarpetIOF5::out_every = 2
+CarpetIOF5::out_vars  = "grid::coordinates"
diff --git a/CarpetDev/CarpetIOF5/par/iof5-uniform-separate-fragment.par b/CarpetDev/CarpetIOF5/par/iof5-uniform-separate-fragment.par
new file mode 100644
index 000000000..9f319c9ed
--- /dev/null
+++ b/CarpetDev/CarpetIOF5/par/iof5-uniform-separate-fragment.par
@@ -0,0 +1,45 @@
+ActiveThorns = "
+        Boundary
+        Carpet
+        CarpetIOF5
+        CarpetLib
+        CartGrid3D
+        CoordBase
+        F5
+        Formaline
+        GSL
+        HDF5
+        InitBase
+        IOUtil
+        LoopControl
+        SymBase
+"
+
+
+
+Cactus::cctk_run_title             = "IOF5"
+Cactus::cctk_full_warnings         = yes
+Cactus::highlight_warning_messages = no
+
+Cactus::cctk_itlast = 2
+
+IO::out_dir = $parfile
+
+Carpet::domain_from_coordbase = yes
+CartGrid3D::type              = "coordbase"
+CoordBase::domainsize         = "minmax"
+CoordBase::xmin               = -1.0
+CoordBase::ymin               = -1.0
+CoordBase::zmin               = -1.0
+CoordBase::xmax               = +1.0
+CoordBase::ymax               = +1.0
+CoordBase::zmax               = +1.0
+CoordBase::dx                 =  0.25
+CoordBase::dy                 =  0.25
+CoordBase::dz                 =  0.25
+
+CarpetIOF5::separate_single_component_tensors = yes
+CarpetIOF5::fragment_contiguous_components    = yes
+
+CarpetIOF5::out_every = 2
+CarpetIOF5::out_vars  = "grid::coordinates"
diff --git a/CarpetDev/CarpetIOF5/par/iof5-uniform-separate-input.par b/CarpetDev/CarpetIOF5/par/iof5-uniform-separate-input.par
new file mode 100644
index 000000000..4cfc8ac32
--- /dev/null
+++ b/CarpetDev/CarpetIOF5/par/iof5-uniform-separate-input.par
@@ -0,0 +1,46 @@
+ActiveThorns = "
+        Boundary
+        Carpet
+        CarpetIOF5
+        CarpetLib
+        CartGrid3D
+        CoordBase
+        F5
+        Formaline
+        GSL
+        HDF5
+        InitBase
+        IOUtil
+        LoopControl
+        SymBase
+"
+
+
+
+Cactus::cctk_run_title             = "IOF5"
+Cactus::cctk_full_warnings         = yes
+Cactus::highlight_warning_messages = no
+
+Cactus::cctk_itlast = 2
+
+IO::out_dir = $parfile
+
+Carpet::domain_from_coordbase = yes
+CartGrid3D::type              = "coordbase"
+CoordBase::domainsize         = "minmax"
+CoordBase::xmin               = +1.0
+CoordBase::ymin               = +1.0
+CoordBase::zmin               = +1.0
+CoordBase::xmax               = +3.0
+CoordBase::ymax               = +3.0
+CoordBase::zmax               = +3.0
+CoordBase::dx                 =  0.25
+CoordBase::dy                 =  0.25
+CoordBase::dz                 =  0.25
+
+IO::filereader_ID_dir   = "iof5-uniform-separate"
+IO::filereader_ID_files = "iof5-uniform-separate"
+IO::filereader_ID_vars  = "grid::coordinates"
+
+CarpetIOF5::out_every = 2
+CarpetIOF5::out_vars  = "grid::coordinates"
diff --git a/CarpetDev/CarpetIOF5/par/iof5-uniform-separate.par b/CarpetDev/CarpetIOF5/par/iof5-uniform-separate.par
new file mode 100644
index 000000000..78d4bd290
--- /dev/null
+++ b/CarpetDev/CarpetIOF5/par/iof5-uniform-separate.par
@@ -0,0 +1,44 @@
+ActiveThorns = "
+        Boundary
+        Carpet
+        CarpetIOF5
+        CarpetLib
+        CartGrid3D
+        CoordBase
+        F5
+        Formaline
+        GSL
+        HDF5
+        InitBase
+        IOUtil
+        LoopControl
+        SymBase
+"
+
+
+
+Cactus::cctk_run_title             = "IOF5"
+Cactus::cctk_full_warnings         = yes
+Cactus::highlight_warning_messages = no
+
+Cactus::cctk_itlast = 2
+
+IO::out_dir = $parfile
+
+Carpet::domain_from_coordbase = yes
+CartGrid3D::type              = "coordbase"
+CoordBase::domainsize         = "minmax"
+CoordBase::xmin               = -1.0
+CoordBase::ymin               = -1.0
+CoordBase::zmin               = -1.0
+CoordBase::xmax               = +1.0
+CoordBase::ymax               = +1.0
+CoordBase::zmax               = +1.0
+CoordBase::dx                 =  0.25
+CoordBase::dy                 =  0.25
+CoordBase::dz                 =  0.25
+
+CarpetIOF5::separate_single_component_tensors = yes
+
+CarpetIOF5::out_every = 2
+CarpetIOF5::out_vars  = "grid::coordinates"
diff --git a/CarpetDev/CarpetIOF5/param.ccl b/CarpetDev/CarpetIOF5/param.ccl
index 6f1b1f9f5..ee44950de 100644
--- a/CarpetDev/CarpetIOF5/param.ccl
+++ b/CarpetDev/CarpetIOF5/param.ccl
@@ -2,19 +2,30 @@
 
 SHARES: IO
 
+USES STRING  filereader_ID_dir AS IO_filereader_ID_dir
+
 USES STRING  out_dir   AS IO_out_dir
 USES INT     out_every AS IO_out_every
 
 USES INT     out_timesteps_per_file
 
+USES STRING  recover_dir AS IO_recover_dir
+
 USES BOOLEAN checkpoint_ID
 USES INT     checkpoint_every
 USES BOOLEAN checkpoint_on_terminate
+USES STRING  checkpoint_dir AS IO_checkpoint_dir
 
 
 
 PRIVATE:
 
+STRING filereader_ID_dir "Input directory (overrides IO::filereader_ID_dir)" STEERABLE=always
+{
+  "^$" :: "Empty: use IO::filereader_ID_dir"
+  ".+" :: "Not empty: directory name"
+} ""
+
 STRING out_dir "Output directory (overrides IO::out_dir)" STEERABLE=always
 {
   "^$" :: "Empty: use IO::out_dir"
@@ -77,6 +88,10 @@ BOOLEAN separate_single_component_tensors "Use a separated representation even f
 {
 } "no"
 
+BOOLEAN fragment_contiguous_components "Use a fragmented representation even for linear representations" STEERABLE=always
+{
+} "no"
+
 
 
 BOOLEAN output_all_timelevels "Output all timelevels (instead of only the current)" STEERABLE=always
@@ -102,6 +117,12 @@ BOOLEAN output_boundary_points "Output outer boundary points" STEERABLE=always
 
 
 
+STRING recover_dir "Recover directory (overrides IO::recover_dir)" STEERABLE=always
+{
+  "^$" :: "Empty: use IO::recover_dir"
+  ".+" :: "Not empty: directory name"
+} ""
+
 BOOLEAN checkpoint "Checkpoint" STEERABLE=always
 {
 } "no"
@@ -110,6 +131,12 @@ BOOLEAN checkpoint_next "Checkpoint at next iteration" STEERABLE=always
 {
 } "no"
 
+STRING checkpoint_dir "Checkpoint directory (overrides IO::checkpoint_dir)" STEERABLE=always
+{
+  "^$" :: "Empty: use IO::checkpoint_dir"
+  ".+" :: "Not empty: directory name"
+} ""
+
 
 
 INT compression_level "Compression level" STEERABLE=always
diff --git a/CarpetDev/CarpetIOF5/src/FIELDS b/CarpetDev/CarpetIOF5/src/FIELDS
new file mode 100644
index 000000000..f502f5fd7
--- /dev/null
+++ b/CarpetDev/CarpetIOF5/src/FIELDS
@@ -0,0 +1,18 @@
+Results from test output comparing different number of processes (np),
+contiguous/separated output, and unfragmented/fragmented output. The
+two rightmost columns describe the group hierarchy of a scalar and a
+vector. The first (leftmost) name of these is the F5 field, the last
+(rightmost) name is the actual dataset.
+
+
+
+                       | scalar                     vector
+-----------------------|-----------------------------------------------------------------
+np=1   contig   unfrag | grid.r/Map_0_Component_0   grid.coordinates/Dx
+np=1   contig   frag   | grid.r/Map_0_Component_0   grid.coordinates/Map_0_Component_0/Dx
+np=1   sep      unfrag | grid.r                     grid.coordinates/Dx
+np=1   sep      frag   | grid.r/Map_0_Component_0   grid.coordinates/Map_0_Component_0/Dx
+np=2   contig   unfrag | grid.r/Map_0_Component_0   grid.coordinates/Map_0_Component_0/Dx
+np=2   contig   frag   | grid.r/Map_0_Component_0   grid.coordinates/Map_0_Component_0/Dx
+np=2   sep      unfrag | grid.r/Map_0_Component_0   grid.coordinates/Map_0_Component_0/Dx
+np=2   sep      frag   | grid.r/Map_0_Component_0   grid.coordinates/Map_0_Component_0/Dx
diff --git a/CarpetDev/CarpetIOF5/src/input.cc b/CarpetDev/CarpetIOF5/src/input.cc
index fcb3f7d12..da0d7dc00 100644
--- a/CarpetDev/CarpetIOF5/src/input.cc
+++ b/CarpetDev/CarpetIOF5/src/input.cc
@@ -17,6 +17,8 @@
 #include <F5/F5iterate.h>
 #include <F5/F5uniform.h>
 
+#include "CactusBase/IOUtil/src/ioutil_CheckpointRecovery.h"
+
 #include <bbox.hh>
 #include <vect.hh>
 
@@ -26,6 +28,29 @@
 
 
 
+// F5 helper
+namespace {
+  char const*
+  ArrayTypeName(ArrayType const array_type)
+  {
+    switch (array_type) {
+    case UnknownArrayType:            return "Unknown";
+    case Contiguous:                  return "Contiguous";
+    case SeparatedCompound:           return "SeparatedCompound";
+    case Constant:                    return "Constant";
+    case FragmentedContiguous:        return "FragmentedContiguous";
+    case FragmentedSeparatedCompound: return "FragmentedSeparatedCompound";
+    case DirectProduct:               return "DirectProduct";
+    case IndexPermutation:            return "IndexPermutation";
+    case UniformSampling:             return "UniformSampling";
+    case FragmentedUniformSampling:   return "FragmentedUniformSampling";
+    default:                          return "(illegal)";
+    }
+  }
+}
+
+
+
 namespace CarpetIOF5 {
   
   using namespace std;
@@ -36,26 +61,37 @@ namespace CarpetIOF5 {
   // Use a class for reading data, so that we have an easy way to pass
   // variables between the various iterators
   class input_iterator_t {
-    cGH *cctkGH;
+    cGH const *cctkGH;
+    
+    vector<bool> const input_var; // whether to input this variable
     
     double time;
     char const *gridname;
     char const *topologyname;
     int index_depth;            // 0=vertex, 1=cell
     int topological_dimension;
+    int reflevel;
     char const *fieldname;
+    ArrayType fieldtype;   // determines how to interpret the fragment
+    int varindex;
     char const *fragmentname;
+    int map, component;
     
     // string chartname;
     
   public:
-    input_iterator_t (cGH *const cctkGH_)
+    input_iterator_t (cGH const *const cctkGH_,
+                      vector<bool> const& input_var_)
       : cctkGH(cctkGH_),
+        input_var(input_var_),
         time(nan),
         gridname(NULL),
         topologyname(NULL), index_depth(-1), topological_dimension(-1),
-        fieldname(NULL),
-        fragmentname(NULL)
+        reflevel(-1),
+        fieldname(NULL), fieldtype(UnknownArrayType),
+        varindex(-1),
+        fragmentname(NULL),
+        map(-1), component(-1)
     {
     }
     
@@ -101,6 +137,7 @@ namespace CarpetIOF5 {
         BEGIN_LOCAL_MAP_LOOP(cctkGH, CCTK_GF) {
           BEGIN_LOCAL_COMPONENT_LOOP(cctkGH, CCTK_GF) {
             DECLARE_CCTK_ARGUMENTS;
+#pragma omp parallel
             CCTK_LOOP3_BOUNDARIES(boundaries, cctkGH,
                                   i,j,k, 1,1,1,
                                   cctk_lsh[0]-1,cctk_lsh[1]-1,cctk_lsh[2]-1)
@@ -130,8 +167,13 @@ namespace CarpetIOF5 {
       
       herr_t herr;
       
-      cout << indent << "topology=\"" << topologyname << "\""
-           << " (" << (index_depth==0 ? "vertex" : "cell") << ")\n";
+      cout << indent
+           << "topology=\"" << topologyname << "\""
+           << " (" << (index_depth==0 ? "vertex" : "cell") << ")\n"
+           << indent
+           << "topological dimension=" << topological_dimension << "\n";
+      
+      
       
       // Ignore topologies that are only an alias for another topology
       H5G_stat_t stat;
@@ -148,6 +190,29 @@ namespace CarpetIOF5 {
         return;
       }
       
+      
+      
+      // Determine refinement level from the topology
+      // (Could also use topology name instead)
+      hsize_t hreffact[FIBER_MAX_RANK];
+      int const iret =
+        F5LAget_dimensions(path->Topology_hid,
+                           FIBER_HDF5_REFINEMENT_INFO, hreffact);
+      assert (iret == dim);
+      hsize_t hreffact2[FIBER_MAX_RANK];
+      void *const pret = F5Tpermute_dimensions(path, dim, hreffact2, hreffact);
+      assert (pret);
+      ivect const reffact = h2v(hreffact2);
+      int rl;
+      for (rl=0; rl<reflevels; ++rl) {
+        if (all(reffact == Carpet::spacereffacts.AT(rl))) break;
+      }
+      assert (rl<reflevels);
+      reflevel = rl;
+      cout << indent << "refinement level " << reflevel << "\n";
+      
+      
+      
       // F5iterate_topology_fields
       //   (path, NULL, field_iterator, this, chartname.c_str(), NULL); 
       F5iterate_topology_fields (path, NULL, field_iterator, this, NULL, NULL); 
@@ -158,15 +223,36 @@ namespace CarpetIOF5 {
       indent_t indent;
       cout << indent << "field=\"" << fieldname << "\"\n";
       
-      if (strcmp(fieldname, "GRID::COORDINATES")!=0 and
-          strcmp(fieldname, "GRID::r")!=0)
-      {
+      interpret_fieldname(cctkGH, fieldname, varindex);
+#warning "TODO: check all variables in the group"
+#warning "TODO: loop over all variables in the group"
+      if (varindex>=0 and input_var.at(varindex)) {
+
+        int major_version, minor_version, release_version;
+        fieldtype =
+          F5Fget_field_enum(path, fieldname,
+                            &major_version, &minor_version, &release_version);
+        cout << indent << "field_type=" << ArrayTypeName(fieldtype) << "\n";
+        
+        int const iret = F5Fopen(path, fieldname);
+        assert(iret);
+        
+        // Do we need to iterate over fragments?
+        int const is_fragmented = F5Fis_fragmented(path, fieldname);
+        cout << indent
+             << (is_fragmented ? "fragmented" : "not fragmented") << "\n";
+        if (is_fragmented) {
+          F5iterate_field_fragments (path, NULL, fragment_iterator, this);
+        } else {
+          read_fragment (path);
+        }
+        
+        F5Fclose(path);
+        
+      } else {
         indent_t indent;
         cout << indent << "ignoring this field\n";
-        return;
       }
-      
-      F5iterate_field_fragments (path, NULL, fragment_iterator, this);
     }
     
     void read_fragment (F5Path *const path)
@@ -176,55 +262,57 @@ namespace CarpetIOF5 {
       int iret;
       void *pret;
       
-      cout << indent << "fragment=\"" << fragmentname << "\"\n";
+      if (fragmentname) {
+        cout << indent << "fragment=\"" << fragmentname << "\"\n";
+      } else {
+        cout << indent << "no fragment\n";
+      }
       
+      int major_version, minor_version, release_version;
+      fieldtype =
+        F5Fget_field_enum(path, fieldname,
+                          &major_version, &minor_version, &release_version);
+      cout << indent << "field_type=" << ArrayTypeName(fieldtype) << "\n";
       
       
-      // Determine refinement level from the topology
-      // (Could also use topology name instead)
-      hsize_t hreffact[FIBER_MAX_RANK];
-      iret = F5LAget_dimensions(path->Topology_hid,
-                                FIBER_HDF5_REFINEMENT_INFO, hreffact);
-      assert (iret == dim);
-      hsize_t hreffact2[FIBER_MAX_RANK];
-      pret = F5Tpermute_dimensions(path, dim, hreffact2, hreffact);
-      assert (pret);
-      ivect const reffact = h2v(hreffact2);
-      int rl;
-      for (rl=0; rl<reflevels; ++rl) {
-        if (all(reffact == Carpet::spacereffacts.AT(rl))) break;
-      }
-      assert (rl<reflevels);
       
       // Determine map and component from fragment name
-      int m, c;
-      interpret_fragmentname (cctkGH, fragmentname, m, c);
-      
-      {
-        indent_t indent;
-        cout << indent << "reading refinement level " << rl << "\n";
-        cout << indent << "reading map " << m << "\n";
+      if (fragmentname) {
+        interpret_fragmentname (cctkGH, fragmentname, map, component);
+      } else {
+        map = 0;
+        component = CCTK_MyProc(cctkGH);
+      }
+      cout << indent << "map " << map << " component " << component << "\n";
+      
+      // This routine has a bug, and does not recognise
+      // "SeparatedCompound" as separated
+      //int const is_separated = F5Fis_separatedcompound(path, fieldname);
+      int const is_separated =
+        fieldtype == FragmentedSeparatedCompound or
+        fieldtype == SeparatedCompound;
+      cout << indent
+           << (is_separated ? "separated" : "contiguous") << "\n";
+      hid_t const type_id = F5Fget_type(path);
+      if (F5Tis_convertible(type_id, H5T_NATIVE_DOUBLE)) {
+        cout << indent << "compound type: scalar\n";
+        read_variable (path, NULL, varindex);
+      } else if (F5Tis_convertible(type_id, F5T_VEC3_DOUBLE)) {
+        cout << indent << "compound type: vector\n";
+        // This assumes a separated vector; we don't support
+        // contiguous vectors
+        assert(is_separated);
+        read_variable (path, "Dx", varindex+0);
+        read_variable (path, "Dy", varindex+1);
+        read_variable (path, "Dz", varindex+2);
+      } else if (F5Tis_convertible(type_id, F5T_METRIC33_DOUBLE)) {
+        cout << indent << "compound type: symmetric tensor\n";
+        // Not yet implemented
+        assert(0);
+      } else {
+        // Unknown tensor type
+        assert(0);
       }
-      
-      ENTER_LEVEL_MODE(cctkGH, rl) {
-        ENTER_SINGLEMAP_MODE(cctkGH, m, CCTK_GF) {
-          
-          if (strcmp(fieldname, "GRID::r") == 0) {
-            
-            read_variable (path, fragmentname, CCTK_VarIndex("grid::r"));
-            
-          } else if (strcmp(fieldname, "GRID::COORDINATES") == 0) {
-            
-            read_variable (path, "Dx", CCTK_VarIndex("grid::x"));
-            read_variable (path, "Dy", CCTK_VarIndex("grid::y"));
-            read_variable (path, "Dz", CCTK_VarIndex("grid::z"));
-            
-          } else {
-            // do nothing
-          }
-          
-        } LEAVE_SINGLEMAP_MODE;
-      } LEAVE_LEVEL_MODE;
     }
     
     void read_variable (F5Path *const path, char const *const name,
@@ -238,36 +326,83 @@ namespace CarpetIOF5 {
       void *pret;
       
       assert (path);
-      assert (name);
       assert (var >= 0);
       
-      cout << indent << "dataset=\"" << name << "\"\n";
+      cout << indent << "dataset=\"" << (name ? name : "(null)") << "\"\n";
       
       assert (var >= 0);
       {
         char *const fullname = CCTK_FullName(var);
-        cout << indent << "variable=" << fullname << "\n";
+        cout << indent << "variable=\"" << fullname << "\"\n";
         free (fullname);
       }
       
       
       
       // Determine fragment properties
-      H5O_info_t info;
-      herr =
-        H5Oget_info_by_name (path->Field_hid, fragmentname, &info, H5P_DEFAULT);
-      assert (not herr);
-      bool const fragment_is_group = info.type == H5O_TYPE_GROUP;
-      hid_t field;
+      
+      // A dataset can either be the same as a field, or can be in a
+      // group containing all fragments, and/or can be in a group
+      // containing all compound (vector) elements.
+      // The possible group hierarchies are (the rightmost element is
+      // the dataset):
+      //    .../FIELD
+      //    .../FIELD/ELEMENT
+      //    .../FIELD/FRAGMENT
+      //    .../FIELD/FRAGMENT/ELEMENT
+      // The field is the Cactus group/variable name, a fragment
+      // describes map and component, the element describes the tensor
+      // element.
+      // We describe this structure with two booleans:
+      //    bool fragment_is_group: the fragment is a group
+      //    bool field_is_dataset:  the field is a dataset
+      // Note that the F5 library will already have openend the field,
+      // either as group or as dataset.
+      
+      // There is a group for all fragments if (a) we expect a
+      // fragment, and (b) there exists a group with this name. (We
+      // probably could examine F5's internal state as well, but
+      // looking for an HDF5 group is the most robust approach.)
+      bool fragment_is_group = false;
+      if (fragmentname) {
+        H5O_info_t info;
+        herr = H5Oget_info_by_name (path->Field_hid, fragmentname,
+                                    &info, H5P_DEFAULT);
+        assert (not herr);
+        fragment_is_group = info.type == H5O_TYPE_GROUP;
+      }
+      cout << indent << "fragment_is_group=" << fragment_is_group << "\n";
+      hid_t fragment;
       if (fragment_is_group) {
-        field = H5Gopen (path->Field_hid, fragmentname, H5P_DEFAULT);
-        assert (field >= 0);
+        fragment = H5Gopen (path->Field_hid, fragmentname, H5P_DEFAULT);
+        assert (fragment >= 0);
       } else {
-        field = path->Field_hid;
+        fragment = path->Field_hid;
       }
       
-      hid_t const fragment = H5Dopen(field, name, H5P_DEFAULT);
-      assert (fragment);
+      // The field consists of a single dataset if we cannot open the
+      // element
+      hid_t element = -1;
+      char const *datasetname = NULL;
+      if (name) {
+        datasetname = name;
+      } else if (fragmentname and not fragment_is_group) {
+        datasetname = fragmentname;
+      }
+      if (datasetname) {
+        H5E_BEGIN_TRY {
+          element = H5Dopen(fragment, datasetname, H5P_DEFAULT);
+        } H5E_END_TRY;
+        assert (element >= 0);
+      }
+      bool const field_is_dataset = element < 0;
+      cout << indent << "field_is_dataset=" << field_is_dataset << "\n";
+      if (field_is_dataset) {
+        // Can't open the fragment -- we assume the field consists of
+        // a single element and has already been opened
+        element = fragment;
+      }
+      assert (element>=0);
       
       // Check index depth
       int index_depth_;
@@ -277,14 +412,17 @@ namespace CarpetIOF5 {
       
       // Read the fragment offset. This is stored with the dataset
       // group.
-      hsize_t hoff[FIBER_MAX_RANK];
-      iret = F5LAget_dimensions(fragment_is_group ? field : fragment,
-                                FIBER_FRAGMENT_OFFSET_ATTRIBUTE, hoff);
-      assert (iret == dim);
-      hsize_t hoff2[FIBER_MAX_RANK];
-      pret = F5Tpermute_dimensions(path, dim, hoff2, hoff);
-      assert (pret);
-      ivect const foff = h2v(hoff2);
+      ivect foff = 0;
+      if (fragmentname) {
+        hsize_t hoff[FIBER_MAX_RANK];
+        iret = F5LAget_dimensions(fragment_is_group ? fragment : element,
+                                  FIBER_FRAGMENT_OFFSET_ATTRIBUTE, hoff);
+        assert (iret == dim);
+        hsize_t hoff2[FIBER_MAX_RANK];
+        pret = F5Tpermute_dimensions(path, dim, hoff2, hoff);
+        assert (pret);
+        foff = h2v(hoff2);
+      }
       assert (all(foff>=0));
       
 #if 0
@@ -301,7 +439,7 @@ namespace CarpetIOF5 {
       ivect const flen = h2v(hlen2);
       assert (all(flen>=0));
 #endif
-      hid_t const space = H5Dget_space(fragment);
+      hid_t const space = H5Dget_space(element);
       assert (space>=0);
       iret = H5Sget_simple_extent_ndims(space);
       assert (iret == dim);
@@ -318,7 +456,8 @@ namespace CarpetIOF5 {
       ibbox const fbox (foff, foff+flen-1, 1);
       {
         indent_t indent;
-        cout << indent << "dataset bbox is " << foff << ":" << foff+flen << "\n";
+        cout << indent
+             << "dataset bbox is " << foff << ":" << foff+flen-1 << "\n";
       }
       
       
@@ -336,111 +475,117 @@ namespace CarpetIOF5 {
       assert (groupdata.stagtype == 0);
       assert (groupdata.dim == dim);
       
-      // TODO: This is too expensive; only traverse Carpet's data
-      // structures instead. (This should be implemented via a
-      // gh::locate_positions, which will turn require support from
-      // the internal tree data structure.)
-      BEGIN_COMPONENT_LOOP(cctkGH, CCTK_GF) {
-        DECLARE_CCTK_ARGUMENTS;
-        
-        ivect lbnd, lsh, lghosts, ughosts;
-        ivect imin, imax, ioff, ilen;
-        for (int d=0; d<dim; ++d) {
-          lbnd[d]    = cctk_lbnd[d];
-          lsh[d]     = cctk_lsh[d];
-          // F5 counts the total overlap, which is the sum of the
-          // ghost zones on this and the adjacent component
-          lghosts[d] = cctk_bbox[2*d  ] ? 0 : 2*cctk_nghostzones[d];
-          ughosts[d] = cctk_bbox[2*d+1] ? 0 : 2*cctk_nghostzones[d];
-          imin[d] = 0;
-          imax[d] = lsh[d];
-          // Do not read in ghost zones
-          imin[d] += lghosts[d] / 2;
-          imax[d] -= ughosts[d] / 2;
-          ioff[d] = lbnd[d] + imin[d];
-          ilen[d] = imax[d] - imin[d];
-        }
-        ibbox const lbox (lbnd, lbnd+lsh-1, 1);
-        ibbox const ibox (ioff, ioff+ilen-1, 1);
-        ibbox const ovlp = ibox & fbox;
-        
-        if (ovlp.empty()) {
-          indent_t indent;
-          cout << indent << "dataset does not intersect component " << component << "; skipping component\n";
-          continue;
-        }
-        
-        indent_t indent;
-        cout << indent << "dataset intersects component " << component << "\n"
-             << indent << "reading bbox " << ovlp.lower() << ":" << ovlp.upper() << " of " << ibox.lower() << ":" << ibox.upper() << "\n";
-        
-        int const proc = vhh.AT(Carpet::map)->processor(reflevel, component);
+      // TODO: This does not work; we need to read the data, and later
+      // distribute it onto the respective Cactus variables
+      ENTER_LEVEL_MODE(cctkGH, reflevel) {
+        ENTER_SINGLEMAP_MODE(cctkGH, map, CCTK_GF) {
+          ENTER_LOCAL_MODE(cctkGH, component, CCTK_GF) {
+            DECLARE_CCTK_ARGUMENTS;
+            
+            ivect lbnd, lsh, lghosts, ughosts;
+            ivect imin, imax, ioff, ilen;
+            for (int d=0; d<dim; ++d) {
+              lbnd[d]    = cctk_lbnd[d];
+              lsh[d]     = cctk_lsh[d];
+              // F5 counts the total overlap, which is the sum of the
+              // ghost zones on this and on the adjacent component
+              lghosts[d] = cctk_bbox[2*d  ] ? 0 : 2*cctk_nghostzones[d];
+              ughosts[d] = cctk_bbox[2*d+1] ? 0 : 2*cctk_nghostzones[d];
+              imin[d] = 0;
+              imax[d] = lsh[d];
+              // Do not read in ghost zones
+              imin[d] += lghosts[d] / 2;
+              imax[d] -= ughosts[d] / 2;
+              ioff[d] = lbnd[d] + imin[d];
+              ilen[d] = imax[d] - imin[d];
+            }
+            ibbox const lbox (lbnd, lbnd+lsh-1, 1);
+            ibbox const ibox (ioff, ioff+ilen-1, 1);
+            ibbox const ovlp = ibox & fbox;
+            
+            if (ovlp.empty()) {
+              indent_t indent;
+              cout << indent << "dataset does not intersect component " << component << "; skipping component\n";
+              continue;
+            }
+            
+            indent_t indent;
+            cout << indent << "dataset intersects component " << component << "\n"
+                 << indent << "reading bbox " << ovlp.lower() << ":" << ovlp.upper() << " of " << ibox.lower() << ":" << ibox.upper() << "\n";
+            
+            #warning "TODO: assert the whole component is read"
+            
+            int const proc = vhh.AT(Carpet::map)->processor(reflevel, component);
 #warning "TODO: send dataset to destination process"
-        assert (proc == dist::rank());
-        
-        int const timelevel = 0;
-        void *const data = CCTK_VarDataPtrI(cctkGH, timelevel, var);
-        assert (data);
-        
-        
-        
-        hvect const hzero(0);
-        hvect mlsh;
-        pret = F5Tpermute_dimensions
-          (path, dim, &mlsh[0], &v2h(lbox.shape())[0]);
-        assert (pret);
-        hid_t const memspace = H5Screate_simple(dim, &mlsh[0], NULL);
-        assert (memspace >= 0);
-        hvect mmin;
-        pret = F5Tpermute_dimensions
-          (path, dim, &mmin[0], &v2h(ovlp.lower()-lbox.lower())[0]);
-        assert (pret);
-        hvect mlen;
-        pret = F5Tpermute_dimensions
-          (path, dim, &mlen[0], &v2h(ovlp.shape())[0]);
-        assert (pret);
-        assert (all(mmin >= hzero));
-        assert (all(mmin+mlen <= mlsh));
-        herr = H5Sselect_hyperslab
-          (memspace, H5S_SELECT_SET, &mmin[0], NULL, &mlen[0], NULL);
-        assert (not herr);
-        // cout << "mlsh=" << mlsh << " mmin=" << mmin << " mlen=" << mlen << "\n";
-        
-        hvect flsh;
-        pret = F5Tpermute_dimensions
-          (path, dim, &flsh[0], &v2h(fbox.shape())[0]);
-        assert (pret);
-        hid_t const filespace = H5Screate_simple(dim, &flsh[0], NULL);
-        assert (filespace >= 0);
-        hvect fmin;
-        pret = F5Tpermute_dimensions
-          (path, dim, &fmin[0], &v2h(ovlp.lower()-fbox.lower())[0]);
-        assert (pret);
-        hvect const flen = mlen;
-        assert (all(fmin >= hzero));
-        assert (all(fmin+flen <= flsh));
-        herr = H5Sselect_hyperslab
-          (filespace, H5S_SELECT_SET, &fmin[0], NULL, &flen[0], NULL);
-        // cout << "flsh=" << flsh << " fmin=" << fmin << " flen=" << flen << "\n";
-        
-        herr = H5Dread (fragment, H5T_NATIVE_DOUBLE, memspace, filespace,
-                        H5P_DEFAULT, data);
-        assert (not herr);
-        
-        herr = H5Sclose(memspace);
-        assert (not herr);
-        herr = H5Sclose(filespace);
-        assert (not herr);
-        
-      } END_COMPONENT_LOOP;
+            assert (proc == dist::rank());
+            
+            int const timelevel = 0;
+            void *const data = CCTK_VarDataPtrI(cctkGH, timelevel, var);
+            assert (data);
+            
+            
+            
+            hvect const hzero(0);
+            hvect mlsh;
+            pret = F5Tpermute_dimensions
+              (path, dim, &mlsh[0], &v2h(lbox.shape())[0]);
+            assert (pret);
+            hid_t const memspace = H5Screate_simple(dim, &mlsh[0], NULL);
+            assert (memspace >= 0);
+            hvect mmin;
+            pret = F5Tpermute_dimensions
+              (path, dim, &mmin[0], &v2h(ovlp.lower()-lbox.lower())[0]);
+            assert (pret);
+            hvect mlen;
+            pret = F5Tpermute_dimensions
+              (path, dim, &mlen[0], &v2h(ovlp.shape())[0]);
+            assert (pret);
+            assert (all(mmin >= hzero));
+            assert (all(mmin+mlen <= mlsh));
+            herr = H5Sselect_hyperslab
+              (memspace, H5S_SELECT_SET, &mmin[0], NULL, &mlen[0], NULL);
+            assert (not herr);
+            // cout << "mlsh=" << mlsh << " mmin=" << mmin << " mlen=" << mlen << "\n";
+            
+            hvect flsh;
+            pret = F5Tpermute_dimensions
+              (path, dim, &flsh[0], &v2h(fbox.shape())[0]);
+            assert (pret);
+            hid_t const filespace = H5Screate_simple(dim, &flsh[0], NULL);
+            assert (filespace >= 0);
+            hvect fmin;
+            pret = F5Tpermute_dimensions
+              (path, dim, &fmin[0], &v2h(ovlp.lower()-fbox.lower())[0]);
+            assert (pret);
+            hvect const flen = mlen;
+            assert (all(fmin >= hzero));
+            assert (all(fmin+flen <= flsh));
+            herr = H5Sselect_hyperslab
+              (filespace, H5S_SELECT_SET, &fmin[0], NULL, &flen[0], NULL);
+            // cout << "flsh=" << flsh << " fmin=" << fmin << " flen=" << flen << "\n";
+            
+            herr = H5Dread (element, H5T_NATIVE_DOUBLE, memspace, filespace,
+                            H5P_DEFAULT, data);
+            assert (not herr);
+            
+            herr = H5Sclose(memspace);
+            assert (not herr);
+            herr = H5Sclose(filespace);
+            assert (not herr);
+            
+          } LEAVE_LOCAL_MODE;
+        } LEAVE_SINGLEMAP_MODE;
+      } LEAVE_LEVEL_MODE;
       
       
       
-      herr = H5Dclose(fragment);
-      assert (not herr);
+      if (not field_is_dataset) {
+        herr = H5Dclose(element);
+        assert (not herr);
+      }
       
       if (fragment_is_group) {
-        herr = H5Gclose(field);
+        herr = H5Gclose(fragment);
         assert (not herr);
       }
     }
@@ -458,17 +603,19 @@ namespace CarpetIOF5 {
     herr_t timeslice_iterator (F5Path *const path, double const time,
                                void *const userdata)
     {
-      input_iterator_t* const iterator = (input_iterator_t*)userdata;
+      input_iterator_t *const iterator = (input_iterator_t*)userdata;
       iterator->time = time;
       iterator->read_timeslice (path);
       return 0;
+      // finish (successfully) after the first timeslice
+      // return 1;
     }
     
     static
     herr_t grid_iterator (F5Path *const path, char const *const gridname,
                           void *const userdata)
     {
-      input_iterator_t* const iterator = (input_iterator_t*)userdata;
+      input_iterator_t *const iterator = (input_iterator_t*)userdata;
       iterator->gridname = gridname;
       iterator->read_grid (path);
       return 0;
@@ -481,11 +628,14 @@ namespace CarpetIOF5 {
                               int const topological_dimension,
                               void *const userdata)
     {
-      input_iterator_t* const iterator = (input_iterator_t*)userdata;
+      input_iterator_t *const iterator = (input_iterator_t*)userdata;
       iterator->topologyname = topologyname;
       iterator->index_depth = index_depth;
       iterator->topological_dimension = topological_dimension;
       iterator->read_topology (path);
+      iterator->topologyname = NULL;
+      iterator->index_depth = -1;
+      iterator->topological_dimension = -1;
       return 0;
     }
     
@@ -493,9 +643,10 @@ namespace CarpetIOF5 {
     herr_t field_iterator (F5Path *const path, char const *const fieldname,
                            void *const userdata)
     {
-      input_iterator_t* const iterator = (input_iterator_t*)userdata;
+      input_iterator_t *const iterator = (input_iterator_t*)userdata;
       iterator->fieldname = fieldname;
       iterator->read_field (path);
+      iterator->fieldname = NULL;
       return 0;
     }
     
@@ -504,9 +655,10 @@ namespace CarpetIOF5 {
                               char const *const fragmentname,
                               void *const userdata)
     {
-      input_iterator_t* const iterator = (input_iterator_t*)userdata;
+      input_iterator_t *const iterator = (input_iterator_t*)userdata;
       iterator->fragmentname = fragmentname;
       iterator->read_fragment (path);
+      iterator->fragmentname = NULL;
       return 0;
     }
     
@@ -514,50 +666,12 @@ namespace CarpetIOF5 {
   
   
   
-  extern "C"
-  void F5_Input (CCTK_ARGUMENTS)
+  void input (cGH const *const cctkGH,
+               hid_t const file,
+               vector<bool> const& input_var)
   {
-    DECLARE_CCTK_ARGUMENTS;
-    DECLARE_CCTK_PARAMETERS;
-    
-    herr_t herr;
-    
-    
-    
-    assert (is_global_mode());
-    CCTK_VInfo (CCTK_THORNSTRING, "F5_Input: iteration=%d", cctk_iteration);
-    
-    
-    
-    // Open file
-    string const basename =
-      generate_basename (cctkGH, CCTK_VarIndex("grid::r"));
-    int const myproc = CCTK_MyProc(cctkGH);
-    int const nprocs = CCTK_nProcs(cctkGH);
-    for (int proc=myproc; ; proc+=nprocs) {
-      string const name =
-        create_filename (cctkGH, basename, proc, false);
-      
-      bool file_exists;
-      H5E_BEGIN_TRY {
-        file_exists = H5Fis_hdf5(name.c_str()) > 0;
-      } H5E_END_TRY;
-      if (not file_exists) break;
-      
-      indent_t indent;
-      cout << indent << "process=" << proc << "\n";
-      
-      hid_t const file = H5Fopen (name.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
-      assert (file >= 0);
-      
-      // Iterate over all time slices
-      input_iterator_t iterator(cctkGH);
-      iterator.iterate (file);
-      
-      // Close file
-      herr = H5Fclose (file);
-      assert (not herr);
-    }
+    input_iterator_t iterator(cctkGH, input_var);
+    iterator.iterate (file);
   }
   
 } // end namespace CarpetIOF5
diff --git a/CarpetDev/CarpetIOF5/src/iof5.cc b/CarpetDev/CarpetIOF5/src/iof5.cc
index 8b93054d6..388708c6e 100644
--- a/CarpetDev/CarpetIOF5/src/iof5.cc
+++ b/CarpetDev/CarpetIOF5/src/iof5.cc
@@ -15,6 +15,7 @@
 #include <hdf5.h>
 #include <iof5.hh>
 
+#include "CactusBase/IOUtil/src/ioGH.h"
 #include "CactusBase/IOUtil/src/ioutil_CheckpointRecovery.h"
 
 
@@ -49,6 +50,9 @@ namespace CarpetIOF5 {
     DECLARE_CCTK_PARAMETERS;
     
     // register I/O method
+    int const ierr = IOUtil_RegisterRecover("CarpetIOF5 recovery", Input);
+    assert(not ierr);
+    
     int const IOMethod = CCTK_RegisterIOMethod ("IOF5");
     CCTK_RegisterIOMethodOutputGH      (IOMethod, OutputGH     );
     CCTK_RegisterIOMethodTimeToOutput  (IOMethod, TimeToOutput );
@@ -182,7 +186,8 @@ namespace CarpetIOF5 {
       string const basename = generate_basename (cctkGH, vindex);
       int const myproc = CCTK_MyProc(cctkGH);
       int const proc = myproc;
-      string const name = create_filename (cctkGH, basename, proc, first_time);
+      string const name =
+        create_filename (cctkGH, basename, proc, io_dir_output, first_time);
       
       indent_t indent;
       cout << indent << "process=" << proc << "\n";
@@ -216,6 +221,7 @@ namespace CarpetIOF5 {
   {
     assert (is_global_mode());
     
+#if 0
     // generate filenames for both the temporary and real checkpoint
     // files
     int const ioproc = CCTK_MyProc(cctkGH);
@@ -226,9 +232,16 @@ namespace CarpetIOF5 {
     char* const tempname =
       IOUtil_AssembleFilename (cctkGH, NULL, ".tmp", ".f5",
                                called_from, ioproc, not parallel_io);
+#endif
+    int const myproc = CCTK_MyProc(cctkGH);
+    int const proc = myproc;
+    string const name =
+      create_filename (cctkGH, "checkpoint", proc, io_dir_checkpoint, true);
+    string const tempname =
+      create_filename (cctkGH, "checkpoint.tmp", proc, io_dir_checkpoint, true);
     
     hid_t const file =
-      H5Fcreate (tempname, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
+      H5Fcreate (tempname.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
     assert (file >= 0);
     
     vector<bool> output_var(CCTK_NumVars());
@@ -266,7 +279,7 @@ namespace CarpetIOF5 {
     // checkpoint files
     // TODO: ensure there were no errors
     CCTK_Barrier (cctkGH);
-    int const ierr = rename (tempname, filename);
+    int const ierr = rename (tempname.c_str(), name.c_str());
     assert (not ierr);
   }
   
@@ -312,4 +325,113 @@ namespace CarpetIOF5 {
     }
   }
   
+  
+  
+  int Input (cGH* const cctkGH,
+             char const* const basefilename, int const called_from)
+  {
+    DECLARE_CCTK_PARAMETERS;
+    
+    herr_t herr;
+    
+    
+    
+    assert (is_level_mode());
+    BEGIN_GLOBAL_MODE(cctkGH) {
+      DECLARE_CCTK_ARGUMENTS;
+      
+      CCTK_VInfo (CCTK_THORNSTRING, "F5::Input: iteration=%d", cctk_iteration);
+      
+      
+      
+      assert (called_from == CP_RECOVER_PARAMETERS or
+              called_from == CP_RECOVER_DATA or
+              called_from == FILEREADER_DATA);
+      bool const in_recovery =
+        called_from == CP_RECOVER_PARAMETERS or
+        called_from == CP_RECOVER_DATA;
+      
+      // We don't know how to do this yet
+      assert (called_from != CP_RECOVER_PARAMETERS);
+      
+      // Determine which variables to read
+      ioGH const* const ioUtilGH =
+        (ioGH const*) CCTK_GHExtension (cctkGH, "IO");
+      vector<bool> input_var(CCTK_NumVars(), true);
+      if (ioUtilGH->do_inVars) {
+        for (int n=0; n<CCTK_NumVars(); ++n) {
+          input_var.at(n) = ioUtilGH->do_inVars[n];
+        }
+      }
+      
+      
+      
+      // Open file
+      // string const basename =
+      //   in_recovery
+      //   ? "checkpoint"
+      //   : generate_basename (cctkGH, CCTK_VarIndex("grid::r"));
+      string const basename = basefilename;
+      
+      // Keep track of which files could be read, and which could not
+      int foundproc = -1, notfoundproc = -1;
+      
+#warning "TODO: Store how many processes contributed to the output, and expect exactly that many files"
+      int const myproc = CCTK_MyProc(cctkGH);
+      int const nprocs = CCTK_nProcs(cctkGH);
+      // Loop over all (possible) files
+      for (int proc=myproc; ; proc+=nprocs) {
+        string const name =
+          create_filename (cctkGH, basename, proc,
+                           in_recovery ? io_dir_recover : io_dir_input, false);
+        
+        bool file_exists;
+        H5E_BEGIN_TRY {
+          file_exists = H5Fis_hdf5(name.c_str()) > 0;
+        } H5E_END_TRY;
+        if (not file_exists) {
+          notfoundproc = proc;
+          break;
+        }
+        foundproc = proc;
+        
+        indent_t indent;
+        cout << indent << "process=" << proc << "\n";
+        
+        hid_t const file = H5Fopen (name.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
+        assert (file >= 0);
+        
+        // Iterate over all time slices
+        input (cctkGH, file, input_var);
+        
+        // Close file
+        herr = H5Fclose (file);
+        assert (not herr);
+      }
+      
+      {
+        int maxfoundproc;
+        MPI_Allreduce(&foundproc, &maxfoundproc, 1, MPI_INT, MPI_MAX,
+                      dist::comm());
+        if (maxfoundproc == -1) {
+          string const name =
+            create_filename (cctkGH, basename, notfoundproc,
+                             in_recovery ? io_dir_recover : io_dir_input,
+                             false);
+          CCTK_VWarn(CCTK_WARN_ALERT, __LINE__, __FILE__, CCTK_THORNSTRING,
+                     "Could not read input file \"%s\"", name.c_str());
+          return 1;
+        }
+        if (notfoundproc <= maxfoundproc) {
+          CCTK_VWarn(CCTK_WARN_ALERT, __LINE__, __FILE__, CCTK_THORNSTRING,
+                     "Could not read file of process %d (but could read file of process %d)",
+                     notfoundproc, maxfoundproc);
+        }
+      }
+      
+    } END_GLOBAL_MODE;
+    
+    return 0;                   // no error
+  }
+  
 } // end namespace CarpetIOF5
diff --git a/CarpetDev/CarpetIOF5/src/iof5.hh b/CarpetDev/CarpetIOF5/src/iof5.hh
index 76837a1e0..7e54eb9c9 100644
--- a/CarpetDev/CarpetIOF5/src/iof5.hh
+++ b/CarpetDev/CarpetIOF5/src/iof5.hh
@@ -63,11 +63,15 @@ namespace CarpetIOF5 {
   
   // Indentation
   class indent_t {
+    static bool const debug = false;
     static int const width = 3;
     static int level;
+  private:
+    indent_t(indent_t const&);
+    indent_t& operator=(indent_t const&);
   public:
-    indent_t() { ++level; }
-    ~indent_t() { --level; }
+    indent_t();
+    ~indent_t();
     ostream& output(ostream& os) const;
   };
   
@@ -183,10 +187,13 @@ namespace CarpetIOF5 {
                      int const variable);
   
   // Create the final file name on a particular processor
+  enum io_dir_t
+    {io_dir_input, io_dir_output, io_dir_recover, io_dir_checkpoint};
   string
   create_filename (cGH const* const cctkGH,
                    string const basename,
                    int const proc,
+                   io_dir_t const io_dir,
                    bool const create_directories);
   
   // Generate a good grid name (simulation name)
@@ -215,6 +222,8 @@ namespace CarpetIOF5 {
   string
   generate_fieldname (cGH const* const cctkGH,
                       int const vi, tensortype_t const tt);
+  void
+  interpret_fieldname(cGH const *const cctkGH, string fieldname, int& vi);
   
   
   
@@ -233,6 +242,10 @@ namespace CarpetIOF5 {
                vector<bool> const& output_var,
                bool const output_everything);
   
+  void input (cGH const* const cctkGH,
+              hid_t const file,
+              vector<bool> const& input_var);
+  
   
   
   // Scheduled routines
@@ -248,7 +261,9 @@ namespace CarpetIOF5 {
   void* SetupGH (tFleshConfig* const fleshconfig,
                  int const convLevel, cGH* const cctkGH);
   
-  // Callbacks for CarpetIOHDF5's I/O method
+  // Callbacks for CarpetIOF5's I/O method
+  int Input (cGH* const cctkGH,
+             char const* const basefilename, int const called_from);
   int OutputGH (cGH const* const cctkGH);
   int TimeToOutput (cGH const* const cctkGH, int const vindex);
   int TriggerOutput (cGH const* const cctkGH, int const vindex);
diff --git a/CarpetDev/CarpetIOF5/src/output.cc b/CarpetDev/CarpetIOF5/src/output.cc
index 36ff01335..d376bbfcf 100644
--- a/CarpetDev/CarpetIOF5/src/output.cc
+++ b/CarpetDev/CarpetIOF5/src/output.cc
@@ -43,7 +43,7 @@ namespace CarpetIOF5 {
     // entries
     cGH* const cctkGH;
     
-    vector<bool> const output_var;    // whether to output this variable
+    vector<bool> const output_var; // whether to output this variable
     bool const output_everything;
     bool const is_multipatch;
     
@@ -304,12 +304,13 @@ namespace CarpetIOF5 {
       DECLARE_CCTK_ARGUMENTS;
       indent_t indent;
       
-      cout << indent << "reflevel=" << reflevel << "\n";
-      
       assert (is_level_mode());
       
       ivect const reffact = spacereffacts.AT(reflevel);
       topologyname = generate_topologyname(cctkGH, group_index, reffact);
+      cout << indent
+           << "reflevel=" << reflevel << " "
+           << "topologyname=" << topologyname << "\n";
       
       // Define grid hierarchy
       map_indices_t const mi(cctkGH, group_index);
@@ -389,13 +390,13 @@ namespace CarpetIOF5 {
       indent_t indent;
       bool error_flag = false;
       
-      cout << indent
-           << "component=" << component << " "
-           << "(local_component=" << local_component << ")\n";
-      
       assert (is_local_mode());
       
       fragmentname = generate_fragmentname(cctkGH, Carpet::map, component);
+      cout << indent
+           << "component=" << component << " "
+           << "(local_component=" << local_component << ") "
+           << "fragmentname=" << fragmentname << "\n";
       
       if (group_type == CCTK_GF) {
         // Define coordinates
@@ -613,6 +614,8 @@ namespace CarpetIOF5 {
         assert (0);
       }
       
+      cout << indent << "fieldname=" << name << "\n";
+      
       // Write data
       assert (type >= 0);
       assert (num_comps > 0);
@@ -692,6 +695,7 @@ namespace CarpetIOF5 {
         // Write single-component tensors into non-separated fractions
         // for convenience (could also use a separated compound
         // instead)
+        // TODO: Extent F5 API to allow writing non-fragmented datasets
         FAILWARN
           (F5Fwrite_fraction (path, name.c_str(),
                               ci.dim,
@@ -702,6 +706,8 @@ namespace CarpetIOF5 {
                               &v2h(ci.lghosts)[0], &v2h(ci.ughosts)[0],
                               fragmentname.c_str(), prop));
       } else {
+        int const full_coverage =
+          will_cover_complete_domain and not fragment_contiguous_components;
         FAILWARN
           (F5FSwrite_fraction (path, name.c_str(),
                                ci.dim,
@@ -711,11 +717,15 @@ namespace CarpetIOF5 {
                                &v2h(ci.ioff)[0],
                                &v2h(ci.lghosts)[0], &v2h(ci.ughosts)[0],
                                fragmentname.c_str(), prop,
-                               will_cover_complete_domain));
+                               full_coverage));
       }
       
       for (int d=0; d<num_comps; ++d) {
-        delete[] data[d];
+        switch (vartype) {
+        case CCTK_VARIABLE_INT:  delete[] (CCTK_INT  const*)data[d]; break;
+        case CCTK_VARIABLE_REAL: delete[] (CCTK_REAL const*)data[d]; break;
+        default: assert(0);
+        }
         data[d] = NULL;
       }
       
diff --git a/CarpetDev/CarpetIOF5/src/poison.cc b/CarpetDev/CarpetIOF5/src/poison.cc
index 2305958e1..8dcc9cbdb 100644
--- a/CarpetDev/CarpetIOF5/src/poison.cc
+++ b/CarpetDev/CarpetIOF5/src/poison.cc
@@ -29,6 +29,7 @@ namespace CarpetIOF5 {
       BEGIN_LOCAL_MAP_LOOP(cctkGH, CCTK_GF) {
         BEGIN_LOCAL_COMPONENT_LOOP(cctkGH, CCTK_GF) {
           DECLARE_CCTK_ARGUMENTS;
+#pragma omp parallel
           CCTK_LOOP3_ALL(F5_Poison, cctkGH, i,j,k) {
             int const ind3d = CCTK_GFINDEX3D(cctkGH, i,j,k);
             x[ind3d] = nan;
@@ -56,6 +57,7 @@ namespace CarpetIOF5 {
       BEGIN_LOCAL_MAP_LOOP(cctkGH, CCTK_GF) {
         BEGIN_LOCAL_COMPONENT_LOOP(cctkGH, CCTK_GF) {
           DECLARE_CCTK_ARGUMENTS;
+#pragma omp parallel
           CCTK_LOOP3_ALL(F5_Check, cctkGH, i,j,k) {
             int const ind3d = CCTK_GFINDEX3D(cctkGH, i,j,k);
             assert (not isnan(x[ind3d]));
diff --git a/CarpetDev/CarpetIOF5/src/util.cc b/CarpetDev/CarpetIOF5/src/util.cc
index 15e01fba7..ce3c83f5e 100644
--- a/CarpetDev/CarpetIOF5/src/util.cc
+++ b/CarpetDev/CarpetIOF5/src/util.cc
@@ -46,6 +46,18 @@ namespace CarpetIOF5 {
   // Indentation
   int indent_t::level = 0;
   
+  indent_t::indent_t()
+  {
+    if (debug) cout << *this << "{{{\n";
+    ++level;
+  }
+  
+  indent_t::~indent_t()
+  {
+    --level;
+    if (debug) cout << *this << "}}}\n";
+  }
+  
   ostream& indent_t::output(ostream& os) const
   {
     return os << string(width*level, ' ');
@@ -139,12 +151,36 @@ namespace CarpetIOF5 {
   create_filename (cGH const *const cctkGH,
                    string const basename,
                    int const proc,
+                   io_dir_t const io_dir,
                    bool const create_directories)
   {
     DECLARE_CCTK_PARAMETERS;
     
-    bool const use_IO_out_dir = strcmp(out_dir, "") == 0;
-    string path = use_IO_out_dir ? IO_out_dir : out_dir;
+    char const *IO_dir = NULL;
+    char const *F5_dir = NULL;
+    switch (io_dir) {
+    case io_dir_input:
+      IO_dir = IO_filereader_ID_dir;
+      F5_dir = filereader_ID_dir;
+      break;
+    case io_dir_output:
+      IO_dir = IO_out_dir;
+      F5_dir = out_dir;
+      break;
+    case io_dir_recover:
+      IO_dir = IO_recover_dir;
+      F5_dir = recover_dir;
+      break;
+    case io_dir_checkpoint:
+      IO_dir = IO_checkpoint_dir;
+      F5_dir = checkpoint_dir;
+      break;
+    default:
+      assert(0);
+    }
+    
+    bool const use_IO_dir = strcmp(F5_dir, "") == 0;
+    string path = use_IO_dir ? IO_dir : F5_dir;
     
     if (create_subdirs) {
       {
@@ -308,12 +344,36 @@ namespace CarpetIOF5 {
     transform (name.begin(), name.end(), name.begin(), ::tolower);
     string const sep = "::";
     size_t const pos = name.find(sep);
-    if (pos != string::npos) {
-      name.replace (pos, sep.size(), ".");
-    }
+    assert (pos != string::npos);
+    name.replace (pos, sep.size(), ".");
     return name;
   }
   
+  void
+  interpret_fieldname(cGH const *const cctkGH, string fieldname, int& vi)
+  {
+    string const sep = ".";
+    size_t const pos = fieldname.find(sep);
+    if (pos == string::npos) {
+      // The field name is not a Cactus group or variable
+      vi = -1;
+      return;
+    }
+    fieldname.replace (pos, sep.size(), "::");
+    
+    vi = CCTK_VarIndex(fieldname.c_str());
+    if (vi < 0) {
+      int const gi = CCTK_GroupIndex(fieldname.c_str());
+      if (gi < 0) {
+        CCTK_VWarn(CCTK_WARN_ALERT, __LINE__, __FILE__, CCTK_THORNSTRING,
+                   "Unknown group variable name \"%s\"", fieldname.c_str());
+        return;
+      }
+      vi = CCTK_FirstVarIndexI(gi);
+      assert (vi>=0);
+    }
+  }
+  
   
   
   char const *const grid_structure = "Grid Structure v5";