Initial revision

darcs-hash:20010301114010-f6438-12fb8a9ffcc80e86c0a97e37b5b0dae0dbc59b79.gz

10 files changed, 1921 insertions, 0 deletions

diff --git a/Carpet/CarpetIOASCII/src/ioascii.cc b/Carpet/CarpetIOASCII/src/ioascii.cc
new file mode 100644
index 000000000..eb28e27d0
--- /dev/null
+++ b/Carpet/CarpetIOASCII/src/ioascii.cc
@@ -0,0 +1,1181 @@
+#include <assert.h>
+#include <limits.h>
+#include <math.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <sys/stat.h>
+#include <sys/types.h>
+
+#include <fstream>
+#include <iomanip>
+#include <ostream>
+#include <sstream>
+#include <string>
+#include <vector>
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+
+#include "CactusBase/IOUtil/src/ioGH.h"
+
+#include "data.hh"
+#include "dist.hh"
+#include "gdata.hh"
+#include "gf.hh"
+#include "ggf.hh"
+#include "vect.hh"
+
+#include "carpet.hh"
+
+#include "ioascii.hh"
+
+extern "C" {
+  static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetIOASCII/src/ioascii.cc,v 1.77 2004/06/23 17:36:41 tradke Exp $";
+  CCTK_FILEVERSION(Carpet_CarpetIOASCII_ioascii_cc);
+}
+
+
+
+// That's a hack
+namespace Carpet {
+  void UnsupportedVarType (const int vindex);
+}
+
+
+
+namespace CarpetIOASCII {
+
+  using namespace std;
+  using namespace Carpet;
+
+  void SetFlag (int index, const char* optstring, void* arg);
+
+
+
+  // Special output routines for complex numbers
+
+#ifdef CCTK_REAL4
+  ostream& operator<< (ostream& os, const CCTK_COMPLEX8& val)
+  {
+    return os << CCTK_Cmplx8Real(val) << " " << CCTK_Cmplx8Imag(val);
+  }
+#endif
+
+#ifdef CCTK_REAL8
+  ostream& operator<< (ostream& os, const CCTK_COMPLEX16& val)
+  {
+    return os << CCTK_Cmplx16Real(val) << " " << CCTK_Cmplx16Imag(val);
+  }
+#endif
+
+#ifdef CCTK_REAL16
+  ostream& operator<< (ostream& os, const CCTK_COMPLEX32& val)
+  {
+    return os << CCTK_Cmplx32Real(val) << " " << CCTK_Cmplx32Imag(val);
+  }
+#endif
+
+
+
+  template<int D,int DD>
+  void WriteASCII (ostream& os,
+		   const gdata<D>* const gfdata,
+		   const bbox<int,D>& gfext,
+		   const int vi,
+		   const int time,
+		   const vect<int,D>& org,
+		   const vect<int,DD>& dirs,
+		   const int rl,
+		   const int ml,
+		   const int m,
+		   const int c,
+		   const int tl,
+		   const CCTK_REAL coord_time,
+		   const vect<CCTK_REAL,D>& coord_lower,
+		   const vect<CCTK_REAL,D>& coord_upper);
+
+
+
+  int CarpetIOASCIIStartup()
+  {
+    IOASCII<0>::Startup();
+    IOASCII<1>::Startup();
+    IOASCII<2>::Startup();
+    IOASCII<3>::Startup();
+    return 0;
+  }
+
+
+
+  void CarpetIOASCIIInit (CCTK_ARGUMENTS)
+  {
+    DECLARE_CCTK_ARGUMENTS;
+
+    for (int d=0; d<4; ++d) {
+      this_iteration[d] = 0;
+      last_output_iteration[d] = 0;
+      last_output_time[d] = cctk_time;
+    }
+  }
+
+
+
+  // Definition of static members
+  template<int outdim> int IOASCII<outdim>::GHExtension;
+  template<int outdim> int IOASCII<outdim>::IOMethod;
+  template<int outdim> vector<bool> IOASCII<outdim>::do_truncate;
+  template<int outdim>
+  vector<vector<vector<int> > > IOASCII<outdim>::last_output;
+
+
+
+  template<int outdim>
+  int IOASCII<outdim>::Startup()
+  {
+    ostringstream msg;
+    msg << "AMR " << outdim << "D ASCII I/O provided by CarpetIOASCII";
+    CCTK_RegisterBanner (msg.str().c_str());
+
+    ostringstream extension_name;
+    extension_name << "CarpetIOASCII_" << outdim << "D";
+    GHExtension = CCTK_RegisterGHExtension(extension_name.str().c_str());
+    CCTK_RegisterGHExtensionSetupGH (GHExtension, SetupGH);
+
+    ostringstream method_name;
+    method_name << "IOASCII_" << outdim << "D";
+    IOMethod = CCTK_RegisterIOMethod (method_name.str().c_str());
+    CCTK_RegisterIOMethodOutputGH (IOMethod, OutputGH);
+    CCTK_RegisterIOMethodOutputVarAs (IOMethod, OutputVarAs);
+    CCTK_RegisterIOMethodTimeToOutput (IOMethod, TimeToOutput);
+    CCTK_RegisterIOMethodTriggerOutput (IOMethod, TriggerOutput);
+
+    return 0;
+  }
+
+
+
+  template<int outdim>
+  void* IOASCII<outdim>
+  ::SetupGH (tFleshConfig* const fc, const int convLevel, cGH* const cgh)
+  {
+    DECLARE_CCTK_PARAMETERS;
+    const void *dummy;
+
+    dummy = &fc;
+    dummy = &convLevel;
+    dummy = &cgh;
+    dummy = &dummy;
+
+    // Truncate all files if this is not a restart
+    do_truncate.resize(CCTK_NumVars(), true);
+
+    // No iterations have yet been output
+    last_output.resize(mglevels);
+    for (int ml=0; ml<mglevels; ++ml) {
+      last_output.at(ml).resize(maxreflevels);
+      for (int rl=0; rl<maxreflevels; ++rl) {
+        last_output.at(ml).at(rl).resize(CCTK_NumVars(), -1);
+      }
+    }
+
+    // We register only once, ergo we get only one handle.  We store
+    // that statically, so there is no need to pass anything to
+    // Cactus.
+    return 0;
+  }
+
+
+
+  template<int outdim>
+  int IOASCII<outdim>
+  ::OutputGH (const cGH* const cgh)
+  {
+    for (int vindex=0; vindex<CCTK_NumVars(); ++vindex) {
+      if (TimeToOutput(cgh, vindex)) {
+	TriggerOutput(cgh, vindex);
+      }
+    }
+    return 0;
+  }
+
+
+
+  template<int outdim>
+  int IOASCII<outdim>
+  ::OutputVarAs (const cGH* const cgh,
+		 const char* const varname, const char* const alias)
+  {
+    DECLARE_CCTK_PARAMETERS;
+
+    assert (is_level_mode());
+
+    const int n = CCTK_VarIndex(varname);
+    if (n<0) {
+      CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+		  "Variable \"%s\" does not exist", varname);
+      return -1;
+    }
+    assert (n>=0 && n<CCTK_NumVars());
+    const int group = CCTK_GroupIndexFromVarI (n);
+    assert (group>=0 && group<(int)Carpet::arrdata.size());
+    const int n0 = CCTK_FirstVarIndexI(group);
+    assert (n0>=0 && n0<CCTK_NumVars());
+    const int var = n - n0;
+    assert (var>=0 && var<CCTK_NumVarsInGroupI(group));
+    const int num_tl = CCTK_NumTimeLevelsFromVarI(n);
+    assert (num_tl>=1);
+
+    // Check for storage
+    if (! CCTK_QueryGroupStorageI(cgh, group)) {
+      CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+		  "Cannot output variable \"%s\" because it has no storage",
+		  varname);
+      return 0;
+    }
+
+    const int grouptype = CCTK_GroupTypeI(group);
+    switch (grouptype) {
+    case CCTK_SCALAR:
+    case CCTK_ARRAY:
+      assert (do_global_mode);
+      break;
+    case CCTK_GF:
+      /* do nothing */
+      break;
+    default:
+      assert (0);
+    }
+    const int rl = grouptype == CCTK_GF ? reflevel : 0;
+
+    const int groupdim = CCTK_GroupDimI(group);
+    if (outdim > groupdim) {
+      CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                  "Cannot produce %dD ASCII output file \"%s\" for variable \"%s\" because it has only %d dimensions", outdim, alias, varname, groupdim);
+      return -1;
+    }
+    assert (outdim <= groupdim);
+
+    // Get grid hierarchy extentsion from IOUtil
+    const ioGH * const iogh = (const ioGH *)CCTK_GHExtension (cgh, "IO");
+    assert (iogh);
+
+    // Create the output directory
+    const char* myoutdir = GetStringParameter("out%dD_dir");
+    if (CCTK_EQUALS(myoutdir, "")) {
+      myoutdir = out_dir;
+    }
+    if (CCTK_MyProc(cgh)==0) {
+      CCTK_CreateDirectory (0755, myoutdir);
+    }
+
+    // Loop over all direction combinations
+    vect<int,outdim> dirs (0);
+
+    bool done;
+    do {
+
+      // Output each combination only once
+      bool ascending = true;
+      for (int d1=0; d1<outdim; ++d1) {
+	for (int d2=d1+1; d2<outdim; ++d2) {
+	  ascending = ascending && dirs[d1] < dirs[d2];
+	}
+      }
+
+      // Skip output if the dimensions are not ascending
+      if (ascending) {
+
+	// If this output is desired at all
+	bool desired;
+	switch (outdim) {
+        case 0:
+	  // Output is always desired (if switched on)
+          desired = true;
+          break;
+	case 1:
+	  switch (dirs[0]) {
+	  case 0:
+	    desired = out1D_x;
+	    break;
+	  case 1:
+	    desired = out1D_y;
+	    break;
+	  case 2:
+	    desired = out1D_z;
+	    break;
+	  default:
+	    assert (0);
+	  }
+	  break;
+	case 2:
+	  if (dirs[0]==0 && dirs[1]==1) {
+	    desired = out2D_xy;
+	  } else if (dirs[0]==0 && dirs[1]==2) {
+	    desired = out2D_xz;
+	  } else if (dirs[0]==1 && dirs[1]==2) {
+	    desired = out2D_yz;
+	  } else {
+	    assert (0);
+	  }
+	  break;
+	case 3:
+	  // Output is always desired (if switched on)
+	  desired = true;
+	  break;
+	default:
+	  assert (0);
+	}
+
+	// Skip output if not desired
+	if (desired) {
+
+	  // Traverse all maps on this refinement and multigrid level
+	  BEGIN_MAP_LOOP(cgh, grouptype) {
+
+            // Find the output offset
+            ivect offset(0);
+            if (grouptype == CCTK_GF) {
+              switch (outdim) {
+              case 0:
+                offset[0] = GetGridOffset
+                  (cgh, 1,
+                   "out%dD_point_xi", /*"out_point_xi"*/ NULL,
+                   "out%dD_point_x",  /*"out_point_x"*/  NULL,
+                   /*out_point_x*/ 0.0);
+                offset[1] = GetGridOffset
+                  (cgh, 2,
+                   "out%dD_point_yi", /*"out_point_yi"*/ NULL,
+                   "out%dD_point_y",  /*"out_point_y"*/  NULL,
+                   /*out_point_y*/ 0.0);
+                offset[2] = GetGridOffset
+                  (cgh, 3,
+                   "out%dD_point_zi", /*"out_point_zi"*/ NULL,
+                   "out%dD_point_z",  /*"out_point_z"*/  NULL,
+                   /*out_point_z*/ 0.0);
+                break;
+              case 1:
+                switch (dirs[0]) {
+                case 0:
+                  offset[1] = GetGridOffset (cgh, 2,
+                                             "out%dD_xline_yi", "out_xline_yi",
+                                             "out%dD_xline_y",  "out_xline_y",
+                                             out_xline_y);
+                  offset[2] = GetGridOffset (cgh, 3,
+                                             "out%dD_xline_zi", "out_xline_zi",
+                                             "out%dD_xline_z",  "out_xline_z",
+                                             out_xline_z);
+                  break;
+                case 1:
+                  offset[0] = GetGridOffset (cgh, 1,
+                                             "out%dD_yline_xi", "out_yline_xi",
+                                             "out%dD_yline_x",  "out_yline_x",
+                                             out_yline_x);
+                  offset[2] = GetGridOffset (cgh, 3,
+                                             "out%dD_yline_zi", "out_yline_zi",
+                                             "out%dD_yline_z",  "out_yline_z",
+                                             out_yline_z);
+                  break;
+                case 2:
+                  offset[0] = GetGridOffset (cgh, 1,
+                                             "out%dD_zline_xi", "out_zline_xi",
+                                             "out%dD_zline_x",  "out_zline_x",
+                                             out_zline_x);
+                  offset[1] = GetGridOffset (cgh, 2,
+                                             "out%dD_zline_yi", "out_zline_yi",
+                                             "out%dD_zline_y",  "out_zline_y",
+                                             out_zline_y);
+                  break;
+                default:
+                  assert (0);
+                }
+                break;
+              case 2:
+                if (dirs[0]==0 && dirs[1]==1) {
+                  offset[2] = GetGridOffset
+                    (cgh, 3,
+                     "out%dD_xyplane_zi", "out_xyplane_zi",
+                     "out%dD_xyplane_z",  "out_xyplane_z",
+                     out_xyplane_z);
+                } else if (dirs[0]==0 && dirs[1]==2) {
+                  offset[1] = GetGridOffset
+                    (cgh, 2,
+                     "out%dD_xzplane_yi", "out_xzplane_yi",
+                     "out%dD_xzplane_y",  "out_xzplane_y",
+                     out_xzplane_y);
+                } else if (dirs[0]==1 && dirs[1]==2) {
+                  offset[0] = GetGridOffset
+                    (cgh, 1,
+                     "out%dD_yzplane_xi", "out_yzplane_xi",
+                     "out%dD_yzplane_x",  "out_yzplane_x",
+                     out_yzplane_x);
+                } else {
+                  assert (0);
+                }
+                break;
+              case 3:
+                // The offset doesn't matter in this case
+                break;
+              default:
+                assert (0);
+              }
+            } // if grouptype is GF
+
+            ofstream file;
+            if (CCTK_MyProc(cgh)==0) {
+
+              // Invent a file name
+              ostringstream filenamebuf;
+              if (new_filename_scheme) {
+                filenamebuf << myoutdir << "/" << alias << ".";
+                if (maps > 1) {
+                  filenamebuf << Carpet::map << ".";
+                }
+                for (int d=0; d<outdim; ++d) {
+                  const char* const coords = "xyz";
+                  filenamebuf << coords[dirs[d]];
+                }
+// The offsets differ per level
+//                 for (int dd=0; dd<groupdim; ++dd) {
+//                   bool print_dir = true;
+//                   for (int d=0; d<outdim; ++d) {
+//                     print_dir = print_dir && dirs[d] != dd;
+//                   }
+//                   if (print_dir) {
+//                     filenamebuf << "." << offset[dd];
+//                   }
+//                 }
+                filenamebuf << ".asc";
+              } else {
+                filenamebuf << myoutdir << "/" << alias << ".";
+                if (maps > 1) {
+                  filenamebuf << Carpet::map << ".";
+                }
+                for (int d=0; d<outdim; ++d) {
+                  assert (dirs[d]>=0 && dirs[d]<3);
+                  const char* const coords = "xyz";
+                  filenamebuf << coords[dirs[d]];
+                }
+                const char* const suffixes = "plpv";
+                filenamebuf << suffixes[outdim];
+              }
+              // we need a persistent temporary here
+              string filenamestr = filenamebuf.str();
+              const char* const filename = filenamestr.c_str();
+
+              // If this is the first time, then write a nice header
+              if (do_truncate.at(n)) {
+                struct stat fileinfo;
+                if (! iogh->recovered
+                    || stat(filename, &fileinfo)!=0) {
+                  file.open (filename, ios::out | ios::trunc);
+                  if (! file.good()) {
+                    CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING,
+                                "Could not open output file \"%s\" for variable \"%s\"",
+                                filename, varname);
+                  }
+                  assert (file.good());
+                  file << "# " << varname;
+                  for (int d=0; d<outdim; ++d) {
+                    file << " " << "xyz"[dirs[d]];
+                  }
+                  file << " (" << alias << ")" << endl;
+                  file << "#" << endl;
+                  assert (file.good());
+                }
+              }
+
+              // Open the file
+              if (! file.is_open()) {
+                file.open (filename, ios::out | ios::app);
+                assert (file.good());
+              }
+
+              file << setprecision(out_precision);
+              assert (file.good());
+
+            } // if on the root processor
+
+            // Traverse and components on this multigrid and
+            // refinement level and map
+            BEGIN_COMPONENT_LOOP(cgh, grouptype) {
+
+              const ggf<dim>* ff = 0;
+
+              assert (var < (int)arrdata.at(group).at(Carpet::map).data.size());
+              ff = (ggf<dim>*)arrdata.at(group).at(Carpet::map).data.at(var);
+
+              const int mintl = output_all_timelevels ? 1-num_tl : 0;
+              const int maxtl = 0;
+              for (int tl=mintl; tl<=maxtl; ++tl) {
+
+                const gdata<dim>* const data
+                  = (*ff) (tl, rl, component, mglevel);
+                ibbox ext = data->extent();
+
+                ivect lo = ext.lower();
+                ivect hi = ext.upper();
+                ivect str = ext.stride();
+
+                // Ignore ghost zones if desired
+                for (int d=0; d<dim; ++d) {
+                  bool output_lower_ghosts
+                    = cgh->cctk_bbox[2*d] ? out3D_outer_ghosts : out3D_ghosts;
+                  bool output_upper_ghosts
+                    = cgh->cctk_bbox[2*d+1] ? out3D_outer_ghosts : out3D_ghosts;
+
+                  if (! output_lower_ghosts) {
+                    lo[d] += cgh->cctk_nghostzones[d] * str[d];
+                  }
+                  if (! output_upper_ghosts) {
+                    hi[d] -= cgh->cctk_nghostzones[d] * str[d];
+                  }
+                }
+                ext = ibbox(lo,hi,str);
+
+                // coordinates
+                const CCTK_REAL coord_time = cgh->cctk_time;
+                rvect global_lower;
+                rvect coord_delta;
+                if (grouptype == CCTK_GF) {
+                  for (int d=0; d<dim; ++d) {
+                    global_lower[d] = cgh->cctk_origin_space[d];
+                    coord_delta[d] = cgh->cctk_delta_space[d] / maxreflevelfact;
+                  }
+                } else {
+                  for (int d=0; d<dim; ++d) {
+                    global_lower[d] = 0.0;
+                    coord_delta[d] = 1.0;
+                  }
+                }
+                const rvect coord_lower
+                  = global_lower + coord_delta * rvect(lo);
+                const rvect coord_upper
+                  = global_lower + coord_delta * rvect(hi);
+
+                ivect offset1;
+                if (grouptype == CCTK_GF) {
+                  const ibbox& baseext = vdd.at(Carpet::map)->bases.at(0).at(mglevel).exterior;
+                  offset1 = baseext.lower() + offset * ext.stride();
+                } else {
+                  offset1 = offset * ext.stride();
+                }
+                for (int d=0; d<outdim; ++d) {
+                  offset1[dirs[d]] = ext.lower()[dirs[d]];
+                }
+
+                WriteASCII (file, data, ext, n, cgh->cctk_iteration,
+                            offset1, dirs,
+                            rl, mglevel, Carpet::map, component, tl,
+                            coord_time, coord_lower, coord_upper);
+
+                // Append EOL after every component
+                if (CCTK_MyProc(cgh)==0) {
+                  if (separate_components) {
+                    assert (file.good());
+                    file << endl;
+                  }
+                }
+                assert (file.good());
+
+              } // for tl
+
+            } END_COMPONENT_LOOP;
+
+            // Append EOL after every complete set of components
+            if (CCTK_MyProc(cgh)==0) {
+              if (separate_grids) {
+                assert (file.good());
+                file << endl;
+              }
+              file.close();
+              assert (file.good());
+            }
+
+            assert (! file.is_open());
+
+          } END_MAP_LOOP;
+
+	} // if (desired)
+
+      } // if (ascending)
+
+      // Next direction combination
+      done = true;
+      for (int d=0; d<outdim; ++d) {
+	++dirs[d];
+	if (dirs[d]<groupdim) {
+	  done = false;
+	  break;
+	}
+	dirs[d] = 0;
+      }
+
+    } while (! done);		// all directions
+
+    // Don't truncate again
+    do_truncate.at(n) = false;
+
+    return 0;
+  }
+
+
+
+  template<int outdim>
+  int IOASCII<outdim>
+  ::TimeToOutput (const cGH* const cctkGH, const int vindex)
+  {
+    DECLARE_CCTK_ARGUMENTS;
+    DECLARE_CCTK_PARAMETERS;
+
+    assert (vindex>=0 && vindex<CCTK_NumVars());
+
+
+
+    const int grouptype = CCTK_GroupTypeFromVarI(vindex);
+    switch (grouptype) {
+    case CCTK_SCALAR:
+    case CCTK_ARRAY:
+      if (! do_global_mode) return 0;
+      break;
+    case CCTK_GF:
+      // do nothing
+      break;
+    default:
+      assert (0);
+    }
+
+
+
+    // check whether to output at this iteration
+    bool output_this_iteration;
+
+    const char* myoutcriterion = GetStringParameter("out%dD_criterion");
+    if (CCTK_EQUALS(myoutcriterion, "default")) {
+      myoutcriterion = out_criterion;
+    }
+
+    if (CCTK_EQUALS (myoutcriterion, "never")) {
+
+      // Never output
+      output_this_iteration = false;
+
+    } else if (CCTK_EQUALS (myoutcriterion, "iteration")) {
+
+      int myoutevery = GetIntParameter("out%dD_every");
+      if (myoutevery == -2) {
+        myoutevery = out_every;
+      }
+      if (myoutevery <= 0) {
+        // output is disabled
+        output_this_iteration = false;
+      } else if (cctk_iteration == this_iteration[outdim]) {
+        // we already decided to output this iteration
+        output_this_iteration = true;
+      } else if (cctk_iteration
+                 >= last_output_iteration[outdim] + myoutevery) {
+        // it is time for the next output
+        output_this_iteration = true;
+        last_output_iteration[outdim] = cctk_iteration;
+        this_iteration[outdim] = cctk_iteration;
+      } else {
+        // we want no output at this iteration
+        output_this_iteration = false;
+      }
+
+    } else if (CCTK_EQUALS (myoutcriterion, "divisor")) {
+
+      int myoutevery = GetIntParameter("out%dD_every");
+      if (myoutevery == -2) {
+        myoutevery = out_every;
+      }
+      if (myoutevery <= 0) {
+        // output is disabled
+        output_this_iteration = false;
+      } else if (cctk_iteration % myoutevery == 0) {
+        // we already decided to output this iteration
+        output_this_iteration = true;
+      } else {
+        // we want no output at this iteration
+        output_this_iteration = false;
+      }
+
+    } else if (CCTK_EQUALS (myoutcriterion, "time")) {
+
+      CCTK_REAL myoutdt = GetRealParameter("out%dD_dt");
+      if (myoutdt == -2) {
+        myoutdt = out_dt;
+      }
+      if (myoutdt < 0) {
+        // output is disabled
+        output_this_iteration = false;
+      } else if (myoutdt == 0) {
+        // output all iterations
+        output_this_iteration = true;
+      } else if (cctk_iteration == this_iteration[outdim]) {
+        // we already decided to output this iteration
+        output_this_iteration = true;
+      } else if (cctk_time / cctk_delta_time
+                 >= (last_output_time[outdim] + myoutdt) / cctk_delta_time - 1.0e-12) {
+        // it is time for the next output
+        output_this_iteration = true;
+        last_output_time[outdim] = cctk_time;
+        this_iteration[outdim] = cctk_iteration;
+      } else {
+        // we want no output at this iteration
+        output_this_iteration = false;
+      }
+
+    } else {
+
+      assert (0);
+
+    } // select output criterion
+
+    if (! output_this_iteration) return 0;
+
+
+
+    // check which variables to output
+    static vector<bool> output_variables;
+    static int output_variables_iteration = -1;
+
+    if (cctk_iteration > output_variables_iteration) {
+      output_variables.resize (CCTK_NumVars());
+
+      const char* const varlist = GetStringParameter("out%dD_vars");
+      if (CCTK_TraverseString (varlist, SetFlag, &output_variables,
+                               CCTK_GROUP_OR_VAR) < 0)
+      {
+        int abort_on_error = output_variables_iteration < 0 &&
+                             strict_io_parameter_check;
+        CCTK_VWarn (abort_on_error ? 0 : 1, __LINE__, __FILE__,CCTK_THORNSTRING,
+                    "error while parsing parameter 'IOASCII::out%dD_vars'",
+                    outdim);
+      }
+
+      output_variables_iteration = cctk_iteration;
+    }
+
+    if (! output_variables.at(vindex)) return 0;
+
+
+
+    if (last_output.at(mglevel).at(reflevel).at(vindex) == cctk_iteration) {
+      // Has already been output during this iteration
+      char* varname = CCTK_FullName(vindex);
+      CCTK_VWarn (5, __LINE__, __FILE__, CCTK_THORNSTRING,
+		  "Skipping output for variable \"%s\", because this variable "
+		  "has already been output during the current iteration -- "
+		  "probably via a trigger during the analysis stage",
+		  varname);
+      free (varname);
+      return 0;
+    }
+
+    assert (last_output.at(mglevel).at(reflevel).at(vindex) < cctk_iteration);
+
+    // Should be output during this iteration
+    return 1;
+  }
+
+
+
+  template<int outdim>
+  int IOASCII<outdim>
+  ::TriggerOutput (const cGH* const cgh, const int vindex)
+  {
+    assert (vindex>=0 && vindex<CCTK_NumVars());
+
+    char* varname = CCTK_FullName(vindex);
+    const int retval = OutputVarAs (cgh, varname, CCTK_VarName(vindex));
+    free (varname);
+
+    last_output.at(mglevel).at(reflevel).at(vindex) = cgh->cctk_iteration;
+
+    return retval;
+  }
+
+
+
+  template<int outdim>
+  int IOASCII<outdim>
+  ::GetGridOffset (const cGH* const cgh, const int dir,
+		   const char* const itempl, const char* const iglobal,
+		   const char* const ctempl, const char* const cglobal,
+		   const CCTK_REAL cfallback)
+  {
+    // First choice: explicit coordinate
+    char cparam[1000];
+    snprintf (cparam, sizeof cparam, ctempl, outdim);
+    const int ncparam = CCTK_ParameterQueryTimesSet (cparam, CCTK_THORNSTRING);
+    assert (ncparam >= 0);
+    if (ncparam > 0) {
+      int ptype;
+      const CCTK_REAL* const pcoord
+	= ((const CCTK_REAL*)CCTK_ParameterGet
+           (cparam, CCTK_THORNSTRING, &ptype));
+      assert (pcoord);
+      const CCTK_REAL coord = *pcoord;
+      assert (ptype == PARAMETER_REAL);
+      return CoordToOffset (cgh, dir, coord, 0);
+    }
+
+    // Second choice: explicit index
+    char iparam[1000];
+    snprintf (iparam, sizeof iparam, itempl, outdim);
+    const int niparam = CCTK_ParameterQueryTimesSet (iparam, CCTK_THORNSTRING);
+    assert (niparam >= 0);
+    if (niparam > 0) {
+      int ptype;
+      const int* const pindex
+	= (const int*)CCTK_ParameterGet (iparam, CCTK_THORNSTRING, &ptype);
+      assert (pindex);
+      const int index = *pindex;
+      assert (ptype == PARAMETER_INT);
+      return index;
+    }
+
+    // Third choice: explicit global coordinate
+    const char* iothorn = CCTK_ImplementationThorn ("IO");
+    assert (iothorn);
+    if (cglobal) {
+      const int ncglobal = CCTK_ParameterQueryTimesSet (cglobal, iothorn);
+      assert (ncglobal >= 0);
+      if (ncglobal > 0) {
+        int ptype;
+        const CCTK_REAL* const pcoord
+          = (const CCTK_REAL*)CCTK_ParameterGet (cglobal, iothorn, &ptype);
+        assert (pcoord);
+        const CCTK_REAL coord = *pcoord;
+        assert (ptype == PARAMETER_REAL);
+        return CoordToOffset (cgh, dir, coord, 0);
+      }
+    }
+
+    // Fourth choice: explicit global index
+    if (iglobal) {
+      const int niglobal = CCTK_ParameterQueryTimesSet (iglobal, iothorn);
+      assert (niglobal >= 0);
+      if (niglobal > 0) {
+        int ptype;
+        const int* const pindex
+          = (const int*)CCTK_ParameterGet (iglobal, iothorn, &ptype);
+        assert (pindex);
+        const int index = *pindex;
+        assert (ptype == PARAMETER_INT);
+        return index;
+      }
+    }
+
+    // Fifth choice: default coordinate
+    return CoordToOffset (cgh, dir, cfallback, 0);
+  }
+
+
+
+  template<int outdim>
+  int IOASCII<outdim>
+  ::CoordToOffset (const cGH* cgh, const int dir, const CCTK_REAL coord,
+		   const int ifallback)
+  {
+    assert (dir>=1 && dir<=dim);
+
+    assert (mglevel!=-1 && reflevel!=-1 && Carpet::map!=-1);
+
+    const CCTK_REAL delta = cgh->cctk_delta_space[dir-1] / cgh->cctk_levfac[dir-1];
+    const CCTK_REAL lower = cgh->cctk_origin_space[dir-1];
+#if 0
+    const int npoints = cgh->cctk_gsh[dir-1];
+    const CCTK_REAL upper = lower + (npoints-1) * delta;
+#endif
+
+    const CCTK_REAL rindex = (coord - lower) / delta;
+    int cindex = (int)floor(rindex + 0.75);
+
+#if 0
+    if (cindex<0 || cindex>=npoints) {
+      cindex = ifallback;
+
+      assert (dir>=1 && dir<=3);
+      CCTK_VWarn (1, __LINE__, __FILE__, CCTK_THORNSTRING,
+                  "The specified coordinate value %g for the %c-direction is not within the grid range [%g,%g] on convergence level %d, refinement level %d, map %d; using %g instead",
+                  coord, "xyz"[dir-1], lower, upper,
+                  mglevel, reflevel, Carpet::map, lower + delta * cindex);
+    }
+
+    assert (cindex>=0 && cindex<npoints);
+#else
+    const void *dummy;
+    dummy = &ifallback;
+    dummy = &dummy;
+#endif
+
+    return cindex;
+  }
+
+
+
+  template<int outdim>
+  const char* IOASCII<outdim>
+  ::GetStringParameter (const char* const parametertemplate)
+  {
+    char parametername[1000];
+    snprintf (parametername, sizeof parametername, parametertemplate, outdim);
+    int ptype;
+    const char* const* const ppval = (const char* const*)CCTK_ParameterGet
+      (parametername, CCTK_THORNSTRING, &ptype);
+    assert (ppval);
+    const char* const pval = *ppval;
+    assert (ptype == PARAMETER_STRING || ptype == PARAMETER_KEYWORD);
+    return pval;
+  }
+
+
+
+  template<int outdim>
+  CCTK_INT IOASCII<outdim>
+  ::GetIntParameter (const char* const parametertemplate)
+  {
+    char parametername[1000];
+    snprintf (parametername, sizeof parametername, parametertemplate, outdim);
+    int ptype;
+    const CCTK_INT* const ppval
+      = (const CCTK_INT*)CCTK_ParameterGet
+      (parametername, CCTK_THORNSTRING, &ptype);
+    assert (ppval);
+    assert (ptype == PARAMETER_INT || ptype == PARAMETER_BOOLEAN);
+    const CCTK_INT pval = *ppval;
+    return pval;
+  }
+
+
+
+  template<int outdim>
+  CCTK_REAL IOASCII<outdim>
+  ::GetRealParameter (const char* const parametertemplate)
+  {
+    char parametername[1000];
+    snprintf (parametername, sizeof parametername, parametertemplate, outdim);
+    int ptype;
+    const CCTK_REAL* const ppval
+      = (const CCTK_REAL*)CCTK_ParameterGet
+      (parametername, CCTK_THORNSTRING, &ptype);
+    assert (ppval);
+    assert (ptype == PARAMETER_REAL);
+    const CCTK_REAL pval = *ppval;
+    return pval;
+  }
+
+
+
+  void SetFlag (int index, const char* optstring, void* arg)
+  {
+    optstring = optstring;
+    vector<bool>& flags = *(vector<bool>*)arg;
+    flags.at(index) = true;
+  }
+
+
+
+  // Output
+  template<int D,int DD>
+  void WriteASCII (ostream& os,
+		   const gdata<D>* const gfdata,
+		   const bbox<int,D>& gfext,
+		   const int vi,
+		   const int time,
+		   const vect<int,D>& org,
+		   const vect<int,DD>& dirs,
+		   const int rl,
+		   const int ml,
+                   const int m,
+		   const int c,
+		   const int tl,
+		   const CCTK_REAL coord_time,
+		   const vect<CCTK_REAL,D>& coord_lower,
+		   const vect<CCTK_REAL,D>& coord_upper)
+  {
+    assert (DD<=D);
+
+    if (gfdata->proc()==0) {
+      // output on processor 0
+
+      int rank;
+      MPI_Comm_rank (dist::comm, &rank);
+      if (rank == 0) {
+
+	assert (os.good());
+
+	os << "# iteration " << time << endl
+	   << "# refinement level " << rl
+           << "   multigrid level " << ml
+	   << "   map " << m
+	   << "   component " << c
+           << "   time level " << tl
+           << endl
+	   << "# column format: it\ttl rl c ml\t";
+	assert (D>=1 && D<=3);
+	const char* const coords = "xyz";
+	for (int d=0; d<D-1; ++d) os << "i" << coords[d] << " "; os << "i" << coords[D-1];
+	os << "\ttime\t";
+	for (int d=0; d<D-1; ++d) os << coords[d] << " "; os << coords[D-1];
+	os << "\tdata" << endl;
+
+	const vect<int,DD> lo = gfext.lower()[dirs];
+	const vect<int,DD> up = gfext.upper()[dirs];
+	const vect<int,DD> str = gfext.stride()[dirs];
+	const bbox<int,DD> ext(lo,up,str);
+
+	// Check whether the output origin is contained in the extent
+	// of the data that should be output
+	ivect org1(org);
+	for (int d=0; d<DD; ++d) org1[dirs[d]] = ext.lower()[d];
+	if (gfext.contains(org1)) {
+
+          typename bbox<int,DD>::iterator it=ext.begin();
+          do {
+
+	    ivect index(org);
+	    for (int d=0; d<DD; ++d) index[dirs[d]] = (*it)[d];
+	    os << time << "\t" << tl << " " << rl << " " << c << " " << ml
+               << "\t";
+	    for (int d=0; d<D-1; ++d) os << index[d] << " "; os << index[D-1];
+	    os << "\t" << coord_time << "\t";
+	    for (int d=0; d<D; ++d) {
+	      assert (gfext.upper()[d] - gfext.lower()[d] >= 0);
+	      if (gfext.upper()[d] - gfext.lower()[d] == 0) {
+                os << coord_lower[d];
+              } else {
+                os << (coord_lower[d] + (index[d] - gfext.lower()[d])
+                              * (coord_upper[d] - coord_lower[d])
+                              / (gfext.upper()[d] - gfext.lower()[d]));
+              }
+              if (d != D-1) os << " ";
+	    }
+	    os << "\t";
+	    switch (CCTK_VarTypeI(vi)) {
+#define TYPECASE(N,T)					\
+	    case N:					\
+	      os << (*(const data<T,D>*)gfdata)[index];	\
+	      break;
+#include "Carpet/Carpet/src/typecase"
+#undef TYPECASE
+	    default:
+	      UnsupportedVarType(vi);
+	    }
+	    os << endl;
+
+            ++it;
+
+	    for (int d=0; d<DD; ++d) {
+	      if ((*it)[d]!=(*ext.end())[d]) break;
+	      os << endl;
+	    }
+
+          } while (it!=ext.end());
+
+	} else {
+
+          os << "#" << endl;
+
+	} // if ! ext contains org
+
+	assert (os.good());
+
+      }
+
+    } else {
+      // copy to processor 0 and output there
+
+      gdata<D>* const tmp = gfdata->make_typed(vi);
+      tmp->allocate(gfdata->extent(), 0);
+      for (comm_state<dim> state; !state.done(); state.step()) {
+        tmp->copy_from (state, gfdata, gfdata->extent());
+      }
+      WriteASCII (os, tmp, gfext, vi, time, org, dirs, rl, ml, m, c, tl,
+		  coord_time, coord_lower, coord_upper);
+      delete tmp;
+
+    }
+  }
+
+
+
+
+
+  // Explicit instantiation for all output dimensions
+  template class IOASCII<0>;
+  template class IOASCII<1>;
+  template class IOASCII<2>;
+  template class IOASCII<3>;
+
+  template
+  void WriteASCII (ostream& os,
+		   const gdata<3>* const gfdata,
+		   const bbox<int,3>& gfext,
+		   const int vi,
+		   const int time,
+		   const vect<int,3>& org,
+		   const vect<int,0>& dirs,
+		   const int rl,
+		   const int ml,
+		   const int m,
+		   const int c,
+		   const int tl,
+		   const CCTK_REAL coord_time,
+		   const vect<CCTK_REAL,3>& coord_lower,
+		   const vect<CCTK_REAL,3>& coord_upper);
+
+  template
+  void WriteASCII (ostream& os,
+		   const gdata<3>* const gfdata,
+		   const bbox<int,3>& gfext,
+		   const int vi,
+		   const int time,
+		   const vect<int,3>& org,
+		   const vect<int,1>& dirs,
+		   const int rl,
+		   const int ml,
+		   const int m,
+		   const int c,
+		   const int tl,
+		   const CCTK_REAL coord_time,
+		   const vect<CCTK_REAL,3>& coord_lower,
+		   const vect<CCTK_REAL,3>& coord_upper);
+
+  template
+  void WriteASCII (ostream& os,
+		   const gdata<3>* const gfdata,
+		   const bbox<int,3>& gfext,
+		   const int vi,
+		   const int time,
+		   const vect<int,3>& org,
+		   const vect<int,2>& dirs,
+		   const int rl,
+		   const int ml,
+		   const int m,
+		   const int c,
+		   const int tl,
+		   const CCTK_REAL coord_time,
+		   const vect<CCTK_REAL,3>& coord_lower,
+		   const vect<CCTK_REAL,3>& coord_upper);
+
+  template
+  void WriteASCII (ostream& os,
+		   const gdata<3>* const gfdata,
+		   const bbox<int,3>& gfext,
+		   const int vi,
+		   const int time,
+		   const vect<int,3>& org,
+		   const vect<int,3>& dirs,
+		   const int rl,
+		   const int ml,
+		   const int m,
+		   const int c,
+		   const int tl,
+		   const CCTK_REAL coord_time,
+		   const vect<CCTK_REAL,3>& coord_lower,
+		   const vect<CCTK_REAL,3>& coord_upper);
+
+} // namespace CarpetIOASCII
diff --git a/Carpet/CarpetIOASCII/src/ioascii.h b/Carpet/CarpetIOASCII/src/ioascii.h
new file mode 100644
index 000000000..60c8b9e02
--- /dev/null
+++ b/Carpet/CarpetIOASCII/src/ioascii.h
@@ -0,0 +1,22 @@
+/* $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetIOASCII/src/ioascii.h,v 1.7 2004/05/21 18:10:37 schnetter Exp $ */
+
+#ifndef CARPETIOASCII_H
+#define CARPETIOASCII_H
+
+#include "cctk_Arguments.h"
+
+#ifdef __cplusplus
+namespace CarpetIOASCII {
+  extern "C" {
+#endif
+    
+    /* Scheduled functions */
+    int CarpetIOASCIIStartup (void);
+    void CarpetIOASCIIInit (CCTK_ARGUMENTS);
+    
+#ifdef __cplusplus
+  } /* extern "C" */
+} /* namespace CarpetIOASCII */
+#endif
+
+#endif /* !defined(CARPETIOASCII_H) */
diff --git a/Carpet/CarpetIOASCII/src/ioascii.hh b/Carpet/CarpetIOASCII/src/ioascii.hh
new file mode 100644
index 000000000..72ea5b97b
--- /dev/null
+++ b/Carpet/CarpetIOASCII/src/ioascii.hh
@@ -0,0 +1,73 @@
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetIOASCII/src/ioascii.hh,v 1.16 2004/04/03 12:38:12 schnetter Exp $
+
+#ifndef CARPETIOASCII_HH
+#define CARPETIOASCII_HH
+
+#include <vector>
+
+#include "cctk.h"
+
+#include "ioascii.h"
+
+
+
+namespace CarpetIOASCII {
+  
+  using namespace std;
+  
+  
+  
+  // Everything is a class template, so that it can easily be
+  // instantiated for all output dimensions.
+  
+  template<int outdim>
+  struct IOASCII {
+    
+    // handle from CCTK_RegisterGHExtension
+    static int GHExtension;
+    
+    // handles from CCTK_RegisterIOMethed
+    static int IOMethod;
+    
+    
+    
+    // Do truncate the output files for a variable
+    static vector<bool> do_truncate;
+    
+    // Last iteration on which a refinement level of a variable was
+    // output (INT_MIN for none)
+    static vector<vector<vector<int> > > last_output; // [ml][rl][var]
+    
+    
+    
+    // scheduled functions
+    static int Startup();
+    
+    
+    
+    // registered functions
+    
+    static void* SetupGH (tFleshConfig* fc, int convLevel, cGH* cgh);
+    
+    static int OutputGH (const cGH* cgh);
+    static int OutputVarAs (const cGH* cgh,
+			    const char* varname, const char* alias);
+    static int TimeToOutput (const cGH* cgh, int vindex);
+    static int TriggerOutput (const cGH* cgh, int vindex);
+    
+    static int GetGridOffset (const cGH* cgh, int dir,
+			      const char* itempl, const char* iglobal,
+			      const char* ctempl, const char* cglobal,
+			      CCTK_REAL cfallback);
+    static int CoordToOffset (const cGH* cgh, int dir, CCTK_REAL coord,
+			      int ifallback);
+    
+    static const char* GetStringParameter (const char* parametertemplate);
+    static CCTK_INT GetIntParameter (const char* parametertemplate);
+    static CCTK_REAL GetRealParameter (const char* parametertemplate);
+    
+  };                            // struct IOASCII
+  
+} // namespace CarpetIOASCII
+
+#endif // !defined(CARPETIOASCII_HH)
diff --git a/Carpet/CarpetIOASCII/src/make.code.defn b/Carpet/CarpetIOASCII/src/make.code.defn
new file mode 100644
index 000000000..e9b78bce0
--- /dev/null
+++ b/Carpet/CarpetIOASCII/src/make.code.defn
@@ -0,0 +1,9 @@
+# Main make.code.defn file for thorn CarpetIOASCII
+# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetIOASCII/src/make.code.defn,v 1.1 2001/03/01 13:40:10 eschnett Exp $
+
+# Source files in this directory
+SRCS = ioascii.cc
+
+# Subdirectories containing source files
+SUBDIRS = 
+
diff --git a/Carpet/CarpetIOASCII/src/make.configuration.defn b/Carpet/CarpetIOASCII/src/make.configuration.defn
new file mode 100644
index 000000000..939aee820
--- /dev/null
+++ b/Carpet/CarpetIOASCII/src/make.configuration.defn
@@ -0,0 +1,4 @@
+# Main make.configuration.defn file for thorn CarpetIOASCII
+# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetIOASCII/src/make.configuration.defn,v 1.1 2002/09/30 15:36:28 schnetter Exp $
+
+ALL_UTILS += carpet2sdf carpet2xgraph
diff --git a/Carpet/CarpetIOASCII/src/make.configuration.deps b/Carpet/CarpetIOASCII/src/make.configuration.deps
new file mode 100644
index 000000000..686ad3112
--- /dev/null
+++ b/Carpet/CarpetIOASCII/src/make.configuration.deps
@@ -0,0 +1,53 @@
+# Main make.configuration.deps file for thorn CarpetIOASCII
+# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetIOASCII/src/make.configuration.deps,v 1.1 2002/09/30 15:36:28 schnetter Exp $
+
+
+
+SDF_INCDIRS := $(HOME)/include
+SDF_LIBDIRS := $(HOME)/lib
+SDF_LIBS := rnpl mfhdf df jpeg z vsso sv m
+
+
+
+# Compile
+
+$(BUILD_DIR)$(DIRSEP)CarpetIOASCII$(DIRSEP)%.o: $(PACKAGE_DIR)$(DIRSEP)Carpet$(DIRSEP)CarpetIOASCII$(DIRSEP)src$(DIRSEP)util$(DIRSEP)%.c
+	@echo "Compiling $<"
+	-$(MKDIR) $(MKDIRFLAGS) $(BUILD_DIR)$(DIRSEP)CarpetIOASCII 2> /dev/null
+	$(CC) $< -DCCODE $(CFLAGS) -I$(CONFIG) -I$(BINDINGS_DIR)$(DIRSEP)include -I$(FLESH_DIR)$(DIRSEP)include -I$(CCTK_HOME)$(DIRSEP)arrangement $(CCOMPILEONLY)$(OPTIONSEP)$@
+
+
+
+# Link
+
+$(UTIL_DIR)$(DIRSEP)%: $(BUILD_DIR)$(DIRSEP)CarpetIOASCII$(DIRSEP)%.o
+	@echo "Creating $* in $(UTIL_DIR) from $<"
+	-$(MKDIR) $(MKDIRFLAGS) $(UTIL_DIR) 2> /dev/null
+	$(LD) $(CREATEEXE)$(OPTIONSEP)$@ $(DEBUG_LD) $(LDFLAGS) $(EXTRAFLAGS) $<
+
+
+
+# Special versions for carpet2sdf:
+
+# Compile
+
+$(BUILD_DIR)$(DIRSEP)CarpetIOASCII$(DIRSEP)carpet2sdf.o: $(PACKAGE_DIR)$(DIRSEP)Carpet$(DIRSEP)CarpetIOASCII$(DIRSEP)src$(DIRSEP)util$(DIRSEP)carpet2sdf.c
+	@echo "Compiling carpet2sdf"
+	-$(MKDIR) $(MKDIRFLAGS) $(BUILD_DIR)$(DIRSEP)CarpetIOASCII 2> /dev/null
+	-$(CC) $< -DCCODE $(CFLAGS) -I$(CONFIG) -I$(BINDINGS_DIR)$(DIRSEP)include -I$(FLESH_DIR)$(DIRSEP)include -I$(CCTK_HOME)$(DIRSEP)arrangements -I$(SDF_INCDIRS:%=-I%) $(CCOMPILEONLY)$(OPTIONSEP)$@
+
+
+
+# Link
+
+$(UTIL_DIR)$(DIRSEP)carpet2sdf: $(BUILD_DIR)$(DIRSEP)CarpetIOASCII$(DIRSEP)carpet2sdf.o
+	@echo "Creating $* in $(UTIL_DIR) from $<"
+	-$(MKDIR) $(MKDIRFLAGS) $(UTIL_DIR) 2> /dev/null
+	-$(LD) $(CREATEEXE)$(OPTIONSEP)$@ $(DEBUG_LD) $(LDFLAGS) $(EXTRAFLAGS) $< $(SDF_LIBDIRS:%=-L%) $(SDF_LIBS:%=-l%)
+	@if [ ! -e $(UTIL_DIR)$(DIRSEP)/carpet2sdf ]; then \
+	echo "*************************************"; \
+	echo "Warning: could not install carpet2sdf"; \
+	echo "*************************************"; \
+	echo "echo \"carpet2sdf could not be installed\"" > $@; \
+	chmod a+x $@; \
+	fi
diff --git a/Carpet/CarpetIOASCII/src/util/Carpet2ygraph.pl b/Carpet/CarpetIOASCII/src/util/Carpet2ygraph.pl
new file mode 100755
index 000000000..0f49aaf01
--- /dev/null
+++ b/Carpet/CarpetIOASCII/src/util/Carpet2ygraph.pl
@@ -0,0 +1,121 @@
+#! /usr/bin/perl -s
+#
+# Blame Ian Hawke for not checking that Scott had already written a C
+# program to do this add so writing a perl version instead
+#
+# Given an output file in CarpetIOASCII 1d format, strip to ygraph format.
+# The arguments should be direction (x=0,y=1,z=2) and filename.
+# Output is to a file. Only the base name should be given.
+# The base will be appended with _<level>.xg
+#
+# Example:
+#
+# Carpet2ygraph.pl 0 alp.xl alp_x
+#
+# produces alp_x_<d>.xg where <d> is the refinement level.
+#
+# The headers of the Carpet files are just copied with the addition of
+# the correct time for ygraph, and the different comment marker.
+#
+# This script is a much altered version of Tom Goodale's convergence
+# testing scripts.
+#
+
+use FileHandle;
+
+if(@ARGV != 3)
+{
+  print "Usage: $0 direction <Inputfile> <Outputfile> \n";
+  exit;
+}
+
+open(CARPETFILE,   "<$ARGV[1]") || die "Unable to find file \"$ARGV[1]\"."; 
+
+#
+# Open the output file for the base grid. 
+#
+
+$file = $ARGV[2]."_0.xg";
+my $fh = new FileHandle(">$file") || die "Unable to open file \"$file\"."; 
+push(@outputfilelist,$fh);
+
+#
+# Find the correct column for the spatial coordinate; requires a magic number
+#
+
+$direction = $ARGV[0]+9; 
+
+$flag = 0;
+$timeset = 0;
+$refinementlevel = 0;
+$componentflag[0] = 0;
+while (<CARPETFILE>)
+{
+    $line = $_;
+
+    if(/^\#/) # The line is a header comment
+    {
+	if ($flag==1) # It's a new level and there is data to output
+        {
+	    $fh = $outputfilelist[$refinementlevel];
+	    print $fh @outputdata;
+	    @outputdata=("\n");
+	    $flag = 0;
+	}
+	if ($line =~ /refinement level ([0-9])/) # Line gives ref. level
+	{
+	    $refinementlevel = $1;
+	    $line =~ /component ([0-9+])/;
+	    $componentflag[$refinementlevel] = $1; # Which component?
+
+	    #
+	    # If no file exists for this refinement level, 
+	    # open and add filehandle to array
+	    #
+
+	    if ($refinementlevel > $#outputfilelist)
+	    { 
+		for ($i = $#outputfileflist+1; $i < $refinementlevel+1; $i++)
+		{
+		    $file = $ARGV[2]."_".$i.".xg";
+		    my $fh = new FileHandle(">$file") || die "Unable to open file \"$file\"."; 
+		    $outputfilelist[$i]=$fh;
+		}
+	    }
+	}
+	# Only output the headers if this is the zero component
+	# FIXME: what happens if component 0 isn't output first?
+	if (0 == $componentflag[$refinementlevel])
+	{
+	    push(@outputdata, ("\"",$line)); # Add ygraph comment marker
+	}
+	else
+	{
+	    $flag = 1;
+	    @outputdata=("");
+	}
+    }
+    else # The line contains real data
+    {
+	@data = split(/ \t]+/,$line);
+	if ($flag== 0) # This is the first line of data
+	{
+	    $flag = 1;
+	    $timeset = $data[8]; # Magic number gives the Cactus time
+	    @outputdata = ("\n\"Time = ",$timeset,@outputdata);
+	}
+	push(@outputdata, $data[$direction], " ", $data[12]);
+    }
+}
+
+#
+# At end of file, so output final data set.
+#
+
+$fh = $outputfilelist[$refinementlevel];
+print $fh @outputdata;
+
+foreach $fh (@outputfilelist)
+{
+    close($fh);
+}
diff --git a/Carpet/CarpetIOASCII/src/util/Carpet2ygraphCat.pl b/Carpet/CarpetIOASCII/src/util/Carpet2ygraphCat.pl
new file mode 100755
index 000000000..669c8e64b
--- /dev/null
+++ b/Carpet/CarpetIOASCII/src/util/Carpet2ygraphCat.pl
@@ -0,0 +1,105 @@
+#! /usr/bin/perl -s
+
+use FileHandle;
+
+if(@ARGV != 3)
+{
+  print "Usage: $0 direction <Inputfile> <Outputfile> \n";
+  exit;
+}
+
+open(CARPETFILE,   "<$ARGV[1]") || die "Unable to find file \"$ARGV[1]\".";
+
+$file = $ARGV[2].".xg";
+my $fh = new FileHandle(">$file") || die "Unable to open file \"$file\".";
+
+my %data;
+my $time = -1;
+my $new = 0;
+my $currentit = -1;
+my $lastit = -1;
+
+my @datatoprint;
+my $nsets = 1;
+my $maxlength = 0;
+my @lengths;
+
+$lengths[0]=0;
+while (<CARPETFILE>)
+{
+    $line = $_;
+    if ($line =~ "iteration") {
+      @itline = split(/ +/,$line);
+      $currentit = @itline[2];
+    }
+    elsif ($line =~ /^#/)
+    {
+	#Do nothing for headers!
+    }
+    elsif ($line =~ /([0-9+-ed.]+)*/)
+    {
+	@dataline = split(/[ \t]+/,$line);
+	if ($currentit != $lastit)
+	{
+	    if ($new)
+	    {
+		push(@datatoprint,"\n\n\"Time = ".$time."\n");
+		my @sortedcoords = sort numerically (keys %data);
+		my $localcoord;
+		foreach $localcoord (@sortedcoords)
+		{
+		    push(@datatoprint, $localcoord." ".$data{$localcoord});
+		}
+		$maxlength = $maxlength > (scalar @sortedcoords) ? $maxlength : (scalar @sortedcoords);
+		$lengths[$nsets-1]=(scalar @sortedcoords);
+		$nsets++;
+		$lengths[$nsets-1]=0;
+		%data=();
+	    }
+	    $new++;
+	    $time = $dataline[8];
+	    $lastit = $currentit;
+	}
+	my $coord = $dataline[9+$ARGV[0]];
+	my $val = $dataline[12];
+	$data{$coord} = $val;
+    }
+}
+push(@datatoprint, "\n\"Time = ",$time,"\n");
+my @sortedcoords = sort numerically (keys %data);
+my $localcoord;
+foreach $localcoord (@sortedcoords)
+{
+    push(@datatoprint, $localcoord." ".$data{$localcoord});
+}
+$maxlength = $maxlength > (scalar @sortedcoords) ? $maxlength : (scalar @sortedcoords);
+$lengths[$nsets-1]=(scalar @sortedcoords);
+$nsets++;
+$lengths[$nsets-1]=0;
+
+my $oldline="";
+$nouts=0;
+my $set=0;
+foreach $line (@datatoprint) {
+  if ($line =~ "Time") {
+    if ($oldline) {
+      for (my $i=$lengths[$set-1]; $i<$maxlength;$i++) {
+	$nouts++;
+	print $fh $oldline;
+      }
+    }
+    $set++;
+    print $fh $line;
+  }
+  else {
+    $nouts++;
+    print $fh $line;
+    $oldline=$line
+  }
+}
+for (my $i=$lengths[$set-1]; $i<$maxlength;$i++) {
+  $nouts++;
+  print $fh $oldline;
+}
+
+sub numerically {$a <=> $b;}
diff --git a/Carpet/CarpetIOASCII/src/util/carpet2sdf.c b/Carpet/CarpetIOASCII/src/util/carpet2sdf.c
new file mode 100644
index 000000000..a8abe71d0
--- /dev/null
+++ b/Carpet/CarpetIOASCII/src/util/carpet2sdf.c
@@ -0,0 +1,167 @@
+

+/*******************************************************************

+ * Program: carpet2sdf

+ * Description: Converts from Carpet ASCII file format to SDF format

+ * Author:  Scott H. Hawley

+ * Date:    June 17, 2002

+ *

+ * Similar to carpet2sdf, except that instead of sending output to

+ * sdtout, it sends it to <infile>.sdf

+ *******************************************************************/

+

+#include <stdio.h>

+#include <stdlib.h>

+#include <string.h>

+#include <bbhutil.h>

+

+

+/*******************************************************************

+ * Function: read_next_set

+ * Supposed to read one time step of data (for one level of refinement)

+ *******************************************************************/

+int read_next_set(FILE *infile, int rlev, int *numelems, 

+		  double *time,

+		  double **coord, double **data, int clip_data, double clip_val)

+{

+   char   in_line[200];

+   int    it,tl,rl,c,ml,ix,iy,iz;

+   double xval, yval, zval;

+   double coordval, dataval;

+   int    hit_first_non_blank = 0;

+   int    stop_reading = 0;

+   int    retval = 0;

+   int    rc = 0;

+

+   *numelems = 0;

+   sprintf(in_line," ");

+

+   /* reads up to first blank line after non-blank lines, or up to

+    * eof */

+   while ((fgets(in_line,200,infile) != NULL) && (!stop_reading) ) {

+      

+      if (in_line[0]=='#') {

+         hit_first_non_blank = 1;

+

+      } else if ((strncmp("\n",in_line,2)==0) ||

+         (in_line[0]==' ')) {

+         if (hit_first_non_blank) {

+            stop_reading = 1;

+	 }

+

+      } else {     /* Assume that we're dealing with a line of numbers */

+         hit_first_non_blank = 1;

+	 retval = sscanf(in_line,"%d %d %d %d %d %d %d %d %lg %lg %lg %lg %lg",

+	   &it,&tl,&rl,&c,&ml,&ix,&iy,&iz,time,&xval,&yval,&zval,&dataval);

+	 coordval = xval;

+   

+	 /* Only add input to arrays if it's for the right ref. level */

+         if ( (retval == 13) && (rl == rlev)) {

+	       (*numelems)++;

+   

+	       if (coord==NULL) {

+	         *coord = (double *)malloc(sizeof(double)* (*numelems));

+	         *data = (double *)malloc(sizeof(double)* (*numelems));

+	       } else {

+	         *coord = (double *)realloc(*coord, sizeof(double)* (*numelems));

+	         *data = (double *)realloc(*data, sizeof(double)* (*numelems));

+	       }

+	       if (coord==NULL) {

+		       fprintf(stderr,"Error: malloc/realloc returned NULL\n");

+		       exit(1);

+	       }

+   

+	       (*coord)[(*numelems)-1] = coordval;

+

+	       /* If desired, "clip" data above a certain value */

+	       if (clip_data && (dataval > clip_val)) {

+	          (*data)[(*numelems)-1] = clip_val;

+	       } else {

+	          (*data)[(*numelems)-1] = dataval;

+	       }

+

+	 }

+

+      }

+

+   }

+   rc = feof(infile);

+

+   return rc;

+}

+

+

+/********************************************************************

+ * Main part of program

+ ********************************************************************/

+int main(int argc, char **argv)

+{

+   char   *infilename;

+   FILE   *infile;

+

+   char   func_name[200];

+   double time;

+   int    shape[1];

+   double bounds[2];

+   int    rank = 1;

+   int    numelems_in_set = 0;

+   int    rlev = 0;

+   double *coord, *data;

+   int clip_data = 0;

+   double clip_val;

+

+

+   /* 

+    * Parse command-line arguments

+    */

+   if (argc <= 1) {

+      fprintf(stderr,"usage: carpet2sdf <infile> [ref_lev] [clip_data_at]\n");

+      exit(1);

+   }

+   if (argc >= 3) {

+      sscanf(argv[2],"%d",&rlev);

+   }

+   if (argc >= 4) {

+      sscanf(argv[3],"%lf",&clip_val);

+      clip_data = 1;

+   }

+

+   /*

+    * Open the input file

+    */

+   infilename = argv[1];

+   infile = fopen(infilename,"r");

+   if (infile == NULL) {

+      fprintf(stderr,"Error opening file '%s'.\n",infilename);

+      exit(1);

+   }

+

+

+   time = 0;

+   coord = NULL;

+   data = NULL;

+   sprintf(func_name,"%s_%d",infilename,rlev);

+

+   /*

+    * Main loop for reading from input file and writing to output file

+    */

+   while ( read_next_set(infile,rlev,&numelems_in_set,&time,&coord,&data,

+			    clip_data, clip_val) == 0 ) {

+      if (numelems_in_set > 0) {

+

+         shape[0] = numelems_in_set;

+         bounds[0] = coord[0];

+         bounds[1] = coord[numelems_in_set-1];

+

+         gft_out_bbox(func_name, time, shape, rank, bounds, data);

+

+         free(coord);

+         free(data);

+         coord = NULL;

+         data = NULL;

+

+      }

+   }

+

+   return 0;

+}

+

diff --git a/Carpet/CarpetIOASCII/src/util/carpet2xgraph.c b/Carpet/CarpetIOASCII/src/util/carpet2xgraph.c
new file mode 100644
index 000000000..61dfa689b
--- /dev/null
+++ b/Carpet/CarpetIOASCII/src/util/carpet2xgraph.c
@@ -0,0 +1,186 @@
+
+/*******************************************************************
+ * Program: carpet2xgraph
+ * Description: Converts from Carpet ASCII file format to Xgraph format
+ * Author:  Scott H. Hawley
+ * Date:    June 17, 2002
+ *
+ * This will select the data for ONE REFINEMENT LEVEL (default 0)
+ * and send the result to stdout.   The "xgraph format" data is
+ * not, however, fully annotated and someone may wish to improve 
+ * this.
+ *
+ * The [clip_data_at] command-line argument was added for use with
+ * 1/r data, to replace the infinities at the origin with whatever
+ * value the user specifies.
+ *******************************************************************/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
+
+/*******************************************************************
+ * Function: read_next_set
+ * Supposed to read one time step of data (for one level of refinement)
+ *******************************************************************/
+int read_next_set(FILE *infile, int rlev, int *numelems, 
+		  double *time,
+		  double **coord, double **data, int clip_data, double clip_val)
+{
+   char   in_line[200];
+   int    it,tl,rl,c,ml,ix,iy,iz;
+   double xval, yval, zval;
+   double coordval, dataval;
+   int    hit_first_non_blank = 0;
+   int    stop_reading = 0;
+   int    retval = 0;
+   int    rc = 0;
+
+   *numelems = 0;
+   sprintf(in_line," ");
+
+   /* reads up to first blank line after non-blank lines, or up to
+    * eof */
+   while ((fgets(in_line,200,infile) != NULL) && (!stop_reading) ) {
+      
+      if (in_line[0]=='#') {
+         hit_first_non_blank = 1;
+
+      } else if ((strncmp("\n",in_line,2)==0) ||
+         (in_line[0]==' ')) {
+         if (hit_first_non_blank) {
+            stop_reading = 1;
+	 }
+
+      } else {     /* Assume that we're dealing with a line of numbers */
+         hit_first_non_blank = 1;
+	 retval = sscanf(in_line,"%d %d %d %d %d %d %d %d %lg %lg %lg %lg %lg",
+	   &it,&tl,&rl,&c,&ml,&ix,&iy,&iz,time,&xval,&yval,&zval,&dataval);
+	 coordval = xval;
+   
+	 /* Only add input to arrays if it's for the right ref. level */
+         if ( (retval == 13) && (rl == rlev)) {
+	       (*numelems)++;
+   
+	       if (coord==NULL) {
+	         *coord = (double *)malloc(sizeof(double)* (*numelems));
+	         *data = (double *)malloc(sizeof(double)* (*numelems));
+	       } else {
+	         *coord = (double *)realloc(*coord, sizeof(double)* (*numelems));
+	         *data = (double *)realloc(*data, sizeof(double)* (*numelems));
+	       }
+	       if (coord==NULL) {
+		       fprintf(stderr,"Error: malloc/realloc returned NULL\n");
+		       exit(1);
+	       }
+   
+	       (*coord)[(*numelems)-1] = coordval;
+
+	       /* If desired, "clip" data above a certain value */
+	       if (clip_data && (dataval > clip_val)) {
+	          (*data)[(*numelems)-1] = clip_val;
+	       } else {
+	          (*data)[(*numelems)-1] = dataval;
+	       }
+
+	 }
+
+      }
+
+   }
+   rc = feof(infile);
+
+   return rc;
+}
+
+
+/********************************************************************
+ * Main part of program
+ ********************************************************************/
+int main(int argc, char **argv)
+{
+   char   *infilename;
+   FILE   *infile;
+
+   char   func_name[200];
+   double time;
+   int    numelems_in_set = 0;
+   int    rlev = 0;
+   double *coord, *data;
+   int clip_data = 0;
+   double clip_val;
+   int    i;
+
+
+   /* 
+    * Parse command-line arguments
+    */
+   if (argc <= 1) {
+      fprintf(stderr," input from stdin...  invoke carpet2xgraph -h for help\n");
+   } else if (strcmp(argv[1],"-h")==0) {
+      fprintf(stderr,"usage: cat file | carpet2xgraph [name_tag] [ref_lev] [clip_data_at]\n");
+      exit(1);
+   }
+
+   if (argc >= 2) {
+     infilename = argv[1];
+   } else {
+     infilename = "stdin";
+   }
+   if (argc >= 3) {
+      sscanf(argv[2],"%d",&rlev);
+   } 
+   if (argc >= 4) {
+      sscanf(argv[3],"%lf",&clip_val);
+      clip_data = 1;
+   }
+
+   /*
+    * Open the input file
+    */
+   infile = stdin;
+
+   /*
+   infile = fopen(infilename,"r");
+   if (infile == NULL) {
+      fprintf(stderr,"Error opening file '%s'.\n",infilename);
+      exit(1);
+   }
+   */
+
+
+   time = 0;
+   coord = NULL;
+   data = NULL;
+   sprintf(func_name,"%s_%d",infilename,rlev);
+
+   printf("\"x-label x\n");
+   printf("\"y-label %s\n\n\n",infilename);
+   printf("\"label = %s\n\n\n",infilename);
+
+   /*
+    * Main loop for reading from input file and writing to output file
+    */
+   while ( read_next_set(infile,rlev,&numelems_in_set,&time,&coord,&data,
+			    clip_data, clip_val) == 0 ) {
+      if (numelems_in_set > 0) {
+
+	 /* xgraph format */
+	 printf("\"Time = %g\n",time);
+	 for (i=0; i<numelems_in_set; i++) {
+		 printf("%g  %17.14g\n",coord[i],data[i]);
+	 }
+	 printf("\n\n");
+
+         free(coord);
+         free(data);
+         coord = NULL;
+         data = NULL;
+
+      }
+   }
+
+   return 0;
+}
+