CarpetIOF5: New thorn

17 files changed, 2970 insertions, 0 deletions

diff --git a/CarpetDev/CarpetIOF5/README b/CarpetDev/CarpetIOF5/README
new file mode 100644
index 000000000..d17fddf2a
--- /dev/null
+++ b/CarpetDev/CarpetIOF5/README
@@ -0,0 +1,9 @@
+Cactus Code Thorn CarpetIOF5
+Author(s)    : Erik Schnetter <schnetter@cct.lsu.edu>
+Maintainer(s): Erik Schnetter <schnetter@cct.lsu.edu>
+Licence      : GPL
+--------------------------------------------------------------------------
+
+1. Purpose
+
+Efficient, binary I/O based on the F5 library
diff --git a/CarpetDev/CarpetIOF5/configuration.ccl b/CarpetDev/CarpetIOF5/configuration.ccl
new file mode 100644
index 000000000..155a0c55b
--- /dev/null
+++ b/CarpetDev/CarpetIOF5/configuration.ccl
@@ -0,0 +1,3 @@
+# Configuration definition for thorn CarpetIOF5
+
+REQUIRES Carpet CarpetLib F5 HDF5 LoopControl
diff --git a/CarpetDev/CarpetIOF5/doc/documentation.tex b/CarpetDev/CarpetIOF5/doc/documentation.tex
new file mode 100644
index 000000000..e73a6e7ed
--- /dev/null
+++ b/CarpetDev/CarpetIOF5/doc/documentation.tex
@@ -0,0 +1,144 @@
+% *======================================================================*
+%  Cactus Thorn template for ThornGuide documentation
+%  Author: Ian Kelley
+%  Date: Sun Jun 02, 2002
+%  $Header$
+%
+%  Thorn documentation in the latex file doc/documentation.tex
+%  will be included in ThornGuides built with the Cactus make system.
+%  The scripts employed by the make system automatically include
+%  pages about variables, parameters and scheduling parsed from the
+%  relevant thorn CCL files.
+%
+%  This template contains guidelines which help to assure that your
+%  documentation will be correctly added to ThornGuides. More
+%  information is available in the Cactus UsersGuide.
+%
+%  Guidelines:
+%   - Do not change anything before the line
+%       % START CACTUS THORNGUIDE",
+%     except for filling in the title, author, date, etc. fields.
+%        - Each of these fields should only be on ONE line.
+%        - Author names should be separated with a \\ or a comma.
+%   - You can define your own macros, but they must appear after
+%     the START CACTUS THORNGUIDE line, and must not redefine standard
+%     latex commands.
+%   - To avoid name clashes with other thorns, 'labels', 'citations',
+%     'references', and 'image' names should conform to the following
+%     convention:
+%       ARRANGEMENT_THORN_LABEL
+%     For example, an image wave.eps in the arrangement CactusWave and
+%     thorn WaveToyC should be renamed to CactusWave_WaveToyC_wave.eps
+%   - Graphics should only be included using the graphicx package.
+%     More specifically, with the "\includegraphics" command.  Do
+%     not specify any graphic file extensions in your .tex file. This
+%     will allow us to create a PDF version of the ThornGuide
+%     via pdflatex.
+%   - References should be included with the latex "\bibitem" command.
+%   - Use \begin{abstract}...\end{abstract} instead of \abstract{...}
+%   - Do not use \appendix, instead include any appendices you need as
+%     standard sections.
+%   - For the benefit of our Perl scripts, and for future extensions,
+%     please use simple latex.
+%
+% *======================================================================*
+%
+% Example of including a graphic image:
+%    \begin{figure}[ht]
+% 	\begin{center}
+%    	   \includegraphics[width=6cm]{MyArrangement_MyThorn_MyFigure}
+% 	\end{center}
+% 	\caption{Illustration of this and that}
+% 	\label{MyArrangement_MyThorn_MyLabel}
+%    \end{figure}
+%
+% Example of using a label:
+%   \label{MyArrangement_MyThorn_MyLabel}
+%
+% Example of a citation:
+%    \cite{MyArrangement_MyThorn_Author99}
+%
+% Example of including a reference
+%   \bibitem{MyArrangement_MyThorn_Author99}
+%   {J. Author, {\em The Title of the Book, Journal, or periodical}, 1 (1999),
+%   1--16. {\tt http://www.nowhere.com/}}
+%
+% *======================================================================*
+
+% If you are using CVS use this line to give version information
+% $Header$
+
+\documentclass{article}
+
+% Use the Cactus ThornGuide style file
+% (Automatically used from Cactus distribution, if you have a
+%  thorn without the Cactus Flesh download this from the Cactus
+%  homepage at www.cactuscode.org)
+\usepackage{../../../../doc/latex/cactus}
+
+\begin{document}
+
+% The author of the documentation
+\author{Erik Schnetter \textless schnetter@cct.lsu.edu\textgreater}
+
+% The title of the document (not necessarily the name of the Thorn)
+\title{CarpetIOF5}
+
+% the date your document was last changed, if your document is in CVS,
+% please use:
+%    \date{$ $Date: 2004-01-07 14:12:39 -0600 (Wed, 07 Jan 2004) $ $}
+\date{July 16 2010}
+
+\maketitle
+
+% Do not delete next line
+% START CACTUS THORNGUIDE
+
+% Add all definitions used in this documentation here
+%   \def\mydef etc
+
+% Add an abstract for this thorn's documentation
+\begin{abstract}
+
+\end{abstract}
+
+% The following sections are suggestive only.
+% Remove them or add your own.
+
+\section{Introduction}
+
+\section{Physical System}
+
+\section{Numerical Implementation}
+
+\section{Using This Thorn}
+
+\subsection{Obtaining This Thorn}
+
+\subsection{Basic Usage}
+
+\subsection{Special Behaviour}
+
+\subsection{Interaction With Other Thorns}
+
+\subsection{Examples}
+
+\subsection{Support and Feedback}
+
+\section{History}
+
+\subsection{Thorn Source Code}
+
+\subsection{Thorn Documentation}
+
+\subsection{Acknowledgements}
+
+
+\begin{thebibliography}{9}
+
+\end{thebibliography}
+
+% Do not delete next line
+% END CACTUS THORNGUIDE
+
+\end{document}
diff --git a/CarpetDev/CarpetIOF5/interface.ccl b/CarpetDev/CarpetIOF5/interface.ccl
new file mode 100644
index 000000000..4c09d8b12
--- /dev/null
+++ b/CarpetDev/CarpetIOF5/interface.ccl
@@ -0,0 +1,46 @@
+# Interface definition for thorn CarpetIOF5
+
+IMPLEMENTS: CarpetIOF5
+
+INHERITS: grid
+
+USES INCLUDE: carpet.hh
+
+
+
+# Return a pointer to an unmodifiable C string
+# which contains a unique ID for this configuration
+CCTK_POINTER_TO_CONST \
+FUNCTION UniqueConfigID (CCTK_POINTER_TO_CONST IN cctkGH)
+USES FUNCTION UniqueConfigID
+
+# Return a pointer to an unmodifiable C string
+# which contains a unique ID for this build
+CCTK_POINTER_TO_CONST \
+FUNCTION UniqueBuildID (CCTK_POINTER_TO_CONST IN cctkGH)
+USES FUNCTION UniqueBuildID
+
+# Return a pointer to an unmodifiable C string
+# which contains a unique ID for this simulation
+CCTK_POINTER_TO_CONST \
+FUNCTION UniqueSimulationID (CCTK_POINTER_TO_CONST IN cctkGH)
+USES FUNCTION UniqueSimulationID
+
+# Return a pointer to an unmodifiable C string
+# which contains a unique ID for this run
+CCTK_POINTER_TO_CONST \
+FUNCTION UniqueRunID (CCTK_POINTER_TO_CONST IN cctkGH)
+USES FUNCTION UniqueRunID
+
+
+
+# Check whether existing output files should be truncated or not
+CCTK_INT FUNCTION IO_TruncateOutputFiles (CCTK_POINTER_TO_CONST IN cctkGH)
+REQUIRES FUNCTION IO_TruncateOutputFiles
+
+
+
+# The setup of the system
+CCTK_INT FUNCTION MultiPatch_GetSystemSpecification \
+    (CCTK_INT OUT maps)
+USES FUNCTION MultiPatch_GetSystemSpecification
diff --git a/CarpetDev/CarpetIOF5/par/iof5-multipatch.par b/CarpetDev/CarpetIOF5/par/iof5-multipatch.par
new file mode 100644
index 000000000..71550e788
--- /dev/null
+++ b/CarpetDev/CarpetIOF5/par/iof5-multipatch.par
@@ -0,0 +1,45 @@
+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"
diff --git a/CarpetDev/CarpetIOF5/par/iof5-refined.par b/CarpetDev/CarpetIOF5/par/iof5-refined.par
new file mode 100644
index 000000000..aa27521b4
--- /dev/null
+++ b/CarpetDev/CarpetIOF5/par/iof5-refined.par
@@ -0,0 +1,48 @@
+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 = 0
+
+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"
diff --git a/CarpetDev/CarpetIOF5/par/iof5-uniform.par b/CarpetDev/CarpetIOF5/par/iof5-uniform.par
new file mode 100644
index 000000000..87e758370
--- /dev/null
+++ b/CarpetDev/CarpetIOF5/par/iof5-uniform.par
@@ -0,0 +1,39 @@
+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
diff --git a/CarpetDev/CarpetIOF5/param.ccl b/CarpetDev/CarpetIOF5/param.ccl
new file mode 100644
index 000000000..f9da760a4
--- /dev/null
+++ b/CarpetDev/CarpetIOF5/param.ccl
@@ -0,0 +1,66 @@
+# Parameter definitions for thorn CarpetIOF5
+
+SHARES: IO
+
+USES STRING out_dir AS IO_out_dir
+
+USES INT    out_timesteps_per_file
+
+
+
+PRIVATE:
+
+STRING out_dir "Output directory (overrides IO::out_dir)" STEERABLE=always
+{
+  "^$" :: "Empty: use IO::out_dir"
+  ".+" :: "Not empty: directory name"
+} ""
+
+KEYWORD file_content "Create one file for every x" STEERABLE=always
+{
+  "group"      :: ""
+  "thorn"      :: ""
+  "everything" :: ""
+} "everything"
+
+INT iteration_digits "Minimum number of digits for iteration number" STEERABLE=always
+{
+  0:* :: ""
+} 10
+
+STRING out_filename "File name (without extension) for file_content='everything'" STEERABLE=always
+{
+  "" :: ""
+} "output"
+
+INT processor_digits "Minimum number of digits for processor number" STEERABLE=always
+{
+  0:* :: ""
+} 6
+
+STRING out_extension "File name extension" STEERABLE=always
+{
+  "" :: "File extension (including a leading dot, if desired)"
+} ".f5"
+
+BOOLEAN create_subdirs "Create subdirectories for the output files to reduce the number of files per directory" STEERABLE=always
+{
+} "no"
+
+BOOLEAN one_dir_per_file "Create one subdirectory per output file to reduce locking overhead" STEERABLE=always
+{
+} "no"
+
+
+
+BOOLEAN output_symmetry_points "Output symmetry and inter-patch boundary points" STEERABLE=always
+{
+} "no"
+
+BOOLEAN output_ghost_points "Output ghost points" STEERABLE=always
+{
+} "no"
+
+BOOLEAN output_boundary_points "Output outer boundary points" STEERABLE=always
+{
+} "yes"
diff --git a/CarpetDev/CarpetIOF5/schedule.ccl b/CarpetDev/CarpetIOF5/schedule.ccl
new file mode 100644
index 000000000..494ce02fe
--- /dev/null
+++ b/CarpetDev/CarpetIOF5/schedule.ccl
@@ -0,0 +1,27 @@
+# Schedule definitions for thorn CarpetIOF5
+
+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/attributes.cc b/CarpetDev/CarpetIOF5/src/attributes.cc
new file mode 100644
index 000000000..9ca1f109c
--- /dev/null
+++ b/CarpetDev/CarpetIOF5/src/attributes.cc
@@ -0,0 +1,184 @@
+#include <cctk.h>
+
+#include <algorithm>
+#include <cassert>
+#include <cstring>
+#include <vector>
+
+#include <hdf5.h>
+
+#include <defs.hh>
+
+#include "iof5.hh"
+
+
+
+namespace CarpetIOF5 {
+  
+  using namespace std;
+  
+  
+  
+  // Write an int attribute
+  bool WriteAttribute (hid_t const group,
+                       char const *const name,
+                       int const ivalue)
+  {
+    bool error_flag = false;
+    
+    hid_t const dataspace = FAILWARN (H5Screate (H5S_SCALAR));
+    hid_t const attribute =
+      FAILWARN (H5Acreate (group, name, H5T_NATIVE_INT, dataspace,
+                           H5P_DEFAULT, H5P_DEFAULT));
+    FAILWARN (H5Awrite (attribute, H5T_NATIVE_INT, &ivalue));
+    FAILWARN (H5Aclose (attribute));
+    FAILWARN (H5Sclose (dataspace));
+    
+    return error_flag;
+  }
+  
+  
+  
+  // Write a double attribute
+  bool WriteAttribute (hid_t const group,
+                       char const *const name,
+                       double const dvalue)
+  {
+    bool error_flag = false;
+    
+    hid_t const dataspace = FAILWARN (H5Screate (H5S_SCALAR));
+    hid_t const attribute =
+      FAILWARN (H5Acreate (group, name, H5T_NATIVE_INT, dataspace,
+                           H5P_DEFAULT, H5P_DEFAULT));
+    FAILWARN (H5Awrite (attribute, H5T_NATIVE_DOUBLE, &dvalue));
+    FAILWARN (H5Aclose (attribute));
+    FAILWARN (H5Sclose (dataspace));
+    
+    return error_flag;
+  }
+  
+  
+  
+  // Write a string attribute
+  bool WriteAttribute (hid_t const group,
+                       char const *const name,
+                       char const *const svalue)
+  {
+    bool error_flag = false;
+    
+    hid_t const datatype = FAILWARN (H5Tcopy (H5T_C_S1));
+    FAILWARN (H5Tset_size (datatype, strlen(svalue) + 1));
+    hid_t const dataspace = FAILWARN (H5Screate (H5S_SCALAR));
+    hid_t const attribute =
+      FAILWARN (H5Acreate (group, name, datatype, dataspace,
+                           H5P_DEFAULT, H5P_DEFAULT));
+    FAILWARN (H5Awrite (attribute, datatype, svalue));
+    FAILWARN (H5Aclose (attribute));
+    FAILWARN (H5Sclose (dataspace));
+    FAILWARN (H5Tclose (datatype));
+    
+    return error_flag;
+  }
+  
+  
+  
+  // Write an array of int attributes
+  bool WriteAttribute (hid_t const group,
+                       char const *const name,
+                       int const *const ivalues,
+                       int const nvalues)
+  {
+    bool error_flag = false;
+    
+    hsize_t const size = nvalues;
+    hid_t const dataspace = FAILWARN (H5Screate_simple (1, &size, NULL));
+    hid_t const attribute =
+      FAILWARN (H5Acreate (group, name, H5T_NATIVE_INT,
+                           dataspace, H5P_DEFAULT, H5P_DEFAULT));
+    FAILWARN (H5Awrite (attribute, H5T_NATIVE_INT, ivalues));
+    FAILWARN (H5Aclose (attribute));
+    FAILWARN (H5Sclose (dataspace));
+    
+    return error_flag;
+  }
+  
+  
+  
+  // Write an array of double attributes
+  bool WriteAttribute (hid_t const group,
+                       char const *const name,
+                       double const *const dvalues,
+                       int const nvalues)
+  {
+    bool error_flag = false;
+    
+    hsize_t const size = nvalues;
+    hid_t const dataspace = FAILWARN (H5Screate_simple (1, &size, NULL));
+    hid_t const attribute =
+      FAILWARN (H5Acreate (group, name, H5T_NATIVE_DOUBLE,
+                           dataspace, H5P_DEFAULT, H5P_DEFAULT));
+    FAILWARN (H5Awrite (attribute, H5T_NATIVE_DOUBLE, dvalues));
+    FAILWARN (H5Aclose (attribute));
+    FAILWARN (H5Sclose (dataspace));
+    
+    return error_flag;
+  }
+  
+  
+  
+  // Write an array of string attributes
+  bool WriteAttribute (hid_t const group,
+                       char const *const name,
+                       char const *const *const svalues,
+                       int const nvalues)
+  {
+    bool error_flag = false;
+    
+    size_t maxstrlen = 0;
+    for (int i=0; i<nvalues; ++i) {
+      maxstrlen = max (maxstrlen, strlen(svalues[i]));
+    }
+    vector<char> svalue (nvalues * (maxstrlen+1));
+    for (int i=0; i<nvalues; ++i) {
+      strncpy (&svalue.at(i*maxstrlen), svalues[i], maxstrlen+1);
+    }
+    
+    hid_t const datatype = FAILWARN (H5Tcopy (H5T_C_S1));
+    FAILWARN (H5Tset_size (datatype, maxstrlen + 1));
+    hsize_t const size = nvalues;
+    hid_t const dataspace = FAILWARN (H5Screate_simple (1, & size, NULL));
+    hid_t const attribute =
+      FAILWARN (H5Acreate (group, name, datatype, dataspace,
+                           H5P_DEFAULT, H5P_DEFAULT));
+    FAILWARN (H5Awrite (attribute, datatype, &svalue.front()));
+    FAILWARN (H5Aclose (attribute));
+    FAILWARN (H5Sclose (dataspace));
+    FAILWARN (H5Tclose (datatype));
+    
+    return error_flag;
+  }
+  
+  
+  
+  // Write a large string attribute
+  bool WriteLargeAttribute (hid_t const group,
+                            char const *const name,
+                            char const *const svalue)
+  {
+    bool error_flag = false;
+    
+    // Create a dataset, since the data may not fit into an attribute
+    hsize_t const size = strlen (svalue) + 1;
+    hid_t const dataspace = FAILWARN (H5Screate_simple (1, & size, NULL));
+    hid_t const dataset =
+      FAILWARN (H5Dcreate (group, name, H5T_NATIVE_CHAR,
+                           dataspace, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT));
+    FAILWARN (H5Dwrite (dataset, H5T_NATIVE_CHAR, H5S_ALL, H5S_ALL,
+                        H5P_DEFAULT, svalue));
+    FAILWARN (H5Dclose (dataset));
+    FAILWARN (H5Sclose (dataspace));
+    
+    return error_flag;
+  }
+
+} // namespace CarpetIOF5
diff --git a/CarpetDev/CarpetIOF5/src/input.cc b/CarpetDev/CarpetIOF5/src/input.cc
new file mode 100644
index 000000000..ecf7dda5c
--- /dev/null
+++ b/CarpetDev/CarpetIOF5/src/input.cc
@@ -0,0 +1,563 @@
+#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>
+
+#include "iof5.hh"
+
+
+
+namespace CarpetIOF5 {
+  
+  using namespace std;
+  using namespace Carpet;
+  
+  
+  
+  // 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;
+    
+    double time;
+    char const *gridname;
+    char const *topologyname;
+    int index_depth;            // 0=vertex, 1=cell
+    int topological_dimension;
+    char const *fieldname;
+    char const *fragmentname;
+    
+    // string chartname;
+    
+  public:
+    input_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)
+    {
+      indent_t indent;
+      cout << indent << "time=" << time << "\n";
+      
+#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.)
+      // F5iterate_grid_atlas?
+      for (int m=0; m<maps; ++m) {
+        indent_t indent;
+        chartname = generate_chartname(cctkGH, m);
+        cout << indent << "chart=\"" << chartname << "\"\n";
+        F5iterate_grids
+          (path, NULL, grid_iterator, this, NULL, chartname.c_str());
+      }
+#endif
+      
+      F5iterate_grids (path, NULL, grid_iterator, this, NULL, NULL);
+      
+      // Synchronise
+      BEGIN_REFLEVEL_LOOP(cctkGH) {
+#if 0
+        int groups[] = {
+          CCTK_GroupIndex("GRID::COORDINATES")
+        };
+        int const ngroups = sizeof groups / sizeof *groups;
+        for (int i=0; i<ngroups; ++i) {
+          assert (groups[i]>=0);
+        }
+        iret = CCTK_SyncGroupsI (cctkGH, ngroups, groups);
+        assert (iret >= 0);
+#endif
+#warning "TODO: don't modify boundaries"
+        BEGIN_LOCAL_MAP_LOOP(cctkGH, CCTK_GF) {
+          BEGIN_LOCAL_COMPONENT_LOOP(cctkGH, CCTK_GF) {
+            DECLARE_CCTK_ARGUMENTS;
+            CCTK_LOOP3_BOUNDARIES(boundaries, cctkGH,
+                                  i,j,k, 1,1,1,
+                                  cctk_lsh[0]-1,cctk_lsh[1]-1,cctk_lsh[2]-1)
+            {
+              int const ind3d = CCTK_GFINDEX3D(cctkGH, i,j,k);
+              x[ind3d] = 0.0;
+              y[ind3d] = 0.0;
+              z[ind3d] = 0.0;
+              r[ind3d] = 0.0;
+            } CCTK_ENDLOOP3_BOUNDARIES(boundaries);
+          } END_LOCAL_COMPONENT_LOOP;
+        } END_LOCAL_MAP_LOOP;
+      } END_REFLEVEL_LOOP;
+    }
+  
+    void read_grid (F5Path *const path)
+    {
+      indent_t indent;
+      cout << indent << "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)
+    {
+      indent_t indent;
+      
+      herr_t herr;
+      
+      cout << indent << "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);
+        indent_t indent;
+        cout << indent << "alias for topology \"" << linkval << "\"\n"
+             << indent << "ignoring this topology\n";
+        return;
+      }
+      
+      // F5iterate_topology_fields
+      //   (path, NULL, field_iterator, this, chartname.c_str(), NULL); 
+      F5iterate_topology_fields (path, NULL, field_iterator, this, NULL, NULL); 
+    }
+    
+    void read_field (F5Path *const path)
+    {
+      indent_t indent;
+      cout << indent << "field=\"" << fieldname << "\"\n";
+      
+      if (strcmp(fieldname, "GRID::COORDINATES")!=0 and
+          strcmp(fieldname, "GRID::r")!=0)
+      {
+        indent_t indent;
+        cout << indent << "ignoring this field\n";
+        return;
+      }
+      
+      F5iterate_field_fragments (path, NULL, fragment_iterator, this);
+    }
+    
+    void read_fragment (F5Path *const path)
+    {
+      indent_t indent;
+      
+      int iret;
+      void *pret;
+      
+      cout << indent << "fragment=\"" << fragmentname << "\"\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 from fragment name
+      // TODO: store map number as attribute?
+      int const m = interpret_fragmentname (cctkGH, fragmentname);
+      
+      {
+        indent_t indent;
+        cout << indent << "reading refinement level " << rl << "\n";
+        cout << indent << "reading map " << m << "\n";
+      }
+      
+      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,
+                        int const var)
+    {
+      indent_t indent;
+      
+      herr_t herr;
+      int ierr;
+      int iret;
+      void *pret;
+      
+      assert (path);
+      assert (name);
+      assert (var >= 0);
+      
+      cout << indent << "dataset=\"" << name << "\"\n";
+      
+      assert (var >= 0);
+      {
+        char *const fullname = CCTK_FullName(var);
+        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;
+      if (fragment_is_group) {
+        field = H5Gopen (path->Field_hid, fragmentname, H5P_DEFAULT);
+        assert (field >= 0);
+      } else {
+        field = path->Field_hid;
+      }
+      
+      hid_t const fragment = H5Dopen(field, name, H5P_DEFAULT);
+      assert (fragment);
+      
+      // Check index depth
+      int index_depth_;
+      iret = F5Tget_index_depth(path, &index_depth_);
+      assert (iret);
+      assert (index_depth_ == index_depth);
+      
+      // 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);
+      assert (all(foff>=0));
+      
+#if 0
+      // Read the fragment size.  This is stored with the field -- why
+      // is this different from the offset?
+      hsize_t hlen[FIBER_MAX_RANK];
+      iret =
+        F5LAget_dimensions(path->Field_hid,
+                           FIBER_FIELD_DATASPACE_DIMENSIONS_ATTRIBUTE, hlen);
+      assert (iret == dim);
+      hsize_t hlen2[FIBER_MAX_RANK];
+      pret = F5Tpermute_dimensions(path, dim, hlen2, hlen);
+      assert (pret);
+      ivect const flen = h2v(hlen2);
+      assert (all(flen>=0));
+#endif
+      hid_t const space = H5Dget_space(fragment);
+      assert (space>=0);
+      iret = H5Sget_simple_extent_ndims(space);
+      assert (iret == dim);
+      hsize_t hlen[dim];
+      iret = H5Sget_simple_extent_dims(space, hlen, NULL);
+      hsize_t hlen2[dim];
+      pret = F5Tpermute_dimensions(path, dim, hlen2, hlen);
+      assert (pret);
+      ivect const flen = h2v(hlen2);
+      assert (all(flen>=0));
+      herr = H5Sclose (space);
+      assert (not herr);
+      
+      ibbox const fbox (foff, foff+flen-1, 1);
+      {
+        indent_t indent;
+        cout << indent << "dataset bbox is " << foff << ":" << foff+flen << "\n";
+      }
+      
+      
+      
+      // Examine Cactus variable
+      int const group = CCTK_GroupIndexFromVarI(var);
+      cGroup groupdata;
+      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 == 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);
+#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;
+      
+      
+      
+      herr = H5Dclose(fragment);
+      assert (not herr);
+      
+      if (fragment_is_group) {
+        herr = H5Gclose(field);
+        assert (not herr);
+      }
+    }
+    
+    
+    
+  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)
+    {
+      input_iterator_t* const iterator = (input_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)
+    {
+      input_iterator_t* const iterator = (input_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)
+    {
+      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);
+      return 0;
+    }
+    
+    static
+    herr_t field_iterator (F5Path *const path, char const *const fieldname,
+                           void *const userdata)
+    {
+      input_iterator_t* const iterator = (input_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)
+    {
+      input_iterator_t* const iterator = (input_iterator_t*)userdata;
+      iterator->fragmentname = fragmentname;
+      iterator->read_fragment (path);
+      return 0;
+    }
+    
+  };                            // class input_iterator_t
+  
+  
+  
+  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"));
+    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);
+    }
+  }
+  
+} // end namespace CarpetIOF5
diff --git a/CarpetDev/CarpetIOF5/src/iof5.cc.old b/CarpetDev/CarpetIOF5/src/iof5.cc.old
new file mode 100644
index 000000000..fe780030b
--- /dev/null
+++ b/CarpetDev/CarpetIOF5/src/iof5.cc.old
@@ -0,0 +1,643 @@
+#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
new file mode 100644
index 000000000..71524fa2c
--- /dev/null
+++ b/CarpetDev/CarpetIOF5/src/iof5.hh
@@ -0,0 +1,214 @@
+#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 <defs.hh>
+#include <vect.hh>
+
+#include <carpet.hh>
+
+
+
+// This requires defining a boolean local variable "error_flag",
+// initialised to false
+#define FAILWARN(_expr)                                                 \
+  failwarn(error_flag, _expr, __LINE__, __FILE__, CCTK_THORNSTRING, #_expr)
+
+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)
+{
+  if (expr < 0) {
+    CCTK_VWarn (CCTK_WARN_ALERT, line, file, thorn,
+                "Expression \"%s\" return %d", msg, (int)expr);
+    error_flag = true;
+  }
+  return expr;
+}
+
+  
+
+namespace CarpetIOF5 {
+  class indent_t;
+};
+ostream& operator<<(ostream& os, CarpetIOF5::indent_t const& indent);
+
+
+
+namespace CarpetIOF5 {
+  
+  using namespace std;
+  using namespace Carpet;
+  
+  
+  
+  // Indentation
+  class indent_t {
+    static int const width = 3;
+    static int level;
+  public:
+    indent_t() { ++level; }
+    ~indent_t() { --level; }
+    ostream& output(ostream& os) const;
+  };
+  
+  
+  
+  // File mode for creating directories
+  int const mode = 0755;
+  
+  // A nan
+  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;
+  
+  
+  
+  // Data types for HDF5 and F5 types
+  typedef vect<hsize_t,dim> hvect;
+  
+  // Conversion operators for these datatypes
+  static inline
+  hvect v2h (ivect const& a)
+  {
+    return hvect(a);
+  }
+  static inline
+  ivect h2v (hvect const& a)
+  {
+    return ivect(a);
+  }
+  static inline
+  ivect h2v (hsize_t const* const a)
+  {
+    ivect res;
+    for (int d=0; d<dim; ++d) {
+      res[d] = a[d];
+    }
+    return res;
+  }
+  
+  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;
+  }
+  
+  
+  
+  // Handle HDF5 attributes more comfortably
+  bool WriteAttribute (hid_t const group,
+                       char const * const name,
+                       int const ivalue);
+  bool WriteAttribute (hid_t const group,
+                       char const * const name,
+                       double const dvalue);
+  bool WriteAttribute (hid_t const group,
+                       char const *const name,
+                       char const *const svalue);
+  bool WriteAttribute (hid_t const group,
+                       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,
+                       int const nvalues);
+  bool WriteAttribute (hid_t const group,
+                       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);
+  
+  
+  
+  // Generate a good file name ("alias") for a variable
+  string
+  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,
+                   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 a good tensor basis name (coordinate system name)
+  string
+  generate_chartname (cGH const *const cctkGH);
+  
+  // Generate a good fragment name (map name)
+  string
+  generate_fragmentname (cGH const *const cctkGH, int const m);
+  int
+  interpret_fragmentname (cGH const *const cctkGH,
+                          char const *const fragmentname);
+  
+  // Generate a good topology name (map and refinement level name)
+  string
+  generate_topologyname (cGH const *const cctkGH,
+                         int const m, ivect const& reffact);
+  
+  
+  
+  void
+  write_metadata (cGH *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);
+  
+} // end namespace CarpetIOF5
diff --git a/CarpetDev/CarpetIOF5/src/make.code.defn b/CarpetDev/CarpetIOF5/src/make.code.defn
new file mode 100644
index 000000000..da381b90b
--- /dev/null
+++ b/CarpetDev/CarpetIOF5/src/make.code.defn
@@ -0,0 +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
+
+# Subdirectories containing source files
+SUBDIRS = 
diff --git a/CarpetDev/CarpetIOF5/src/output.cc b/CarpetDev/CarpetIOF5/src/output.cc
new file mode 100644
index 000000000..0966e5dfb
--- /dev/null
+++ b/CarpetDev/CarpetIOF5/src/output.cc
@@ -0,0 +1,573 @@
+#include <cctk.h>
+#include <cctk_Arguments.h>
+#include <cctk_Parameters.h>
+#include <cctk_Version.h>
+
+#include <cassert>
+#include <cstdlib>
+#include <cstring>
+#include <iomanip>
+#include <iostream>
+#include <limits>
+#include <sstream>
+#include <string>
+#include <vector>
+
+#include <hdf5.h>
+#include <F5/F5F.h>
+#include <F5/F5R.h>
+#include <F5/F5X.h>
+#include <F5/F5iterate.h>
+#include <F5/F5uniform.h>
+
+#include "CactusBase/IOUtil/src/ioutil_CheckpointRecovery.h"
+
+#include <bbox.hh>
+#include <vect.hh>
+
+#include <carpet.hh>
+
+#include "iof5.hh"
+
+
+
+namespace CarpetIOF5 {
+  
+  using namespace std;
+  using namespace Carpet;
+  
+  
+  
+  class output_iterator_t {
+    cGH *const cctkGH;
+    
+    bool const is_multipatch;
+    
+    string gridname;
+    string chartname;
+    string fragmentname;
+    string topologyname;
+    
+  public:
+    output_iterator_t (cGH *const cctkGH_)
+      : cctkGH(cctkGH_),
+        is_multipatch
+        (CCTK_IsFunctionAliased("MultiPatch_GetSystemSpecification"))
+    {
+    }
+    
+    void iterate (hid_t const file)
+    {
+      int const rl = 0;
+      ENTER_LEVEL_MODE(cctkGH, rl) {
+        output_simulation (file);
+      } LEAVE_LEVEL_MODE;
+    }
+    
+    
+    
+  private:
+    
+    void output_simulation (hid_t const file)
+    {
+      assert (is_level_mode() and reflevel==0);
+      DECLARE_CCTK_ARGUMENTS;
+      indent_t indent;
+      
+      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);
+      
+      // 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;
+      
+      herr_t herr;
+      
+      assert (is_singlemap_mode() and reflevel==0);
+      
+      cout << indent << "map=" << Carpet::map << "\n";
+      
+      fragmentname = generate_fragmentname(cctkGH, Carpet::map);
+      
+      // 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) {
+              output_reflevel (file);
+            } LEAVE_SINGLEMAP_MODE;
+          }
+        } END_REFLEVEL_LOOP;
+      } END_GLOBAL_MODE;
+    }
+    
+    
+    
+    void output_reflevel (hid_t const file)
+    {
+      DECLARE_CCTK_ARGUMENTS;
+      indent_t indent;
+      
+      cout << indent << "reflevel=" << reflevel << "\n";
+      
+      ivect const reffact = spacereffacts.AT(reflevel);
+      topologyname = generate_topologyname(cctkGH, Carpet::map, reffact);
+      
+      // Define grid hierarchy
+      int const indexdepth = 0; // vertices
+      F5Path *const path =
+        F5Rcreate_coordinate_topology (file, &cctk_time,
+                                       gridname.c_str(), chartname.c_str(),
+                                       topologyname.c_str(),
+                                       indexdepth,
+                                       dim, dim, &v2h(reffact)[0]);
+      
+      // 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];
+      }
+      
+      if (not is_multipatch) {
+        // 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) {
+        output_component (path);
+      } END_LOCAL_COMPONENT_LOOP;
+      
+      // Close topology
+      F5close (path);
+    }
+  
+  
+  
+    void output_component (F5Path *const path)
+    {
+      DECLARE_CCTK_ARGUMENTS;
+      indent_t indent;
+      
+      cout << indent << "component=" << 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];
+      }
+      
+      // Define coordinates
+      if (not is_multipatch) {
+        // (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());
+      } else {
+        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"));
+    }
+    
+    
+    
+    void output_variable (F5Path *const path, int const var,
+                          bool const write_positions = false)
+    {
+      DECLARE_CCTK_ARGUMENTS;
+      DECLARE_CCTK_PARAMETERS;
+      indent_t indent;
+      
+      int ierr;
+      
+      assert (var >= 0);
+      {
+        char *const fullname = CCTK_FullName(var);
+        cout << indent << "variable=" << fullname << "\n";
+        free (fullname);
+      }
+      
+      int const group = CCTK_GroupIndexFromVarI(var);
+      int const v0 = CCTK_FirstVarIndexI(group);
+      
+      cGroup groupdata;
+      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;
+      
+      // 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;
+        }
+        ioff[d] = lbnd[d] + imin[d];
+        ilen[d] = imax[d] - imin[d];
+      }
+      
+#warning "TODO: Do not output symmetry zones (unless requested by the user)"
+#warning "TODO: Do not output buffer zones (is that easily possible?)"
+      
+      
+      
+      enum tensortype_t {
+        tt_scalar, tt_vector, tt_error
+      };
+      
+      tensortype_t tensortype = tt_error;
+      
+      int const coordinates_group = CCTK_GroupIndex ("grid::coordinates");
+      assert (coordinates_group >= 0);
+      
+      if (group == coordinates_group) {
+        
+        // Special case
+        if (var >= v0 and var < v0+dim) {
+          tensortype = tt_vector;
+        } else if (var == v0+dim) {
+          tensortype = tt_scalar;
+        } else {
+          assert (0);
+        }
+        
+      } else {
+        
+        // TODO: use the tensor type instead
+        switch (groupdata.numvars) {
+        case 1:
+          tensortype = tt_scalar; break;
+        case 3:
+          tensortype = tt_vector; break;
+        default:
+          // fallback: output group as scalars
+          tensortype = tt_scalar; break;
+        }
+        
+      }
+      
+      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);
+        return;
+      }
+      
+      
+      
+      switch (tensortype) {
+        
+      case tt_scalar: {
+        // Scalar field
+        CCTK_REAL const *const rdata =
+          (CCTK_REAL const*)CCTK_VarDataPtrI (cctkGH, timelevel, var);
+        assert (rdata);
+        
+        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];
+            }
+          }
+        }
+        void const *const data = &idata.front();
+        
+        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
+        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);
+        break;
+      }
+        
+      case tt_vector: {
+        // Vector field
+        CCTK_REAL const* rdata[cctk_dim];
+        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();
+          
+        }
+        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);
+        break;
+      }
+        
+      default:
+        assert (0);
+      }
+    }
+    
+  };                            // class output_iterator_t
+  
+  
+  
+  void write_metadata (cGH *const cctkGH, hid_t const file)
+  {
+    DECLARE_CCTK_PARAMETERS;
+    
+    assert (cctkGH);
+    assert (file >= 0);
+    
+    herr_t herr;
+    
+    
+    CCTK_INFO ("Writing simulation metadata...");
+    
+    // NOTE: We could write the metadata into only one of the process
+    // files (to save space), or write it into a separate metadata
+    // file
+    
+    // Create a group to hold all metadata
+    hid_t const group =
+      H5Gcreate (file, metadata_group, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+    assert (group >= 0);
+    
+    // General information
+    WriteAttribute (group, "Cactus version", CCTK_FullVersion());
+    
+    // Unique identifiers
+    if (CCTK_IsFunctionAliased ("UniqueConfigID")) {
+      WriteAttribute
+        (group, "config id", (char const*) UniqueConfigID(cctkGH));
+    }
+    if (CCTK_IsFunctionAliased ("UniqueBuildID")) {
+      WriteAttribute
+        (group, "build id", (char const*) UniqueBuildID(cctkGH));
+    }
+    if (CCTK_IsFunctionAliased ("UniqueSimulationID")) {
+      WriteAttribute
+        (group, "simulation id", (char const*) UniqueSimulationID(cctkGH));
+    }
+    if (CCTK_IsFunctionAliased ("UniqueRunID")) {
+      WriteAttribute
+        (group, "run id", (char const*) UniqueRunID(cctkGH));
+    }
+    
+    // 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);
+    
+    // Write parameters into a separate dataset (they may be too large
+    // for an attribute)
+    int const get_all = 1;
+    char *const parameters = IOUtil_GetAllParameters (cctkGH, get_all);
+    assert (parameters);
+    WriteLargeAttribute (group, all_parameters, parameters);
+    free (parameters);
+    
+    // Grid structure
+    string const gs = serialise_grid_structure (cctkGH);
+    WriteLargeAttribute (group, grid_structure, gs.c_str());
+    
+    herr = H5Gclose (group);
+    assert (not herr);
+  }
+  
+  
+  
+  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);
+    
+    
+    
+    // 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);
+    iterator.iterate (file);
+    
+    // Close file
+    herr = H5Fclose (file);
+    assert (not herr);
+  }
+  
+} // end namespace CarpetIOF5
diff --git a/CarpetDev/CarpetIOF5/src/poison.cc b/CarpetDev/CarpetIOF5/src/poison.cc
new file mode 100644
index 000000000..6029b55a8
--- /dev/null
+++ b/CarpetDev/CarpetIOF5/src/poison.cc
@@ -0,0 +1,77 @@
+#include <cctk.h>
+#include <cctk_Arguments.h>
+#include <cctk_Parameters.h>
+
+#include <carpet.hh>
+#include <loopcontrol.h>
+
+#include "iof5.hh"
+
+
+
+namespace CarpetIOF5 {
+  
+  using namespace std;
+  using namespace Carpet;
+  
+  
+  
+  extern "C"
+  void F5_Poison (CCTK_ARGUMENTS)
+  {
+    DECLARE_CCTK_PARAMETERS;
+    
+    assert (is_global_mode());
+    CCTK_VInfo (CCTK_THORNSTRING,
+                "F5_Poison: iteration=%d", cctkGH->cctk_iteration);
+    
+    BEGIN_REFLEVEL_LOOP(cctkGH) {
+      BEGIN_LOCAL_MAP_LOOP(cctkGH, CCTK_GF) {
+        BEGIN_LOCAL_COMPONENT_LOOP(cctkGH, CCTK_GF) {
+          DECLARE_CCTK_ARGUMENTS;
+          LC_LOOP3(F5_Poison, i,j,k,
+                   0,0,0, cctk_lsh[0],cctk_lsh[1],cctk_lsh[2],
+                   cctk_lsh[0],cctk_lsh[1],cctk_lsh[2])
+          {
+            int const ind3d = CCTK_GFINDEX3D(cctkGH, i,j,k);
+            x[ind3d] = nan;
+            y[ind3d] = nan;
+            z[ind3d] = nan;
+            r[ind3d] = nan;
+          } LC_ENDLOOP3(F5_Poison);
+        } END_LOCAL_COMPONENT_LOOP;
+      } END_LOCAL_MAP_LOOP;
+    } END_REFLEVEL_LOOP;
+  }
+  
+  
+  
+  extern "C"
+  void F5_Check (CCTK_ARGUMENTS)
+  {
+    DECLARE_CCTK_PARAMETERS;
+    
+    assert (is_global_mode());
+    CCTK_VInfo (CCTK_THORNSTRING,
+                "F5_Check: iteration=%d", cctkGH->cctk_iteration);
+    
+    BEGIN_REFLEVEL_LOOP(cctkGH) {
+      BEGIN_LOCAL_MAP_LOOP(cctkGH, CCTK_GF) {
+        BEGIN_LOCAL_COMPONENT_LOOP(cctkGH, CCTK_GF) {
+          DECLARE_CCTK_ARGUMENTS;
+          LC_LOOP3(F5_Check, i,j,k,
+                   0,0,0, cctk_lsh[0],cctk_lsh[1],cctk_lsh[2],
+                   cctk_lsh[0],cctk_lsh[1],cctk_lsh[2])
+          {
+            int const ind3d = CCTK_GFINDEX3D(cctkGH, i,j,k);
+            assert (not isnan(x[ind3d]));
+            assert (not isnan(y[ind3d]));
+            assert (not isnan(z[ind3d]));
+            assert (not isnan(r[ind3d]));
+          } LC_ENDLOOP3(F5_Check);
+        } END_LOCAL_COMPONENT_LOOP;
+      } END_LOCAL_MAP_LOOP;
+    } END_REFLEVEL_LOOP;
+  }
+  
+} // end namespace CarpetIOF5
diff --git a/CarpetDev/CarpetIOF5/src/util.cc b/CarpetDev/CarpetIOF5/src/util.cc
new file mode 100644
index 000000000..57ccfd6db
--- /dev/null
+++ b/CarpetDev/CarpetIOF5/src/util.cc
@@ -0,0 +1,282 @@
+#include <cctk.h>
+#include <cctk_Arguments.h>
+#include <cctk_Parameters.h>
+
+#include <cassert>
+#include <cstdio>
+#include <cstdlib>
+#include <cstring>
+#include <iomanip>
+#include <iostream>
+#include <limits>
+#include <sstream>
+#include <string>
+#include <vector>
+
+#include <hdf5.h>
+#include <F5/F5F.h>
+#include <F5/F5R.h>
+#include <F5/F5iterate.h>
+#include <F5/F5uniform.h>
+
+#include <bbox.hh>
+#include <defs.hh>
+#include <vect.hh>
+
+#include <carpet.hh>
+
+#include "iof5.hh"
+
+
+
+ostream& operator<<(ostream& os, CarpetIOF5::indent_t const& indent)
+{ return indent.output(os); }
+
+
+
+namespace CarpetIOF5 {
+  
+  using namespace std;
+  using namespace Carpet;
+  
+  
+  
+  // Indentation
+  int indent_t::level = 0;
+  
+  ostream& indent_t::output(ostream& os) const
+  {
+    return os << string(width*level, ' ');
+  }
+  
+  
+  
+  // Generate a good file name ("alias") for a variable
+  string
+  generate_basename (cGH const *const cctkGH,
+                     int const variable)
+  {
+    DECLARE_CCTK_PARAMETERS;
+    
+    ostringstream filename_buf;
+    
+    if (CCTK_EQUALS (file_content, "variable")) {
+      assert (variable >= 0);
+      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")) {
+      assert (variable >= 0);
+      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")) {
+      assert (variable >= 0);
+      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,
+                   int const proc,
+                   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;
+    
+    if (create_subdirs) {
+      {
+        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);
+        }
+      }
+      {
+        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) {
+      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();
+  }
+  
+  
+  
+  // 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 *const suffix = ".par";
+    if (strlen(parfilename) >= strlen(suffix)) {
+      char *const p = parfilename + strlen(parfilename) - strlen(suffix);
+      if (strcmp(p, suffix) == 0) {
+        *p = '\0';
+      }
+    }
+    char *const s = strrchr(parfilename, '/');
+    char *const p = s ? s+1 : parfilename;
+    
+    return p;
+  }
+  
+  // Generate a good tensor basis name (chart name)
+  string
+  generate_chartname (cGH const *const cctkGH)
+  {
+    // Assume a global tensor basis, where all maps use the same chart
+    return FIBER_HDF5_DEFAULT_CHART;
+    
+#if 0
+    // TODO: Add some interesting information here, e.g. the name of
+    // the multi-patch system, or a good name of the patch?
+    ostringstream buf;
+    buf << "Map_" << m;
+    return buf.str();
+#endif
+  }
+  
+  // Generate a good fragment name (map name)
+  string
+  generate_fragmentname (cGH const *const cctkGH, int const m)
+  {
+    ostringstream buf;
+    buf << "Map_" << m;
+    return buf.str();
+  }
+  
+  int
+  interpret_fragmentname (cGH const *const cctkGH,
+                          char const *const fragmentname)
+  {
+    int m = -1;
+    sscanf (fragmentname, "Map_%d", &m);
+    assert (m>=0 and m<Carpet::maps);
+    return m;
+  }
+  
+  // Generate a good topology name (map and refinement level name)
+  string
+  generate_topologyname (cGH const *const cctkGH,
+                         int const m, ivect const& reffact)
+  {
+    ostringstream buf;
+    // buf << "Map_" << m << "_VertexLevel_";
+    buf << "VertexLevel_";
+    for (int d=0; d<dim; ++d) {
+      if (d>0) buf << "x";
+      buf << reffact[d];
+    }
+    return buf.str();
+  }
+  
+  
+  
+  char const *const grid_structure = "Grid Structure v5";
+  
+  string
+  serialise_grid_structure (cGH const *const cctkGH)
+  {
+    ostringstream buf;
+    buf << setprecision(17);
+    buf << grid_structure << "{";
+    buf << "maps=" << maps << ",";
+    buf << "[";
+    for (int m=0; m<maps; ++m) {
+      buf << "[" << m << "]={";
+      buf << "superregions:" << vhh.AT(m)->superregions << ",";
+      // We could omit the regions
+      buf << "regions:" << vhh.AT(m)->regions.AT(0) << ",";
+      buf << "ghost_widths:" << vdd.AT(m)->ghost_widths << ",";
+      buf << "buffer_widths:" << vdd.AT(m)->buffer_widths << ",";
+      buf << "prolongation_orders_space:"
+          << vdd.AT(m)->prolongation_orders_space << ",";
+      buf << "},";
+    }
+    buf << "],";
+    buf << "times:" << *tt << ",";
+    buf << "}.";
+    return buf.str();
+  }
+  
+} // end namespace CarpetIOF5