CarpetIOF5: Provide standard Cactus output interface

Provide standard Cactus output interface.
Begin to implement checkpointing (support for dim!=3 variables still missing).

15 files changed, 1213 insertions, 1011 deletions

diff --git a/CarpetDev/CarpetIOF5/interface.ccl b/CarpetDev/CarpetIOF5/interface.ccl
index 4c09d8b12..fd263229d 100644
--- a/CarpetDev/CarpetIOF5/interface.ccl
+++ b/CarpetDev/CarpetIOF5/interface.ccl
@@ -44,3 +44,8 @@ REQUIRES FUNCTION IO_TruncateOutputFiles
 CCTK_INT FUNCTION MultiPatch_GetSystemSpecification \
     (CCTK_INT OUT maps)
 USES FUNCTION MultiPatch_GetSystemSpecification
+
+
+
+CCTK_INT this_iteration        TYPE=scalar
+CCTK_INT next_output_iteration TYPE=scalar
diff --git a/CarpetDev/CarpetIOF5/par/iof5-checkpoint.par b/CarpetDev/CarpetIOF5/par/iof5-checkpoint.par
new file mode 100644
index 000000000..3306b73e6
--- /dev/null
+++ b/CarpetDev/CarpetIOF5/par/iof5-checkpoint.par
@@ -0,0 +1,51 @@
+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"
+
+CarpetIOF5::checkpoint = yes
+IO::checkpoint_every   = 1024
diff --git a/CarpetDev/CarpetIOF5/par/iof5-multipatch.par b/CarpetDev/CarpetIOF5/par/iof5-multipatch.par
index 71550e788..9af0afa2f 100644
--- a/CarpetDev/CarpetIOF5/par/iof5-multipatch.par
+++ b/CarpetDev/CarpetIOF5/par/iof5-multipatch.par
@@ -43,3 +43,6 @@ CarpetRegrid2::radius_1[1]    =  0.5
 CarpetRegrid2::radius_1[2]    =  0.2
 
 InitBase::initial_data_setup_method = "init_all_levels"
+
+CarpetIOF5::out_every = 1
+CarpetIOF5::out_vars  = "grid::coordinates"
diff --git a/CarpetDev/CarpetIOF5/par/iof5-refined-cell.par b/CarpetDev/CarpetIOF5/par/iof5-refined-cell.par
new file mode 100644
index 000000000..75dc0c17b
--- /dev/null
+++ b/CarpetDev/CarpetIOF5/par/iof5-refined-cell.par
@@ -0,0 +1,58 @@
+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"
+
+CarpetIOF5::out_every = 256
+CarpetIOF5::out_vars  = "grid::coordinates"
diff --git a/CarpetDev/CarpetIOF5/par/iof5-refined.par b/CarpetDev/CarpetIOF5/par/iof5-refined.par
index aa27521b4..3fdc219f5 100644
--- a/CarpetDev/CarpetIOF5/par/iof5-refined.par
+++ b/CarpetDev/CarpetIOF5/par/iof5-refined.par
@@ -22,7 +22,7 @@ Cactus::cctk_run_title             = "IOF5"
 Cactus::cctk_full_warnings         = yes
 Cactus::highlight_warning_messages = no
 
-Cactus::cctk_itlast = 0
+Cactus::cctk_itlast = 1024
 
 IO::out_dir = $parfile
 
@@ -46,3 +46,6 @@ CarpetRegrid2::radius_1[1]    =  0.5
 CarpetRegrid2::radius_1[2]    =  0.2
 
 InitBase::initial_data_setup_method = "init_all_levels"
+
+CarpetIOF5::out_every = 256
+CarpetIOF5::out_vars  = "grid::coordinates"
diff --git a/CarpetDev/CarpetIOF5/par/iof5-uniform.par b/CarpetDev/CarpetIOF5/par/iof5-uniform.par
index 87e758370..102096f6b 100644
--- a/CarpetDev/CarpetIOF5/par/iof5-uniform.par
+++ b/CarpetDev/CarpetIOF5/par/iof5-uniform.par
@@ -37,3 +37,6 @@ CoordBase::zmax               = +1.0
 CoordBase::dx                 =  0.25
 CoordBase::dy                 =  0.25
 CoordBase::dz                 =  0.25
+
+CarpetIOF5::out_every = 2
+CarpetIOF5::out_vars  = "grid::coordinates"
diff --git a/CarpetDev/CarpetIOF5/param.ccl b/CarpetDev/CarpetIOF5/param.ccl
index f9da760a4..6f1b1f9f5 100644
--- a/CarpetDev/CarpetIOF5/param.ccl
+++ b/CarpetDev/CarpetIOF5/param.ccl
@@ -2,9 +2,14 @@
 
 SHARES: IO
 
-USES STRING out_dir AS IO_out_dir
+USES STRING  out_dir   AS IO_out_dir
+USES INT     out_every AS IO_out_every
 
-USES INT    out_timesteps_per_file
+USES INT     out_timesteps_per_file
+
+USES BOOLEAN checkpoint_ID
+USES INT     checkpoint_every
+USES BOOLEAN checkpoint_on_terminate
 
 
 
@@ -16,6 +21,20 @@ STRING out_dir "Output directory (overrides IO::out_dir)" STEERABLE=always
   ".+" :: "Not empty: directory name"
 } ""
 
+STRING out_vars "Variables to output in F5 format" STEERABLE=always
+{
+  "" :: "List of group or variable names"
+} ""
+
+INT out_every "How often to do CarpetIOF5 output, overrides out_every" STEERABLE=always
+{
+  1:*  :: "Output every so many time steps"
+  -1:0 :: "No output"
+  -2   :: "Use 'IO::out_every'"
+} -2
+
+
+
 KEYWORD file_content "Create one file for every x" STEERABLE=always
 {
   "group"      :: ""
@@ -30,8 +49,9 @@ INT iteration_digits "Minimum number of digits for iteration number" STEERABLE=a
 
 STRING out_filename "File name (without extension) for file_content='everything'" STEERABLE=always
 {
-  "" :: ""
-} "output"
+  ""   :: "use the parameter file name"
+  ".+" :: "use this file name"
+} ""
 
 INT processor_digits "Minimum number of digits for processor number" STEERABLE=always
 {
@@ -53,6 +73,21 @@ BOOLEAN one_dir_per_file "Create one subdirectory per output file to reduce lock
 
 
 
+BOOLEAN separate_single_component_tensors "Use a separated representation even for single-component tensors" STEERABLE=always
+{
+} "no"
+
+
+
+BOOLEAN output_all_timelevels "Output all timelevels (instead of only the current)" STEERABLE=always
+{
+} "no"
+
+INT minimum_component_overlap "Minimum overlap between components, even when ghosts are not output" STEERABLE=always
+{
+  0:* :: ""
+} 1
+
 BOOLEAN output_symmetry_points "Output symmetry and inter-patch boundary points" STEERABLE=always
 {
 } "no"
@@ -64,3 +99,25 @@ BOOLEAN output_ghost_points "Output ghost points" STEERABLE=always
 BOOLEAN output_boundary_points "Output outer boundary points" STEERABLE=always
 {
 } "yes"
+
+
+
+BOOLEAN checkpoint "Checkpoint" STEERABLE=always
+{
+} "no"
+
+BOOLEAN checkpoint_next "Checkpoint at next iteration" STEERABLE=always
+{
+} "no"
+
+
+
+INT compression_level "Compression level" STEERABLE=always
+{
+  -1  :: "no compression"
+  0:9 :: "higher numbers compress better"
+} 1
+
+BOOLEAN use_checksums "Include a checksum for the data" STEERABLE=always
+{
+} "yes"
diff --git a/CarpetDev/CarpetIOF5/schedule.ccl b/CarpetDev/CarpetIOF5/schedule.ccl
index 494ce02fe..59b86f1cf 100644
--- a/CarpetDev/CarpetIOF5/schedule.ccl
+++ b/CarpetDev/CarpetIOF5/schedule.ccl
@@ -1,27 +1,68 @@
 # Schedule definitions for thorn CarpetIOF5
 
-SCHEDULE F5_Output AT initial
+STORAGE: this_iteration next_output_iteration
+
+
+
+# Initialisation
+
+SCHEDULE CarpetIOF5_Startup AT startup
 {
   LANG: C
-  OPTIONS: global-late
-} "Create an output file"
+} "Register I/O method"
 
+SCHEDULE CarpetIOF5_Init AT basegrid
+{
+  LANG: C
+  OPTIONS: global
+} "Initialise I/O data structures"
 
 
-SCHEDULE F5_Poison AT initial AFTER F5_Output BEFORE F5_Input
+
+# Checkpointing
+
+schedule CarpetIOF5_InitialDataCheckpoint at CPINITIAL
 {
   LANG: C
-  OPTIONS: global-late
-} "Poison all variables"
+  OPTIONS: global
+} "Initial data checkpoint routine"
 
-SCHEDULE F5_Input AT initial AFTER F5_Output
+schedule CarpetIOF5_EvolutionCheckpoint at CHECKPOINT
 {
   LANG: C
-  OPTIONS: global-late
-} "Read from file"
+  OPTIONS: global
+} "Evolution checkpoint routine"
 
-SCHEDULE F5_Check AT initial AFTER F5_Input
+schedule CarpetIOF5_TerminationCheckpoint at TERMINATE
 {
   LANG: C
-  OPTIONS: global-late
-} "Check all variables for poison"
+  OPTIONS: global
+} "Termination checkpoint routine"
+
+
+
+#SCHEDULE F5_Output AT initial
+#{
+#  LANG: C
+#  OPTIONS: global-late
+#} "Create an output file"
+#
+#
+#
+#SCHEDULE F5_Poison AT initial AFTER F5_Output BEFORE F5_Input
+#{
+#  LANG: C
+#  OPTIONS: global-late
+#} "Poison all variables"
+#
+#SCHEDULE F5_Input AT initial AFTER F5_Output
+#{
+#  LANG: C
+#  OPTIONS: global-late
+#} "Read from file"
+#
+#SCHEDULE F5_Check AT initial AFTER F5_Input
+#{
+#  LANG: C
+#  OPTIONS: global-late
+#} "Check all variables for poison"
diff --git a/CarpetDev/CarpetIOF5/src/input.cc b/CarpetDev/CarpetIOF5/src/input.cc
index ecf7dda5c..fcb3f7d12 100644
--- a/CarpetDev/CarpetIOF5/src/input.cc
+++ b/CarpetDev/CarpetIOF5/src/input.cc
@@ -70,8 +70,8 @@ namespace CarpetIOF5 {
       
 #if 0
       // TODO: Loop over all charts that exist for the grid, or at
-      // least remember how many maps there are.  (This is also
-      // written as part of the grid structure.)
+      // least remember how many maps there are. (This is also written
+      // as part of the grid structure.)
       // F5iterate_grid_atlas?
       for (int m=0; m<maps; ++m) {
         indent_t indent;
@@ -196,9 +196,9 @@ namespace CarpetIOF5 {
       }
       assert (rl<reflevels);
       
-      // Determine map from fragment name
-      // TODO: store map number as attribute?
-      int const m = interpret_fragmentname (cctkGH, fragmentname);
+      // Determine map and component from fragment name
+      int m, c;
+      interpret_fragmentname (cctkGH, fragmentname, m, c);
       
       {
         indent_t indent;
@@ -275,7 +275,7 @@ namespace CarpetIOF5 {
       assert (iret);
       assert (index_depth_ == index_depth);
       
-      // Read the fragment offset.  This is stored with the dataset
+      // 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,
@@ -288,7 +288,7 @@ namespace CarpetIOF5 {
       assert (all(foff>=0));
       
 #if 0
-      // Read the fragment size.  This is stored with the field -- why
+      // Read the fragment size. This is stored with the field -- why
       // is this different from the offset?
       hsize_t hlen[FIBER_MAX_RANK];
       iret =
@@ -337,7 +337,7 @@ namespace CarpetIOF5 {
       assert (groupdata.dim == dim);
       
       // TODO: This is too expensive; only traverse Carpet's data
-      // structures instead.  (This should be implemented via a
+      // 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) {
diff --git a/CarpetDev/CarpetIOF5/src/iof5.cc b/CarpetDev/CarpetIOF5/src/iof5.cc
new file mode 100644
index 000000000..8b93054d6
--- /dev/null
+++ b/CarpetDev/CarpetIOF5/src/iof5.cc
@@ -0,0 +1,315 @@
+#include <cctk.h>
+#include <cctk_Arguments.h>
+#include <cctk_Parameters.h>
+#include <util_Table.h>
+
+#include <cassert>
+#include <cstdlib>
+#include <cstring>
+#include <iomanip>
+#include <iostream>
+#include <limits>
+#include <sstream>
+#include <string>
+
+#include <hdf5.h>
+#include <iof5.hh>
+
+#include "CactusBase/IOUtil/src/ioutil_CheckpointRecovery.h"
+
+
+
+namespace CarpetIOF5 {
+  
+  using namespace std;
+  using namespace Carpet;
+  
+  
+  
+  // Checkpointing state
+  static int last_checkpoint_iteration = -1;
+  
+  
+  
+  // Scheduled startup routine
+  int CarpetIOF5_Startup ()
+  {
+    CCTK_RegisterBanner ("AMR F5 I/O provided by CarpetIOF5");
+    
+    int const GHExtension = CCTK_RegisterGHExtension (CCTK_THORNSTRING);
+    CCTK_RegisterGHExtensionSetupGH (GHExtension, SetupGH);
+    
+    return 0;
+  }
+  
+  // Registered GH extension setup routine
+  void* SetupGH (tFleshConfig* const fleshconfig,
+                 int const convLevel, cGH* const cctkGH)
+  {
+    DECLARE_CCTK_PARAMETERS;
+    
+    // register I/O method
+    int const IOMethod = CCTK_RegisterIOMethod ("IOF5");
+    CCTK_RegisterIOMethodOutputGH      (IOMethod, OutputGH     );
+    CCTK_RegisterIOMethodTimeToOutput  (IOMethod, TimeToOutput );
+    CCTK_RegisterIOMethodTriggerOutput (IOMethod, TriggerOutput);
+    CCTK_RegisterIOMethodOutputVarAs   (IOMethod, OutputVarAs  );
+    
+    // there no actual extension data structure
+    return NULL;
+  }
+  
+  // Scheduled initialisation routine
+  void CarpetIOF5_Init (CCTK_ARGUMENTS)
+  {
+    DECLARE_CCTK_ARGUMENTS;
+    
+    *this_iteration        = -1;
+    *next_output_iteration =  0;
+    
+    last_checkpoint_iteration = cctk_iteration;
+  }
+  
+  
+  
+  // Callbacks for CarpetIOHDF5's I/O method
+
+  struct callback_arg_t {
+    cGH const* cctkGH;
+  };
+  void do_output (int const vindex,
+                  char const* const optstring,
+                  void* const callback_arg);
+  
+  int OutputGH (cGH const* const cctkGH)
+  {
+    DECLARE_CCTK_ARGUMENTS;
+    DECLARE_CCTK_PARAMETERS;
+    
+    static Carpet::Timer timer ("F5::OutputGH");
+    timer.start();
+    
+    callback_arg_t callback_arg;
+    callback_arg.cctkGH = cctkGH;
+    CCTK_TraverseString (out_vars, do_output, &callback_arg, CCTK_GROUP_OR_VAR);
+    
+    timer.stop(0);
+    
+    return 0;
+  }
+  
+  void do_output (int const vindex,
+                  char const* const optstring,
+                  void* const callback_arg_)
+  {
+    callback_arg_t& callback_arg =
+      * static_cast<callback_arg_t*>(callback_arg_);
+    cGH const* const cctkGH = callback_arg.cctkGH;
+    if (TimeToOutput (cctkGH, vindex)) {
+      TriggerOutput (cctkGH, vindex);
+    }
+  }
+  
+  int TimeToOutput (cGH const* const cctkGH, int const vindex)
+  {
+    DECLARE_CCTK_ARGUMENTS;
+    DECLARE_CCTK_PARAMETERS;
+    
+    int const numvars = CCTK_NumVars();
+    assert (vindex>=0 and vindex<numvars);
+    
+    // Output only in global mode
+    if (not do_global_mode) return 0;
+    
+    // whether to output at this iteration
+    bool output_this_iteration = false;
+    
+    int const my_out_every = out_every == -2 ? IO_out_every : out_every;
+    if (my_out_every > 0) {
+      if (*this_iteration == cctk_iteration) {
+        // we already decided to output this iteration
+        output_this_iteration = true;
+      } else if (cctk_iteration >= *next_output_iteration) {
+        // it is time for the next output
+        output_this_iteration = true;
+        *this_iteration = cctk_iteration;
+        *next_output_iteration = cctk_iteration + my_out_every;
+      }
+    }
+    
+    return output_this_iteration ? 1 : 0;
+  }
+  
+  int TriggerOutput (cGH const* const cctkGH, int const vindex)
+  {
+    DECLARE_CCTK_PARAMETERS;
+    
+    char* const fullname = CCTK_FullName(vindex);
+    int const gindex = CCTK_GroupIndexFromVarI(vindex);
+    char* const groupname = CCTK_GroupName(gindex);
+    for (char* p=groupname; *p; ++p) *p=tolower(*p);
+    int const retval = OutputVarAs (cctkGH, fullname, groupname);
+    free (groupname);
+    free (fullname);
+    
+    return retval;
+  }
+  
+  int OutputVarAs (cGH const* const cctkGH,
+                   const char* const varname, const char* const alias)
+  {
+    DECLARE_CCTK_PARAMETERS;
+    
+    assert (is_level_mode());
+    BEGIN_GLOBAL_MODE(cctkGH) {
+      DECLARE_CCTK_ARGUMENTS;
+      
+      CCTK_VInfo (CCTK_THORNSTRING,
+                  "F5::OutputVarAs: iteration=%d", cctk_iteration);
+      
+      
+      
+      // We don't know how to open multiple files yet
+      assert (CCTK_EQUALS (file_content, "everything"));
+      
+      // Open file
+      static bool first_time = true;
+      
+      // The file name doesn't matter since we currently write
+      // everything into a single file
+      int const vindex = CCTK_VarIndex (varname);
+      assert (vindex >= 0);
+      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);
+      
+      indent_t indent;
+      cout << indent << "process=" << proc << "\n";
+      
+      bool const truncate_file = first_time and IO_TruncateOutputFiles(cctkGH);
+      hid_t const file =
+        truncate_file ?
+        H5Fcreate (name.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT) :
+        H5Fopen   (name.c_str(), H5F_ACC_RDWR , H5P_DEFAULT);
+      assert (file >= 0);
+      first_time = false;
+      
+      vector<bool> output_var(CCTK_NumVars());
+      output_var.at(vindex) = true;
+      output (cctkGH, file, output_var, false);
+      
+      // Close file
+      herr_t const herr = H5Fclose (file);
+      assert (not herr);
+      
+    } END_GLOBAL_MODE;
+    
+    return 0;                   // no error
+  }
+  
+  
+  
+  // Checkpointing
+  
+  void Checkpoint (cGH const* const cctkGH, int const called_from)
+  {
+    assert (is_global_mode());
+    
+    // generate filenames for both the temporary and real checkpoint
+    // files
+    int const ioproc = CCTK_MyProc(cctkGH);
+    int const parallel_io = 1;
+    char* const filename =
+      IOUtil_AssembleFilename (cctkGH, NULL, "", ".f5",
+                               called_from, ioproc, not parallel_io);
+    char* const tempname =
+      IOUtil_AssembleFilename (cctkGH, NULL, ".tmp", ".f5",
+                               called_from, ioproc, not parallel_io);
+    
+    hid_t const file =
+      H5Fcreate (tempname, H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
+    assert (file >= 0);
+    
+    vector<bool> output_var(CCTK_NumVars());
+    for (int gindex=0; gindex<CCTK_NumGroups(); ++gindex) {
+      // only checkpoint groups with storage
+      if (CCTK_QueryGroupStorageI(cctkGH, gindex) <= 0) continue;
+      if (CCTK_NumVarsInGroupI(gindex) == 0) continue;
+      
+      // do not checkpoint groups with a "checkpoint=no" tag
+      cGroup gdata;
+      CCTK_GroupData (gindex, &gdata);
+      int const len =
+        Util_TableGetString (gdata.tagstable, 0, NULL, "checkpoint");
+      if (len > 0) {
+        char buf[1000];
+        Util_TableGetString (gdata.tagstable, sizeof buf, buf, "checkpoint");
+        if (CCTK_EQUALS(buf, "no")) continue;
+        assert (CCTK_EQUALS(buf, "yes"));
+      }
+      
+      int const first_vindex = CCTK_FirstVarIndexI (gindex);
+      int const num_vars = CCTK_NumVarsInGroupI (gindex);
+      for (int vindex=first_vindex; vindex<first_vindex+num_vars; ++vindex) {
+        output_var.at(vindex) = true;
+      }
+    }
+    
+    output (cctkGH, file, output_var, true);
+    
+    // Close file
+    herr_t const herr = H5Fclose (file);
+    assert (not herr);
+    
+    // Wait until all files have been written, then rename the
+    // checkpoint files
+    // TODO: ensure there were no errors
+    CCTK_Barrier (cctkGH);
+    int const ierr = rename (tempname, filename);
+    assert (not ierr);
+  }
+  
+  void CarpetIOF5_InitialDataCheckpoint (CCTK_ARGUMENTS)
+  {
+    DECLARE_CCTK_ARGUMENTS;
+    DECLARE_CCTK_PARAMETERS;
+    
+    if (checkpoint and checkpoint_ID) {
+      Checkpoint (cctkGH, CP_INITIAL_DATA);
+    }
+  }
+  
+  void CarpetIOF5_EvolutionCheckpoint (CCTK_ARGUMENTS)
+  {
+    DECLARE_CCTK_ARGUMENTS;
+    DECLARE_CCTK_PARAMETERS;
+    
+    bool const checkpoint_by_iteration =
+      checkpoint_every > 0 and
+      cctk_iteration >= last_checkpoint_iteration + checkpoint_every;
+    bool const do_checkpoint =
+      checkpoint and (checkpoint_by_iteration or checkpoint_next);
+    
+    if (do_checkpoint) {
+      Checkpoint (cctkGH, CP_EVOLUTION_DATA);
+    }
+  }
+  
+  void CarpetIOF5_TerminationCheckpoint (CCTK_ARGUMENTS)
+  {
+    DECLARE_CCTK_ARGUMENTS;
+    DECLARE_CCTK_PARAMETERS;
+    
+    if (checkpoint and checkpoint_on_terminate) {
+      
+      bool const did_checkpoint = cctk_iteration == last_checkpoint_iteration;
+      bool const do_checkpoint = not did_checkpoint;
+      
+      if (do_checkpoint) {
+        Checkpoint (cctkGH, CP_EVOLUTION_DATA);
+      }
+    }
+  }
+  
+} // end namespace CarpetIOF5
diff --git a/CarpetDev/CarpetIOF5/src/iof5.cc.old b/CarpetDev/CarpetIOF5/src/iof5.cc.old
deleted file mode 100644
index fe780030b..000000000
--- a/CarpetDev/CarpetIOF5/src/iof5.cc.old
+++ /dev/null
@@ -1,643 +0,0 @@
-#include <cctk.h>
-#include <cctk_Arguments.h>
-#include <cctk_Parameters.h>
-
-#include <cassert>
-#include <cstdlib>
-#include <cstring>
-#include <iomanip>
-#include <iostream>
-#include <limits>
-#include <sstream>
-#include <string>
-
-#include <hdf5.h>
-#include <F5/F5F.h>
-#include <F5/F5R.h>
-#include <F5/F5iterate.h>
-#include <F5/F5uniform.h>
-
-#include <bbox.hh>
-#include <vect.hh>
-
-#include <carpet.hh>
-
-
-
-namespace CarpetIOF5 {
-  
-  using namespace std;
-  using namespace Carpet;
-  
-  
-  
-  // File mode for creating directories
-  int const mode = 0755;
-  
-  // A nan
-  CCTK_REAL const nan = numeric_limits<CCTK_REAL>::quiet_NaN();
-  
-  // Indentation
-  inline string indent (int const n)
-  {
-    int const width = 3;
-    return string(width*n, ' ');
-  }
-  
-  
-  
-  typedef vect<hsize_t,dim> hvect;
-  typedef vect<float,dim> fvect;
-  
-  
-  
-  static inline
-  hvect v2h (ivect const& a)
-  {
-    return hvect(a);
-  }
-  
-  static inline
-  F5_vec3_point_t v2p (rvect const& a)
-  {
-    F5_vec3_point_t res;
-    res.x = a[0];
-    res.y = a[1];
-    res.z = a[2];
-    return res;
-  }
-  
-  static inline
-  F5_vec3_float_t v2f (rvect const& a)
-  {
-    F5_vec3_float_t res;
-    res.x = a[0];
-    res.y = a[1];
-    res.z = a[2];
-    return res;
-  }
-  
-  static inline
-  F5_vec3_double_t v2d (rvect const& a)
-  {
-    F5_vec3_double_t res;
-    res.x = a[0];
-    res.y = a[1];
-    res.z = a[2];
-    return res;
-  }
-  
-  
-  
-  // Generate a good grid name (simulation name)
-  string
-  generate_gridname (cGH const *const cctkGH)
-  {
-#if 0
-    char const *const gridname = (char const*) UniqueSimulationID(cctkGH);
-    assert (gridname);
-    return gridname;
-#endif
-    // Use the parameter file name, since the simulation ID is too
-    // long and doesn't look nice
-    char parfilename[10000];
-    CCTK_ParameterFilename (sizeof parfilename, parfilename);
-    char *const p = strstr(parfilename, ".par");
-    if (p) *p = '\0';
-    char *const s = strrchr(parfilename, '/');
-    char *const gridname = s ? s+1 : parfilename;
-    return gridname;
-  }
-  
-  
-  
-  // Generate a good file name ("alias") for a variable
-  string
-  generate_basename (cGH const *const cctkGH,
-                     int const variable)
-  {
-    DECLARE_CCTK_PARAMETERS;
-    
-    assert (variable >= 0);
-    
-    ostringstream filename_buf;
-    
-    if (CCTK_EQUALS (file_content, "variable")) {
-      char *const varname = CCTK_FullName(variable);
-      assert (varname);
-      for (char *p = varname; *p; ++p) {
-        *p = tolower(*p);
-      }
-      filename_buf << varname;
-      free (varname);
-    }
-    else if (CCTK_EQUALS (file_content, "group")) {
-      char *const groupname = CCTK_GroupNameFromVarI(variable);
-      assert (groupname);
-      for (char *p = groupname; *p; ++p) {
-        *p = tolower(*p);
-      }
-      filename_buf << groupname;
-      free (groupname);
-    }
-    else if (CCTK_EQUALS (file_content, "thorn")) {
-      char const *const impname = CCTK_ImpFromVarI(variable);
-      char *const thornname = strdup(impname);
-      assert (thornname);
-      char *const colon = strchr(thornname, ':');
-      assert (colon);
-      *colon = '\0';
-      for (char *p = thornname; *p; ++p) {
-        *p = tolower(*p);
-      }
-      filename_buf << thornname;
-      free (thornname);
-    }
-    else if (CCTK_EQUALS (file_content, "everything")) {
-      filename_buf << out_filename;
-    }
-    else {
-      assert (0);
-    }
-    
-    if (out_timesteps_per_file > 0) {
-      int const iteration = (cctkGH->cctk_iteration
-                             / out_timesteps_per_file * out_timesteps_per_file);
-      filename_buf << ".it"
-                   << setw (iteration_digits) << setfill ('0') << iteration;
-    }
-    
-    return filename_buf.str();
-  }
-  
-  
-  
-  // Create the final file name on a particular processor
-  string
-  create_filename (cGH const *const cctkGH,
-                   string const basename,
-                   bool const create_directories)
-  {
-    DECLARE_CCTK_PARAMETERS;
-    
-    int const proc   = CCTK_MyProc(cctkGH);
-    int const nprocs = CCTK_nProcs(cctkGH);
-    
-    bool const use_IO_out_dir = strcmp(out_dir, "") == 0;
-    string path = use_IO_out_dir ? IO_out_dir : out_dir;
-    
-    if (create_subdirs) {
-      if (nprocs >= 10000) {
-        ostringstream buf;
-        buf << path + "/"
-            << basename
-            << ".p" << setw (max (0, processor_digits - 4)) << setfill ('0')
-            << proc / 10000
-            << "nnnn/";
-        path = buf.str();
-        if (create_directories) {
-          check (CCTK_CreateDirectory (mode, path.c_str()) >= 0);
-        }
-      }
-      if (nprocs >= 100) {
-        ostringstream buf;
-        buf << path + "/"
-            << basename
-            << ".p" << setw (max (0, processor_digits - 2)) << setfill ('0')
-            << proc / 100
-            << "nn/";
-        path = buf.str();
-        if (create_directories) {
-          check (CCTK_CreateDirectory (mode, path.c_str()) >= 0);
-        }
-      }
-    }
-    if (one_dir_per_file and nprocs >= 1) {
-      ostringstream buf;
-      buf << path
-          << basename
-          << ".p" << setw (processor_digits) << setfill ('0')
-          << proc
-          << "/";
-      path = buf.str();
-      if (create_directories) {
-        check (CCTK_CreateDirectory (mode, path.c_str()) >= 0);
-      }
-    }
-    ostringstream buf;
-    buf << path + "/"
-        << basename
-        << ".p" << setw (processor_digits) << setfill ('0') << proc
-        << out_extension;
-    return buf.str();
-  }
-  
-  
-  
-  extern "C"
-  void F5_Output (CCTK_ARGUMENTS)
-  {
-    DECLARE_CCTK_ARGUMENTS;
-    DECLARE_CCTK_PARAMETERS;
-    
-    herr_t herr;
-    
-    
-    
-    assert (is_global_mode());
-    CCTK_VInfo (CCTK_THORNSTRING, "F5_Output: iteration=%d", cctk_iteration);
-    
-    string const gridname = generate_gridname(cctkGH);
-    
-    
-    
-    // Open file
-    static bool first_time = true;
-    
-    string const basename =
-      generate_basename (cctkGH, CCTK_VarIndex("grid::r"));
-    string const name =
-      create_filename (cctkGH, basename, first_time);
-    
-    bool const truncate_file = first_time and IO_TruncateOutputFiles(cctkGH);
-    hid_t const file =
-      truncate_file ?
-      H5Fcreate (name.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT) :
-      H5Fopen (name.c_str(), H5F_ACC_RDWR, H5P_DEFAULT);
-    assert (file >= 0);
-    first_time = false;
-    
-    
-    
-    BEGIN_REFLEVEL_LOOP (cctkGH) {
-      DECLARE_CCTK_ARGUMENTS;
-      
-#if 0
-      // Write data
-#warning "TODO: Save the coordinates in double precision"
-      F5Path *const path =
-        F5Fwrite_uniform_cartesian3D
-        (file, cctk_time, gridname, &v2p(lower), &v2f(delta), &v2h(lsh)[0],
-         "radius" /* group name */, H5T_NATIVE_DOUBLE, r, NULL, F5P_DEFAULT);
-#endif
-      
-      // Define grid hierarchy
-      ivect const reffact = spacereffacts.AT(reflevel);
-      F5Path *const path = F5Rcreate_vertex_refinement3D
-        (file, cctk_time, gridname.c_str(), &v2h(reffact)[0], NULL);
-      
-      // Define iteration
-      F5Rset_timestep (path, cctk_iteration);
-      
-#if 0
-      F5Path *refpath = NULL;
-      static CCTK_REAL last_root_time = -1;
-      if (reflevel == 0) {
-        last_root_time = cctk_time;
-      } else {
-        refpath = F5Rcreate_relative_vertex_Qrefinement3D
-          (file, cctk_time, gridname, &v2h(reffact)[0],
-           last_root_time, &v2h(ivect(1))[0]);
-      }
-#endif
-      
-      
-      
-      BEGIN_LOCAL_MAP_LOOP (cctkGH, CCTK_GF) {
-        DECLARE_CCTK_ARGUMENTS;
-        
-        // Determine level coordinates
-        ivect gsh;
-        rvect origin, delta;
-        rvect lower, upper;
-        for (int d=0; d<dim; ++d) {
-          gsh[d]    = cctk_gsh[d];
-          origin[d] = CCTK_ORIGIN_SPACE(d);
-          delta[d]  = CCTK_DELTA_SPACE(d);
-          lower[d]  = origin[d];
-          upper[d]  = lower[d] + (gsh[d]-1) * delta[d];
-        }
-        
-        // Define level topology
-        if (reflevel == 0) {
-          // Choose (arbitrarily) the root level as default topology,
-          // for readers which don't understand AMR
-          F5Rlink_default_vertex_topology (path, &v2h(reffact)[0]);
-        }
-        
-        // Define level geometry
-        F5Fwrite_linear (path, FIBER_HDF5_POSITIONS_STRING,
-                         dim, &v2h(gsh)[0],
-                         F5T_COORD3_DOUBLE, &v2d(lower), &v2d(delta));
-        F5Fset_range (path, &v2d(lower), &v2d(upper));
-        
-        
-        
-        BEGIN_LOCAL_COMPONENT_LOOP (cctkGH, CCTK_GF) {
-          DECLARE_CCTK_ARGUMENTS;
-          
-          // Determine component coordinates
-          ivect lbnd, lsh, lghosts, ughosts;
-          rvect clower, cupper;
-#if 0
-          F5_refinement3D_point_t iorigin, idelta;
-#endif
-          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];
-            clower[d] = origin[d] + cctk_lbnd[d] * delta[d];
-            cupper[d] = clower[d] + (lsh[d]-1) * delta[d];
-#if 0
-            iorigin.d[d].num   =
-              cctk_levoff[d] + cctk_lbnd[d] * cctk_levoffdenom[d]; 
-            iorigin.d[d].denom = cctk_levoffdenom[d] * reffact[d];
-            idelta .d[d].num   = 1; 
-            idelta .d[d].denom = reffact[d];
-#endif
-          }
-          ostringstream buf2;
-          buf2 << "component=" << component;
-          string const cname = buf2.str();
-        
-#if 0
-          if (reflevel > 0) {
-            F5Fwrite_linear_fraction (refpath, FIBER_HDF5_POSITIONS_STRING,
-                                      cctk_dim, &v2h(gsh)[0], &v2h(lsh)[0],
-                                      F5T_COORD3_INT_FRACTION32,
-                                      &iorigin, &idelta, &v2h(lbnd)[0],
-                                      cname.c_str());
-          }
-#endif
-          
-          // Define coordinates
-          // (This is redundant, since the level's overall bounding
-          // box was already defined above, but it provides the
-          // individual components' bounding boxes.)
-          F5Fwrite_linear_fraction (path, FIBER_HDF5_POSITIONS_STRING,
-                                    cctk_dim, &v2h(gsh)[0], &v2h(lsh)[0],
-                                    F5T_COORD3_DOUBLE,
-                                    &clower, &delta, &v2h(lbnd)[0],
-                                    cname.c_str());
-          
-          // Write data
-          // scalar field
-          F5Fwrite_fraction (path, "radius" /* group name */,
-                             cctk_dim, &v2h(gsh)[0], &v2h(lsh)[0],
-                             H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE,
-                             r,
-                             &v2h(lbnd)[0], &v2h(lghosts)[0], &v2h(ughosts)[0],
-                             cname.c_str(), H5P_DEFAULT);
-          
-          // vector field
-          void const* data[dim];
-          data[0] = x;
-          data[1] = y;
-          data[2] = z;
-          F5FSwrite_fraction (path, "coordinates" /* group name */,
-                              cctk_dim, &v2h(gsh)[0], &v2h(lsh)[0],
-                              F5T_VEC3_DOUBLE, F5T_VEC3_DOUBLE,
-                              data,
-                              &v2h(lbnd)[0], &v2h(lghosts)[0], &v2h(ughosts)[0],
-                              cname.c_str(), H5P_DEFAULT, 0);
-          
-        } END_LOCAL_COMPONENT_LOOP;
-        
-      } END_LOCAL_MAP_LOOP;
-      
-      // Close topology
-      F5close (path);
-      
-    } END_REFLEVEL_LOOP;
-    
-    // Close file
-    herr = H5Fclose (file);
-    assert (not herr);
-  }
-  
-  
-  
-  // Use a class for reading data, so that we have an easy way to pass
-  // variables between the various iterators
-  class iterator_t {
-    cGH *cctkGH;
-    double time;
-    char const *gridname;
-    char const *topologyname;
-    int index_depth;            // 0=vertex, 1=cell
-    int topological_dimension;
-    char const *fieldname;
-    char const *fragmentname;
-    
-  public:
-    iterator_t (cGH *const cctkGH_)
-      : cctkGH(cctkGH_),
-        time(nan),
-        gridname(NULL),
-        topologyname(NULL), index_depth(-1), topological_dimension(-1),
-        fieldname(NULL),
-        fragmentname(NULL)
-    {
-    }
-    
-    
-    
-  private:
-    
-    void read_timeslice (F5Path *const path)
-    {
-      cout << indent(1) << "time=" << time << "\n";
-      F5iterate_grids (path, NULL, grid_iterator, this, NULL, NULL);
-    }
-  
-    void read_grid (F5Path *const path)
-    {
-      cout << indent(2) << "grid=\"" << gridname << "\"\n";
-      // F5iterate_vertex_fields (path, NULL, field_iterator, this, NULL, NULL);
-      F5iterate_topologies (path, NULL, topology_iterator, this);
-    }
-    
-    void read_topology (F5Path *const path)
-    {
-      herr_t herr;
-      
-      cout << indent(3) << "topology=\"" << topologyname << "\""
-           << " (" << (index_depth==0 ? "vertex" : "cell") << ")\n";
-      
-      // Ignore topologies that are only an alias for another topology
-      H5G_stat_t stat;
-      herr = H5Gget_objinfo (path->Grid_hid, topologyname, false, &stat);
-      assert (not herr);
-      if (stat.type == H5G_LINK) {
-        char linkval[100000];
-        herr = H5Gget_linkval
-          (path->Grid_hid, topologyname, sizeof linkval, linkval);
-        assert (not herr);
-        cout << indent(4) << "alias for topology \"" << linkval << "\"\n"
-             << indent(4) << "ignoring this topology\n";
-        return;
-      }
-      
-      F5iterate_topology_fields (path, NULL, field_iterator, this, NULL, NULL); 
-    }
-    
-    void read_field (F5Path *const path)
-    {
-      cout << indent(4) << "field=\"" << fieldname << "\"\n";
-      F5iterate_field_fragments (path, NULL, fragment_iterator, this);
-    }
-    
-    void read_fragment (F5Path *const path)
-    {
-      herr_t herr;
-      
-      cout << indent(5) << "fragment=\"" << fragmentname << "\"\n";
-      
-      if (strcmp(fieldname, "radius") == 0) {
-        
-        hid_t const group = path->Field_hid;
-        read_variable (group, fragmentname, CCTK_VarIndex("grid::r"));
-        
-      } else if (strcmp(fieldname, "coordinates") == 0) {
-        
-        hid_t const group =
-          H5Gopen (path->Field_hid, fragmentname, H5P_DEFAULT);
-        assert (group);
-        read_variable (group, "x", CCTK_VarIndex("grid::x"));
-        read_variable (group, "y", CCTK_VarIndex("grid::y"));
-        read_variable (group, "z", CCTK_VarIndex("grid::z"));
-        herr = H5Gclose (group);
-        assert (not herr);
-        
-      } else {
-        // do nothing
-      }
-    }
-    
-    void read_variable (hid_t const group, char const *const name,
-                        int const var)
-    {
-      assert (group >= 0);
-      assert (name);
-      assert (var >= 0);
-      
-      cout << indent(6) << "dataset=\"" << name << "\"\n";
-      //  hid_t const fragment = H5Dopen(group, name);
-      //  assert (fragment);
-      //  
-      //  hid_t const space = H5Dget_space(fragment);
-      //  assert (space);
-      //  int const rank = H5Sget_simple_extent_dims(space, NULL, NULL);
-      //  assert (rank>=0);
-      //  
-      //  assert (rank == cctk_dim);
-      //  
-      //  H5Sclose(space_id);
-    }
-    
-    
-    
-  public:
-    
-    void iterate (hid_t const object)
-    {
-      F5iterate_timeslices (object, NULL, timeslice_iterator, this);
-    }
-    
-    static
-    herr_t timeslice_iterator (F5Path *const path, double const time,
-                               void *const userdata)
-    {
-      iterator_t* const iterator = (iterator_t*)userdata;
-      iterator->time = time;
-      iterator->read_timeslice (path);
-      return 0;
-    }
-    
-    static
-    herr_t grid_iterator (F5Path *const path, char const *const gridname,
-                          void *const userdata)
-    {
-      iterator_t* const iterator = (iterator_t*)userdata;
-      iterator->gridname = gridname;
-      iterator->read_grid (path);
-      return 0;
-    }
-    
-    static
-    herr_t topology_iterator (F5Path *const path,
-                              char const *const topologyname,
-                              int const index_depth,
-                              int const topological_dimension,
-                              void *const userdata)
-    {
-      iterator_t* const iterator = (iterator_t*)userdata;
-      iterator->topologyname = topologyname;
-      iterator->index_depth = index_depth;
-      iterator->topological_dimension = topological_dimension;
-      iterator->read_topology (path);
-      return 0;
-    }
-    
-    static
-    herr_t field_iterator (F5Path *const path, char const *const fieldname,
-                           void *const userdata)
-    {
-      iterator_t* const iterator = (iterator_t*)userdata;
-      iterator->fieldname = fieldname;
-      iterator->read_field (path);
-      return 0;
-    }
-    
-    static
-    herr_t fragment_iterator (F5Path *const path,
-                              char const *const fragmentname,
-                              void *const userdata)
-    {
-      iterator_t* const iterator = (iterator_t*)userdata;
-      iterator->fragmentname = fragmentname;
-      iterator->read_fragment (path);
-      return 0;
-    }
-    
-  };
-  
-  
-  
-  extern "C"
-  void F5_Input (CCTK_ARGUMENTS)
-  {
-    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"));
-    string const name =
-      create_filename (cctkGH, basename, false);
-    
-    hid_t const file = H5Fopen (name.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
-    assert (file >= 0);
-    
-    // Iterate over all time slices
-    iterator_t iterator(cctkGH);
-    iterator.iterate (file);
-    
-    // Close file
-    herr = H5Fclose (file);
-    assert (not herr);
-  }
-  
-} // end namespace CarpetIOF5
diff --git a/CarpetDev/CarpetIOF5/src/iof5.hh b/CarpetDev/CarpetIOF5/src/iof5.hh
index 71524fa2c..76837a1e0 100644
--- a/CarpetDev/CarpetIOF5/src/iof5.hh
+++ b/CarpetDev/CarpetIOF5/src/iof5.hh
@@ -10,6 +10,7 @@
 #include <limits>
 #include <sstream>
 #include <string>
+#include <vector>
 
 #include <hdf5.h>
 #include <F5/F5F.h>
@@ -33,8 +34,8 @@
 template<typename T>
 static
 T failwarn (bool& error_flag, T const expr,
-            int const line, char const *const file, char const *const thorn,
-            char const *const msg)
+            int const line, char const* const file, char const* const thorn,
+            char const* const msg)
 {
   if (expr < 0) {
     CCTK_VWarn (CCTK_WARN_ALERT, line, file, thorn,
@@ -79,9 +80,16 @@ namespace CarpetIOF5 {
   CCTK_REAL const nan = numeric_limits<CCTK_REAL>::quiet_NaN();
   
   // Special group and attribute names
-  char const *const metadata_group = "Parameters and Global Attributes";
-  char const *const all_parameters = "All Parameters";
-  extern char const *const grid_structure;
+  char const* const metadata_group = "Parameters and Global Attributes";
+  char const* const all_parameters = "All Parameters";
+  extern char const* const grid_structure;
+  
+  
+  
+  // Tensor types
+  enum tensortype_t {
+    tt_scalar, tt_vector, tt_symtensor, tt_tensor, tt_error
+  };
   
   
   
@@ -143,72 +151,108 @@ namespace CarpetIOF5 {
   
   // Handle HDF5 attributes more comfortably
   bool WriteAttribute (hid_t const group,
-                       char const * const name,
+                       char const*  const name,
                        int const ivalue);
   bool WriteAttribute (hid_t const group,
-                       char const * const name,
+                       char const*  const name,
                        double const dvalue);
   bool WriteAttribute (hid_t const group,
-                       char const *const name,
-                       char const *const svalue);
+                       char const* const name,
+                       char const* const svalue);
   bool WriteAttribute (hid_t const group,
-                       char const *const name,
-                       int const *const ivalues,
+                       char const* const name,
+                       int const* const ivalues,
                        int const nvalues);
   bool WriteAttribute (hid_t const group,
-                       char const *const name,
-                       double const *const dvalues,
+                       char const* const name,
+                       double const* const dvalues,
                        int const nvalues);
   bool WriteAttribute (hid_t const group,
-                       char const *const name,
-                       char const *const *const svalues,
+                       char const* const name,
+                       char const* const* const svalues,
                        int const nvalues);
   bool WriteLargeAttribute (hid_t const group,
-                            char const *const name,
-                            char const *const svalue);
+                            char const* const name,
+                            char const* const svalue);
   
   
   
   // Generate a good file name ("alias") for a variable
   string
-  generate_basename (cGH const *const cctkGH,
+  generate_basename (cGH const* const cctkGH,
                      int const variable);
   
   // Create the final file name on a particular processor
   string
-  create_filename (cGH const *const cctkGH,
+  create_filename (cGH const* const cctkGH,
                    string const basename,
                    int const proc,
                    bool const create_directories);
   
   // Generate a good grid name (simulation name)
   string
-  generate_gridname (cGH const *const cctkGH);
+  generate_gridname (cGH const* const cctkGH);
   
-  // Generate a good tensor basis name (coordinate system name)
+  // Generate a good topology name (refinement level name)
   string
-  generate_chartname (cGH const *const cctkGH);
+  generate_topologyname (cGH const* const cctkGH,
+                         int const gi,
+                         ivect const& reffact);
   
-  // Generate a good fragment name (map name)
+  // Generate a good chart name (tensor basis name)
   string
-  generate_fragmentname (cGH const *const cctkGH, int const m);
-  int
-  interpret_fragmentname (cGH const *const cctkGH,
-                          char const *const fragmentname);
+  generate_chartname (cGH const* const cctkGH);
   
-  // Generate a good topology name (map and refinement level name)
+  // Generate a good fragment name (map and component name)
   string
-  generate_topologyname (cGH const *const cctkGH,
-                         int const m, ivect const& reffact);
+  generate_fragmentname (cGH const* const cctkGH, int const m, int const c);
+  void
+  interpret_fragmentname (cGH const* const cctkGH,
+                          char const* const fragmentname,
+                          int& m, int& c);
+  
+  // Generate a good field name (group or variable name)
+  string
+  generate_fieldname (cGH const* const cctkGH,
+                      int const vi, tensortype_t const tt);
   
   
   
   void
-  write_metadata (cGH *const cctkGH, hid_t const file);
+  write_metadata (cGH const* const cctkGH, hid_t const file);
   
   // Handle Carpet's grid structure (this should move to Carpet and/or
   // CarpetLib)
   string
-  serialise_grid_structure (cGH const *const cctkGH);
+  serialise_grid_structure (cGH const* const cctkGH);
+  
+  
+  
+  void output (cGH const* const cctkGH,
+               hid_t const file,
+               vector<bool> const& output_var,
+               bool const output_everything);
+  
+  
+  
+  // Scheduled routines
+  extern "C" {
+    int CarpetIOF5_Startup ();
+    void CarpetIOF5_Init (CCTK_ARGUMENTS);
+    void CarpetIOF5_InitialDataCheckpoint (CCTK_ARGUMENTS);
+    void CarpetIOF5_EvolutionCheckpoint (CCTK_ARGUMENTS);
+    void CarpetIOF5_TerminationCheckpoint (CCTK_ARGUMENTS);
+  }
+  
+  // Registered GH extension setup routine
+  void* SetupGH (tFleshConfig* const fleshconfig,
+                 int const convLevel, cGH* const cctkGH);
+  
+  // Callbacks for CarpetIOHDF5's I/O method
+  int OutputGH (cGH const* const cctkGH);
+  int TimeToOutput (cGH const* const cctkGH, int const vindex);
+  int TriggerOutput (cGH const* const cctkGH, int const vindex);
+  int OutputVarAs (cGH const* const cctkGH,
+                   const char* const varname, const char* const alias);
   
 } // end namespace CarpetIOF5
diff --git a/CarpetDev/CarpetIOF5/src/make.code.defn b/CarpetDev/CarpetIOF5/src/make.code.defn
index da381b90b..4603ab0bc 100644
--- a/CarpetDev/CarpetIOF5/src/make.code.defn
+++ b/CarpetDev/CarpetIOF5/src/make.code.defn
@@ -1,7 +1,7 @@
 # Main make.code.defn file for thorn CarpetIOF5
 
 # Source files in this directory
-SRCS = attributes.cc input.cc output.cc poison.cc util.cc
+SRCS = attributes.cc input.cc iof5.cc output.cc poison.cc util.cc
 
 # Subdirectories containing source files
 SUBDIRS = 
diff --git a/CarpetDev/CarpetIOF5/src/output.cc b/CarpetDev/CarpetIOF5/src/output.cc
index 123c40524..b35d43775 100644
--- a/CarpetDev/CarpetIOF5/src/output.cc
+++ b/CarpetDev/CarpetIOF5/src/output.cc
@@ -39,18 +39,133 @@ namespace CarpetIOF5 {
   
   
   class output_iterator_t {
-    cGH *const cctkGH;
+    // Can't be cGH const, since the mode loops change change its
+    // entries
+    cGH* const cctkGH;
     
+    vector<bool> const output_var;    // whether to output this variable
+    bool const output_everything;
     bool const is_multipatch;
     
+    int group_type;             // CCTK_GF or CCTK_ARRAY
+    int group_index;            // if group_type != CCTK_GF; else -1
+    vector<int> vindices;       // variable indices to output
+    
     string gridname;
     string chartname;
-    string fragmentname;
     string topologyname;
+    string fragmentname;
+    
+    
+    
+    struct map_indices_t {
+      int dim;
+      ivect gsh;
+      rvect origin, delta;
+      rvect lower, upper;
+      
+      map_indices_t (cGH const* const cctkGH, int const gindex)
+      {
+        DECLARE_CCTK_ARGUMENTS;
+        
+        if (gindex == -1) {
+          // grid function
+          dim = ::dim;
+          for (int d=0; d<(::dim); ++d) {
+            gsh[d]    = cctk_gsh[d];
+            origin[d] = CCTK_ORIGIN_SPACE(d);
+            delta[d]  = CCTK_DELTA_SPACE(d);
+          }
+        } else {
+          // grid array
+          cGroupDynamicData dyndata;
+          int const ierr = CCTK_GroupDynamicData (cctkGH, gindex, &dyndata);
+          assert (not ierr);
+          // HDF5 and F5 can't handle dim=0
+          dim = max(dyndata.dim, 1);
+          for (int d=0; d<dyndata.dim; ++d) {
+            gsh[d]    = dyndata.gsh[d];
+            origin[d] = 0.0;
+            delta[d]  = 1.0;
+          }
+          for (int d=dyndata.dim; d<(::dim); ++d) {
+            gsh[d]    = 1;
+            origin[d] = 0.0;
+            delta[d]  = 1.0;
+          }
+        }
+        for (int d=0; d<(::dim); ++d) {
+          lower[d] = origin[d];
+          upper[d] = lower[d] + (gsh[d]-1) * delta[d];
+        }
+      }
+    };
+    
+    struct component_indices_t: map_indices_t {
+      // elements >=dim remain undefined
+      ivect lbnd, lsh, lghosts, ughosts;
+      ivect imin, imax, ioff, ilen;
+      rvect clower, cupper;
+      
+      component_indices_t (cGH const* const cctkGH, int const gindex)
+        : map_indices_t(cctkGH, gindex)
+      {
+        DECLARE_CCTK_ARGUMENTS;
+        DECLARE_CCTK_PARAMETERS;
+        
+        if (gindex == -1) {
+          // grid function
+          for (int d=0; d<(::dim); ++d) {
+            lbnd[d]    = cctk_lbnd[d];
+            lsh[d]     = cctk_lsh[d];
+            lghosts[d] = cctk_bbox[2*d  ] ? 0 : cctk_nghostzones[d];
+            ughosts[d] = cctk_bbox[2*d+1] ? 0 : cctk_nghostzones[d];
+          }
+        } else {
+          // grid array
+          cGroupDynamicData dyndata;
+          int const ierr = CCTK_GroupDynamicData (cctkGH, gindex, &dyndata);
+          assert (not ierr);
+          for (int d=0; d<dyndata.dim; ++d) {
+            lbnd[d]    = dyndata.lbnd[d];
+            lsh[d]     = dyndata.lsh[d];
+            lghosts[d] = dyndata.bbox[2*d  ] ? 0 : dyndata.nghostzones[d];
+            ughosts[d] = dyndata.bbox[2*d+1] ? 0 : dyndata.nghostzones[d];
+          }
+          for (int d=dyndata.dim; d<(::dim); ++d) {
+            lbnd[d]    = 0;
+            lsh[d]     = 1;
+            lghosts[d] = 0;
+            ughosts[d] = 0;
+          }
+        }
+        for (int d=0; d<(::dim); ++d) {
+          imin[d] = 0;
+          imax[d] = lsh[d];
+          if (not output_ghost_points) {
+            int const overlap = min(ughosts[d], minimum_component_overlap);
+            imin[d] += lghosts[d];
+            imax[d] -= ughosts[d] - overlap;
+            lghosts[d] = 0;
+            ughosts[d] = overlap;
+          }
+          ioff[d] = lbnd[d] + imin[d];
+          ilen[d] = imax[d] - imin[d];
+          clower[d] = lower[d] + ioff[d] * delta[d];
+          cupper[d] = clower[d] + (ilen[d]-1) * delta[d];
+        }
+      }
+    };
+    
+    
     
   public:
-    output_iterator_t (cGH *const cctkGH_)
+    output_iterator_t (cGH* const cctkGH_,
+                       vector<bool> const& output_var_,
+                       bool const output_everything_)
       : cctkGH(cctkGH_),
+        output_var(output_var_),
+        output_everything(output_everything_),
         is_multipatch
         (CCTK_IsFunctionAliased("MultiPatch_GetSystemSpecification"))
     {
@@ -58,10 +173,40 @@ namespace CarpetIOF5 {
     
     void iterate (hid_t const file)
     {
-      int const rl = 0;
-      ENTER_LEVEL_MODE(cctkGH, rl) {
+      // Iterate over the variables in groups, first all grid
+      // functions, then all non-GF groups
+      group_type = CCTK_GF;
+      group_index = -1;
+      vindices.clear();
+      for (int vindex=0; vindex<CCTK_NumVars(); ++vindex) {
+        if (output_var.at(vindex)) {
+          int const gindex = CCTK_GroupIndexFromVarI(vindex);
+          if (CCTK_GroupTypeI(gindex) == CCTK_GF) {
+            vindices.push_back (vindex);
+          }
+        }
+      }
+      if (not vindices.empty()) {
         output_simulation (file);
-      } LEAVE_LEVEL_MODE;
+      }
+      
+      group_type = CCTK_ARRAY;
+      for (group_index=0; group_index<CCTK_NumGroups(); ++group_index) {
+        if (CCTK_GroupTypeI(group_index) != CCTK_GF) {
+          vindices.clear();
+          int const first_vindex = CCTK_FirstVarIndexI (group_index);
+          int const num_vars = CCTK_NumVarsInGroupI (group_index);
+          for (int vindex=first_vindex; vindex<first_vindex+num_vars; ++vindex)
+          {
+            if (output_var.at(vindex)) {
+              vindices.push_back (vindex);
+            }
+          }
+          if (not vindices.empty()) {
+            output_simulation (file);
+          }
+        }
+      }
     }
     
     
@@ -70,65 +215,86 @@ namespace CarpetIOF5 {
     
     void output_simulation (hid_t const file)
     {
-      assert (is_level_mode() and reflevel==0);
-      DECLARE_CCTK_ARGUMENTS;
+      DECLARE_CCTK_PARAMETERS;
       indent_t indent;
+      bool error_flag = false;
       
       gridname = generate_gridname(cctkGH);
-      chartname = generate_chartname(cctkGH);
-      
       cout << indent << "simulation=" << gridname << "\n";
       
-      ivect const reffact = spacereffacts.AT(reflevel);
-      F5Path *const globalpath = F5Rcreate_vertex_refinement3D
-        (file, cctk_time, gridname.c_str(), &v2h(reffact)[0],
-         chartname.c_str());
-      
-      // Choose (arbitrarily) the root level as default topology, for
-      // readers which don't understand AMR
-      F5Rlink_default_vertex_topology (globalpath, &v2h(reffact)[0]);
-      
-      // Define iteration
-      F5Rset_timestep (globalpath, cctk_iteration);
-      
-      // Attach Cactus/Carpet metadata
-      write_metadata (cctkGH, globalpath->Grid_hid);
-      
-      // Close topology
-      F5close (globalpath);
+      assert (is_global_mode());
       
-      // Iterate over all maps (on the root level)
-      BEGIN_MAP_LOOP(cctkGH, CCTK_GF) {
-        output_map (file);
-      } END_MAP_LOOP;
-    }
-    
-    
-    
-    void output_map (hid_t const file)
-    {
-      DECLARE_CCTK_ARGUMENTS;
-      indent_t indent;
-      
-      assert (is_singlemap_mode() and reflevel==0);
-      
-      cout << indent << "map=" << Carpet::map << "\n";
-      
-      fragmentname = generate_fragmentname(cctkGH, Carpet::map);
+      chartname = generate_chartname(cctkGH);
       
-      // Iterate over all refinement levels of this map
-      int const m = Carpet::map;
-      BEGIN_GLOBAL_MODE(cctkGH) {
-        BEGIN_REFLEVEL_LOOP(cctkGH) {
-          if (vhh.AT(m)->local_components(reflevel) > 0) {
-            // Continue only if this process owns a component of this
-            // map (on this level)
-            ENTER_SINGLEMAP_MODE(cctkGH, m, CCTK_GF) {
+      int const max_rl = group_type == CCTK_GF ? reflevels : 1;
+      for (int rl=0; rl<max_rl; ++rl) {
+        // Continue only if this level exists at this iteration
+        assert (maxtimereflevelfact % timereffacts.AT(rl) == 0);
+        int const do_every =
+          group_type == CCTK_GF ? maxtimereflevelfact / timereffacts.AT(rl) : 1;
+        if (cctkGH->cctk_iteration % do_every == 0) {
+          ENTER_LEVEL_MODE(cctkGH, rl) {
+            DECLARE_CCTK_ARGUMENTS;
+            
+            assert (timelevel == 0);
+            int const max_tl =
+              output_everything or output_all_timelevels ?
+              (group_type == CCTK_GF ?
+               timelevels :
+               CCTK_NumTimeLevelsI(group_index)) :
+              1;
+            for (timelevel=0; timelevel<max_tl; ++timelevel) {
+              cctkGH->cctk_time = tt->get_time (mglevel, reflevel, timelevel);
+              
+              // Choose (arbitrarily) the root level as default
+              // topology, for readers which don't understand AMR
+              if (reflevel == 0) {
+                ivect const reffact = spacereffacts.AT(reflevel);
+                F5Path *const globalpath = F5Rcreate_vertex_refinement3D
+                  (file, cctk_time, gridname.c_str(), &v2h(reffact)[0],
+                   chartname.c_str());
+                assert (globalpath);
+                FAILWARN (F5Rlink_default_vertex_topology (globalpath,
+                                                           &v2h(reffact)[0]));
+                
+                // Define iteration
+                FAILWARN (F5Rset_timestep (globalpath, cctk_iteration));
+                
+                // Attach Cactus/Carpet metadata
+                if (timelevel == 0) {
+                  // hid_t const metadata_group = globalpath->Grid_hid;
+                  ostringstream pathname;
+                  pathname << FIBER_CONTENT_GRIDS << "/" << gridname;
+                  hid_t group;
+                  group = H5Gopen (globalpath->ContentsGroup_hid,
+                                   pathname.str().c_str(),
+                                   H5P_DEFAULT);
+                  if (group < 0) {
+                    group = H5Gcreate (globalpath->ContentsGroup_hid,
+                                       pathname.str().c_str(),
+                                       H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+                  }
+                  assert (group >= 0);
+                  write_metadata (cctkGH, group);
+                  herr_t const herr = H5Gclose (group);
+                  assert (not herr);
+                }
+                
+                // Close topology
+                F5close (globalpath);
+                
+              } // if reflevel==0
+              
               output_reflevel (file);
-            } LEAVE_SINGLEMAP_MODE;
-          }
-        } END_REFLEVEL_LOOP;
-      } END_GLOBAL_MODE;
+              
+            } // for timelevel
+            timelevel = 0;
+            cctkGH->cctk_time = tt->get_time (mglevel, reflevel, timelevel);
+            
+          } LEAVE_LEVEL_MODE;
+        } // if do_every
+      }   // for rl
+      
     }
     
     
@@ -140,98 +306,120 @@ namespace CarpetIOF5 {
       
       cout << indent << "reflevel=" << reflevel << "\n";
       
+      assert (is_level_mode());
+      
       ivect const reffact = spacereffacts.AT(reflevel);
-      topologyname = generate_topologyname(cctkGH, Carpet::map, reffact);
+      topologyname = generate_topologyname(cctkGH, group_index, reffact);
       
       // Define grid hierarchy
-      int const indexdepth = 0; // vertices
+      map_indices_t const mi(cctkGH, group_index);
+      int const indexdepth = vhh.at(0)->refcent == vertex_centered ? 0 : 1;
       F5Path *const path =
         F5Rcreate_coordinate_topology (file, &cctk_time,
                                        gridname.c_str(), chartname.c_str(),
                                        topologyname.c_str(),
                                        indexdepth,
-                                       dim, dim, &v2h(reffact)[0]);
+                                       mi.dim, mi.dim, &v2h(reffact)[0]);
+      assert (path);
       
-      // Determine level coordinates
-      ivect gsh;
-      rvect origin, delta;
-      rvect lower, upper;
-      for (int d=0; d<dim; ++d) {
-        gsh[d]    = cctk_gsh[d];
-        origin[d] = CCTK_ORIGIN_SPACE(d);
-        delta[d]  = CCTK_DELTA_SPACE(d);
-        lower[d]  = origin[d];
-        upper[d]  = lower[d] + (gsh[d]-1) * delta[d];
-      }
+      BEGIN_LOCAL_MAP_LOOP (cctkGH, CCTK_GF) {
+        output_map (path);
+      } END_LOCAL_MAP_LOOP;
+      
+      // Close topology
+      F5close (path);
+    }
+    
+    
+    
+    void output_map (F5Path *const path)
+    {
+      DECLARE_CCTK_ARGUMENTS;
+      indent_t indent;
+      bool error_flag = false;
+      
+      cout << indent << "map=" << Carpet::map << "\n";
+      
+      assert (is_singlemap_mode());
       
       if (not is_multipatch) {
         // Define level geometry
-        F5_vec3_double_t const vlower = v2d(lower);
-        F5_vec3_double_t const vupper = v2d(upper);
-        F5_vec3_double_t const vdelta = v2d(delta);
-        F5Fwrite_linear (path, FIBER_HDF5_POSITIONS_STRING,
-                         dim, &v2h(gsh)[0],
-                         F5T_COORD3_DOUBLE,
-                         &vlower, &vdelta);
-        F5Fset_range (path, &vlower, &vupper);
+        map_indices_t const mi(cctkGH, group_index);
+        F5_vec3_double_t const vlower = v2d(mi.lower);
+        F5_vec3_double_t const vupper = v2d(mi.upper);
+        F5_vec3_double_t const vdelta = v2d(mi.delta);
+        static vector<ChartDomain_IDs*> charts;
+        if (charts.size() < dim+1) {
+          charts.resize(dim+1, NULL);
+          charts.at(3) = F5B_standard_cartesian_chart3D();
+        }
+        if (not charts.at(mi.dim)) {
+          assert (mi.dim != 0);
+          char const* coordnames[] = {"x", "y", "z"};
+          ostringstream chartnamebuf;
+          chartnamebuf << "Cartesian " << mi.dim << "D";
+          charts.at(mi.dim) =
+            F5B_new_global_float_chart(coordnames,
+                                       mi.dim, chartnamebuf.str().c_str(),
+                                       F5_FORTRAN_ORDER);
+          assert (charts.at(mi.dim));
+        }
+        hid_t const type = charts.at(mi.dim)->DoublePrecision.Point_hid_t;
+        FAILWARN (F5Fwrite_linear (path, FIBER_HDF5_POSITIONS_STRING,
+                                   mi.dim, &v2h(mi.gsh)[0],
+                                   type,
+                                   &vlower, &vdelta));
+#warning "TODO: path and chart don't match"
+        FAILWARN (F5Fset_range (path, &vlower, &vupper));
       }
       
       BEGIN_LOCAL_COMPONENT_LOOP (cctkGH, CCTK_GF) {
         output_component (path);
       } END_LOCAL_COMPONENT_LOOP;
-      
-      // Close topology
-      F5close (path);
     }
-  
-  
-  
+    
+    
+    
     void output_component (F5Path *const path)
     {
       DECLARE_CCTK_ARGUMENTS;
+      DECLARE_CCTK_PARAMETERS;
       indent_t indent;
+      bool error_flag = false;
       
-      cout << indent << "component=" << component << "\n";
+      cout << indent
+           << "component=" << component << " "
+           << "(local_component=" << local_component << ")\n";
       
-      // Determine component coordinates
-      ivect gsh;
-      ivect lbnd, lsh, lghosts, ughosts;
-      rvect origin, delta;
-      rvect clower, cupper;
-      for (int d=0; d<dim; ++d) {
-        gsh[d]     = cctk_gsh[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];
-        origin[d]  = CCTK_ORIGIN_SPACE(d);
-        delta[d]   = CCTK_DELTA_SPACE(d);
-        clower[d]  = origin[d] + cctk_lbnd[d] * delta[d];
-        cupper[d]  = clower[d] + (lsh[d]-1) * delta[d];
-      }
+      assert (is_local_mode());
+      
+      fragmentname = generate_fragmentname(cctkGH, Carpet::map, component);
       
       // Define coordinates
-      if (not is_multipatch) {
+      // TODO: also define and use is_cartesian for each map
+      if (not is_multipatch and group_type == CCTK_GF) {
         // (This is redundant, since the level's overall bounding box
         // was already defined above, but it provides the individual
         // components' bounding boxes.)
-        F5Fwrite_linear_fraction (path, FIBER_HDF5_POSITIONS_STRING,
-                                  cctk_dim, &v2h(gsh)[0], &v2h(lsh)[0],
-                                  F5T_COORD3_DOUBLE,
-                                  &clower, &delta, &v2h(lbnd)[0],
-                                  fragmentname.c_str());
+        component_indices_t const ci(cctkGH, group_index);
+        FAILWARN (F5Fwrite_linear_fraction (path, FIBER_HDF5_POSITIONS_STRING,
+                                            ci.dim,
+                                            &v2h(ci.gsh)[0], &v2h(ci.ilen)[0],
+                                            F5T_COORD3_DOUBLE,
+                                            &ci.clower, &ci.delta,
+                                            &v2h(ci.ioff)[0],
+                                            fragmentname.c_str()));
       } else {
+        // Output coordinates
         output_variable (path, CCTK_VarIndex("grid::x"), true);
       }
       
-      // Output some variables
-      output_variable (path, CCTK_VarIndex("grid::r"));
-      
-      output_variable (path, CCTK_VarIndex("grid::x"));
-      output_variable (path, CCTK_VarIndex("grid::y"));
-      output_variable (path, CCTK_VarIndex("grid::z"));
+      // Output variables
+      for (vector<int>::const_iterator
+             vi = vindices.begin(); vi != vindices.end(); ++vi)
+      {
+        output_variable (path, *vi);
+      }
     }
     
     
@@ -242,13 +430,19 @@ namespace CarpetIOF5 {
       DECLARE_CCTK_ARGUMENTS;
       DECLARE_CCTK_PARAMETERS;
       indent_t indent;
+      bool error_flag = false;
       
       int ierr;
       
       assert (var >= 0);
-      {
+      if (write_positions) {
+        cout << indent << "positions\n";
+      } else {
         char *const fullname = CCTK_FullName(var);
-        cout << indent << "variable=" << fullname << "\n";
+        char *const groupname = CCTK_GroupNameFromVarI(var);
+        cout << indent
+             << "variable=" << fullname << " (group=" << groupname << ")\n";
+        free (groupname);
         free (fullname);
       }
       
@@ -259,49 +453,36 @@ namespace CarpetIOF5 {
       ierr = CCTK_GroupData (group, &groupdata);
       assert (not ierr);
       
-      // TODO
-      assert (groupdata.grouptype == CCTK_GF);
-      assert (groupdata.vartype == CCTK_VARIABLE_REAL);
-      assert (groupdata.disttype == CCTK_DISTRIB_DEFAULT);
-      assert (groupdata.stagtype == 0);
-      assert (groupdata.dim == cctk_dim);
-      
-      int const timelevel = 0;
+      assert ((groupdata.grouptype == CCTK_GF) == (group_type == CCTK_GF));
       
-      // Determine component coordinates
-      int const dim = cctk_dim;
-      ivect gsh;
-      ivect lbnd, lsh, lghosts, ughosts;
-      ivect imin, imax, ioff, ilen;
-      ivect const izero(0);
-      for (int d=0; d<dim; ++d) {
-        gsh[d]  = cctk_gsh[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];
-        if (not output_ghost_points) {
-          imin[d] += lghosts[d] / 2;
-          imax[d] -= ughosts[d] / 2;
-          lghosts[d] = 0;
-          ughosts[d] = 0;
-        }
-        ioff[d] = lbnd[d] + imin[d];
-        ilen[d] = imax[d] - imin[d];
+      // Output distrib=constant variables only on process 0
+      switch (groupdata.disttype) {
+      case CCTK_DISTRIB_CONSTANT:
+        if (CCTK_MyProc(cctkGH) != 0) return;
+        break;
+      case CCTK_DISTRIB_DEFAULT:
+        // do nothing
+        break;
+      default:
+        assert (0);
       }
       
+      assert (groupdata.stagtype == 0);
+      assert (groupdata.dim == dim);
+      
 #warning "TODO: Do not output symmetry zones (unless requested by the user)"
 #warning "TODO: Do not output buffer zones (is that easily possible?)"
       
+      int const will_cover_complete_domain = not is_multipatch and reflevel==0;
+      
+      cGroupDynamicData dyndata;
+      ierr = CCTK_GroupDynamicData (cctkGH, group, &dyndata);
+      assert (not ierr);
+      
+      // Only output active timelevels
+      if (timelevel >= dyndata.activetimelevels) return;
       
       
-      enum tensortype_t {
-        tt_scalar, tt_vector, tt_error
-      };
       
       tensortype_t tensortype = tt_error;
       
@@ -327,6 +508,10 @@ namespace CarpetIOF5 {
           tensortype = tt_scalar; break;
         case 3:
           tensortype = tt_vector; break;
+        case 6:
+          tensortype = tt_symtensor; break;
+        case 9:
+          tensortype = tt_tensor; break;
         default:
           // fallback: output group as scalars
           tensortype = tt_scalar; break;
@@ -337,129 +522,174 @@ namespace CarpetIOF5 {
       if (tensortype != tt_scalar and var != v0) {
         // TODO: don't do this; instead, keep track of which groups
         // have been output
-        char *const fullname = CCTK_FullName(var);
-        indent_t indent;
-        cout << indent << "skipping output since it is not the first variable in its group\n";
-        free (fullname);
+        indent_t indent2;
+        cout << indent2 << "skipping output since it is not the first variable in its group\n";
         return;
       }
       
       
       
+      hid_t type = -1;
+      int num_comps = 0;
+      string name;
+      
       switch (tensortype) {
         
       case tt_scalar: {
         // Scalar field
-        CCTK_REAL const *const rdata =
-          (CCTK_REAL const*)CCTK_VarDataPtrI (cctkGH, timelevel, var);
-        assert (rdata);
+        assert (not write_positions);
+        switch (groupdata.vartype) {
+        case CCTK_VARIABLE_INT:  type = H5T_NATIVE_INT;    break;
+        case CCTK_VARIABLE_REAL: type = H5T_NATIVE_DOUBLE; break;
+        default: assert(0);
+        }
+        num_comps = 1;
+        name = generate_fieldname (cctkGH, var, tensortype);
+        break;
+      }
         
-        int const vartype = CCTK_VarTypeI(var);
-        assert (vartype == CCTK_VARIABLE_REAL);
-        vector<CCTK_REAL> idata (prod(ilen));
-        for (int k=0; k<ilen[2]; ++k) {
-          for (int j=0; j<ilen[1]; ++j) {
-            for (int i=0; i<ilen[0]; ++i) {
-              int const isrc =
-                CCTK_GFINDEX3D(cctkGH, imin[0]+i,imin[1]+j,imin[2]+k);
-              int const idst = i + ilen[0] * (j + ilen[1] * k);
-              idata[idst] = rdata[isrc];
+      case tt_vector: {
+        // Vector field, or positions
+        switch (groupdata.vartype) {
+        case CCTK_VARIABLE_REAL:
+          type = write_positions ? F5T_COORD3_DOUBLE : F5T_VEC3_DOUBLE;
+          break;
+        default: assert(0);
+        }
+        num_comps = dim;
+        name =
+          write_positions ?
+          FIBER_HDF5_POSITIONS_STRING :
+          generate_fieldname (cctkGH, var, tensortype);
+        if (write_positions) {
+          htri_t const exists =
+            H5Lexists (path->Representation_hid, name.c_str(), H5P_DEFAULT);
+          assert (exists >= 0);
+          if (exists) {
+            string const fragmentpath = name + "/" + fragmentname;
+            htri_t const exists2 =
+              H5Lexists (path->Representation_hid, fragmentpath.c_str(),
+                         H5P_DEFAULT);
+            assert (exists2 >= 0);
+            if (exists2) {
+              // Write positions only once
+              indent_t indent2;
+              cout << indent2 << "skipping output since the positions have already been output\n";
+              return;
             }
           }
         }
-        void const *const data = &idata.front();
+        break;
+      }
         
-        char *const groupname = CCTK_GroupName(group);
-        char *const fullname = CCTK_FullName(var);
-        // Use the variable name instead of the group name if we may
-        // output several variables per group
+      case tt_symtensor: {
+        // Symmetric tensor field
         assert (not write_positions);
-        char const *const name = groupdata.numvars == 1 ? groupname : fullname;
-#if 0
-        F5Fwrite_fraction (path, name,
-                           dim,
-                           is_multipatch ? NULL : &v2h(gsh)[0], &v2h(lsh)[0],
-                           H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE,
-                           data,
-                           &v2h(lbnd)[0], &v2h(lghosts)[0], &v2h(ughosts)[0],
-                           fragmentname.c_str(), H5P_DEFAULT);
-#endif
-        F5Fwrite_fraction (path, name,
-                           dim,
-                           is_multipatch ? NULL : &v2h(gsh)[0], &v2h(ilen)[0],
-                           H5T_NATIVE_DOUBLE, H5T_NATIVE_DOUBLE,
-                           data,
-                           &v2h(ioff)[0], &v2h(izero)[0], &v2h(izero)[0],
-                           fragmentname.c_str(), H5P_DEFAULT);
-        free (groupname);
-        free (fullname);
+        switch (groupdata.vartype) {
+        case CCTK_VARIABLE_REAL: type = F5T_METRIC33_DOUBLE; break;
+        default: assert(0);
+        }
+        num_comps = dim*(dim+1)/2;
+        name = generate_fieldname (cctkGH, var, tensortype);
         break;
       }
         
-      case tt_vector: {
-        // Vector field
-        CCTK_REAL const* rdata[cctk_dim];
-        vector<vector<CCTK_REAL> > idata(cctk_dim);
-        void const* data[cctk_dim];
-        for (int d=0; d<cctk_dim; ++d) {
-          rdata[d] =
-            (CCTK_REAL const*)CCTK_VarDataPtrI (cctkGH, timelevel, v0+d);
-          assert (rdata[d]);
-          
-          int const vartype = CCTK_VarTypeI(var);
-          assert (vartype == CCTK_VARIABLE_REAL);
-          idata[d].resize (prod(ilen));
-          for (int k=0; k<ilen[2]; ++k) {
-            for (int j=0; j<ilen[1]; ++j) {
-              for (int i=0; i<ilen[0]; ++i) {
-                int const isrc =
-                  CCTK_GFINDEX3D(cctkGH, imin[0]+i,imin[1]+j,imin[2]+k);
-                int const idst = i + ilen[0] * (j + ilen[1] * k);
-                idata[d][idst] = rdata[d][isrc];
-              }
-            }
-          }
-          data[d] = &idata[d].front();
-          
+      case tt_tensor: {
+        // Non-symmetric tensor field
+        assert (not write_positions);
+        switch (groupdata.vartype) {
+        case CCTK_VARIABLE_REAL: type = F5T_BIVEC3_DOUBLE; break;
+        default: assert(0);
         }
-        char *const groupname = CCTK_GroupName(group);
-        char const *const name =
-          write_positions ? FIBER_HDF5_POSITIONS_STRING : groupname;
-        hid_t const type =
-          write_positions ? F5T_COORD3_DOUBLE : F5T_VEC3_DOUBLE;
-        int const will_cover_complete_domain = 0;
-#if 0
-        F5FSwrite_fraction (path, name,
-                            dim,
-                            is_multipatch ? NULL : &v2h(gsh)[0], &v2h(lsh)[0],
-                            type, type,
-                            data,
-                            &v2h(lbnd)[0], &v2h(lghosts)[0], &v2h(ughosts)[0],
-                            fragmentname.c_str(), H5P_DEFAULT,
-                            will_cover_complete_domain);
-#endif
-        F5FSwrite_fraction (path, name,
-                            dim,
-                            is_multipatch ? NULL : &v2h(gsh)[0], &v2h(ilen)[0],
-                            type, type,
-                            data,
-                            &v2h(ioff)[0], &v2h(izero)[0], &v2h(izero)[0],
-                            fragmentname.c_str(), H5P_DEFAULT,
-                            will_cover_complete_domain);
-        free (groupname);
+        num_comps = dim*dim;
+        name = generate_fieldname (cctkGH, var, tensortype);
         break;
       }
         
       default:
         assert (0);
       }
+      
+      // Write data
+      assert (type >= 0);
+      assert (num_comps > 0);
+      assert (name.size() > 0);
+      
+      component_indices_t const ci(cctkGH, group_index);
+      
+      CCTK_REAL const* rdata[num_comps];
+      vector<vector<CCTK_REAL> > idata(num_comps);
+      void const* data[num_comps];
+      for (int d=0; d<num_comps; ++d) {
+        rdata[d] = (CCTK_REAL const*)CCTK_VarDataPtrI(cctkGH, timelevel, var+d);
+        assert (rdata[d]);
+        
+        int const vartype = CCTK_VarTypeI(var);
+        assert (vartype == CCTK_VARIABLE_REAL);
+        idata[d].resize (prod(ci.ilen));
+        for (int k=0; k<ci.ilen[2]; ++k) {
+          for (int j=0; j<ci.ilen[1]; ++j) {
+            for (int i=0; i<ci.ilen[0]; ++i) {
+              int const isrc =
+                ci.imin[0]+i + ci.lsh[0] *
+                (ci.imin[1]+j + ci.lsh[1] * (ci.imin[2]+k));
+              int const idst = i + ci.ilen[0] * (j + ci.ilen[1] * k);
+              idata[d][idst] = rdata[d][isrc];
+            }
+          }
+        }
+        data[d] = &idata[d].front();
+      }
+      
+      // Dataset properties
+      hid_t const prop = H5Pcreate (H5P_DATASET_CREATE);
+      assert (prop >= 0);
+      // Enable compression if requested
+      if (compression_level >= 0) {
+        FAILWARN (H5Pset_chunk (prop, dim, &v2h(ci.ilen).reverse()[0]));
+        FAILWARN (H5Pset_deflate (prop, compression_level));
+      }
+      // Enable checksums if requested
+      if (use_checksums) {
+        FAILWARN (H5Pset_chunk (prop, dim, &v2h(ci.ilen).reverse()[0]));
+        FAILWARN (H5Pset_fletcher32 (prop));
+      }
+      
+      if (num_comps == 1 and not separate_single_component_tensors) {
+        // Write single-component tensors into non-separated fractions
+        // for convenience (could also use a separated compound
+        // instead)
+        FAILWARN
+          (F5Fwrite_fraction (path, name.c_str(),
+                              dim,
+                              is_multipatch ? NULL : &v2h(ci.gsh)[0],
+                              &v2h(ci.ilen)[0],
+                              type, type, data[0],
+                              &v2h(ci.ioff)[0],
+                              &v2h(ci.lghosts)[0], &v2h(ci.ughosts)[0],
+                              fragmentname.c_str(), prop));
+      } else {
+        FAILWARN
+          (F5FSwrite_fraction (path, name.c_str(),
+                               dim,
+                               is_multipatch ? NULL : &v2h(ci.gsh)[0],
+                               &v2h(ci.ilen)[0],
+                               type, type, data,
+                               &v2h(ci.ioff)[0],
+                               &v2h(ci.lghosts)[0], &v2h(ci.ughosts)[0],
+                               fragmentname.c_str(), prop,
+                               will_cover_complete_domain));
+      }
+      
+      FAILWARN (H5Pclose (prop));
+      
     }
     
   };                            // class output_iterator_t
   
   
   
-  void write_metadata (cGH *const cctkGH, hid_t const file)
+  void write_metadata (cGH const* const cctkGH, hid_t const file)
   {
     DECLARE_CCTK_PARAMETERS;
     
@@ -475,6 +705,12 @@ namespace CarpetIOF5 {
     // files (to save space), or write it into a separate metadata
     // file
     
+    // Create metadata only once
+    // TODO: update metadata at some point
+    htri_t const exists = H5Lexists (file, metadata_group, H5P_DEFAULT);
+    assert (exists >= 0);
+    if (exists) return;
+    
     // Create a group to hold all metadata
     hid_t const group =
       H5Gcreate (file, metadata_group, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
@@ -501,10 +737,14 @@ namespace CarpetIOF5 {
         (group, "run id", (char const*) UniqueRunID(cctkGH));
     }
     
+    // Don't write this attribute; the number of files may change
+    // after recovering, or after recombining files
+#if 0
     // Number of I/O processes (i.e. the number of output files)
     int const nprocs = CCTK_nProcs(cctkGH);
     int const nioprocs = nprocs;
     WriteAttribute (group, "nioprocs", nioprocs);
+#endif
     
     // Write parameters into a separate dataset (they may be too large
     // for an attribute)
@@ -524,53 +764,14 @@ namespace CarpetIOF5 {
   
   
   
-  extern "C"
-  void F5_Output (CCTK_ARGUMENTS)
+  void output (cGH const* const cctkGH,
+               hid_t const file,
+               vector<bool> const& output_var,
+               bool const output_everything)
   {
-    DECLARE_CCTK_ARGUMENTS;
-    DECLARE_CCTK_PARAMETERS;
-    
-    herr_t herr;
-    
-    
-    
-    assert (is_global_mode());
-    CCTK_VInfo (CCTK_THORNSTRING, "F5_Output: iteration=%d", cctk_iteration);
-    
-    
-    
-    // We don't know how to open multiple files yet
-    assert (strcmp (file_content, "everything") == 0);
-    
-    // Open file
-    static bool first_time = true;
-    
-    string const basename =
-      generate_basename (cctkGH, CCTK_VarIndex("grid::r"));
-    int const myproc = CCTK_MyProc(cctkGH);
-    int const proc = myproc;
-    string const name =
-      create_filename (cctkGH, basename, proc, first_time);
-    
-    indent_t indent;
-    cout << indent << "process=" << proc << "\n";
-    
-    bool const truncate_file = first_time and IO_TruncateOutputFiles(cctkGH);
-    hid_t const file =
-      truncate_file ?
-      H5Fcreate (name.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT) :
-      H5Fopen (name.c_str(), H5F_ACC_RDWR, H5P_DEFAULT);
-    assert (file >= 0);
-    first_time = false;
-    
-    // write_metadata (cctkGH, file);
-    
-    output_iterator_t iterator(cctkGH);
+    output_iterator_t
+      iterator(const_cast<cGH*>(cctkGH), output_var, output_everything);
     iterator.iterate (file);
-    
-    // Close file
-    herr = H5Fclose (file);
-    assert (not herr);
   }
   
 } // end namespace CarpetIOF5
diff --git a/CarpetDev/CarpetIOF5/src/util.cc b/CarpetDev/CarpetIOF5/src/util.cc
index 57ccfd6db..15e01fba7 100644
--- a/CarpetDev/CarpetIOF5/src/util.cc
+++ b/CarpetDev/CarpetIOF5/src/util.cc
@@ -2,7 +2,9 @@
 #include <cctk_Arguments.h>
 #include <cctk_Parameters.h>
 
+#include <algorithm>
 #include <cassert>
+#include <cctype>
 #include <cstdio>
 #include <cstdlib>
 #include <cstring>
@@ -95,7 +97,26 @@ namespace CarpetIOF5 {
       free (thornname);
     }
     else if (CCTK_EQUALS (file_content, "everything")) {
-      filename_buf << out_filename;
+      if (CCTK_EQUALS(out_filename, "")) {
+        // Obtain the parameter file name
+        char buf[10000];
+        int const ilen = CCTK_ParameterFilename (sizeof buf, buf);
+        assert (ilen < int(sizeof buf) - 1);
+        char* parfilename = buf;
+        // Remove directory prefix, if any
+        char* const slash = strrchr(parfilename, '/');
+        if (slash) {
+          parfilename = &slash[1];
+        }
+        // Remove suffix, if it is there
+        char* const suffix = strrchr(parfilename, '.');
+        if (suffix and strcmp(suffix, ".par")==0) {
+          suffix[0] = '\0';
+        }
+        filename_buf << parfilename;
+      } else {
+        filename_buf << out_filename;
+      }
     }
     else {
       assert (0);
@@ -200,7 +221,33 @@ namespace CarpetIOF5 {
     return p;
   }
   
-  // Generate a good tensor basis name (chart name)
+  // Generate a good topology name (refinement level name)
+  string
+  generate_topologyname (cGH const *const cctkGH,
+                         int const gi,
+                         ivect const& reffact)
+  {
+    ostringstream buf;
+    if (gi == -1) {
+      // grid function
+      buf << "VertexLevel_";
+      for (int d=0; d<dim; ++d) {
+        if (d>0) buf << "x";
+        buf << reffact[d];
+      }
+    } else {
+      // grid array
+      char *const groupname = CCTK_GroupName(gi);
+      for (char *p=groupname; *p; ++p) {
+        *p = tolower(*p);
+      }
+      buf << "Group_" << groupname;
+      free (groupname);
+    }
+    return buf.str();
+  }
+  
+  // Generate a good chart name (tensor basis name)
   string
   generate_chartname (cGH const *const cctkGH)
   {
@@ -216,38 +263,55 @@ namespace CarpetIOF5 {
 #endif
   }
   
-  // Generate a good fragment name (map name)
+  // Generate a good fragment name (map and component name)
+  // (We assume that fragment names need to be unique only within a
+  // topology, not across topologies.)
   string
-  generate_fragmentname (cGH const *const cctkGH, int const m)
+  generate_fragmentname (cGH const *const cctkGH, int const m, int const c)
   {
     ostringstream buf;
-    buf << "Map_" << m;
+    buf << "Map_" << m << "_Component_" << c;
     return buf.str();
   }
   
-  int
+  void
   interpret_fragmentname (cGH const *const cctkGH,
-                          char const *const fragmentname)
+                          char const *const fragmentname,
+                          int& m, int& c)
   {
-    int m = -1;
-    sscanf (fragmentname, "Map_%d", &m);
+    m = -1;
+    c = -1;
+    sscanf (fragmentname, "Map_%d_Component_%d", &m, &c);
     assert (m>=0 and m<Carpet::maps);
-    return m;
+    assert (c>=0);
   }
   
-  // Generate a good topology name (map and refinement level name)
+  // Generate a good field name (group or variable name)
   string
-  generate_topologyname (cGH const *const cctkGH,
-                         int const m, ivect const& reffact)
+  generate_fieldname (cGH const *const cctkGH,
+                      int const vi, tensortype_t const tt)
   {
-    ostringstream buf;
-    // buf << "Map_" << m << "_VertexLevel_";
-    buf << "VertexLevel_";
-    for (int d=0; d<dim; ++d) {
-      if (d>0) buf << "x";
-      buf << reffact[d];
+    int const gi = CCTK_GroupIndexFromVarI(vi);
+    int const numvars = CCTK_NumVarsInGroupI(gi);
+    string name;
+    // Use the variable name instead of the group name if we may
+    // output several variables per group
+    if (tt == tt_scalar and numvars >1) {
+      char *const fullname = CCTK_FullName(vi);
+      name = fullname;
+      free (fullname);
+    } else {
+      char *const groupname = CCTK_GroupName(gi);
+      name = groupname;
+      free (groupname);
     }
-    return buf.str();
+    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(), ".");
+    }
+    return name;
   }