The hopefully final renaming of the FlexIO thorn.

git-svn-id: http://svn.cactuscode.org/arrangements/CactusPUGHIO/IOFlexIO/trunk@2 ebee0441-1374-4afa-a3b5-247f3ba15b9a

25 files changed, 7075 insertions, 0 deletions

diff --git a/README b/README
new file mode 100644
index 0000000..9b79a05
--- /dev/null
+++ b/README
@@ -0,0 +1,20 @@
+Cactus Code Thorn IOFlexIO
+Authors    : John Shalf
+Managed by : Thomas Radke <tradke@aei-potsdam.mpg.de>
+Version    : ...
+CVS info   : $Header$
+--------------------------------------------------------------------------
+
+1. Purpose of the thorn
+
+This thorn provides John Shalf's IEEEIO (FlexIO) library as a thorn to Cactus.
+
+2. Dependencies of the thorn
+
+no dependencies
+
+3. Thorn distribution
+
+This thorn is available to ...
+
+4. Additional information
diff --git a/interface.ccl b/interface.ccl
new file mode 100644
index 0000000..b552ebd
--- /dev/null
+++ b/interface.ccl
@@ -0,0 +1,3 @@
+# Interface definition for thorn IOFlexIO
+
+implements: IOFlexIO
diff --git a/param.ccl b/param.ccl
new file mode 100644
index 0000000..b9f6d66
--- /dev/null
+++ b/param.ccl
@@ -0,0 +1,160 @@
+# Parameter definitions for thorn IOFlexIO
+
+
+#############################################################################
+### declare IOFlexIO parameters
+#############################################################################
+private:
+
+##########################
+# What variables to output
+##########################
+STRING output2D "Variables to output in 2D FlexIO file format"
+{
+  :: "default is to output nothing"
+} ""
+STRING output3D "Variables to output in 3D FlexIO file format"
+{
+  :: "default is to output nothing"
+} ""
+################
+# Various things
+################
+LOGICAL checkpoint_FlexIO "Do checkpointing with IOFlexIO"
+{
+} "no"
+LOGICAL reuse_fh "Reuse file handles by closing files after each write operation."
+{
+} "no"
+
+
+#############################################################################
+### import IOUtil parameters
+#############################################################################
+friend: IO
+
+####################
+# Output directories
+####################
+EXTEND STRING IO_outdir "Name of IO output directory"
+{
+} "."
+EXTEND STRING IO_outdir2D "Name of IO 2D output directory, overrides IO_outdir"
+{
+} "IO_outdir"
+EXTEND STRING IO_outdir3D "Name of IO 3D output directory, overrides IO_outdir"
+{
+} "IO_outdir"
+
+
+########################
+# How often to do output
+########################
+EXTEND INTEGER IO_every "How often to do IO output"
+{
+  -1:* ::
+} -1
+EXTEND INTEGER IO_2Devery "How often to do 2D output, overrides IO_every"
+{
+  -1:* ::		
+} -1
+EXTEND INTEGER IO_3Devery "How often to do 3D output, overrides IO_every"
+{
+  -1:* ::
+} -1
+
+
+################
+# various things
+################
+EXTEND LOGICAL IO_verbose "Give extended screen output in IO?"
+{
+} "no"
+EXTEND LOGICAL IO_datestamp "Write date as attribute to IO 3D output file?"
+{
+} "yes"
+EXTEND LOGICAL IO_parameters "Write parameters to IO 3D output file?"
+{
+} "yes"
+EXTEND LOGICAL IO_structures "Write structures to IO 3D output file?"
+{
+} "yes"
+
+
+#######################
+# Specific to 3D output
+#######################
+EXTEND KEYWORD iomode "Which mode for 3D IO"
+{
+  "proc"    :: "every processor writes its share of data into a separate file"
+  "np"      :: "data is collected and written by every N'th processor into a separate file, where N is specified by the parameter ioproc_every"
+  "onefile" :: "all output is written into a single file by processor 0"
+} "proc"
+
+EXTEND INTEGER ioproc_every "Do IO on every N processors."
+{
+  1:*       :: "Must be a positive integer"
+} 8
+EXTEND LOGICAL onefileperslice "Write one file per time slice, as opposed to all data in one file"
+{
+} "no"
+EXTEND LOGICAL unchunked "Don't write data in chunks."
+{
+} "no"
+
+
+##############################################
+# Downsampling parameters (only for 3D output)
+##############################################
+EXTEND INTEGER downsample_x "Factor by which to downsample output in x direction. Point (0,0,0) is always included."
+{
+  1:*       :: "Must be a positive integer"
+} 1
+EXTEND INTEGER downsample_y "Factor by which to downsample output in y direction. Point (0,0,0) is always included."
+{
+  1:*       :: "Must be a positive integer"
+} 1
+EXTEND INTEGER downsample_z "Factor by which to downsample output in z direction. Point (0,0,0) is always included."
+{
+  1:*       :: "Must be a positive integer"
+} 1
+EXTEND LOGICAL out_single "Output 3D data in single precision ? This parameter is ignored for Cactus compiled with SINGLE_PRECISION."
+{
+} "no"
+
+
+###################################
+# Checkpointing/recovery parameters
+###################################
+EXTEND LOGICAL checkpoint_ID "Checkpoint initial data ?"
+{
+} "no"
+EXTEND LOGICAL checkpoint_keep_all "Keep all checkpoint files ?"
+{
+} "no"
+EXTEND LOGICAL recover "Recover from a checkpoint file ?"
+{
+} "no"
+EXTEND INTEGER checkpoint_every "Checkpoint every x iterations"
+{
+  -1:* :: "negative values disable checkpointing"
+} -1
+EXTEND INTEGER checkpoint_keep "How many checkpoint files to keep"
+{
+  1:* :: "1 overwrites the latest checkpoint file"
+} 1
+EXTEND STRING checkpoint_file "File name for regular checkpoint"
+{
+} "checkpoint"
+EXTEND STRING checkpoint_ID_file "File name for initial data checkpoint"
+{
+} "checkpointID"
+EXTEND STRING recover_file "File name of recovery file"
+{
+} "checkpoint"
+EXTEND STRING checkpoint_dir "Output directory for checkpoint files"
+{
+} "."
+EXTEND STRING recovery_dir "Directory to look for the recovery file"
+{
+} "."
diff --git a/schedule.ccl b/schedule.ccl
new file mode 100644
index 0000000..ddf2d19
--- /dev/null
+++ b/schedule.ccl
@@ -0,0 +1,31 @@
+# Schedule definitions for thorn IOFlexIO
+
+########################################################################
+### register IOFlexIO routines
+########################################################################
+schedule IOFlexIO_Startup at STARTUP
+{
+  LANG:C
+} "IOFlexIO startup routine"
+
+
+########################################################################
+### register checkpointing routines
+########################################################################
+if (checkpoint_FlexIO && checkpoint_ID) {
+  schedule IOFlexIO_InitialDataDumpGH at CCTK_CPINITIAL {
+    LANG:C
+  } "Initial data checkpoint routine"
+}
+
+if (checkpoint_FlexIO && checkpoint_every > 0) {
+  schedule IOFlexIO_ConditionallyDumpGH at CCTK_CHECKPOINT {
+    LANG:C
+  } "Regular checkpoint routine"
+}
+
+if (checkpoint_FlexIO) {
+  schedule IOFlexIO_TerminationDumpGH at CCTK_TERMINATE {
+    LANG:C
+  } "Termination checkpoint routine"
+}
diff --git a/src/DumpGH.c b/src/DumpGH.c
new file mode 100644
index 0000000..fd366ba
--- /dev/null
+++ b/src/DumpGH.c
@@ -0,0 +1,379 @@
+ /*@@
+   @file      DumpGH.c
+   @date      Wed Jun 10 14:13:35 1998
+   @author    Paul Walker
+   @desc 
+   DumpGH dumps an entire grid hierarchy (except for 1D arrays) to a 
+   checkpoint file. This file also contains the different wrappers for
+   DumpGH: in case of initial data, termination or regular checkpointing
+   @enddesc 
+   @version $Id$
+ @@*/
+
+static char *rcsid = "$Id$";
+
+#include <stdio.h>
+#include <string.h>
+#include <stdlib.h>
+
+#include "cctk.h"
+#include "flesh.h"
+#include "declare_parameters.h"
+#include "Groups.h"
+#include "Comm.h"
+#include "GHExtensions.h"
+#include "CactusBase/IOUtil/src/ioGH.h"
+#include "ioFlexGH.h"
+
+
+/* we define the values for cp_type; cp_type will be set by the 
+   the wrappers indicating what type of checkpointiung is going to happen;
+   That way we allow for more flex. treatment, eg. directory/filenames 
+   If you put more than 5 CP_*s here, enlarge the made_cpdir array. */
+
+#define CP_RECOVER_DATA   0
+#define CP_EVOLUTION_DATA 1
+#define CP_INITIAL_DATA   2
+
+/* these have to be consistent with the defaults in param.ccl */
+#define CP_DEFAULT_FNAME   "checkpoint"
+#define CP_DEFAULT_IDFNAME "checkpointID"
+
+/* cp_type indicates the type of checkpoint (cp); default: recovering data 
+   The reason why RECOVER is default is, that RecoverGH doesn't share the scope
+   of the cp_type variable the way the routines in this file do. Fix for this ? */
+static int cp_type=CP_RECOVER_DATA;
+
+void IOFlexIO_DumpGH (cGH *GH, int called_from);
+
+
+
+ /*@@
+   @routine    IOFlexIO_ConditionallyDumpGH
+   @date       Fri Aug 21 14:38:25 1998
+   @author    Gabrielle Allen
+   @desc 
+     registered as CCTK_CHECKPOINT, checks if it's time for 
+     checkpointing, sets the checkpoint type and calls DumpGH 
+   @enddesc 
+   @calls    IOFlexIO_DumpGH
+   @calledby rfrTraverse
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+void IOFlexIO_ConditionallyDumpGH (cGH *GH)
+{
+  DECLARE_PARAMETERS
+
+  if (checkpoint_every > 0 &&
+      GH->iteration % checkpoint_every == 0) {
+    if (IO_verbose) {
+      printf ("------------------------------------------------------------\n");
+      printf ("Dumping periodic checkpoint file at iteration %d\n",
+              GH->iteration);
+    }
+    IOFlexIO_DumpGH (GH, CP_EVOLUTION_DATA);
+  }
+}
+
+
+ /*@@
+   @routine    IOFlexIO_TerminationDumpGH
+   @date       Fri Aug 21 14:40:21 1998
+   @author     Gabrielle Allen
+   @desc 
+     This routine is registered as CCTK_TERMINATE, it checks if a
+     termination signal is raised (that is: global var cactus_terminate
+     set to TERMINATION_RAISED_BRDCAST by ./steppers/TerminationStepper.c)
+     and checkpoints all available grid hierarchies.
+   @enddesc 
+   @calls     IOFlexIO_DumpGH
+   @calledby  rfrTraverse
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+void IOFlexIO_TerminationDumpGH (cGH *GH)
+{
+/*** FIXME: no cactus_terminate flag anymore ? ***/
+#if	0
+#ifdef CACTUSBASE_PUGH
+  if (cactus_terminate == TERMINATION_RAISED_BRDCAST) {
+    IOFlexIO_DumpGH (GH, CP_EVOLUTION_DATA); 
+  }
+#endif /* CACTUSBASE_PUGH */
+#else
+  CCTK_WARN (1, "TerminationDumpGH() not yet implemented !");
+#endif
+}
+
+ /*@@
+   @routine    IOFlexIO_InitialDataDumpGH
+   @date       Fri Aug 21 14:46:28 1998
+   @author     Gerd Lanfermann
+   @desc 
+     This routine dumps the initial data of the GHs, registered 
+     as CCTK_CPINITIAL; sets cp_type to CP_INITIAL_DATA and calls
+     DumpGH; BoxInBox needs special treatment since for now.
+   @enddesc 
+   @calls     IOFlexIO_DumpGH
+   @calledby  rfrTraverse
+   @history 
+ 
+   @endhistory 
+
+@@*/
+
+void IOFlexIO_InitialDataDumpGH (cGH *GH)
+{
+  cGH *helperGH;
+  static int maxlev = 0;
+  int l;
+
+/*** FIXME ***/
+#if	0
+#ifdef THORN_BOXINBOX
+  /* This will go once we reorganize the CINitial calling structure! */
+  /* or have a function to deregister a rfr routine */
+  /* BoxinBox does its own: get maxlevel (once!) and do dumpGH for all GH*/
+  /* when this is called by the highest box */
+  if (Contains("boxinbox","yes")) {
+    if (maxlev==0) while (GetGHbyLevel(0,maxlev+1,0)) maxlev++;
+    if (maxlev==0) return; 
+    
+    if (GH->level==maxlev) 
+      for (l=0;l<=maxlev;l++) {
+	printf("Dumping BoxinBox level %d of %d \n",l,maxlev);
+	IOFlexIO_DumpGH(GetGHbyLevel(0,l,0), CP_INITIAL_DATA);
+      }
+    return;
+  }
+#endif
+#endif
+  IOFlexIO_DumpGH (GH, CP_INITIAL_DATA);
+}
+
+
+void IOFlexIO_IEEEIOStructDump (cGH *GH, IOFile iof)
+{
+  CCTK_INT4 i_temp;
+  CCTK_REAL d_temp;
+  ioGH *ioUtilGH;
+
+  /* Get the handle for IOUtil extensions */
+  ioUtilGH = (ioGH *) GH->extensions [CCTK_GetGHExtensionHandle ("IO")];
+
+  i_temp = GH->iteration;
+  CACTUS_IEEEIO_ERROR (IOwriteAttribute (iof, "GH$iteration", FLEXIO_INT4,
+                       1, &i_temp));
+
+  i_temp = ioUtilGH->ioproc_every;
+  CACTUS_IEEEIO_ERROR (IOwriteAttribute (iof, "GH$ioproc_every", FLEXIO_INT4,
+                       1, &i_temp));
+
+  i_temp = CCTK_GetnProcs (GH);
+  CACTUS_IEEEIO_ERROR (IOwriteAttribute (iof, "GH$nprocs", FLEXIO_INT4,
+                       1, &i_temp));
+
+  d_temp = GH->time;
+  CACTUS_IEEEIO_ERROR (IOwriteAttribute (iof, "GH$time", FLEXIO_REAL,
+                       1, &d_temp));
+}
+
+
+ /*@@
+   @routine    IOFlexIO_DumpGH
+   @date       Fri Aug 21 15:13:13 1998
+   @author     Paul Walker
+   @desc 
+     The heart of checkpointing. Called by the different wrappers, this
+     routines get the appropriate filename by IOUtil_PrepareFilename() and 
+     then dumps away using the dump routines from IO ...
+   @enddesc 
+   @calls IO_DumpGF IO_DumpArray IO_DumpScalar
+   @history 
+   @hdate Oct 4 1998 @hauthor Gabrielle Allen
+   @hdesc Removed code which checked and reset unchunked, assume that this 
+          is done when the variable is first set.
+   @hdate Nov 4 1998 @hauthor Gabrielle Allen
+   @hdesc Do a forced synchronization before checkpointing
+   @hdate Apr 16 1999 @hauthor Thomas Radke
+   @hdesc Removed forced sync before checkpointing (just do a normal sync)
+          Introduced a ring buffer for keeping last <checkpoint_keep> files
+   @endhistory 
+
+@@*/
+
+void IOFlexIO_DumpGH (cGH *GH, int called_from)
+{
+  DECLARE_PARAMETERS
+  char dumpfname [1024], tmpfname [1024];
+  IOFile iof;
+  int index, wrotech;
+  int timelevel, current_timelevel;
+  int forceSync;
+  int old_downsample_x, old_downsample_y, old_downsample_z;
+  int old_out_single;
+  ioGH *ioUtilGH;
+  cTimer total_time, write_time;
+  static char **dumpfnames = NULL;  /* dump filename ring buffer */
+  static int findex = 0;            /* index into ring buffer */
+
+
+  /* allocate the ring buffer for filenames */
+  if (! dumpfnames)
+    dumpfnames = (char **) calloc (checkpoint_keep, sizeof (char *));
+
+  /* initialize timers */
+  CactusResetTimer (&total_time);
+  CactusResetTimer (&write_time);
+
+  CactusStartTimer (&total_time);
+
+  /* Get the handle for IOUtil xtensions */
+  ioUtilGH = (ioGH *) GH->extensions [CCTK_GetGHExtensionHandle ("IO")];
+
+  /* disable downsampling after saving original downsampling params */
+  old_downsample_x = ioUtilGH->downsample_x;
+  old_downsample_y = ioUtilGH->downsample_y;
+  old_downsample_z = ioUtilGH->downsample_z;
+  ioUtilGH->downsample_x = ioUtilGH->downsample_y = ioUtilGH->downsample_z = 1;
+
+  /* disable output in single precision */
+  old_out_single = ioUtilGH->out_single;
+  ioUtilGH->out_single = 0;
+
+  /* sync all groups */
+  for (index = 0; index < CCTK_GetNumGroups (); index++) {
+    char *gname = CCTK_GetGroupName (index);
+
+    CCTK_SyncGroup (GH, gname);
+    free (gname);
+  }
+
+  /* get the base filename ... */
+  IOUtil_PrepareFilename (GH, NULL, dumpfname, called_from,
+                          CCTK_GetMyProc (GH) / ioUtilGH->ioproc_every,
+                          ioUtilGH->unchunked);
+
+  /* ... and append the extension */
+  sprintf (tmpfname,  "%s.chkpt_tmp.ieee", dumpfname);
+  sprintf (dumpfname, "%s.chkpt.ieee",     dumpfname);
+
+  /* Now open the file */
+  if (CCTK_GetMyProc (GH) == ioUtilGH->ioproc) {
+    if (IO_verbose)
+      printf ("Creating file %s\n", tmpfname);
+    iof = IEEEopen (tmpfname, "w");
+    if (! IOisValid (iof)) {
+      CCTK_WARN (1, "Can't open checkpoint file. Checkpointing is skipped");
+      return;
+    }
+  } else
+    iof = (IOFile) NULL;
+
+  /* Great; Now start dumping away! */
+  CactusResetTimer (&write_time);
+  CactusStartTimer (&write_time);
+
+  if (IO_verbose)
+    printf ("Dumping Params ...\n -- ");
+
+  /* first the parameters ... */
+  if (iof) {
+/*** FIXME ***/
+#if	0
+    if (IO_parameters)
+      IOFlexIO_IEEEIOparamDump (iof);
+#endif
+    if (IO_structures)
+      IOFlexIO_IEEEIOStructDump (GH, iof);
+  }
+
+  wrotech = 0;
+ 
+  /* ... now the variables */
+  for (index = 0; index < CCTK_GetNumVars (); index++) {
+
+    /* let only variables pass which have storage assigned */
+    if (! CCTK_QueryGroupStorage_ByIndex (GH, 
+                  CCTK_GetGroupIndexFromVar_ByIndex (index)))
+      continue;
+
+    if (IO_verbose) {
+      char *varname = CCTK_GetVarName (index);
+
+      printf (" %s", varname);
+      wrotech += strlen (varname) + 1;
+      if (wrotech > 65) {
+        printf ("\n    ");
+        wrotech = 0;
+      }
+    }
+
+    /* get the current timelevel */
+    current_timelevel = CCTK_GetNumTimeLevels (index) - 1;
+    if (current_timelevel > 0)
+      current_timelevel--;
+
+    /* now dump all timelevels up to the current */
+    for (timelevel = 0; timelevel <= current_timelevel; timelevel++)
+      IOFlexIO_DumpVar (GH, index, timelevel, iof);
+  } /* end of loop over all variables */
+
+  CactusStopTimer (&write_time);
+
+  CACTUS_IEEEIO_ERROR (IOclose (iof));
+  
+  if (IO_verbose)
+    printf ("\nTime to write: %lf sec (%lf sec per Grid Func)\n",
+/*** FIXME: choose right component of basic[] ***/
+            write_time.total.basic [0],
+            write_time.total.basic [0] / CCTK_GetNumVars ());
+
+#ifndef WIN32
+  /* PW: This should work on WIN32 also, since it
+     is in stdio.h. Someone should check this
+     if/when we go into WIN32 production.
+  */
+  if (CCTK_GetMyProc (GH) == ioUtilGH->ioproc) {
+    if (rename (tmpfname, dumpfname)) {
+      char msg [512];
+
+      sprintf (msg, "Could not rename temporary checkpoint file %s to %s", tmpfname, dumpfname);
+      CCTK_WARN (0, msg);
+      /* should not return here */
+    }
+  }
+#endif
+
+  /* delete the oldest dumpfile if checkpoint_keep_all isn't set
+     and put the new filename into the ring buffer */
+  if (dumpfnames [findex]) {
+    if (! checkpoint_keep_all)
+      remove (dumpfnames [findex]);
+    free (dumpfnames [findex]);
+  }
+  dumpfnames [findex] = strdup (dumpfname);
+  findex = (findex + 1) % checkpoint_keep;
+
+  /* restore output precision flag */
+  ioUtilGH->out_single = old_out_single;
+
+  /* restore original downsampling params */
+  ioUtilGH->downsample_x = old_downsample_x;
+  ioUtilGH->downsample_y = old_downsample_y;
+  ioUtilGH->downsample_z = old_downsample_z;
+
+  CactusStopTimer (&total_time);
+  if (IO_verbose) {
+/*** FIXME: choose right component of basic[] ***/
+    printf ("Time to checkpoint: %lf sec\n", total_time.total.basic [0]);
+    printf ("------------------------------------------------------------\n");
+  }
+}
diff --git a/src/DumpVar.c b/src/DumpVar.c
new file mode 100644
index 0000000..77325b8
--- /dev/null
+++ b/src/DumpVar.c
@@ -0,0 +1,941 @@
+/*@@
+   @file      DumpVar.c
+   @desc 
+   Do the actual writing of a 3D grid function, for output or for checkpointing
+   @enddesc 
+   @history
+   @hauthor Gabrielle Allen @hdate Oct 17 1998
+   @hdesc Split into subroutines and tidied.
+   @hauthor Gabrielle Allen @hdate 19 Oct 1998      
+   @hdesc Changed names ready for thorn_IO
+   @hauthor Thomas Radke @hdate 16 Mar 1999
+   @hdesc Converted to Cactus 4.0
+          introduced separate downsample values for each dimension
+   @hauthor Thomas Radke @hdate 17 Mar 1999
+   @hdesc Changed naming: IEEEIO -> IOFlexIO
+   @hauthor Thomas Radke @hdate 17 Mar 1999
+   @hdesc fixed errors for MPI
+   @hauthor Thomas Radke @hdate 12 Apr 1999
+   @hdesc Store full variable name in IOFlexIO_AddCommonAttributes()
+   @hauthor Thomas Radke @hdate 16 Apr 1999
+   @hdesc Added groupname, grouptype, ntimelevels, and current timelevel
+          as common attributes
+   @hendhistory
+   @version $Id$
+ @@*/
+
+static char *rcsid = "$Id$";
+
+#include <stdio.h>
+#include <stdlib.h>
+#ifdef SGI
+#include <time.h>
+#endif
+
+#include "cctk.h"
+#include "flesh.h"
+#include "Groups.h"
+#include "GroupsOnGH.h"
+#include "Comm.h"
+#include "WarnLevel.h"
+#include "GHExtensions.h"
+#include "declare_parameters.h"
+#ifdef CACTUSBASE_PUGH
+#include "CactusBase/pugh/src/include/pugh.h"
+#endif
+#include "CactusBase/IOUtil/src/ioGH.h"
+#include "ioFlexGH.h"
+
+
+#define IOTAGBASE 20000 /* This may break on more than 2000 processors */
+
+
+static char *char_time_date = NULL;
+
+/* local function prototypes */
+void IOFlexIO_DumpScalar (cGH *GH, int index, int timelevel, IOFile iof,
+                        int ioflex_type);
+void IOFlexIO_DumpGF (cGH *GH, int index, int timelevel, IOFile iof,
+                    int ioflex_type, int mpi_type, int element_size);
+void IOFlexIO_DumpArray (cGH *GH, int index, int timelevel, IOFile iof,
+                       int ioflex_type, int mpi_type, int element_size);
+void IOFlexIO_getDumpData (cGH *GH, int index, int timelevel, void **outme,
+                         int *free_outme, CCTK_INT4 bnd [9], int element_size);
+void IOFlexIO_eachProcDump (cGH *GH, int index, int timelevel, void *outme,
+                          CCTK_INT4 bnd [9], IOFile iof, int ioflex_type);
+void IOFlexIO_procDump (IOFile iof, cGH *GH, int index, void *outme,
+                      CCTK_INT4 bnd [9], int ioflex_type);
+void IOFlexIO_AddChunkAttributes (cGH *GH, int index, CCTK_INT4 chunk_origin [3],
+                                IOFile iof);
+void IOFlexIO_AddCommonAttributes (cGH *GH, int index, int timelevel,
+                                 CCTK_INT4 gsz [3], IOFile iof);
+
+
+
+/*@@
+  @routine    IOFlexIO_DumpVar
+  @date       16 Apr 1999
+  @author     Thomas Radke
+  @desc 
+              Generic dump routine, just calls the appropriate dump routine
+  @enddesc 
+  @calls
+  @calledby   IOFlexIO_DumpGH IOFlexIO_Write3D
+
+  @var        GH
+  @vdesc      Pointer to CCTK grid hierarchy
+  @vtype      cGH
+  @vio        in
+  @endvar
+  @var        index
+  @vdesc      global index of the variable to be dumped
+  @vtype      int
+  @vio        in
+  @endvar
+  @var        timelevel
+  @vdesc      the timelevel to store
+  @vtype      int
+  @vio        in
+  @endvar
+  @var        iof
+  @vdesc      the IEEEIO file to dump to
+  @vtype      IOFile
+  @vio        in
+  @endvar
+  @history 
+  @endhistory 
+
+@@*/
+void IOFlexIO_DumpVar (cGH *GH, int index, int timelevel, IOFile iof)
+{
+  char msg [80];
+  ioGH *ioUtilGH;
+  int variable_type;
+  int ioflex_type, mpi_type = 0, element_size;
+
+  /* Get the handle for IOUtil extensions */
+  ioUtilGH = (ioGH *) GH->extensions [CCTK_GetGHExtensionHandle ("IO")];
+
+  variable_type = CCTK_GetVarVType (index);
+
+  switch (variable_type) {
+    case CCTK_VARIABLE_REAL:
+      ioflex_type = ioUtilGH->out_single ? FLEXIO_REAL4 : FLEXIO_REAL;
+      element_size = ioUtilGH->out_single ?
+                     sizeof (CCTK_REAL4) : sizeof (CCTK_REAL);
+#ifdef MPI
+      mpi_type = ioUtilGH->out_single ? PUGH_MPI_REAL4 : PUGH_MPI_REAL;
+#endif
+      break;
+
+    case CCTK_VARIABLE_INTEGER:
+      ioflex_type = FLEXIO_INT;
+      element_size = sizeof (CCTK_INT);
+#ifdef MPI
+      mpi_type = PUGH_MPI_INT;
+#endif
+      break;
+
+    case CCTK_VARIABLE_CHAR:
+      ioflex_type = FLEXIO_CHAR;
+      element_size = sizeof (CCTK_CHAR);
+#ifdef MPI
+      mpi_type = PUGH_MPI_CHAR;
+#endif
+      break;
+
+    case CCTK_VARIABLE_COMPLEX:
+      CCTK_WARN (1,
+               "Don't know how to dump complex datatype in IOFlexIO_DumpScalar");
+      return;
+
+    default:
+      sprintf (msg, "Unknown variable type %d in IOFlexIO_DumpVar", variable_type);
+      CCTK_WARN (1, msg);
+      return;
+  }
+
+  switch (CCTK_GetVarGType (index)) {
+    case GROUP_SCALAR:
+      IOFlexIO_DumpScalar (GH, index, timelevel, iof, ioflex_type);
+      break;
+
+    case GROUP_ARRAY:
+      IOFlexIO_DumpArray (GH, index, timelevel, iof, ioflex_type, mpi_type,
+                        element_size);
+      break;
+
+    case GROUP_GF:
+      IOFlexIO_DumpGF (GH, index, timelevel, iof, ioflex_type, mpi_type,
+                     element_size);
+      break;
+
+    default:
+      sprintf (msg, "Invalid group type %d in IOFlexIO_DumpVar",
+               CCTK_GetVarGType (index));
+      CCTK_WARN (1, msg);
+      return;
+  }
+}
+
+
+/*@@
+  @routine    IOFlexIO_DumpScalar
+  @date       16 Apr 1999
+  @author     Thomas Radke
+  @desc 
+              Dumps a SCALAR type variable into a IEEEIO file
+  @enddesc 
+  @calls
+  @calledby   IOFlexIO_DumpVar
+
+  @var        GH
+  @vdesc      Pointer to CCTK grid hierarchy
+  @vtype      cGH
+  @vio        in
+  @endvar
+  @var        index
+  @vdesc      global index of the variable to be dumped
+  @vtype      int
+  @vio        in
+  @endvar
+  @var        timelevel
+  @vdesc      the timelevel to store
+  @vtype      int
+  @vio        in
+  @endvar
+  @var        iof
+  @vdesc      the IEEEIO file to dump to
+  @vtype      IOFile
+  @vio        in
+  @endvar
+  @var        ioflex_type
+  @vdesc      the IOFlexIO datatype to store
+  @vtype      int
+  @vio        in
+  @endvar
+  @history 
+  @endhistory 
+
+@@*/
+void IOFlexIO_DumpScalar (cGH *GH, int index, int timelevel, IOFile iof,
+                        int ioflex_type)
+{
+  ioGH *ioUtilGH;
+  CCTK_INT gsz [] = {0, 0, 0};    /* not needed here ? */
+  int dim [] = {1};                   /* size of dimensions for scalar */
+
+
+  /* Get the handle for IOUtil extensions */
+  ioUtilGH = (ioGH *) GH->extensions [CCTK_GetGHExtensionHandle ("IO")];
+
+  if (CCTK_GetMyProc (GH) != ioUtilGH->ioproc)
+    return;
+
+  /* first dump the data then add the attributes */
+  CACTUS_IEEEIO_ERROR (IOwrite (iof, ioflex_type, 1, dim,
+                       CCTK_GetVarDataPtr_ByIndex (GH, timelevel, index)));
+  IOFlexIO_AddCommonAttributes (GH, index, timelevel, gsz, iof);
+}
+
+
+/*@@
+  @routine    IOFlexIO_DumpArray
+  @date       16 Apr 1999
+  @author     Thomas Radke
+  @desc 
+              Dumps a SCALAR type variable into a IEEEIO file
+  @enddesc 
+  @calls
+  @calledby   IOFlexIO_DumpVar
+
+  @var        GH
+  @vdesc      Pointer to CCTK grid hierarchy
+  @vtype      cGH
+  @vio        in
+  @endvar
+  @var        index
+  @vdesc      global index of the variable to be dumped
+  @vtype      int
+  @vio        in
+  @endvar
+  @var        timelevel
+  @vdesc      the timelevel to store
+  @vtype      int
+  @vio        in
+  @endvar
+  @var        iof
+  @vdesc      the IEEEIO file to dump to
+  @vtype      IOFile
+  @vio        in
+  @endvar
+  @var        ioflex_type
+  @vdesc      the IOFlexIO datatype to store
+  @vtype      int
+  @vio        in
+  @endvar
+  @var        mpi_type
+  @vdesc      the MPI datatype to communicate
+  @vtype      int
+  @vio        in
+  @endvar
+  @var        element_size
+  @vdesc      the size of an element of variable
+  @vtype      int
+  @vio        in
+  @endvar
+  @history 
+  @endhistory 
+
+@@*/
+void IOFlexIO_DumpArray (cGH *GH, int index, int timelevel, IOFile iof,
+                       int hdf5io_type, int mpi_type, int element_size)
+{
+  CCTK_WARN (1, "Dumping of arrays not yet supported");
+}
+
+
+/*@@
+  @routine    IOFlexIO_DumpGF
+  @date       July 1998
+  @author     Paul Walker
+  @desc 
+  
+  @enddesc 
+  @calls     
+  @calledby   IOFlexIO_DumpVar
+
+  @var        GH
+  @vdesc      Pointer to CCTK grid hierarchy
+  @vtype      cGH
+  @vio        in
+  @endvar
+  @var        index
+  @vdesc      global index of the variable to be dumped
+  @vtype      int
+  @vio        in
+  @endvar
+  @var        timelevel
+  @vdesc      the timelevel to store
+  @vtype      int
+  @vio        in
+  @endvar
+  @var        iof
+  @vdesc      the IEEEIO file to dump to
+  @vtype      IOFile
+  @vio        in
+  @endvar
+  @var        ioflex_type
+  @vdesc      the IOFlexIO datatype to store
+  @vtype      int
+  @vio        in
+  @endvar
+  @var        mpi_type
+  @vdesc      the MPI datatype to communicate
+  @vtype      int
+  @vio        in
+  @endvar
+  @var        element_size
+  @vdesc      the size of an element of variable
+  @vtype      int
+  @vio        in
+  @endvar
+
+  @history 
+  @hdate Wed Sep  2 10:14:17 1998 @hauthor Tom Goodale
+  @hdesc Merged in Szu-Wen's Panda changes 
+  @hdate 26 Sep 1998
+  @hauthor Gabrielle Allen
+  @hdesc Moved into separate file, tidied up, changed parameter names.
+         Removed zoom_active flag
+  @hauthor    Thomas Radke
+  @hdesc      New parameter timelevel specifies the timelevel to dump
+  @endhistory 
+
+@@*/
+void IOFlexIO_DumpGF (cGH *GH, int index, int timelevel, IOFile iof,
+                    int ioflex_type, int mpi_type, int element_size)
+{
+  DECLARE_PARAMETERS
+  ioGH *ioUtilGH;
+  pGH *pughGH;
+  int myproc;
+  int nprocs;
+  CCTK_INT4 bnd [9]; /* Lower bounds, size and global size of data we send */
+  void *outme;          /* pointer to the data to be dumped */
+  int free_outme;       /* indicates whether data buffer needs freeing */
+#ifdef SGI
+  /* Get the date if we can */
+  time_t t = time (NULL);
+  char_time_date = asctime (localtime (&t));
+#endif
+#ifdef MPI
+  int i;
+  MPI_Status ms;
+#endif
+
+  /* Get the handles for pugh and IO extensions */
+  ioUtilGH = (ioGH *) GH->extensions [CCTK_GetGHExtensionHandle ("IO")];
+  pughGH = (pGH *) GH->extensions [CCTK_GetGHExtensionHandle ("PUGH")];
+
+  myproc = CCTK_GetMyProc (GH);
+  nprocs = CCTK_GetnProcs (GH);
+
+  /* If I'm not vertex centered, bail (for now) */
+  if (pughGH->stagger != PUGH_VERTEXCTR) {
+    CCTK_WARN (1, "3D Output for staggered variables not available");
+    return;
+  }
+
+  /* Get the pointer to the data we want to output (includes downsampling) */
+  IOFlexIO_getDumpData (GH, index, timelevel, &outme, &free_outme, bnd,
+                      element_size);
+
+  if (ioUtilGH->ioproc_every == 1) {
+  
+    /* 
+     * OUTPUT ON EVERY PROCESSOR
+     */
+
+    IOFlexIO_eachProcDump (GH, index, timelevel, outme, bnd, iof, ioflex_type);
+
+  }
+
+#ifdef MPI
+  else if (myproc != ioUtilGH->ioproc) {
+
+    int outgoing;
+
+    outgoing = bnd [3] * bnd [4] * bnd [5];
+      
+    /* 
+     * OUTPUT ON A SUBSET OF PROCESSORS
+     */
+
+    /* Send the geometry and data from the processors which aren't 
+     * outputing to the ones that are 
+     */
+
+    /* GEOMETRY SEND */
+    CACTUS_MPI_ERROR (MPI_Send (bnd, 6, PUGH_MPI_INT4, ioUtilGH->ioproc,
+                      2*myproc+IOTAGBASE+1, pughGH->PUGH_COMM_WORLD));
+    /* DATA SEND */
+    CACTUS_MPI_ERROR (MPI_Send (outme, outgoing, mpi_type, ioUtilGH->ioproc,
+                      2*myproc+IOTAGBASE, pughGH->PUGH_COMM_WORLD));
+    if (IO_verbose) {
+      printf ("Sent %d data points\n", outgoing);
+      fflush (stdout);
+    }
+
+  } else if (myproc == ioUtilGH->ioproc) {
+
+    int lsz [3];
+
+    /* First calculate the local size and reserve the chunk */
+    if (ioUtilGH->ioproc_every >= nprocs) {
+
+      /* One output file ... the local chunk size is the global size */
+      lsz [0] = bnd [6]; lsz [1] = bnd [7]; lsz [2] = bnd [8];
+    }
+
+    /* Dump data held on output processor */
+
+    if (ioUtilGH->unchunked)
+      IOreserveChunk (iof, ioflex_type, 3, lsz);
+
+    /* first the data then the attributes */
+    IOFlexIO_procDump (iof, GH, index, outme, bnd, ioflex_type);
+    IOFlexIO_AddCommonAttributes (GH, index, timelevel, &bnd [6], iof);
+
+    /* Dump data from all other processors */
+
+    for (i = myproc+1; i < myproc+ioUtilGH->ioproc_every && i < nprocs; i++) {
+
+      /* Allocate temp to the right size (based on GH->ub [i] et..) */
+      void *tmpd;
+      int incoming;
+
+      /* Receive bounds information */
+      CACTUS_MPI_ERROR (MPI_Recv (bnd, 6, PUGH_MPI_INT4, i, 2*i+IOTAGBASE+1,
+                        pughGH->PUGH_COMM_WORLD, &ms));
+
+      incoming = bnd [3] * bnd [4] * bnd [5];
+
+      tmpd = malloc (incoming * element_size);
+      CACTUS_MPI_ERROR (MPI_Recv (tmpd, incoming, mpi_type, i,
+                        2*i+IOTAGBASE, pughGH->PUGH_COMM_WORLD, &ms));
+
+      IOFlexIO_procDump (iof, GH, index, tmpd, bnd, ioflex_type);
+      free (tmpd);
+      
+    } /* End loop over processors */
+  } /* End myproc = 0 */
+
+  CACTUS_MPI_ERROR (MPI_Barrier (pughGH->PUGH_COMM_WORLD));
+#endif /* MPI */
+
+  if (free_outme)
+    free (outme);
+}
+
+
+ /*@@
+   @routine    IOFlexIO_AddCommonAttributes
+   @date       July 1998
+   @author     Paul Walker
+   @desc 
+   Add "Common" attributes, these are the GF name, the current date,
+   simulation time, origin, bounding box, gridspacings (both downsampled
+   and evolution), global grid size, number of processors and iteration
+   number. Note that the datestamp should be turned of if you are byte 
+   comparing two output files.
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+   @hdate Wed Sep  2 10:15:22 1998 @hauthor Tom Goodale
+   @hdesc Abstracted WRITE_ATTRIBUTE to merge in Szu-Wen's Panda changes. 
+   @hdate Apr 16 1999 @hauthor Thomas Radke
+   @hdesc Added attributes groupname, grouptype, ntimelevels,
+          and current timelevel to be stored
+   @hdate May 02 1999 @hauthor Thomas Radke
+   @hdesc Made chunked attribute nioprocs common so that it can be found by
+          the recombiner even for unchunked datasets (eg. SCALARs)
+   @hdate May 05 1999 @hauthor Thomas Radke
+   @hdesc Added "unchunked" attribute to distinguish between chunked/unchunked
+          output files
+   @endhistory 
+
+@@*/
+void IOFlexIO_AddCommonAttributes (cGH *GH, int index, int timelevel,
+                                 CCTK_INT4 gsz [3], IOFile iof)
+{
+  DECLARE_PARAMETERS
+  CCTK_REAL d3_to_IO [3]; /* buffer for writing doubles to IEEEIO */
+  CCTK_INT4 i_to_IO;      /* buffer for writing an int to IEEEIO */
+  char *name, *gname;
+  pGH *pughGH;
+  ioGH *ioUtilGH;
+
+  /* Get the handles for pugh and IO extensions */
+  pughGH = (pGH *) GH->extensions [CCTK_GetGHExtensionHandle ("PUGH")];
+  ioUtilGH = (ioGH *) GH->extensions [CCTK_GetGHExtensionHandle ("IO")];
+
+  name = CCTK_GetFullName (index);
+  CACTUS_IEEEIO_ERROR (IOwriteAttribute (iof, "name", BYTE,
+                       strlen (name) + 1, name));
+  free (name);
+
+  gname = CCTK_GetGroupNameFromVar_ByIndex (index);
+  CACTUS_IEEEIO_ERROR (IOwriteAttribute (iof, "groupname", BYTE,
+                       strlen (gname) + 1, gname));
+
+  i_to_IO = CCTK_GetVarGType (index);
+  CACTUS_IEEEIO_ERROR (IOwriteAttribute (iof, "grouptype", FLEXIO_INT4,
+                       1, &i_to_IO));
+
+  i_to_IO = CCTK_GetNumTimeLevels (index);
+  CACTUS_IEEEIO_ERROR (IOwriteAttribute (iof, "ntimelevels", FLEXIO_INT4,
+                       1, &i_to_IO));
+
+  i_to_IO = timelevel;
+  CACTUS_IEEEIO_ERROR (IOwriteAttribute (iof, "timelevel", FLEXIO_INT4,
+                       1, &i_to_IO));
+
+  if (char_time_date && IO_datestamp)
+    CACTUS_IEEEIO_ERROR (IOwriteAttribute (iof, "date", BYTE,
+                         strlen (char_time_date) + 1, char_time_date));
+
+  CACTUS_IEEEIO_ERROR (IOwriteAttribute (iof, "time", FLEXIO_REAL, 1,&GH->time));
+
+  /* FIXME : Origin (d3_to_IO) */
+  d3_to_IO [0] = 0.0;
+  d3_to_IO [1] = 0.0;
+  d3_to_IO [2] = 0.0;
+  CACTUS_IEEEIO_ERROR (IOwriteAttribute (iof, "origin", FLEXIO_REAL,3,d3_to_IO));
+  CACTUS_IEEEIO_ERROR (IOwriteAttribute (iof, "min_ext",FLEXIO_REAL,3,d3_to_IO));
+
+  d3_to_IO [0] = GH->delta_space [0] * ioUtilGH->downsample_x;
+  d3_to_IO [1] = GH->delta_space [1] * ioUtilGH->downsample_y;
+  d3_to_IO [2] = GH->delta_space [2] * ioUtilGH->downsample_z;
+  CACTUS_IEEEIO_ERROR (IOwriteAttribute (iof, "delta", FLEXIO_REAL, 3,d3_to_IO));
+
+  if (ioUtilGH->downsample_x > 1 ||
+      ioUtilGH->downsample_y > 1 ||
+      ioUtilGH->downsample_z > 1) {
+    d3_to_IO [0] = GH->delta_space [0];
+    d3_to_IO [1] = GH->delta_space [1];
+    d3_to_IO [2] = GH->delta_space [2];
+    CACTUS_IEEEIO_ERROR (IOwriteAttribute (iof, "evolution_delta", FLEXIO_REAL,
+                         3, d3_to_IO));
+  }
+
+  /* FIXME : Origin (d3_to_IO) */
+  d3_to_IO [0] = 0.0 + gsz [0] * GH->delta_space [0] * ioUtilGH->downsample_x;
+  d3_to_IO [1] = 0.0 + gsz [1] * GH->delta_space [1] * ioUtilGH->downsample_y;
+  d3_to_IO [2] = 0.0 + gsz [2] * GH->delta_space [2] * ioUtilGH->downsample_z;
+  CACTUS_IEEEIO_ERROR (IOwriteAttribute (iof, "max_ext",FLEXIO_REAL,3,d3_to_IO));
+
+  CACTUS_IEEEIO_ERROR (IOwriteAttribute (iof, "global_size", FLEXIO_INT4,
+                       3, gsz));
+
+  i_to_IO = CCTK_GetnProcs (GH);
+  CACTUS_IEEEIO_ERROR (IOwriteAttribute (iof, "nprocs", FLEXIO_INT4,
+                       1, &i_to_IO));
+
+  i_to_IO = ioUtilGH->ioproc_every;
+  CACTUS_IEEEIO_ERROR (IOwriteAttribute (iof, "ioproc_every", FLEXIO_INT4,
+                       1, &i_to_IO));
+
+  i_to_IO = ioUtilGH->unchunked;
+  CACTUS_IEEEIO_ERROR (IOwriteAttribute (iof, "unchunked", FLEXIO_INT4,
+                       1, &i_to_IO));
+
+  i_to_IO = GH->iteration;
+  CACTUS_IEEEIO_ERROR (IOwriteAttribute (iof, "iteration", FLEXIO_INT4,
+                       1, &i_to_IO));
+}
+
+/*@@
+  @routine IOFlexIO_AddChunkAttributes
+  @author  Paul Walker 
+  @date    Feb 1997
+  @desc
+  This routine adds chunk attributes to a data set. That is,
+  is adds attributes to the chunks of data passed to the
+  io processors. 
+  @enddesc
+  @history 
+  @hauthor Gabrielle Allen
+  @hdate Jul 10 1998
+  @hdesc Added the name of the grid function 
+  @hdate Wed Sep  2 10:08:31 1998 @hauthor Tom Goodale
+  @hdesc Abstracted WRITE_ATTRIBUTE to merge in Szu-Wen's Panda calls.
+  @endhistory 
+  @@*/
+
+void IOFlexIO_AddChunkAttributes (cGH *GH, int index, CCTK_INT4 chunk_origin [3],
+                                IOFile iof)
+{
+  char *name;
+  CCTK_INT4 i_to_IO;
+
+  /* there is nothing to do for a serial run */
+  if (CCTK_GetnProcs (GH) == 1)
+    return;
+
+  CACTUS_IEEEIO_ERROR (IOwriteAttribute (iof, "chunk_origin", FLEXIO_INT4,
+                       3, chunk_origin));
+
+  i_to_IO = GH->iteration;
+  CACTUS_IEEEIO_ERROR (IOwriteAttribute (iof, "chunk_dataset", FLEXIO_INT4,
+                       1, &i_to_IO));
+
+  name = CCTK_GetFullName (index);
+  CACTUS_IEEEIO_ERROR (IOwriteAttribute (iof, "name", BYTE,
+                       strlen (name)+1, name));
+}
+
+
+/*@@
+  @routine IOFlexIO_getDumpData
+  @author  Paul Walker 
+  @date    Feb 1997
+  @desc
+  Bounds and data to be output, takes into account downsampling
+  @enddesc
+  @history 
+  @hauthor Gabrielle Allen @hdate Oct 5 1998
+  @hdesc Made into subroutine
+  @endhistory 
+  @var     GH
+  @vdesc   Identifies grid hierachy
+  @vtype   pGH *
+  @vio     in
+  @vcomment 
+  @endvar 
+  @var     GF
+  @vdesc   Identifies grid function
+  @vtype   pGF *
+  @vio     in
+  @vcomment 
+  @endvar 
+  @var     outme
+  @vdesc   data segment to output
+  @vtype   void **
+  @vio     out
+  @vcomment This is just GF->data if there is no downsampling
+  @endvar 
+  @var     free_outme
+  @vdesc   Specifies whether or not to free the storage associated with outme
+  @vtype   int *
+  @vio     out
+  @vcomment 1: free storage; 0: don't free storage
+  @endvar 
+  @var     bnd
+  @vdesc   bounds, local size and global shape of the output
+  @vtype   CCTK_INT4 [9]
+  @vio     out
+  @vcomment bnd [0..2] lower bounds and bnd[3],bnd[4],bnd[5] 
+            bnd [3..5] local size of data to send
+            bnd [6..8] global shape
+  @endvar 
+  @var     element_size
+  @vdesc   the size of an element of variable
+  @vtype   int
+  @vio     in
+  @endvar
+  @@*/
+
+void IOFlexIO_getDumpData (cGH *GH, int index, int timelevel, void **outme,
+                         int *free_outme, CCTK_INT4 bnd [9], int element_size) 
+{
+  DECLARE_PARAMETERS
+  int i;
+  int myproc;
+  ioGH *ioUtilGH;
+  pGH *pughGH;
+  CCTK_REAL4 *single_ptr;
+  CCTK_REAL *real_ptr;
+  CCTK_CHAR *char_ptr;
+  CCTK_INT *int_ptr;
+  void *data = CCTK_GetVarDataPtr_ByIndex (GH, timelevel, index);
+
+  ioUtilGH = (ioGH *) GH->extensions [CCTK_GetGHExtensionHandle ("IO")];
+  pughGH = (pGH *) GH->extensions [CCTK_GetGHExtensionHandle ("PUGH")];
+
+  myproc = CCTK_GetMyProc (GH);
+
+  if (ioUtilGH->downsample_x == 1 &&
+      ioUtilGH->downsample_y == 1 &&
+      ioUtilGH->downsample_z == 1) {
+
+    if (ioUtilGH->out_single) {
+      single_ptr = (CCTK_REAL4 *) malloc (pughGH->npoints*sizeof (CCTK_REAL4));
+
+      for (i = 0; i < pughGH->npoints; i++)
+        single_ptr [i] = (CCTK_REAL4) ((CCTK_REAL *) data) [i];
+
+      *outme = single_ptr;
+      *free_outme = 1;
+    } else {
+      *outme = data;
+      *free_outme = 0;
+    }
+
+    for (i = 0; i < 3; i++) {
+      bnd [i] = GH->lower_bound[i];     /* the bounds */
+      bnd [i+3] = GH->local_shape[i];   /* the sizes */
+      bnd [i+6] = GH->global_shape[i];  /* the global space */
+    }
+
+  } else {
+
+    int start [3], end [3];
+    int i, j, k, l;
+
+    /* Downsampling code ... */
+    bnd [6] = GH->global_shape[0] / ioUtilGH->downsample_x;
+    if (GH->global_shape[0] % ioUtilGH->downsample_x)
+      bnd [6]++;
+    bnd [7] = GH->global_shape[1] / ioUtilGH->downsample_y;
+    if (GH->global_shape[1] % ioUtilGH->downsample_y)
+      bnd [7]++;
+    bnd [8] = GH->global_shape[2] / ioUtilGH->downsample_z;
+    if (GH->global_shape[2] % ioUtilGH->downsample_z)
+      bnd [8]++;
+
+    if (IO_verbose) 
+      printf ("Downsampled sizes (%d, %d, %d) -> (%d, %d, %d)\n",
+              GH->global_shape[0], GH->global_shape[1], GH->global_shape[2],
+              (int) bnd [6], (int) bnd [7], (int) bnd [8]);
+    
+    /* Now figure out the local downsampling */
+    /* The local starts are the lb modded into the downsample */
+    for (i = 0; i < 3; i++) {
+      int downsample;
+
+      if (i == 0)
+        downsample = ioUtilGH->downsample_x;
+      else if (i == 1)
+        downsample = ioUtilGH->downsample_y;
+      else
+        downsample = ioUtilGH->downsample_z;
+
+      bnd [i]   = GH->lower_bound[i] / downsample;
+      start [i] = bnd [i] * downsample;
+      if (start [i] <
+          GH->lower_bound[i] + pughGH->ownership [PUGH_VERTEXCTR][i][0]) {
+        start [i] += downsample;
+        bnd [i] ++;
+      }
+      end [i]   = ((GH->lower_bound [i] +
+                   pughGH->ownership [PUGH_VERTEXCTR][i][1] - 1) / downsample)
+                   * downsample;
+      bnd [i+3] = (end [i] - start [i]) / downsample + 1;
+    }
+
+    if (IO_verbose) {
+      printf ("Downsample ranges (%d, %d, %d) -> (%d, %d, %d)\n",
+        start [0], start [1], start [2],
+        end [0], end [1], end [2]);
+      printf ("Local size/bound  (%d, %d, %d)  (%d, %d, %d)\n",
+        (int) bnd [3], (int) bnd [4], (int) bnd [5],
+        (int) bnd [0], (int) bnd [1], (int) bnd [2]);
+    }
+
+    /* compute local ranges */
+    for (i = 0; i < 3; i++) {
+      start [i] -= GH->lower_bound [i];
+      end [i] -= GH->lower_bound [i];
+    }
+
+    *outme = malloc (bnd [3] * bnd [4] * bnd [5] * element_size);
+    *free_outme = 1;
+
+    /* I hate it to repeat the loops for each case label
+       but that way produces much more efficient code */
+    l = 0;
+    switch (CCTK_GetVarVType (index)) {
+      case CCTK_VARIABLE_CHAR:
+        char_ptr = (CCTK_CHAR *) *outme;
+        for (k = start [2]; k <= end [2]; k += ioUtilGH->downsample_z)
+          for (j = start [1]; j <= end [1]; j += ioUtilGH->downsample_y)
+            for (i = start [0]; i <= end [0]; i += ioUtilGH->downsample_x)
+              char_ptr [l++] = ((CCTK_CHAR *) data) [DI (pughGH, i, j, k)];
+        break;
+
+      case CCTK_VARIABLE_INTEGER:
+        int_ptr = (CCTK_INT *) *outme;
+        for (k = start [2]; k <= end [2]; k += ioUtilGH->downsample_z)
+          for (j = start [1]; j <= end [1]; j += ioUtilGH->downsample_y)
+            for (i = start [0]; i <= end [0]; i += ioUtilGH->downsample_x)
+              int_ptr [l++] = ((CCTK_INT *) data) [DI (pughGH, i, j, k)];
+        break;
+
+      case CCTK_VARIABLE_REAL:
+        if (ioUtilGH->out_single)
+          single_ptr = (CCTK_REAL4 *) *outme;
+        else
+          real_ptr = (CCTK_REAL *) *outme;
+        for (k = start [2]; k <= end [2]; k += ioUtilGH->downsample_z)
+          for (j = start [1]; j <= end [1]; j += ioUtilGH->downsample_y)
+            for (i = start [0]; i <= end [0]; i += ioUtilGH->downsample_x)
+              if (ioUtilGH->out_single)
+                single_ptr [l++] = (CCTK_REAL4)
+                                  (((CCTK_REAL *) data) [DI (pughGH, i, j, k)]);
+              else
+                real_ptr [l++] = ((CCTK_REAL *) data) [DI (pughGH, i, j, k)];
+        break;
+
+      default:
+        CCTK_WARN (1, "Unsupported variable type in IOFlexIO_getDumpData");
+        return;
+    }
+  }
+
+  if (IO_verbose) {
+    printf ("Global size: %d %d %d\n",
+            (int) bnd [6], (int) bnd [7], (int) bnd [8]);
+    printf ("Lower bound: %d %d %d\n",
+            (int) bnd [0], (int) bnd [1], (int) bnd [2]);
+    printf ("Chunk size : %d %d %d\n",
+            (int) bnd [3], (int) bnd [4], (int) bnd [5]);
+  }
+}
+
+
+/*@@
+  @routine IOFlexIO_eachProcDump
+  @author  Paul Walker 
+  @date    Feb 1997
+  @desc
+  Dump data from each processor
+  @enddesc
+  @history 
+  @hauthor Gabrielle Allen @hdate Oct 5 1998
+  @hdesc Made into subroutine
+  @endhistory 
+  @@*/
+
+void IOFlexIO_eachProcDump (cGH *GH, int index, int timelevel, void *outme,
+                          CCTK_INT4 bnd [9], IOFile iof, int ioflex_type) 
+{
+  int chunk_dims [3];
+  ioGH *ioUtilGH;
+
+  ioUtilGH = (ioGH *) GH->extensions [CCTK_GetGHExtensionHandle ("IO")];
+
+  /* So set up the local shape. */
+  chunk_dims [0] = bnd [3];
+  chunk_dims [1] = bnd [4];
+  chunk_dims [2] = bnd [5];
+  
+  /* Dump the data */
+  CACTUS_IEEEIO_ERROR (IOwrite (iof, ioflex_type, 3, chunk_dims, outme));
+  
+  /* Add attributes for global space */
+  IOFlexIO_AddCommonAttributes (GH, index, timelevel, &bnd [6], iof);
+  
+  /* Add chunk attributes ... lower bound of this chunk and
+   * number of processors participating. 
+   */
+  CACTUS_IEEEIO_ERROR (IOwriteAttribute (iof, "subchunk_lb", FLEXIO_INT4,
+                       3, &bnd [0]));
+  CACTUS_IEEEIO_ERROR (IOwriteAttribute (iof, "global_size", FLEXIO_INT4,
+                       3, &bnd [6]));
+  
+  IOFlexIO_AddChunkAttributes (GH, index, bnd, iof);
+}
+
+
+/*@@
+  @routine IOFlexIO_procDump
+  @author  Paul Walker 
+  @date    Feb 1997
+  @desc
+  Dump data 
+  @enddesc
+  @history 
+  @hauthor Gabrielle Allen @hdate Oct 5 1998
+  @hdesc Made into subroutine
+  @endhistory 
+  @@*/
+
+void IOFlexIO_procDump (IOFile iof, cGH *GH, int index, void *outme,
+                      CCTK_INT4 bnd [9], int ioflex_type)
+{
+  int chunk_dims [3], chunk_origin [3];      /* Chunk dim and origin */
+  ioGH *ioUtilGH;
+
+  ioUtilGH = (ioGH *) GH->extensions [CCTK_GetGHExtensionHandle ("IO")];
+
+  chunk_origin [0] = bnd [0];
+  chunk_origin [1] = bnd [1];
+  chunk_origin [2] = bnd [2];
+  chunk_dims [0] = bnd [3];
+  chunk_dims [1] = bnd [4];
+  chunk_dims [2] = bnd [5];
+
+  if (ioUtilGH->unchunked) {
+    
+    /* Unchunked data */
+    CACTUS_IEEEIO_ERROR (IOwriteChunk (iof, chunk_dims, chunk_origin, outme));
+    
+  } else {
+    
+    /* Chunked data */
+    CACTUS_IEEEIO_ERROR (IOwrite (iof, ioflex_type, 3, chunk_dims, outme));
+    
+    /* Write chunk attributes */
+    IOFlexIO_AddChunkAttributes (GH, index, bnd, iof);
+
+    /* Add chunk attributes ... lower bound of this chunk and
+     * number of processors participating. 
+     */
+    CACTUS_IEEEIO_ERROR (IOwriteAttribute (iof, "subchunk_lb", FLEXIO_INT4,
+                         3, &bnd [0]));
+    CACTUS_IEEEIO_ERROR (IOwriteAttribute (iof, "global_size", FLEXIO_INT4,
+                         3, &bnd [6]));
+  }
+}
diff --git a/src/GHExtension.c b/src/GHExtension.c
new file mode 100644
index 0000000..824a7a2
--- /dev/null
+++ b/src/GHExtension.c
@@ -0,0 +1,91 @@
+ /*@@
+   @file      GHExtension.c
+   @date      Fri May 21 1999
+   @author    Thomas Radke
+   @desc 
+   IOFlexIO GH extension stuff.
+   @enddesc 
+   @history
+   @hauthor Thomas Radke @hdate May 21 1999
+   @hdesc Just copied from thorn IO.
+   @endhistory
+ @@*/
+
+/*#define DEBUG_IO*/
+
+#include <stdlib.h>
+#include <string.h>
+#include <stdio.h>
+
+#include "flesh.h"
+#include "Groups.h"
+#include "Comm.h"
+#include "Misc.h"
+#include "GHExtensions.h"
+#include "declare_parameters.h"
+#ifdef	CACTUSBASE_PUGH
+#include "CactusBase/pugh/src/include/pugh.h"
+#endif
+#include "CactusBase/IOUtil/src/ioGH.h"
+#include "ioFlexGH.h"
+
+
+void *IOFlexIO_SetupGH (tFleshConfig *config, int convergence_level, cGH *GH)
+{
+  int i, numvars;
+  flexioGH *newGH;
+
+  numvars = CCTK_GetNumVars ();
+
+  newGH = (flexioGH *) malloc (sizeof (flexioGH));
+  newGH->IO_2Dnum = (int *) malloc (numvars * sizeof (int));
+  newGH->IO_3Dnum = (int *) malloc (numvars * sizeof (int));
+  newGH->IO_2Dlast = (int *) malloc (numvars * sizeof (int));
+  newGH->IO_3Dlast = (int *) malloc (numvars * sizeof (int));
+
+  newGH->IEEEfile_2D = (IOFile **) malloc (numvars * sizeof (IOFile *));
+  for (i = 0; i < numvars; i++)
+    newGH->IEEEfile_2D [i] = (IOFile *) malloc (3 * sizeof (IOFile));
+
+  newGH->IEEEfname_3D = (char **) malloc (numvars * sizeof (char **));
+
+  for (i = 0; i < numvars; i++)
+    newGH->IEEEfname_3D [i] = (char *) malloc (512 * sizeof (char));
+  newGH->IEEEfile_3D = (IOFile *) malloc (numvars * sizeof (IOFile));
+
+  return (newGH);
+}
+
+int IOFlexIO_InitGH (cGH *GH)
+{
+  DECLARE_PARAMETERS
+  int i;
+  ioGH *ioUtilGH;
+  flexioGH *myGH;
+
+  /* get the handles for IOUtil and IOFlexIO extensions */
+  ioUtilGH = (ioGH *) GH->extensions [CCTK_GetGHExtensionHandle ("IO")];
+  myGH = (flexioGH *) GH->extensions [CCTK_GetGHExtensionHandle ("IOFlexIO")];
+
+  InitIONum (myGH->IO_2Dnum, output2D);
+  InitIONum (myGH->IO_3Dnum, output3D);
+
+  for (i=0; i<CCTK_GetNumVars(); i++)
+    myGH->IO_2Dlast [i] = myGH->IO_3Dlast [i] = -1;
+
+  myGH->reuse_fh = reuse_fh;
+
+  /* Only have reuse for chunked data */
+  if (myGH->reuse_fh && ! ioUtilGH->unchunked) {
+    CCTK_Warn (2, "Cannot reuse handles with unchunked data. "
+                  "Ignoring parameter 'reuse_fh'");
+    myGH->reuse_fh = 0;
+  }
+
+  return (0);
+}
+
+int IOFlexIO_rfrTraverseGH (cGH *GH, int rfrpoint)
+{
+  return 0;
+}
diff --git a/src/IEEEIO/Arch.h b/src/IEEEIO/Arch.h
new file mode 100644
index 0000000..e96147c
--- /dev/null
+++ b/src/IEEEIO/Arch.h
@@ -0,0 +1,72 @@
+#ifndef _MACHDEFS_H_
+#define _MACHDEFS_H_
+
+#ifdef ANSI
+#define PROTO(x) x
+#else
+#define PROTO(x) ()
+#endif
+
+#if defined(SGI) || defined(CM5) || defined (DEC) 
+#define F77NAME(a,b,c) a
+#else
+#ifdef HP
+#define  F77NAME(a,b,c) b
+#else
+#ifdef CRAY
+#define F77NAME(a,b,c) c
+#endif
+#endif
+#endif
+
+#ifndef F77NAME
+#define F77NAME(a,b,c) a
+#endif
+
+
+#ifndef SGI
+typedef long Long8;
+#else /* default SGI behavior.  long longs=8 bytes long=4 bytes*/
+typedef long long Long8;
+#endif
+
+#ifdef _CONFIG_H_
+/* config.h gives us INTEGER4 */
+typedef CCTK_INT4 Long; /* for now Long *MUST* be 32bit integer for PC compatability */
+#else
+#ifdef T3E
+typedef short Int; /* T3E uses 8-byte integers as default */
+typedef short Long; /* this is wierd, but its a T3E... */
+#else
+typedef int Int;  /* every other sane design uses 4-byte ints */
+typedef int Long; /* for now Long *MUST* be 32bit integer for PC compatability */
+#endif
+/* if we run into problems later we'll just need have a separate rev of the file format for PC's */
+#endif
+
+#ifdef WIN32 /* this aint a happenin thang yet... */
+/* #include <sys/types.h> ... bahh!!! Win32 doesn't have this! */
+typedef unsigned int IOFile; /* cast to integer for pointers :( */
+union Integer8 {
+  long l;
+  int i[2];
+  char c[8];
+}; /* what can be said about the byte order though? */
+#else /* its a Mac or Unix box probably */
+#include <sys/types.h> 
+typedef caddr_t IOFile; /* use caddr_t to store pointers */
+#endif
+
+#ifdef HP
+#define NO_RECURSIVE_INLINE
+#endif
+
+#ifndef CONST_BROKEN
+#define CONST const
+#else
+#ifdef CONST
+#undef CONST
+#endif
+#endif
+
+#endif
diff --git a/src/IEEEIO/FlexArrayTmpl.H b/src/IEEEIO/FlexArrayTmpl.H
new file mode 100644
index 0000000..b965e19
--- /dev/null
+++ b/src/IEEEIO/FlexArrayTmpl.H
@@ -0,0 +1,180 @@
+// FlexArrayTmpl.hh
+#ifndef __FLEXARRAY_TMPL_HH_
+#define __FLEXARRAY_TMPL_HH_
+
+#ifdef SGI
+#include <new.h>
+#endif
+
+#include "Arch.h"
+//#define DOUBLEWORD 8 // sgi doubleword == -align64 at a minimum
+#define DOUBLEWORD 1
+/*-------------
+ Template: FlexArray
+ Purpose: Simplifies dynamically allocated array management since it contains
+ array size information and manages array resizing.  However, since speed-critical
+ operations like indexing (the [] operator) are inlined, it adds no speed
+ penalty for the most common speed-critical operations.
+ Method Summary:
+ constructor: The initial "size" is 0 (the size of the array that is accessible
+ using the indexing operator [].  However, the internally allocated size is
+ a minimum of 8 elements.
+ [] : Indexing operator.  Acts just like [] would for a generic array.  Since
+ it is inlined, there is actually no penalty in performance compared to
+ using a generic array.
+ append : Add an element to the end of the array, resizing the array if necessary.
+ The arrays are allocated in blocks so it doesn't need to realloc on every append.
+ For arrays < 1024 elements long, it doubles the array size on each realloc.  When
+ larger than 1024 elements, it adds 1024 elements to the array on every realloc.
+ setSize():   Sets the internal and externally availible array sizes to the
+ specified size; shrinking a larger internal array if necessary.
+ getSize():   Return the current externally availible array size.  The internal
+ size may be larger.	
+---------------*/
+template<class T>
+class FlexArray {
+  int size,maxsize;
+  T *data;
+  // Internal reallocation procedure.
+  void grow(int targetSize=1);
+public:
+	// constructor
+  FlexArray(int sz=0);
+  // copy constructor
+  FlexArray(FlexArray &src);
+  FlexArray &operator=(FlexArray<T> &src);
+  // destructor
+  ~FlexArray();
+  inline int getSize(){return size;}
+  int setSize(int s);
+  int append(T value);
+  int append(T *values,int nvals);
+#ifdef FLEXARRAYDEBUG
+  T &operator[](int index){
+    if(index>=size || index<0) { // range check
+      printf("FlexArray:: Bad index. Arraysize=%u Index=%u\n",
+	     size,index);
+      index=0;
+    }
+    return (data[index]); 
+  }
+#else
+  inline T &operator[](int index){ return (data[index]); }
+#endif
+  // purge? which forces array delete
+  inline void purge(){ delete[] data; data=0; maxsize=size=0; }
+  //inline operator T*() { return data; }  // typecasting operator
+  inline T *getData() { return data; } // sometimes you need direct access to the contained array.
+};
+
+template<class T>
+FlexArray<T>::FlexArray(int sz):size(sz){
+  /*
+  if(sz>DOUBLEWORD) maxsize=sz;
+  else maxsize=DOUBLEWORD;  always allocate in doubleword chunks
+  data = new T[maxsize];
+  */
+  maxsize=sz; // problem here is that virtually guarauntees reallocation
+  if(maxsize>0) data = new T[maxsize];
+  else data = 0;
+}
+
+template<class T>
+FlexArray<T>::~FlexArray(){
+  if(maxsize>0 && data)
+    delete[] data;
+  size=maxsize=0;
+  data=0;
+}
+
+
+template<class T>
+void FlexArray<T>::grow(int targetSize){
+  do{
+    if(!maxsize) maxsize=1;
+    //if(maxsize<1024)
+      maxsize<<=1; // double in size
+      //else
+      // maxsize+=1024; // increase by 1k element block
+  }while(maxsize<targetSize);
+  //printf("FlexArray::grow(%u)->%u\n",targetSize,maxsize);
+  T *newdata=new T[maxsize];
+  // has to do its own copy because realloc is not always compatible with
+  // new/delete on all machines
+  if(data){
+    for(int i=0;i<size;i++)
+      newdata[i]=data[i];
+    delete[] data;
+  }
+  data=newdata;
+}
+
+template<class T>
+FlexArray<T>::FlexArray(FlexArray<T> &src):data(0),size(0),maxsize(0){
+  setSize(src.getSize());
+  //puts("\tFlexArray<T>copy constructor***********");
+  for(int i=0;i<src.getSize();i++)
+    data[i]=src.data[i];
+    //this->data[i]=src->data[i];
+
+    //    (*this)[i]=src[i]; // a full stupid copy!!
+}
+
+template<class T>
+FlexArray<T> &FlexArray<T>::operator=(FlexArray<T> &src){
+  if(this == &src) return *this;
+  setSize(src.getSize());
+  //puts("\t\tFlexArray<T>copy operator");
+  for(int i=0;i<src.getSize();i++)
+    data[i]=src[i]; // a full stupid copy!!
+  return *this;
+}
+
+template<class T>
+int FlexArray<T>::setSize(int s){
+  if(s<size){
+    size=s;
+    return size; // can't allow shrinkage below current size
+  }
+  if(s==maxsize){
+    size=maxsize;
+    return size;
+  }
+  // otherwise if maxsize>s>size>
+  maxsize=s;
+  T *newdata=new T[maxsize];
+  if(data){
+    for(int i=0;i<size && i<maxsize;i++)
+      newdata[i]=data[i];
+    delete[] data;
+  }
+  data=newdata;
+  size=maxsize;
+  return size;
+}
+
+template<class T>
+int FlexArray<T>::append(T value){
+  if(size>=maxsize) // failsafe check
+    grow(size+1); // this array only increases in size
+  data[size++]=value;
+  return size;
+}
+
+template<class T>
+int FlexArray<T>::append(T *values,int nvals){
+  if(size+nvals >= maxsize)
+    grow(size+nvals);//this array only increases in size
+  for(int i=0;i<nvals;i++)
+    data[size++]=values[i];
+  return size;
+}
+
+/*
+  1) Refcount must be fixed
+  2) Need to have more flexibility in how the block allocation
+     is defined (eg. memory allocation blocksize)  Use obstacks?
+     This would require a more global memory pool
+*/
+
+#endif // __FLEXARRAY_TMPL_HH_
diff --git a/src/IEEEIO/IEEEIO.cc b/src/IEEEIO/IEEEIO.cc
new file mode 100644
index 0000000..5bcbbbc
--- /dev/null
+++ b/src/IEEEIO/IEEEIO.cc
@@ -0,0 +1,1363 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include "IEEEIO.hh"
+
+#ifdef WIN32
+// Are we Microsoft VC++ 5.0 or 6.0?
+#if defined(_MSC_VER) && ( (_MSC_VER == 1100) || (_MSC_VER == 1200) ) // yes we are
+#include <fcntl.h>
+#define O_RDONLY _O_RDONLY|_O_BINARY
+#define O_RDWR _O_RDWR|_O_BINARY
+#define O_WRONLY _O_WRONLY|_O_BINARY
+#define O_CREAT _O_CREAT|_O_BINARY
+#define O_TRUNC _O_TRUNC|_O_BINARY
+#define O_APPEND _O_APPEND|_O_BINARY
+#else // not an MSC compiler (use the old defines)
+#define O_RDONLY _O_RDONLY
+#define O_RDWR _O_RDWR
+#define O_WRONLY _O_WRONLY
+#define O_CREAT _O_CREAT
+#define O_TRUNC _O_TRUNC
+#define O_APPEND _O_APPEND
+#endif // not an MSC compiler
+#endif // WIN32
+
+#ifdef FFIO
+#include <ffio.h>
+#define open(x,y,z) ffopens(x,y,z, 0, 0, NULL, "bufa.bufsize=256.num_buffers=4")
+#define close(x) ffclose(x)
+#endif
+
+#ifdef T3E
+int IEEEIO::ffioIsTrue(){
+#ifdef FFIO
+  return 1;
+#else
+  return 0;
+#endif
+}
+#endif
+
+
+// int IEEEIO::nextRecord() (stores state in current_record)
+// also garauntees integrity of datasetnumber and 
+// the current_rec and current_dat
+
+IEEEIO::AttribRef::AttribRef(IEEEIO::AttribRef &src){ 
+  //puts("AttribRef::copy constructor"); 
+  rec=src.rec;
+  offset=src.offset;
+  strcpy(name,src.name);
+}
+
+IEEEIO::AttribRef::AttribRef(){
+  //puts("empty attribref constructor");
+}
+/*
+IEEEIO::AttribRef::~AttribRef(){
+  //puts("delete an attrib ref");
+}
+IEEEIO::RecRef::~RecRef(){
+  //puts("delete a recref");
+}
+*/
+
+IEEEIO::AttribRef &IEEEIO::AttribRef::operator=(IEEEIO::AttribRef &src){
+  if(this != &src) {
+    //puts("\tAttribRef::copy operator");
+    rec=src.rec;
+    offset=src.offset;
+    //name=src.name; // uses the flexarray copy constructor
+    strcpy(name,src.name);
+  }
+  return *this;
+}
+
+IEEEIO::DataRef::DataRef():end_offset(0){
+  //puts("Empty DataRef Constructor!!");
+}
+
+
+IEEEIO::DataRef::DataRef(IEEEIO::DataRef &src){
+  //puts("DataRef:: copy constructor");
+  rec=src.rec;
+  offset=src.offset;
+  end_offset=src.end_offset;
+  annotations = src.annotations;
+  attributes = src.attributes; 
+  //puts("DataRef:: copy constructor DONE");
+}
+/*
+IEEEIO::DataRef::~DataRef(){
+  // puts("dataref destructor");
+}*/
+
+IEEEIO::DataRef &IEEEIO::DataRef::operator=(IEEEIO::DataRef &src){
+  if(this != &src){
+    rec=src.rec;
+    offset=src.offset;
+    end_offset=src.end_offset;
+    annotations = src.annotations;
+    attributes = src.attributes;
+  }
+  return *this;
+}
+
+void IEEEIO::initPosition(){
+  // init
+#ifdef FFIO
+  actual_position = ::ffseek(fid,0,L_INCR);
+#else 
+  actual_position = ::lseek(fid,0,L_INCR);
+#endif
+  virtual_position = actual_position + writebuffercount;
+  // get actual filesize
+  
+#ifdef FFIO
+  file_length = ::ffseek(fid,0,L_XTND);
+  ::ffseek(fid,actual_position,L_SET); // seek back to original pos
+#else
+  file_length = ::lseek(fid,0,L_XTND);
+  ::lseek(fid,actual_position,L_SET); // seek back to original pos
+#endif
+  if(virtual_position > file_length) 
+    file_length = virtual_position;
+  // printf("Init position: a=%ld v=%ld f=%ld\n",
+  //	 actual_position,virtual_position,file_length);
+}
+
+
+//long IEEEIO::getPosition(int immediate) {
+//#ifdef FFIO
+//  if(!writebuffer || immediate)
+//    return ::ffseek(fid,0,L_INCR);
+//  else 
+//    return ::ffseek(fid,0,L_INCR)+writebuffercount;
+//#else
+//  if(!writebuffer || immediate)
+//    return ::lseek(fid,0,L_INCR);
+//  else 
+//    return ::lseek(fid,0,L_INCR)+writebuffercount;
+//#endif
+//}
+
+void IEEEIO::restart(){
+  // this->initPosition();
+  //puts("restart seek");
+  IEEEIO::lseek(FileHdrSize,L_SET); // T3E Kludge
+  //puts("restart");
+  current_rec_offset=current_dat_offset=FileHdrSize; // T3E Kludge
+  current_rec.recordsize=0;
+  current_dat.datasize=0;
+  current_dat.rank=0;
+  datasetnumber=0; // initial
+  hasread=0; // and it hasn't been read yet...
+}
+
+int IEEEIO::nextRecord(){
+  RecordHdr rec;
+  // should have it check the datasetnumber before seeking
+  // this will require consistancy checks for the datasetnumber
+  // elsewhere in the code
+  long src=current_rec_offset;
+  long dst=src+current_rec.recordsize;
+  current_rec_offset=IEEEIO::lseek(dst,L_SET); // seek from start
+  // or use rec.read(rec,fid)
+  // or IEEEIO::RecordHdr::read(rec,fid)
+  if(RecordHdr::read(this,rec) <= 0)
+    return 0;
+  current_rec_offset=dst+RecordHdrSize;
+  current_rec=rec;
+  return 1;
+}
+
+int IEEEIO::nextDataRecord(){
+  hasread=0; // reset the hasread field
+  if(current_rec.recordtype==DataRecord && 
+     datasetnumber!=current_rec.sequenceID){
+    IEEEIO::lseek(current_rec_offset,L_SET); // seek to record start
+  }
+  else {
+    do { // scan through records for the next DataRecord
+      if(!nextRecord()) 
+	return 0; // reached end of file
+    } while(current_rec.recordtype!=DataRecord);
+  }
+  // dat_offset is same as current_rec_offset for the data record
+  // so nextRecord() can be used for scanning for annotations
+  current_dat_offset=current_rec_offset;
+  DataRecordHdr::read(this,current_dat);
+  datasetnumber=current_rec.sequenceID;
+  return 1;
+}
+
+//long IEEEIO::getLength(int immediate){
+//#ifdef FFIO
+//  long currentpos=::ffseek(fid,0,L_INCR); // find currentposition
+//  long auspos = ::ffseek(fid,0,L_XTND); // seek to end to find length
+//#else
+//  long currentpos=::lseek(fid,0,L_INCR); // find currentposition
+//  long auspos = ::lseek(fid,0,L_XTND); // seek to end to find length
+//#endif
+//  if(writebuffer && writebuffercount && !immediate){
+//    // must compute real end of file
+//    register long pos1=currentpos+writebuffercount;
+//    if(pos1>auspos) auspos=pos1;
+//  }
+//#ifdef FFIO
+//  ::ffseek(fid,currentpos,L_SET); // seek back to original position
+//#else
+//  ::lseek(fid,currentpos,L_SET); // seek back to original position
+//#endif
+//  return auspos; // return end position
+//}
+
+void IEEEIO::byteswapBuffer(void *buf,long nelements,int elementsize){
+  char *buffer=(char *)buf; // treat as a character buffer
+  if(elementsize<=1) return;
+  for(long i=0;i<nelements;i++,buffer+=elementsize){
+    register int s,d;
+    register char c;
+    // do the swap thang on each element
+    for(s=0,d=elementsize-1;s<d;s++,d--){
+      c=buffer[s];
+      buffer[s]=buffer[d];
+      buffer[d]=c;
+    }
+  }
+}
+
+void IEEEIO::byteswapBuffer(void *source,void *dest,long nelements,int elementsize){
+  if(elementsize<=1) return;
+  // Lets optimize for integers
+  switch(elementsize){
+  case 4:
+    {
+      Int *src=(Int*)source;
+      Int *dst=(Int*)dest;
+      for(long i=0;i<nelements;i++){
+	union IntChar {
+	  int i;
+	  char a[4];
+	};
+	register IntChar ts,td;
+	// do the swap thang on each element
+	td.i=ts.i=src[i];
+	td.a[0]=ts.a[3];
+	td.a[1]=ts.a[2];
+	td.a[2]=ts.a[1];
+	td.a[3]=ts.a[0];
+	dst[i]=td.i;
+      }
+    }
+    break;
+#ifndef WIN32 // Win32 can't have 8byte integers
+  case 8:
+    {
+      Long *src=(Long*)source;
+      Long *dst=(Long*)dest;
+      for(long i=0;i<nelements;i++){
+	union LongChar {
+	  Long l;
+	  char a[8];
+	};
+	register LongChar ts,td;
+	// do the swap thang on each element
+	td.l=ts.l=src[i];
+	td.a[0]=ts.a[7];
+	td.a[1]=ts.a[6];
+	td.a[2]=ts.a[5];
+	td.a[3]=ts.a[4];
+	td.a[4]=ts.a[3];
+	td.a[5]=ts.a[2];
+	td.a[6]=ts.a[1];
+	td.a[7]=ts.a[0];
+	dst[i]=td.l;
+      }
+    }
+    break;
+#endif
+  default:
+    {
+      char *src=(char *)source; // treat as a character buffer
+      char *dst=(char *)dest;
+      for(long i=0;i<nelements;i++,src+=elementsize,dst+=elementsize){
+	register int s,d;
+	// do the swap thang on each element
+	for(s=0,d=elementsize-1;s<d;s++,d--){
+	  dst[d]=src[s];
+	  dst[s]=src[d];
+	}
+      }
+    }
+    break;
+  }
+}
+
+int IEEEIO::writeFileHeader(){
+  // long pos=getPosition(); // use virtual_position
+  // this->flush();
+  //puts("WriteHeader");
+  for(char c=0;c<4;c++) // make sure magic is correct
+    file_header.magic.c[c]=c+1;
+  if(swapbytes)
+    file_header.byteorder.i=IEEE_REVERSE_MAGIC;
+  else 
+    file_header.byteorder.i=IEEE_NATIVE_MAGIC;
+  file_header.majorversion = IEEE_MAJORVERSION;
+  file_header.minorversion = IEEE_MINORVERSION;
+  IEEEIO::lseek(0,L_SET); // will force a flush (if needed)
+  //puts("done writeheader seek");
+  if(swapbytes) file_header.byteswap();
+  // T3E Kludge
+  int sz=FileHdr::write(this,file_header);
+  if(swapbytes) file_header.byteswap();
+  // if(pos>=FileHdrSize)  // T3E Kludge
+  //  IEEEIO::lseek(pos,L_SET);
+  return sz;
+}
+
+// reads the start of the file (including the magic number)
+// and determines if this is 
+//       1) A valid IEEEIO file
+//       2) what the byte order is
+//       3) How many datasets are contained
+int IEEEIO::readFileHeader(){
+  long pos=getPosition(); // use virtual_position
+  IEEEIO::lseek(0,L_SET);
+  // T3E Kludge
+  int sz=FileHdr::read(this,file_header);
+  IEEEIO::lseek(pos,L_SET); // return to original position
+  for(char c=0;c<4;c++){
+    if((c+1)!=file_header.magic.c[c]){
+      file_header.ndatasets=0;
+      fprintf(stderr,"IEEEIO: File %s has the wrong magic number\n",filename);
+      return -1; // bad magic
+    }
+  }
+  // now compare byte order
+  if(file_header.byteorder.i==IEEE_NATIVE_MAGIC)
+    swapbytes=0;
+  else{
+    swapbytes=1;
+    //puts("byteswapping is On !!!!!!!!!!!!!!!");
+  }
+  if(swapbytes) file_header.byteswap();
+  ndatasets=file_header.ndatasets;
+  return 1;
+}
+
+void IEEEIO::rebuildFileHeader(){
+  int lndatasets;
+  long pos=getPosition(); // save file position (use virtual position
+  restart();
+  // steal seek loop from the nextDataRecord()
+  for(lndatasets=0;
+      nextDataRecord()>0;
+      lndatasets++){
+    // now check the length of the record
+    // for each record.
+  }
+  file_header.ndatasets=lndatasets;
+  writeFileHeader();
+  IEEEIO::lseek(pos,L_SET);
+}
+
+void IEEEIO::appendRecordTable(){
+  //printf("appending records==================\n");
+  while(nextRecord()>0){
+    //printf("appending record\n");
+    switch(current_rec.recordtype){
+    case DataRecord:
+      //fprintf(stderr,"DataRecord");
+      (rec_table[current_rec.sequenceID]).rec=current_rec;
+      // use virtual_position
+      (rec_table[current_rec.sequenceID]).offset=getPosition()-RecordHdrSize; 
+      break;
+    case AnnotationRecord:{
+      //fprintf(stderr,"\tAnnotationRecord\n");
+      RecRef ref;
+      ref.rec=current_rec;
+      ref.offset=getPosition()-	RecordHdrSize; // T3E Kludge
+      (rec_table[current_rec.sequenceID]).annotations.append(ref);
+    }
+    break;
+    case AttributeRecord:{
+      //fprintf(stderr,"\tAttributeRecord\n");
+      AttribRef ref,*refp;
+      AttributeRecordHdr attribhdr;
+      ref.rec=current_rec;
+      ref.offset=getPosition()- RecordHdrSize; // T3E Kludge
+      (rec_table[current_rec.sequenceID]).attributes.append(ref);
+      int idx=(rec_table[current_rec.sequenceID]).attributes.getSize()-1;
+      refp = &((rec_table[current_rec.sequenceID]).attributes[idx]);
+      AttributeRecordHdr::read(this,attribhdr);
+      IEEEIO::read(refp->name,attribhdr.namesize);
+    }
+    break;
+    default:
+      fprintf(stderr,"\tIEEEIO::Error UNKNOWN RECORD TYPE [%d]... recovering.",
+	      (int)(current_rec.recordtype));
+      return;
+    }
+    (rec_table[current_rec.sequenceID]).end_offset = 
+      current_rec_offset + current_rec.recordsize;
+  }
+}
+
+void IEEEIO::buildRecordTable(){
+  if(file_header.ndatasets<0)
+    return; // failure
+  restart();
+  rec_table.setSize(file_header.ndatasets);
+  appendRecordTable();
+}
+
+void IEEEIO::openFile(CONST char *fname,IObase::AccessMode access,int swbytes){
+  switch(access){
+  case IObase::Read:
+    fid=open(fname,O_RDONLY,0644);
+    if(fid<0){
+      fprintf(stderr,"IEEEIO: Failed to open %s for reading\n",fname);
+      return;
+    }
+    initPosition();
+    if(readFileHeader()<=0){
+      close(fid);
+      fprintf(stderr,"IEEEIO: File %s is empty (opened for reading)\n",fname);
+      fid=-1; // invalid file
+      return;
+    }
+    if(file_header.ndatasets<0){
+      // must recover from crash
+      close(fid);
+      fid=open(fname,O_RDWR,0644);
+      rebuildFileHeader();
+      close(fid);
+      fid=open(fname,O_RDONLY,0644);
+      readFileHeader(); // reread..
+    }
+    buildRecordTable();
+    restart(); // go to start of file
+    ndatasets = file_header.ndatasets;
+    nextRecord(); // prime the recordnumber
+    break;
+  case IObase::Write: // truncates
+    fid=open(fname,O_WRONLY|O_CREAT|O_TRUNC,0644);
+    if(fid<0){
+      fprintf(stderr,"IEEEIO: Failed to create %s for writing\n",fname);
+      return;
+    }
+    initPosition();
+    file_header.ndatasets=-1;
+    ndatasets=0;
+    writeFileHeader();
+    restart();
+    break;
+  case IObase::Append:
+  default:
+    fid=open(fname,O_RDWR,0644); // PW: We *don't* want O_APPEND here.
+				// O_Append moves the file ptr to end
+				// of file before a write; this means
+				// when we do the writeHeader() below
+				// this is stuck at the end of this file
+				// which makes it puke.
+    if(fid<0){
+      fprintf(stderr,"IEEEIO: Failed to open %s for append\n",fname);
+      return;
+    }
+    initPosition();
+    if(readFileHeader()<=0){
+      close(fid);
+      fprintf(stderr,"IEEEIO: File %s is empty (opened for append)\n",fname);
+      fid=-1; // invalid file
+      return;
+    }
+    if(file_header.ndatasets<0){
+      // must recover from crash
+      close(fid);
+      fid=open(fname,O_RDWR,0644);
+      initPosition();
+      rebuildFileHeader();
+      close(fid);
+      fid=open(fname,O_RDWR,0644);	// PW: We want to open RDWR again
+      initPosition();
+      readFileHeader(); // reread..
+    }
+    buildRecordTable();
+    restart(); // go to start of file
+    ndatasets = file_header.ndatasets;
+    // now we sync up in the way we would for a write
+    file_header.ndatasets=-1;
+    writeFileHeader(); // write header with -1 to indicate we are writing
+    // could refcount by increasing the - of numbers in header 
+    // (needs extra flag for synced)
+    restart();
+    nextRecord(); // prime the recordnumber
+    //seek(ndatasets);
+    break;
+  }
+}
+
+IEEEIO::IEEEIO(IEEEIO *file): // read only dup
+  IObase(file->filename,IObase::Read),fid(dup(file->fid)),
+  swapbytes(file->swapbytes),datasetnumber(0),
+  ndatasets(file->ndatasets),hasread(0),writebuffer(0),
+  writebuffersize(0),writebuffercount(0),savedposition(-1),
+  file_length(-1),actual_position(-1),virtual_position(-1)
+{
+  if(file->masterfile) masterfile=file->masterfile;
+  else masterfile=file;
+  initPosition(); // initial positional pointers
+  ndatasets = file_header.ndatasets = file->ndatasets;
+  buildRecordTable();
+  restart(); // go to start of file
+  nextRecord(); // prime the recordnumber 
+}
+
+// IEEEIO::AccesMode has scope commented out to satisfy Microsoft compiler
+IEEEIO::IEEEIO(CONST char *fname,/*IEEEIO::*/AccessMode access,int swbytes):
+  IObase(fname,mode(access)),fid(-1),swapbytes(swbytes),datasetnumber(0),
+  ndatasets(0),hasread(0),masterfile(0),writebuffer(0),
+  writebuffersize(0),writebuffercount(0),savedposition(-1),
+  file_length(-1),actual_position(-1),virtual_position(-1)
+{
+  if(access==IEEEIO::SharedRead){
+    masterfile=this; // we are multi-reading
+    fid=open(fname,O_RDONLY,0644);
+    if(fid<0){
+      fprintf(stderr,"IEEEIO: Failed to open %s for reading\n",fname);
+      return;
+    }
+    this->initPosition();
+    if(readFileHeader()<=0){
+      close(fid);
+      fprintf(stderr,"IEEEIO: File %s is empty (opened for reading)\n",fname);
+      fid=-1; // invalid file
+      return;
+    }
+    buildRecordTable();
+    restart(); // go to start of file
+    ndatasets = file_header.ndatasets;
+    nextRecord(); // prime the recordnumber
+  }
+  else
+    openFile(fname,mode(access),swbytes);
+}
+
+IEEEIO::IEEEIO(CONST char *fname,IObase::AccessMode access,int swbytes):
+  IObase(fname,access),fid(-1),swapbytes(swbytes),datasetnumber(0),
+  ndatasets(0),hasread(0),masterfile(0),writebuffer(0),
+  writebuffersize(0),writebuffercount(0),savedposition(-1)
+{
+  //  long fpos;
+  openFile(fname,access,swbytes);
+}
+
+IEEEIO::~IEEEIO(){
+  if(fid<0 && savedposition>=0) resume();
+  // resume IO only if it is paused so that
+  // we can do the final writes to the file before
+  // closing
+  //puts("do bufferoff");
+  if(writebuffer) bufferOff(); // automatically flushes buffer
+  //puts("now rewrite header");
+  if(fid>=0){
+    if(accessmode!=IObase::Read){
+      file_header.ndatasets=ndatasets;
+      writeFileHeader();
+    }
+    close(fid);
+  }
+  fid=-1;
+}
+
+int IEEEIO::isValid(){
+  if(fid>=0) return 1;
+  else return 0;
+}
+
+int IEEEIO::write(IObase::DataType typeID,int rank,CONST int *dims,void *data){
+  int i;
+  RecordHdr rec;
+  DataRecordHdr hdr;
+  // make sure its the EOF
+  if(accessmode==IObase::Read) return 0;
+  hasread=0; // reset hasread; (JMS: changed from local Thu Mar 12 13:53:13 CST 1998)
+  if(rec_table.getSize()>0){
+    long endpos = rec_table[rec_table.getSize()-1].end_offset;
+    IEEEIO::lseek(endpos,L_SET);  // don't know if this is costly
+  }
+  else 
+    IEEEIO::lseek(0,L_XTND); // seek to end of file
+  // if 0-length seeks are costly, then alternative
+  // logic can be constructed to ensure writes to end of file
+  for(i=0,hdr.datasize=1;i<rank;i++) 
+    hdr.datasize*=dims[i];
+
+  if(chunkdims.getSize()<rank)
+    chunkdims.setSize(rank);
+  for(i=0;i<rank;i++) chunkdims[i]=dims[i];
+  
+  hdr.datasize*=sizeOf(typeID);
+  hdr.numbertype=typeID;
+  hdr.rank=rank;
+  // if last annotation slot is filled, it is a pointer
+  // to another block of 8 annotation pointers.
+  rec.recordtype = DataRecord;
+  rec.recordsize = hdr.datasize + 
+    DataRecordHdrSize +
+    sizeof(Int) * rank;
+  rec.sequenceID = datasetnumber = ndatasets++;
+  // need to byteswap the header and records if swapbytes==True
+  current_dat=hdr; // first copy native info to current
+  current_rec=rec;
+  RecordHdr::write(this,rec);
+  current_dat_offset=current_rec_offset=getPosition();
+  DataRecordHdr::write(this,hdr);
+  // write the dims.... (byteswap if necessary)
+  if(swapbytes) byteswapBuffer(chunkdims.getData(),rank,sizeof(Int));
+  IEEEIO::write(chunkdims.getData(),sizeof(Int)*rank);
+  if(swapbytes) byteswapBuffer(chunkdims.getData(),rank,sizeof(Int)); // swap back
+  if(swapbytes) byteswapBuffer(data,hdr.datasize/sizeOf(typeID),sizeOf(typeID));
+  int sz = IEEEIO::write(data,hdr.datasize);
+  // yeah! double bytswap seem stupid, but consider the alternatives
+  // write one byte at a time (swapping as we go) == systemcall overhead
+  // copy to a temporary buffer and swap bytes == malloc, copy, and free!
+  // strangely its faster to swap twice
+  if(swapbytes) byteswapBuffer(data,hdr.datasize/sizeOf(typeID),sizeOf(typeID));
+  DataRef dref;
+  dref.rec=current_rec;
+  // Must find exact size!!!
+  dref.offset = current_rec_offset - RecordHdrSize;
+  dref.end_offset = current_rec_offset + current_rec.recordsize;
+  if(accessmode!=IObase::Write) // do not store if in write-only mode
+    rec_table.append(dref);
+  return sz;
+}
+
+//setCurrentRec(offset);
+// storing in 8*1024 byte blocks (always aligned to chunk boundaries)
+// will result in more efficient seeking behavior due to 
+// disk block size... thus file is always aligned to unix
+// filesystem blocks. (for now, we'll use die method maximo-stupido...
+int IEEEIO::readInfo(IObase::DataType &typeID,int &rank,int *dims,int maxdims){
+  if(accessmode!=IObase::Read && accessmode!=IObase::Append)
+    return 0;
+  //int sz;
+  if(hasread) datasetnumber++; // increment to the next one if this has been read
+  if(datasetnumber>=rec_table.getSize()) return 0; // end of file
+  // read the record + header
+  hasread=1; // we've read this, so next time we hit it, we'll increment again
+  IEEEIO::lseek(rec_table[datasetnumber].offset,L_SET);
+  // T3E Kludge
+  RecordHdr::read(this,current_rec);
+  DataRecordHdr::read(this,current_dat);
+  rank=current_dat.rank;
+  if(chunkdims.getSize()<rank) chunkdims.setSize(rank);
+  typeID=Int2DataType(current_dat.numbertype);
+  IEEEIO::read(chunkdims.getData(),sizeof(Int)*rank);
+  if(swapbytes) byteswapBuffer(chunkdims.getData(),(rank>maxdims)?maxdims:rank,sizeof(Int));
+  
+  for(int i=0;i<maxdims && i<rank;i++) dims[i]=chunkdims[i];
+  streaming=0;
+  return 1;
+}
+
+/*
+  It does not appear that the bytswapping is really configured here...
+ */
+int IEEEIO::read(void *data){
+  if(accessmode!=IObase::Read && accessmode!=IObase::Append)
+    return 0; // can't readif write only
+  if(datasetnumber>=rec_table.getSize()) return 0; // seek past end.
+  long datapos=rec_table[datasetnumber].offset+RecordHdrSize+DataRecordHdrSize+
+    sizeof(Int)*current_dat.rank;
+  // should make certain current file position is correct
+  // if(getPosition() != datapos) (redundant call to lseek())
+  IEEEIO::lseek(datapos,L_SET); // seek to position
+  int sz=IEEEIO::read(data,current_dat.datasize); // read nbytes
+  int typelen = sizeOf(Int2DataType(current_dat.numbertype));// compute elem sz
+  if(swapbytes) byteswapBuffer(data,current_dat.datasize/typelen,typelen);
+  return sz;
+}
+
+int IEEEIO::seek(int idx){ 
+  if(accessmode!=IObase::Read &&
+     accessmode!=IObase::Append) // can't seek unless readable
+    return -1; // failed to seek
+  if (idx >= ndatasets || idx < 0) {
+    return -1;	// do a bounds check (probably should just clip)
+  }
+  // bound the index
+  if((1+idx)>rec_table.getSize()) idx=rec_table.getSize();
+  if(idx<0) idx=0;
+  index = idx;
+  IEEEIO::lseek(rec_table[index].offset,L_SET);
+  // T3E Kludge
+  RecordHdr::read(this,current_rec);
+  current_rec_offset=getPosition();
+  datasetnumber=current_rec.sequenceID; // Changed:  had -1
+  hasread=0; // new file position
+  return current_rec.sequenceID;
+}
+
+int IEEEIO::nDatasets() {
+  // can work across threads or across processes since
+  // it gathers information directly from the file.
+  // it does not need to share info with masterfile
+  // in fact you can delete the masterfile
+  if(masterfile){ // we are passive read-only. update nDatasets
+    // seek-scan to find datasets?
+    int oldlength = (rec_table.getSize())?(rec_table[rec_table.getSize()-1].end_offset):0;
+    //printf("oldlength=%u newlength=%u\n",oldlength,(unsigned int)(getLength()));
+    if(getLength()>oldlength){
+      // lets scan to find the end
+      int lndatasets;
+      puts("We should clearly not be here!!!");
+      restart(); // go to beginning
+      // steal seek loop from the nextDataRecord()
+      for(lndatasets=0;
+	  nextDataRecord()>0;
+	  lndatasets++){
+	//printf("scan dataset[%u]\n",lndatasets);
+      }
+      file_header.ndatasets = lndatasets;
+      rec_table.setSize(lndatasets);
+      // printf("counted %u datasets. old was %u datasets\n",lndatasets,ndatasets);
+      seek(ndatasets); // seek to current last
+      ndatasets = lndatasets; // then reset the last
+      appendRecordTable(); // append new stuff to record table
+    }
+  }
+  return ndatasets; 
+}
+
+int IEEEIO::writeAnnotation(CONST char *annotation){
+  int stringlen;
+  RecRef aref;
+  RecordHdr rec;
+  if(accessmode==IObase::Read || !annotation) return 0;
+  stringlen=strlen(annotation)+1;
+
+  if(rec_table.getSize()>0){
+    long endpos = rec_table[rec_table.getSize()-1].end_offset;
+    IEEEIO::lseek(endpos,L_SET);  // don't know if this is costly
+  }
+  else 
+    IEEEIO::lseek(0,L_XTND); // seek to end of file
+
+  rec.recordtype=AnnotationRecord;
+  rec.recordsize=stringlen;
+  if(datasetnumber>=0) rec.sequenceID=datasetnumber;
+  else rec.sequenceID=current_rec.sequenceID; // a kludge for error immunity
+  current_rec = rec;
+  // T3E Kludge
+  current_rec_offset = getPosition() + RecordHdrSize;
+  // T3E Kludge
+  RecordHdr::write(this,rec);
+  // if(swapbytes) byteswap(rec);
+  int sz=IEEEIO::write(annotation,stringlen);
+
+  aref.rec=current_rec;
+  // T3E Kludge
+  aref.offset=current_rec_offset-RecordHdrSize;
+  if(accessmode!=IObase::Write){ // don't write to index cache if in write-only mode
+    rec_table[datasetnumber].annotations.append(aref); 
+    rec_table[datasetnumber].end_offset=getPosition();
+  }
+  return sz;
+}
+/*
+	If objects could have all attributes implicitly availible for querying
+	when passed to subroutines.  Basic object methods for data movement like
+	linearize() or object.linearized[index] object.linearized.size.
+	Then for objects that a receiver can't deal with either
+		1: have compiletime "gatekeepers" that restrict allowable types
+			that it can receive
+		2: Do a treesearch through object space to find a convertor sequence
+			that can convert the current type into the one that the object
+			accepts.  If no such object can be found, indicate runtime
+			object adaptor failure.
+*/
+int IEEEIO::readAnnotationInfo(int number,int &length){
+  if(datasetnumber<0 ||
+     number>=rec_table[datasetnumber].annotations.getSize()){ 
+    length=0;
+    return -1;
+  }
+  length=rec_table[datasetnumber].annotations[number].rec.recordsize;
+  return length;
+}
+
+int IEEEIO::readAnnotation(int number,char *annotation,int maxsize){
+  // returns actual size of annotation or -1 if error
+  if(datasetnumber<0 || accessmode==IObase::Write ||
+     number>=rec_table[datasetnumber].annotations.getSize()){ 
+    return -1;
+  }
+  IEEEIO::lseek(rec_table[datasetnumber].annotations[number].offset,L_SET);
+  // T3E Kludge
+  RecordHdr::read(this,current_rec);
+  IEEEIO::read(annotation,(maxsize<rec_table[datasetnumber].annotations[number].rec.recordsize)?
+       maxsize:(rec_table[datasetnumber].annotations[number].rec.recordsize));
+  annotation[maxsize-1]='\0';
+  return rec_table[datasetnumber].annotations[number].rec.recordsize;
+}
+
+int IEEEIO::nAnnotations(){
+  if(datasetnumber<0 || accessmode==IObase::Write) 
+    return -1;
+  return rec_table[datasetnumber].annotations.getSize();
+}
+
+// for attributes
+int IEEEIO::writeAttribute(CONST char *name,IObase::DataType typeID,Long length,void *data){
+  int stringlen;
+  AttribRef aref;
+  RecordHdr rec;
+  AttributeRecordHdr attrib;
+
+  if(datasetnumber<0){
+    fprintf(stderr,"IEEEIO::writeAttribute():  Error, cannot write attribute before any datasets have been written!\n");
+    return 0;
+  }
+  if(accessmode==IObase::Read) return 0;
+  stringlen=strlen(name)+1;
+
+  if(rec_table.getSize()>0){
+    long endpos = rec_table[rec_table.getSize()-1].end_offset;
+    IEEEIO::lseek(endpos,L_SET);  // don't know if this is costly
+  }
+  else 
+    IEEEIO::lseek(0,L_XTND); // seek to end of file
+  attrib.datasize=length*sizeOf(typeID);
+  attrib.namesize=stringlen;
+  attrib.numbertype=typeID;
+  rec.recordtype=AttributeRecord;
+  rec.recordsize=attrib.datasize+attrib.namesize+AttributeRecordHdrSize;
+  if(datasetnumber>=0) rec.sequenceID=datasetnumber;
+  else rec.sequenceID=current_rec.sequenceID; // a kludge for error immunity
+  current_rec=rec;
+  current_rec_offset=getPosition() + RecordHdrSize;
+  aref.rec=current_rec;
+  aref.offset=current_rec_offset-RecordHdrSize;
+  // no copy constructor for aref, so must do manually
+  if(accessmode != IObase::Write){
+    int lastindex;
+    rec_table[datasetnumber].attributes.append(aref);
+    lastindex=rec_table[datasetnumber].attributes.getSize()-1;
+    strcpy(rec_table[datasetnumber].attributes[lastindex].name,name);
+  }
+  // T3E Kludge
+  RecordHdr::write(this,rec);
+  AttributeRecordHdr::write(this,attrib);
+  //printf("\tnow write the string\n");
+  IEEEIO::write(name,stringlen);
+  //printf("\tdone IEEEIO::stringlen=%u typeid=%u\n",attrib.namesize,attrib.numbertype);
+  // data is not byte-reversed...
+  int sz=0;
+  if(swapbytes) byteswapBuffer(data,length,sizeOf(typeID));
+  IEEEIO::write(data,length*sizeOf(typeID));
+  if(swapbytes) byteswapBuffer(data,length,sizeOf(typeID)); // swap back
+  if(accessmode != IObase::Write)
+    rec_table[datasetnumber].end_offset=getPosition();
+  // doesn't appear to store attributes properly
+  return sz;
+}
+
+int IEEEIO::readAttributeInfo(int number,char *name,IObase::DataType &typeID,
+			      Long &nelem,int maxnamelen){
+  if(accessmode==IObase::Write) return -1;
+  FlexArray<AttribRef> *attribs= &(rec_table[datasetnumber].attributes);
+  if(number>=(*attribs).getSize()) return -1; // > number of attributes
+  AttribRef *attrib=&((*attribs)[number]);
+  if(strlen(attrib->name)>maxnamelen){
+    strncpy(name,attrib->name,maxnamelen); // don't we want to copy the other way?
+    name[maxnamelen-1]='\0';
+  }
+  else strcpy(name,attrib->name);
+  name[maxnamelen-1]='\0'; // make certain it is null capped
+  AttributeRecordHdr attribhdr;
+  IEEEIO::lseek((*attribs)[number].offset,L_SET);
+  RecordHdr::read(this,current_rec);
+  AttributeRecordHdr::read(this,attribhdr);
+//printf("IEEEIO:attribrechdr = ds,nt,ns %d,%d,%d\n",
+//	 attribhdr.datasize,
+//	 attribhdr.numbertype,
+//	 attribhdr.namesize);
+  typeID = Int2DataType(attribhdr.numbertype);
+  nelem = attribhdr.datasize/sizeOf(typeID);
+  return -1;
+}
+
+int IEEEIO::readAttributeInfo(CONST char *name,IObase::DataType &typeID,Long &nelem){
+  // by name
+  if(accessmode==IObase::Write) return -1;
+  FlexArray<AttribRef> *attribs= &(rec_table[datasetnumber].attributes);
+  for(int i=0;i<attribs->getSize();i++){
+    if(!strcmp(name,(*attribs)[i].name)){
+      // must read the record to get this info
+      AttributeRecordHdr attribhdr;
+      IEEEIO::lseek((*attribs)[i].offset,L_SET);
+      RecordHdr::read(this,current_rec);
+      AttributeRecordHdr::read(this,attribhdr);
+      typeID = Int2DataType(attribhdr.numbertype);
+      nelem = attribhdr.datasize/sizeOf(typeID);
+      return i;
+    }
+  }
+  return -1; // Attribute not found
+}
+
+int IEEEIO::readAttribute(int number,void *data){
+  if(accessmode==IObase::Write) return -1;
+  FlexArray<AttribRef> *attribs= &(rec_table[datasetnumber].attributes);
+  if(number>=(*attribs).getSize()) return -1; // > number of attributes
+  AttributeRecordHdr attribhdr;
+  IEEEIO::lseek((*attribs)[number].offset,L_SET);
+  RecordHdr::read(this,current_rec);
+  AttributeRecordHdr::read(this,attribhdr);
+  IObase::DataType typeID = Int2DataType(attribhdr.numbertype);
+  if(typeID==IObase::Error) return -1; // read failed due to bad datatype
+  long nelem = attribhdr.datasize/sizeOf(typeID);
+  IEEEIO::lseek(attribhdr.namesize,L_INCR);
+  int sz=IEEEIO::read(data,attribhdr.datasize)/sizeOf(typeID);
+  if(swapbytes) byteswapBuffer(data,sz,sizeOf(typeID));
+  if(typeID==IObase::String){
+    char *cdata=(char *)data;
+    cdata[sz]='\0'; /* Null Terminate String data */
+  }
+  return sz; // returns number of elements read
+}
+
+int IEEEIO::nAttributes(){
+  if(datasetnumber<0 || accessmode == IObase::Write) return -1;
+  return rec_table[datasetnumber].attributes.getSize();
+}
+
+//================Chunking Interface-----------------------
+void  IEEEIO::clearChunk(int nbytes){
+#ifdef T3E
+#define DISKBLOCKSIZE 128*1024
+#else
+#define DISKBLOCKSIZE 8192
+#endif
+// This is a wasteful feature that will be eliminated in later releases
+// This is currently done because I needed the disk space to be zero'ed 
+// prior to writing for debugging purposes...  Could have a writeStream 
+// class.
+  int nwritten;
+  char dummy[DISKBLOCKSIZE];
+  for(int i=0;i<DISKBLOCKSIZE;i++) 
+    dummy[i]=0; // bzero it (reserve space with zero'ed data)
+  while(nbytes>DISKBLOCKSIZE){
+    nwritten=IEEEIO::write(dummy,DISKBLOCKSIZE);
+    if(nwritten<=0) return; // this aint working
+    nbytes-=nwritten;
+  }
+  if(nbytes>0) IEEEIO::write(dummy,nbytes); // write the remaining bytes
+#undef DISKBLOCKSIZE
+}
+
+int IEEEIO::reserveChunk(IObase::DataType typeID,int rank,CONST int *dims){
+  reserveStream(typeID,rank,dims); // sets up the chunking state machine
+  // clearing does not appear to be needed... just need header
+  //clearChunk(IObase::nBytes(typeID,rank,dims)); // this pre-clears the reserved area to 0. (inefficient)
+  return 1;
+}
+
+// Chunking should probably include striding....
+// But thats a pain since you'll need to read + write to
+// Accomplish that (yuck!).  Lots of seeking.
+// Lets leave striding out for now.
+
+// should also check for contiguous data (eg. slicing) and optimize for that layout
+// for now, the streaming interface best serves contiguous data.
+int IEEEIO::writeChunk(CONST int *dims,CONST int *origin,void *data){
+  int i,sz;
+  // make sure its the EOF
+  if(accessmode==IObase::Read) {
+    fprintf(stderr,"IEEEIO::writeChunk(): Error!  Access is ReadOnly\n");
+    return 0;
+  }
+  // should have a "chunk reserved" flag set.
+  if(current_reserved_chunk!=datasetnumber){
+    fprintf(stderr,"IEEEIO::writeChunk(): Error! You forgot to reserve space for the chunk using IO::reserveChunk()\n");
+    return 0;
+  }
+  
+  // now we need to seek to the position of the data
+  int rank = current_dat.rank;
+  IObase::DataType typeID = IObase::Int2DataType(current_dat.numbertype);
+  int typesize=sizeOf(typeID);
+  // long basefileoffset = rec_table[datasetnumber].offset+
+  long basefileoffset = stream_offset +
+    RecordHdrSize+DataRecordHdrSize+sizeof(Int)*current_dat.rank; // T3E Kludge
+  long chunkcolumnsize=dims[0]*typesize; // stride between columns in chunk
+  long filecolumnsize=chunkdims[0]*typesize; // stride between columns on disk (full array size)
+  long originoffset;
+  long accumdims;
+  for(originoffset=0,accumdims=typesize,i=0;i<rank;i++){
+    originoffset += (origin[i]*accumdims);
+    accumdims *= chunkdims[i];
+  }
+  for(i=0;i<rank;i++){
+    if((origin[i]+dims[i])>chunkdims[i]){
+      fprintf(stderr,"IEEEIO::writeChunk(): ERROR!!  specified dims and origin exceed reserved block size\n");
+      fprintf(stderr,"\tfor dimension %u origin=%u and dims=%u\n",i,origin[i],
+	      dims[i]);
+      fprintf(stderr,"\torigin+dims=%u whereas the maximum must be less than %u\n",origin[i]+dims[i],chunkdims[i]);
+    }
+  }
+  originoffset+=basefileoffset; // for absolute seeking
+  long maxindex,minindex;
+  minindex=basefileoffset;
+  maxindex=basefileoffset+current_dat.datasize;
+  int ncolumns; // computed in loop below
+  for(ncolumns=1,i=1;i<rank;i++) ncolumns*=dims[i];
+  if(swapbytes) 
+    byteswapBuffer(data,IObase::nElements(rank,dims),typesize);
+  for(i=0;i<ncolumns;i++){ // read the columns
+    if(originoffset<minindex){
+      fprintf(stderr,"WriteChunk() inconsistency. Writing less than min index\n");
+      fprintf(stderr,"\tCol[%u]: Requested %u, but minidex= %u\n",
+	      i,(unsigned int)originoffset,(unsigned int)minindex);
+    }
+    if(originoffset>maxindex){
+      fprintf(stderr,"WriteChunk() inconsistency. Writing greater than maximum index\n");
+      fprintf(stderr,"\tCol[%u]: Requested %u, but maxindex= %u\n",
+	      i,(unsigned int)originoffset,(unsigned int)maxindex);
+    }
+    if((originoffset+chunkcolumnsize)>maxindex){
+      fprintf(stderr,"WriteChunk() inconsistency. This write will overrun the reserved data\n");
+      fprintf(stderr,"\tCol[%u]: Requested %u, maxindex %u, and the write of %u will run to %u\n",
+	      i,(unsigned int)originoffset,(unsigned int)maxindex,(unsigned int)chunkcolumnsize,(unsigned int)(originoffset+chunkcolumnsize));
+    }
+    IEEEIO::lseek(originoffset,L_SET);
+    sz=IEEEIO::write(((char*)data)+i*chunkcolumnsize,(int)chunkcolumnsize);
+    originoffset+=filecolumnsize;
+    for(long j=1,idx=dims[1],planesize=filecolumnsize;
+	j<(rank-1) && !((i+1)%idx);
+	idx*=dims[++j]){
+      long extraoffset=planesize*(chunkdims.getData()[j]-dims[j]);
+      originoffset+=extraoffset;
+      planesize*=dims[j];
+    }
+  }
+  // now swap the data back to native order...
+  if(swapbytes) 
+    byteswapBuffer(data,IObase::nElements(rank,dims),typesize);
+
+  return sz;
+}
+int IEEEIO::readChunk(CONST int *dims,CONST int *origin,void *data){ 
+  int sz,i;
+  if(accessmode!=IObase::Read && accessmode!=IObase::Append)
+    return 0;
+  // gonna have to stride through this sucker... (yuck)
+  // use the same chunkdims for this.  Must set during readinfo...
+  // now we need to seek to the position of the data
+  int rank = current_dat.rank;
+  IObase::DataType typeID = IObase::Int2DataType(current_dat.numbertype);
+  int typesize=sizeOf(typeID);
+  long basefileoffset = rec_table[datasetnumber].offset+
+    // sizeof(RecordHdr)+sizeof(DataRecordHdr)+sizeof(Int)*current_dat.rank;
+    RecordHdrSize+DataRecordHdrSize+sizeof(Int)*current_dat.rank; // T3E Kludge
+  long chunkcolumnsize=dims[0]*typesize; // stride between columns in chunk
+  long filecolumnsize=chunkdims[0]*typesize; // stride between columns on disk (full array size)
+  long originoffset;
+  long accumdims;
+  // compute the offset into the data on disk required by the chunk origin (initial offset)
+  for(originoffset=0,accumdims=typesize,i=0;i<rank;i++){
+    originoffset += (origin[i]*accumdims);
+    accumdims *= chunkdims[i];
+  }
+  originoffset+=basefileoffset; // for absolute seeking
+  long ncolumns; // computed in loop below
+  for(ncolumns=1,i=1;i<rank;i++) ncolumns*=dims[i];
+  for(i=0;i<ncolumns;i++){ // read the columns
+    IEEEIO::lseek(originoffset,L_SET);
+    sz=IEEEIO::read(((char*)data)+i*chunkcolumnsize,(int)chunkcolumnsize); 
+    originoffset+=filecolumnsize;
+    for(long j=1,idx=dims[1],planesize=filecolumnsize;
+	j<(rank-1) && !((i+1)%idx);
+	idx*=dims[++j]){
+      long extraoffset=planesize*(chunkdims.getData()[j]-dims[j]);
+      originoffset+=extraoffset;
+      planesize*=dims[j];
+    }
+  }
+  // now swap the data back to native order...
+  if(swapbytes)
+    byteswapBuffer(data,IObase::nElements(rank,dims),typesize);
+  return sz;
+}
+// Nearly identical to reserveChunk except that we don't need to pre-clear
+// the space to 0 for streamed data.
+int IEEEIO::reserveStream(IObase::DataType typeID,int rank,CONST int *dims){
+  int i;
+  RecordHdr rec;
+  DataRecordHdr hdr;
+  // make sure its the EOF
+  if(accessmode==IObase::Read) return 0;
+  
+  if(rec_table.getSize()>0){
+    long endpos = rec_table[rec_table.getSize()-1].end_offset;
+    IEEEIO::lseek(endpos,L_SET);  // don't know if this is costly
+  }
+  else {
+    IEEEIO::lseek(0,L_XTND); // seek to end of file
+  }
+  // if 0-length seeks are costly, then alternative
+  // logic can be constructed to ensure writes to end of file
+  if(chunkdims.getSize() != rank)
+    chunkdims.setSize(rank); // KLUDGE!
+  for(i=0,hdr.datasize=1;i<rank;i++) {
+    hdr.datasize*=dims[i];
+    chunkdims[i]=dims[i];
+  }
+  cur_type_size=sizeOf(typeID); // KLUDGE: for writeStream()
+  hdr.datasize*=sizeOf(typeID);
+  hdr.numbertype=typeID;
+  hdr.rank=rank;
+  // if last annotation slot is filled, it is a pointer
+  // to another block of 8 annotation pointers.
+  rec.recordtype = DataRecord;
+  rec.recordsize = hdr.datasize + 
+    DataRecordHdrSize +  // T3E Kludge
+    sizeof(Int) * rank;
+  rec.sequenceID = datasetnumber = ndatasets++;
+  // need to byteswap the header and records if swapbytes==True
+  current_dat=hdr; // first copy native info to current
+  current_rec=rec;
+  RecordHdr::write(this,rec);
+  current_dat_offset=current_rec_offset=getPosition();
+  DataRecordHdr::write(this,hdr);
+  if(swapbytes) byteswapBuffer(chunkdims.getData(),rank,sizeof(Int));
+  IEEEIO::write(chunkdims.getData(),sizeof(Int)*rank); // for better for worse...
+  if(swapbytes) byteswapBuffer(chunkdims.getData(),rank,sizeof(Int)); // swap back
+  // OK, now writing 8k at a time, reserve the space with 0's
+  current_reserved_chunk=datasetnumber;
+  DataRef dref;
+  dref.rec=current_rec;
+  dref.offset=current_rec_offset-RecordHdrSize;
+  dref.end_offset=current_rec_offset+current_rec.recordsize;
+  if(accessmode!=IObase::Write)
+    rec_table.append(dref);
+  // else (JMS debugging change 5/12/98)
+  stream_offset = current_rec_offset-RecordHdrSize;
+  streaming=0;
+  return 1;
+}
+
+int IEEEIO::writeStream(void *data,int length){
+  int len;
+  if(!streaming){// seek to starting offset
+    long basefileoffset = stream_offset +
+      RecordHdrSize+DataRecordHdrSize+sizeof(Int)*current_dat.rank;
+    IEEEIO::lseek(basefileoffset,L_SET);
+    streaming=1;
+  }
+  // need to do a bounds check on the stream... for now it'll play dumb...
+  long nbytes=cur_type_size * length;
+  len = IEEEIO::write((char *)data,nbytes); // record file position??
+  if(len<0)
+    return -1; // write failure
+  return len;
+}
+
+int IEEEIO::readStream(void *data,int length){
+  int len;
+  if(!streaming){// seek to starting offset
+    long basefileoffset = rec_table[datasetnumber].offset+
+      RecordHdrSize+DataRecordHdrSize+sizeof(Int)*current_dat.rank;
+    IEEEIO::lseek(basefileoffset,L_SET);
+    streaming=1;
+  }
+  // need to do a bounds check on the stream... for now it'll play dumb...
+  int typesize= sizeOf(Int2DataType(current_dat.numbertype));
+  long nbytes = typesize * length;
+  len = IEEEIO::read((char *)data,nbytes);
+  if(len<0)
+    return -1; // read failure
+  else
+    if(swapbytes) byteswapBuffer(data,length,typesize);
+  return len;
+}
+
+void IEEEIO::bufferOn(long size){
+  if(writebuffer) {
+    if(size==0)
+      bufferOff();
+    return;
+  }
+  else {
+#if !defined(WIN32) && !defined(T3E)// if not Windows, use fstat() to choose optimal blocksiz
+    if(size<0){
+      struct stat mystat;
+      fstat(fid,&mystat);
+      size=mystat.st_blksize;
+    }
+#endif 
+    // if windows then take a wild guess
+#ifdef T3E
+    size=(size<=64)?8*1024*1024:size; // 8Meg is a good guess for T3E
+#else
+    size=(size<=64)?8*1024:size; //(8k is a good guess for workstations/WinTel)
+#endif
+    writebuffer = new char[size];
+  }
+  writebuffersize=size;
+  writebuffercount=0;
+}
+
+void IEEEIO::bufferOff(){
+  if(!writebuffer) return;
+  if(fid<0 && writebuffercount>0){
+    fprintf(stderr,"IEEEIO::bufferOff() ERROR!!!!");
+    fprintf(stderr,"\tFile %s is either paused or invalid, so I can't deallocate the write buffer safely\n",filename);
+    return;
+  }
+  //printf("bufferoff on a=%ld v=%ld f=%ld b=%d\n",
+  //	 actual_position,virtual_position,file_length,writebuffercount);
+  this->flush();
+  delete writebuffer;
+  writebuffer=0;
+}
+
+int IEEEIO::pause(){
+  if(fid<0) return 0; // fail
+  this->flush(); // flush the buffer to pause
+  savedposition=getPosition();
+  close(fid);
+  fid=-1;
+  return 1; // success
+}
+
+int IEEEIO::resume(){
+  if(fid>=0 || savedposition<0) return 0; // fail
+  switch(IObase::accessmode){
+  case IObase::Read:
+    fid=open(filename,O_RDONLY,0644);
+    break;
+  case IObase::Write:
+    fid=open(filename,O_WRONLY,0644);
+    break;
+  case IObase::Append:
+    fid=open(filename,O_RDWR,0644);
+    break;
+  }
+  if(fid<0){
+    printf("IEEEIO::resume failed for file %s\n",filename);
+    return -1;
+  }
+  ::lseek(fid,savedposition,L_SET);  // this is OK..
+  savedposition=-1;
+  return 1; // success
+}
+
+//*********************Fortran Interface**********************
+Long8 f_ieee_open (char *file,char *accessname,int flen,int alen){
+  // would have used tolower(), but it doesn't exist everywhere.... :(
+  IObase::AccessMode mode;
+  if(*accessname=='R' || *accessname=='r')
+    mode=IObase::Read;
+  else if(*accessname=='W' || *accessname=='w' ||
+	  *accessname=='C' || *accessname=='c')
+    mode=IObase::Write;
+  else if(*accessname=='A' || *accessname=='a')
+    mode=IObase::Append;
+  else {
+    fprintf(stderr,"IEEEopen(): Error unknown option [%s] to open file %s\n",
+	    accessname,file);
+    return 0;
+  }
+  IObase *fid=new IEEEIO(file,mode);
+  if(fid->isValid()) 
+    return (Long8)fid;
+  else
+    delete fid; // file open failed
+  return 0;
+}
+
+Long8 f_ieee_openr(char *filename,int namelen){
+  filename[namelen]='\0';
+  return (Long8)(new IEEEIO(filename,IObase::Read));
+}
+
+Long8 f_ieee_openw(char *filename,int namelen){
+  filename[namelen]='\0';
+  return (Long8)(new IEEEIO(filename,IObase::Write));
+}
+
+Long8 f_ieee_opena(char *filename,int namelen){
+  filename[namelen]='\0';
+  return (Long8)(new IEEEIO(filename,IObase::Append));
+}
+
+void ieee_bufon(Long8 *fileID,int bufsize){
+  IEEEIO *f=(IEEEIO*)(*fileID);
+  if(bufsize<0) f->bufferOn(); // use default size
+  else f->bufferOn(bufsize); // or specify
+}
+
+void ieee_bufoff(Long8 *fileID){
+  IEEEIO *f=(IEEEIO*)(*fileID);
+  f->bufferOff();
+}
+
+IOFile IEEEopen(char *filename,char *accessname){
+  // Parse all of the ansi stdio access option strings
+  IObase::AccessMode mode;
+  if(!strcmp(accessname,"read") ||
+     !strcmp(accessname,"r") ||
+     !strcmp(accessname,"rb"))
+    mode=IObase::Read;
+  else if(*accessname=='a')
+    mode=IObase::Append;
+  else if(!strcmp(accessname,"write") || 
+	  !strcmp(accessname,"create") ||
+	  !strcmp(accessname,"w") ||
+	  !strcmp(accessname,"wb") ||
+	  !strcmp(accessname,"w+") ||
+	  !strcmp(accessname,"w+b") ||
+	  !strcmp(accessname,"wb+"))
+    mode=IObase::Write;
+  else{
+    fprintf(stderr,"IEEEopen(): Error unknown option [%s] to open file %s\n",
+	    accessname,filename);
+    return 0;
+  }
+  IObase *fid=new IEEEIO(filename,mode);
+  if(fid->isValid()) 
+    return (IOFile)fid;
+  else
+    delete fid; // file open failed
+  return 0; // unknown option
+}
+
+IOFile IEEEopenRead(char *filename){
+  return (IOFile)(new IEEEIO(filename,IObase::Read));
+}
+
+IOFile IEEEopenWrite(char *filename){
+  return (IOFile)(new IEEEIO(filename,IObase::Create));
+}
+
+IOFile IEEEopenAppend(char *filename){
+  return (IOFile)(new IEEEIO(filename,IObase::Append));
+}
+
+void IEEEbufferOn(IOFile fileID,int bufsize){
+  IEEEIO *f=(IEEEIO*)fileID;
+  if(bufsize<0) f->bufferOn(); // use default size
+  else f->bufferOn(bufsize); // or specify
+}
+
+void IEEEbufferOff(IOFile fileID){
+  IEEEIO *f=(IEEEIO*)fileID;
+  f->bufferOff();
+}
+
diff --git a/src/IEEEIO/IEEEIO.h b/src/IEEEIO/IEEEIO.h
new file mode 100644
index 0000000..e17da17
--- /dev/null
+++ b/src/IEEEIO/IEEEIO.h
@@ -0,0 +1,13 @@
+#ifndef __IEEEIO_H_
+#define __IEEEIO_H_
+
+#include "IO.h"
+#include "Arch.h"
+
+IOFile IEEEopen PROTO((char *filename,char *accessname));
+IOFile IEEEopenRead PROTO((char *filename));
+IOFile IEEEopenWrite PROTO((char *filename));
+IOFile IEEEopenAppend PROTO((char *filename));
+void IEEEbufferOn PROTO((IOFile fileID,int bufsize));
+void IEEEbufferOff PROTO((IOFile fileID));
+#endif
diff --git a/src/IEEEIO/IEEEIO.hh b/src/IEEEIO/IEEEIO.hh
new file mode 100644
index 0000000..4e79647
--- /dev/null
+++ b/src/IEEEIO/IEEEIO.hh
@@ -0,0 +1,598 @@
+#ifndef __IEEEIO_HH_
+#define __IEEEIO_HH_
+#ifndef WIN32
+#include <stdio.h>
+#include <unistd.h>
+#include <sys/file.h>
+#else
+#include <io.h>
+#endif
+
+#include "IO.hh"
+#include "FlexArrayTmpl.H"
+#include "Arch.h"
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <fcntl.h>
+#include <string.h>
+#if defined(WIN32) || defined(SOLARIS)
+// I'm a WIN32 box and I ignore standards
+// or I'm a Sun box and I hate BSD.
+#define bcopy(src,dest,len) memcpy(dest,src,len)
+#else
+#include <strings.h> // available non-ANSI/Posix compliant 
+#endif
+
+#ifdef	FFIO
+#include <ffio.h>
+#endif
+
+// Win32/Solaris lseek crap...
+#if defined(WIN32) || defined(SOLARIS)
+#define L_INCR SEEK_CUR
+#ifndef L_SET // for that wanky Solaris compiler
+#define L_SET  SEEK_SET
+#endif
+#define L_XTND SEEK_END
+#endif
+// T3E lseek crap... 
+#ifdef T3E // I'm a T3E and I ignore ANSI/Posix standards
+#define L_INCR SEEK_CUR
+#define L_SET  SEEK_SET
+#define L_XTND SEEK_END
+#endif
+
+#define IEEE_NATIVE_MAGIC 0x01020304
+#define IEEE_REVERSE_MAGIC 0x04030201
+#define IEEE_ALTERNATING_MAGIC 0x02010403
+
+#define IEEE_BIG_ENDIAN IEEE_NATIVE_MAGIC
+#define IEEE_LITTLE_ENDIAN IEEE_REVERSE_MAGIC
+#define IEEE_VAX_ENDIAN IEEE_ALTERNATING_MAGIC
+
+#define IEEE_MAJORVERSION 1
+#define IEEE_MINORVERSION 0
+
+class IEEEIO : public IObase {
+public:
+  enum RecordType {DataRecord=1,AnnotationRecord,AttributeRecord}; 
+  union btorder {
+    char c[4]; // write as Int, read as Int
+    Int  i;
+  };
+#define FileHdrSize (4+4+4+2+2) // T3E Kludge
+  struct FileHdr {
+    // do magic twice so we can figure out which byte-order
+    // the current machine is in addition (without hardcoding) in
+    // addition to which byte order the file is in
+    btorder magic;
+    btorder byteorder;
+    Long ndatasets;
+    // T3E will be stupid for this, but thats OK
+    short majorversion;
+    short minorversion;
+    // ** Inline these because they are 
+    // ** actually kludges for the T3E
+    static int write(IEEEIO *ieeeio,FileHdr &rec){
+      int r=0;
+#ifdef T3E // data structures are not packed contiguously in memory
+      r+=ieeeio->write((rec.magic.c),4);
+      r+=ieeeio->write((rec.byteorder.c),4);
+      r+=ieeeio->write(&(rec.ndatasets),4);
+      // this may be incorrect, but it will work for now
+      // It assumes the low order bytes will be first to
+      // be written.
+      r+=ieeeio->write(&(rec.majorversion),2);
+      r+=ieeeio->write(&(rec.minorversion),2);
+#else
+      // problem if this includes pointers to static methods
+      // will require an extra degree of inheritance. (inherit
+      // from a base structure.
+      r=ieeeio->write(&rec,FileHdrSize);
+#endif
+      return r;
+    }
+    static int read(IEEEIO *ieeeio,FileHdr &rec){
+      int r=0;
+#ifdef T3E // data structures are not packed contiguously in memory
+      char buffer[FileHdrSize];
+      r=ieeeio->read(buffer,FileHdrSize); // read full buffer in and then copy to struct
+      ::bcopy(buffer,(rec.magic.c),4);
+      ::bcopy(buffer+4,(rec.byteorder.c),4);
+      ::bcopy(buffer+8,&(rec.ndatasets),4);
+      // May be invalid if low order bytes are not the
+      // proper ones to load. (eg. may fail on big-endian systems)
+      ::bcopy(buffer+12,&(rec.majorversion),2);
+      ::bcopy(buffer+14,&(rec.minorversion),2);
+#else
+      r=ieeeio->read(&rec,FileHdrSize);
+#endif
+      return r;
+    }
+    static void byteswap(FileHdr &hdr){
+      IEEEIO::byteswapBuffer(&(hdr.ndatasets),1,sizeof(hdr.ndatasets));
+      // because the t3e doesn't have a 2-byte datatype, must do this manually
+      IEEEIO::byteswapBuffer(&(hdr.majorversion),1,2);
+      IEEEIO::byteswapBuffer(&(hdr.minorversion),1,2);
+    }
+    // May be a problem to have a nonstatic class member here because
+    // this will create a vtable entry.
+    void byteswap(){
+      IEEEIO::byteswapBuffer(&ndatasets,1,sizeof(ndatasets));
+      // because the t3e doesn't have a 2-byte datatype, must do this manually
+      IEEEIO::byteswapBuffer(&majorversion,1,2);
+      IEEEIO::byteswapBuffer(&minorversion,1,2);
+    }
+  };
+#define RecordHdrSize (4+4+4) // T3E Kludge
+  struct RecordHdr {
+    Int recordtype;
+    Long recordsize;
+    Int sequenceID;
+    // ** Inline these because they are 
+    // ** actually kludges for the T3E
+    static int write(IEEEIO *ieeeio,RecordHdr &rec){
+      int r=0;
+      // copy into buffer
+      if(ieeeio->swapbytes){
+	char buffer[RecordHdrSize];
+	IEEEIO::byteswapBuffer(&(rec.recordtype),buffer,1,4);
+	IEEEIO::byteswapBuffer(&(rec.recordsize),buffer+4,1,4);
+	IEEEIO::byteswapBuffer(&(rec.sequenceID),buffer+8,1,4);
+	r=ieeeio->write(buffer,RecordHdrSize);
+      }
+      else {
+#ifdef T3E // data structures are not packed contiguously in memory
+	r+=ieeeio->write(&(rec.recordtype),4);
+	r+=ieeeio->write(&(rec.recordsize),4);
+	r+=ieeeio->write(&(rec.sequenceID),4);
+#else
+	r=ieeeio->write(&rec,RecordHdrSize);
+	//printf("Wrote rechdr %u,%u,%u\n",rec.recordtype,rec.recordsize,rec.sequenceID);
+#endif
+      }
+      return r;
+    }
+    static int read(IEEEIO *ieeeio,RecordHdr &rec){
+      char buffer[RecordHdrSize];
+      int r=ieeeio->read(buffer,RecordHdrSize);
+      if(ieeeio->swapbytes)
+	IEEEIO::byteswapBuffer(buffer,3,4);//3 elements of 4 bytes 
+#ifdef T3E // data structures are not packed contiguously in memory
+      // do it the stupid way first 
+      // would be smarter to pack buffers really
+      ::bcopy(buffer,&(rec.recordtype),4);
+      ::bcopy(buffer+4,&(rec.recordsize),4);
+      ::bcopy(buffer+8,&(rec.sequenceID),4);
+#else
+      ::bcopy(buffer,(char*)&rec,RecordHdrSize);
+#endif
+      return r;
+    }
+  };
+  
+#define DataRecordHdrSize (4+4+4) // T3E Kludge   
+  struct DataRecordHdr {
+    Long datasize;
+    Int numbertype;
+    Int rank;
+    // ** Inline these because they are 
+    // ** actually kludges for the T3E
+    static int write(IEEEIO *ieeeio,DataRecordHdr &rec){
+      int r=0;
+      // copy into buffer
+      if(ieeeio->swapbytes){
+	char buffer[DataRecordHdrSize];
+	IEEEIO::byteswapBuffer(&(rec.datasize),buffer,1,4);
+	IEEEIO::byteswapBuffer(&(rec.numbertype),buffer+4,1,4);
+	IEEEIO::byteswapBuffer(&(rec.rank),buffer+8,1,4);
+	r=ieeeio->write(buffer,DataRecordHdrSize);
+      }
+      else {
+#ifdef T3E
+	r+=ieeeio->write(&(rec.datasize),4);
+	r+=ieeeio->write(&(rec.numbertype),4);
+	r+=ieeeio->write(&(rec.rank),4);
+#else
+	r=ieeeio->write(&rec,DataRecordHdrSize);
+#endif
+      }
+      return r;
+    }
+    static int read(IEEEIO *ieeeio,DataRecordHdr &rec){	
+      char buffer[DataRecordHdrSize];
+      int r=ieeeio->read(buffer,DataRecordHdrSize);
+      if(ieeeio->swapbytes)
+	IEEEIO::byteswapBuffer(buffer,3,4);//3 elements of 4 bytes 
+#ifdef T3E
+      // do it the stupid way first 
+      // would be smarter to pack buffers really	
+      ::bcopy(buffer,&(rec.datasize),4);
+      ::bcopy(buffer+4,&(rec.numbertype),4);
+      ::bcopy(buffer+8,&(rec.rank),4);
+#else
+      ::bcopy(buffer,(char*)&rec,DataRecordHdrSize);
+#endif
+      return r;
+    }
+  };
+#define AttributeRecordHdrSize (4+4+4) // T3E Kludge  
+  struct AttributeRecordHdr {
+    Long datasize;
+    Int numbertype;
+    Int namesize; // attributes are named
+    // ** Inline these because they are 
+    // ** actually kludges for the T3E
+    static int write(IEEEIO *ieeeio,AttributeRecordHdr &rec){	
+      int r=0;
+      // copy into buffer
+      if(ieeeio->swapbytes){ // pack a buffer
+	char buffer[AttributeRecordHdrSize];
+	IEEEIO::byteswapBuffer(&(rec.datasize),buffer,1,4);
+	IEEEIO::byteswapBuffer(&(rec.numbertype),buffer+4,1,4);
+	IEEEIO::byteswapBuffer(&(rec.namesize),buffer+8,1,4);
+	r=ieeeio->write(buffer,RecordHdrSize);
+      }
+      else {
+#ifdef T3E
+	r+=ieeeio->write(&(rec.datasize),4);
+	r+=ieeeio->write(&(rec.numbertype),4);
+	r+=ieeeio->write(&(rec.namesize),4);
+#else
+	//printf("AttribRec::writing ds,nt,ns %d,%d,%d\n",rec.datasize,rec.numbertype,rec.namesize);
+	r = ieeeio->write(&rec,AttributeRecordHdrSize);
+	//printf("Wrote attribhdr %d,%d,%d\n",rec.datasize,rec.numbertype,rec.namesize);
+	//printf("\tr=%u\n",r);
+#endif
+      }
+      return r;
+    }
+    static int read(IEEEIO *ieeeio,AttributeRecordHdr &rec){
+      int r=0;
+      char buffer[AttributeRecordHdrSize];
+      r=ieeeio->read(buffer,AttributeRecordHdrSize);
+      if(ieeeio->swapbytes)
+	IEEEIO::byteswapBuffer(buffer,3,4);
+#ifdef T3E
+      ::bcopy(buffer,&(rec.datasize),4);
+      ::bcopy(buffer+4,&(rec.numbertype),4);
+      ::bcopy(buffer+8,&(rec.namesize),4);
+#else
+      ::bcopy(buffer,(char*)&rec,AttributeRecordHdrSize);
+      //printf("AttribRec::reading ds,nt,ns %d,%d,%d\n",rec.datasize,rec.numbertype,rec.namesize);
+#endif
+      return r;
+    }
+  };
+
+  struct RecRef {
+    RecordHdr rec;
+    long offset;
+    //virtual ~RecRef();
+    // ** Inline these because they are 
+    // ** actually kludges for the T3E
+  };
+  
+  struct AttribRef : public RecRef {
+    //FlexArray<char> name;
+    char name[64];
+    AttribRef();
+    AttribRef(AttribRef &src);
+    //virtual ~AttribRef();
+    AttribRef &operator=(AttribRef &src);
+    // ** Inline these because they are 
+    // ** actually kludges for the T3E
+  };
+  
+  struct DataRef : public RecRef {
+    FlexArray<RecRef> annotations;
+    FlexArray<AttribRef> attributes;
+    long end_offset;
+    DataRef();
+    //virtual ~DataRef();
+    DataRef(DataRef &src);
+    DataRef &operator=(DataRef &src);
+  };
+private:
+  friend class AttribRef;
+  friend class DataRef;
+  friend class FileHdr;
+  friend class RecordHdr;
+  friend class DataRecordHdr;
+  friend class AttributeRecordHdr;
+  // convenience utilities
+#ifdef T3E
+  int ffioIsTrue();
+#endif
+  void initPosition();
+  inline long getPosition(){return this->lseek(0,L_INCR);}
+  void restart();
+  int nextRecord();
+  int nextDataRecord();
+  inline long getLength(){ return file_length; }
+  static void byteswapBuffer(void *buf,long nelements,int elementsize);
+  static void byteswapBuffer(void *source,void *dest,long nelements,int elementsize);
+  int writeFileHeader();
+  int readFileHeader();
+  void rebuildFileHeader();
+  void appendRecordTable();
+  void buildRecordTable();
+  void clearChunk(int nbytes);
+  void openFile(CONST char *fname,IObase::AccessMode access,int swbytes);
+  // annotations don't need a header
+  // for c, size includes null terminator.
+  // must decrement size for f77
+  //inline int write(int fid,void *data,int size);
+protected:
+  FlexArray<DataRef> rec_table;
+  FlexArray<Int> chunkdims;
+  FileHdr file_header;
+  RecordHdr current_rec,current_data_rec;
+  DataRecordHdr current_dat;
+  long current_rec_offset; // offset past the record header
+  long stream_offset; // kludge to allow index cache to be turned off
+  long current_dat_offset; // to datarec
+  int fid,datasetnumber,ndatasets,swapbytes,current_reserved_chunk;
+  int hasread,streaming,cur_type_size;
+  IEEEIO *masterfile;
+  int writebuffersize,writebuffercount;
+  char *writebuffer;
+  long savedposition;
+  long virtual_position,actual_position;
+  long file_length;
+protected:
+  inline int write(const void *data,size_t size);
+  inline void flush();
+  inline int lseek(long len, int direction);
+  inline int read(void *data,size_t size);
+public:
+  //------------------------core stuff.....................
+  enum AccessMode {Read=0,Write=1,Create=1,Append=2,SharedRead=4};
+  IObase::AccessMode mode(IEEEIO::AccessMode amode){
+    switch(amode){
+    case Read:
+      return IObase::Read;
+    case Write:
+      return IObase::Write;
+    case Append:
+      return IObase::Append;
+    default:
+      return IObase::Read;
+    }
+  }
+  // the IEEEIO:: commented out to satisfy Microsoft Compiler
+  IEEEIO(CONST char *fname, /*IEEEIO::*/AccessMode access,int swbytes=0);
+  IEEEIO(CONST char *fname,IObase::AccessMode access,int swbytes=0);
+  IEEEIO(IEEEIO *file); // read-only dup an existing open file
+  virtual ~IEEEIO();
+  virtual int isValid();
+  // could use overloading to differentiate type here... (but I'm going simple)
+  virtual int write(IObase::DataType typeID,int rank,CONST int *dims,void *data);
+  virtual int readInfo(IObase::DataType &typeID,int &rank,int *dims,int maxdims=3);
+  virtual int read(void *data);
+  //virtual int readChunk(int dims,int origin,int stride,void *data)=0;
+  virtual int seek(int dataset_index);
+  virtual int nDatasets();
+  
+  virtual int writeAnnotation(CONST char *annotation);
+  virtual int readAnnotationInfo(int number,int &length); // returns length (-1 if none left)
+  virtual int readAnnotation(int number,char *annotation,int maxsize=128);
+  virtual int nAnnotations();
+
+  virtual int writeAttribute(CONST char *name,IObase::DataType typeID,Long length,void *data);
+  // returns number
+  virtual int readAttributeInfo(int number,char *name,IObase::DataType &typeID,Long &nelem,int maxnamelen=128);
+  virtual int readAttributeInfo(CONST char *name,IObase::DataType &typeID,Long &nelem); // returns number
+  virtual int readAttribute(int number,void *data);
+  // virtual Long readAttribute(CONST char *name,void *data);
+  virtual int nAttributes();
+
+  //-----------------Chunking Utilities..................
+  virtual int reserveChunk(IObase::DataType typeID,int rank,CONST int *dims);
+  virtual int writeChunk(CONST int *chunkdims,CONST int *chunkorigin,void *data);
+  virtual int readChunk(CONST int *chunkdims,CONST int *chunkorigin,void *data);
+  // Streaming interface is for support of PANDA, Sockets etc..
+  virtual int reserveStream(IObase::DataType typeID,int rank,CONST int *dims);
+  virtual int writeStream(void *data,int length);
+  virtual int readStream(void *data,int length);
+  virtual int pause();
+  virtual int resume();
+  void bufferOn(long size=-1);
+  void bufferOff();
+};
+/*
+int bcopy(char *src,char *dst,int len){
+  // copying bytes on the DEC Alpha processor can be costly due
+  // to the penalty of non-aligned memory accesses.  This tries
+  // to align the copys.
+
+  You can use the program copyperf to see how bcopy() improves performance
+  You must build copyperf separately using
+     gmake copyperf
+  
+  The following will test copy a 64meg buffer 1 time.
+     copyperf 64m
+  The following will test copy a 32k buffer 1000 times
+     copyperf 32k 1000
+  The performance differences are 4-to-1 or greater in favor of bcopy().
+}
+*/
+int IEEEIO::write(const void *data,size_t size){
+  //puts("***********write()");
+#ifdef T3E
+  static int ffio_inline_true=0;
+  if(ffio_inline_true >= 0){
+    // do this test only once on first write
+    // this just gives someone a clue that they've done the
+    // wrong thing when they compiled
+#ifdef FFIO
+    if(!ffioIsTrue()){
+      fprintf(stderr,"Error!!! Inconsistent FFIO flags on a T3E!\n");
+      fprintf(stderr,"Your libraries have been compiled without -DFFIO, but your executable has been compled and linked with them with -DFFIO.\n");
+      fprintf(stderr,"\t*****Please either rebuild the IEEEIO library with -DFFIO or *remove* -DFFIO from your own build.******\n");
+    }
+#else
+    if(ffioIsTrue()){
+      fprintf(stderr,"Error!!! Inconsistent FFIO flags on a T3E!\n");
+      fprintf(stderr,"Your libraries have been compiled with -DFFIO, but your executable has been compled and linked with them without -DFFIO\n");
+      fprintf(stderr,"\t*****Please recompile with -DFFIO******\n");
+    }
+#endif
+    ffio_inline_true = -1;
+  }
+#endif
+  if(!writebuffer){
+    register long r=0;
+#ifdef FFIO
+    r = ::ffwrite(fid,data,size);
+#else
+    r = ::write(fid,data,size);
+#endif
+    if(r>0) actual_position += r;
+    virtual_position = actual_position;
+    // if(r>0) virtual_position = actual_position = (actual_position+r);
+    if(virtual_position>file_length) file_length = virtual_position;
+    return r;
+  }
+  // otherwise, do a buffered write
+  int retval=0;
+  char *datap = (char *)data;
+  if(writebuffercount>0){
+    // copy as much as you can to the writebuffer
+    // should use bcopy... its a whole lot faster!!!
+    long copysize = (size>(writebuffersize-writebuffercount))?(writebuffersize-writebuffercount):size;
+    bcopy(datap,writebuffer+writebuffercount,copysize);
+    writebuffercount += copysize;
+    virtual_position += copysize;
+    retval += copysize;
+    if(writebuffercount>=writebuffersize)
+      this->flush(); // dump immediately (updates both virtual and actual pos)
+    size-=copysize; datap+=copysize;
+    if(size<=0) {
+      // virtual_position = actual_position + writebuffercount;
+      if(virtual_position>file_length) file_length=virtual_position;
+      return retval;
+    }
+  }
+  // at this point size>0 and writebuffercount=0
+  if(size<=0 || writebuffercount!=0){
+    printf("***Write Failed: size=%u writebuffercount=%u\n",
+	   (unsigned int)size,writebuffercount);
+  }
+  if(size<writebuffersize){// store to buffer (size>0 is implicit in loop)
+    bcopy(datap,writebuffer+writebuffercount,size);
+    writebuffercount += size;
+    virtual_position += size;
+    retval += size;
+  }
+  else {
+    register long r=0;
+    // its larger than buffer, so bypass it and
+    // write it directly to the disk.
+#ifdef FFIO
+    r = ::ffwrite(fid,datap,size);
+#else
+    r = ::write(fid,datap,size);
+#endif
+    actual_position += r;
+    virtual_position = actual_position; // essentially equiv to flush
+    if(r>0) retval+=r;
+  }
+  // double-check virtual position
+  // virtual_position = actual_position + writebuffercount;
+  if(virtual_position>file_length) file_length = virtual_position;
+  return retval;
+}
+void IEEEIO::flush(){
+  register long r=0;
+  if(writebuffercount){
+#ifdef  FFIO
+    r = ::ffwrite(fid,writebuffer,writebuffercount);
+#else
+    r = ::write(fid,writebuffer,writebuffercount);
+#endif
+  }
+  writebuffercount=0;
+  if((r+actual_position)!=virtual_position){
+    fprintf(stderr,"IEEEIO::flush() : inconsistent file positions! r=%ld v=%ld a=%ld\n",
+	 r,virtual_position,actual_position);
+  }
+  actual_position = virtual_position;
+  if(virtual_position>file_length) file_length=virtual_position;
+}
+int IEEEIO::lseek(long len, int direction){
+  // compute new pos
+  register long npos=0;
+  // can seek beyond end of file
+  switch(direction){
+  case L_XTND: /* should do special case for 0 len */
+    npos = file_length + len;
+    break;
+  case L_SET:
+    npos = len;
+    break;
+  case L_INCR:;
+    npos = virtual_position + len;
+    break;
+  }
+  if(npos>file_length)
+    file_length = npos;
+  if(virtual_position != npos) {
+    if(writebuffer)
+      this->flush(); // must flush current buffer contents
+    // printf("Actually seeking from %ld to %ld\n",virtual_position,npos);
+    virtual_position = actual_position = npos;
+#ifdef  FFIO
+    // if(!writebuffer) 
+    return ::ffseek(fid,npos,L_SET);
+#else
+    // if(!writebuffer) 
+    return ::lseek(fid,npos,L_SET);
+#endif    
+  }
+  return npos; // otherwise no change required
+}
+/* Will eventually want to buffer the reads as well... */
+int IEEEIO::read(void *data,size_t size){
+  register int r=0;
+  if(writebuffercount>0) 
+    flush();
+#ifdef  FFIO
+  r = ::ffread(fid,data,size);
+#else
+  r = ::read(fid,data,size);
+#endif
+  if(r>0) virtual_position += r;
+  actual_position=virtual_position;
+  // can't read past end of file of course so no need to update file_length;
+  return r;
+}
+
+
+//==========F77 Interface=============
+#define f_ieee_open F77NAME(ieee_open_,ieee_open,IEEE_OPEN)
+#define f_ieee_openr F77NAME(ieee_openr_,ieee_openr,IEEE_OPENR)
+#define f_ieee_openw F77NAME(ieee_openw_,ieee_openw,IEEE_OPENW)
+#define f_ieee_opena F77NAME(ieee_opena_,ieee_opena,IEEE_OPENA)
+#define f_ieee_bufon F77NAME(ieee_bufon_,ieee_bufon,IEEE_BUFON)
+#define f_ieee_bufoff F77NAME(ieee_bufoff_,ieee_bufoff,IEEE_BUFOFF)
+extern "C"{
+#ifdef CRAY // Note: This isn't really implemented yet...
+#include <fortran.h>
+  Long8 f_ieee_open (_fcd fcfilename,_fcd fcaccessmode);
+  Long8 f_ieee_openr (_fcd fcfilename);
+  Long8 f_ieee_openw (_fcd fcfilename);
+  Long8 f_ieee_opena (_fcd fcfilename);
+#else
+  Long8 f_ieee_open (char *filename,char *accessmode,int namelen,int accesslength);
+  Long8 f_ieee_openr (char *filename,int namelen);
+  Long8 f_ieee_openw (char *filename,int namelen); // actually IObase::Create
+  Long8 f_ieee_opena (char *filename,int namelen);
+#endif
+  void f_ieee_bufon (Long8 *fileID,int *bufsize);
+  void f_ieee_bufoff (Long8 *fileID);
+  //=====================ANSI C interface
+#include "IEEEIO.h" // ANSI C interface
+	  }
+
+#endif
diff --git a/src/IEEEIO/IO.cc b/src/IEEEIO/IO.cc
new file mode 100644
index 0000000..e347f2a
--- /dev/null
+++ b/src/IEEEIO/IO.cc
@@ -0,0 +1,371 @@
+
+#include <stdio.h>
+#include <string.h>
+#include <limits.h>
+#include "IO.hh"
+
+IObase::IObase(CONST char *fname,AccessMode access):index(0),nitems(0),accessmode(access){
+  strcpy(filename,fname);
+}
+
+IObase::DataType IObase::Int2DataType(int dt){
+  if(dt==Byte)
+    return Byte;
+  else if(dt==Int16)
+    return Int16;
+  else if(dt==Int32)
+    return Int32;
+  else if(dt==Int64)
+    return Int64;
+  else if(dt==Float32)
+    return Float32;
+  else if(dt==Float64)
+    return Float64;
+  else if(dt==uInt8)
+    return uInt8;
+  else if(dt==uInt16)
+    return uInt16;
+  else if(dt==uInt32)
+    return uInt32;
+  else if(dt==uInt64)
+    return uInt64;
+  else if(dt==Char8)
+    return Char8;
+  else if(dt==Char16)
+    return Char16;
+  else return Error;
+}
+
+int IObase::sizeOf(IObase::DataType dt){
+  switch(dt){
+  case Int8:
+  case uInt8:
+  case Char8:
+    return 1;
+  case Int16:
+  case uInt16:
+  case Char16:
+    return 2;
+  case Int32:
+  case uInt32:
+  case Float32:
+    return 4;
+  case Int64:
+  case uInt64:
+  case Float64:
+    return 8;
+  default:
+    return 0;
+  }
+}
+
+int IObase::nElements(int rank,CONST int *dims){
+  int nelem,i;
+  for(i=0,nelem=1;i<rank;i++) nelem*=dims[i];
+  return nelem;
+}
+
+int IObase::nBytes(IObase::DataType dt, int rank,CONST int *dims){
+  return nElements(rank,dims)*sizeOf(dt);
+}
+
+int f_io_close(Long8 *deviceID){
+  IObase *dev=(IObase*)(*deviceID);
+  delete dev;
+  *deviceID=0;
+  return 1;
+}
+
+int f_io_isvalid(Long8 *deviceID){
+  if(!deviceID || !*deviceID) return 0;
+  IObase *dev=(IObase*)(*deviceID);
+  return dev->isValid();
+}
+
+int f_io_sizeof(int *type){
+  return IObase::sizeOf(IObase::Int2DataType(*type));
+}
+
+int f_io_nelements(int *rank,int *dims){
+  return IObase::nElements(*rank,dims);
+}
+
+int f_io_nbytes(int *type,int *rank,int *dims){
+  return IObase::nBytes(IObase::Int2DataType(*type),*rank,dims);
+}
+
+int f_io_write(Long8 *deviceID,int *typeID,int *rank,int *dims,void *data){
+  IObase *dev=(IObase*)(*deviceID);
+  dev->write(dev->Int2DataType(*typeID),*rank,dims,data);
+  return 1;
+}
+
+int f_io_readinfo(Long8 *deviceID,int *typeID,int *rank,int *dims,int *maxdims){
+  IObase *dev=(IObase*)(*deviceID);
+  IObase::DataType tid;
+  if(*maxdims<0) *maxdims=0;
+  int r=dev->readInfo(tid,*rank,dims,*maxdims);
+  //int r=dev->readInfo(tid,*rank,dims);
+  *typeID=tid;
+  return r;
+}
+
+int f_io_numdata(Long8 *deviceID){ 
+  IObase *dev=(IObase*)(*deviceID);
+  return dev->nDatasets();
+}
+int f_io_read(Long8 *deviceID,void *data){
+  IObase *dev=(IObase*)(*deviceID);
+  return dev->read(data);
+}
+
+int f_io_seek(Long8 *deviceID,int *dataset_index){
+  IObase *dev=(IObase*)(*deviceID);
+  return dev->seek(*dataset_index);
+}
+//--------------Annotations
+int f_io_writenote(Long8 *deviceID,char *annotation,int size){
+  IObase *dev=(IObase*)(*deviceID);
+  annotation[size]='\0'; // Yeah! its unsafe I know..!! (but it works!)
+  return dev->writeAnnotation(annotation);
+}
+
+int f_io_readnote(Long8 *deviceID,int *index,char *annotation,int maxsize){ 
+  IObase *dev=(IObase*)(*deviceID);
+  return dev->readAnnotation(*index,annotation,maxsize);
+}
+
+int f_io_noteinfo(Long8 *deviceID,int *index,int *length){ 
+  IObase *dev=(IObase*)(*deviceID);
+  return dev->readAnnotationInfo(*index,*length);
+}
+
+int f_io_numnote(Long8 *deviceID){ 
+  IObase *dev=(IObase*)(*deviceID);
+  return dev->nAnnotations();
+}
+
+//---------------Attributes
+int f_io_writeatt(Long8 *deviceID,char *name,
+		 int *datatype,Long *nelements,void *data,int namesize){
+  IObase *dev=(IObase*)(*deviceID);
+  name[namesize]='\0'; // cap name (to be certain).. unsafe but it works
+  // should copy into a flexarray which will be destroyed on return
+  IObase::DataType typeID = IObase::Int2DataType(*datatype);
+  return dev->writeAttribute(name,typeID,*nelements,data);
+}
+
+int f_io_attinfo(Long8 *deviceID,char *name,
+		int *datatype,Long *nelements,int namesize){
+  IObase *dev=(IObase*)(*deviceID);
+  IObase::DataType typeID;
+  if(namesize>0)
+    name[namesize]='\0';
+  int i=dev->readAttributeInfo(name,typeID,*nelements);
+  *datatype=typeID;
+  return i;
+}
+
+int f_io_iattinfo(Long8 *deviceID,int *index,char *name,
+		 int *datatype,Long *nelements,int namesize){
+  int i;
+  IObase *dev=(IObase*)(*deviceID);
+  IObase::DataType typeID;
+  for(i=0;i<namesize;i++) name[i]='\0';
+  //printf("io_iattinfo(): Namesize=%u Name=[%s]\n",namesize,name);
+  i=dev->readAttributeInfo(*index,name,typeID,*nelements,namesize);
+  // need to zero the array
+  //printf("io_iattinfo(): Newname=[%s]\n",name); 
+  *datatype=typeID;
+  //printf("io_iattinfo(): attribs are index=%u type=%u nelements=%u namesize=%u\n",
+  //	 *index,*datatype,*nelements,namesize);
+  return i;
+}
+
+int f_io_readatt(Long8 *deviceID,int *number,void *data){
+  IObase *dev=(IObase*)(*deviceID);
+  return dev->readAttribute(*number,data);
+}
+
+int f_io_numatt(Long8 *deviceID){
+  IObase *dev=(IObase*)(*deviceID);
+  return dev->nAttributes();
+}
+
+//==========F77 Chunking interface--------------------
+int f_io_reserveck(Long8 *deviceID,int *typeID,int *rank,int *dims){
+  IObase *dev=(IObase*)(*deviceID);
+  return dev->reserveChunk(dev->Int2DataType(*typeID),*rank,dims);
+}
+
+int f_io_writeck(Long8 *deviceID,int *chunkdims,int *chunkorigin,void *data){
+  IObase *dev=(IObase*)(*deviceID);
+  return dev->writeChunk(chunkdims,chunkorigin,data);
+}
+
+int f_io_readck(Long8 *deviceID,int *chunkdims,int *chunkorigin,void *data){
+  IObase *dev=(IObase*)(*deviceID);
+  return dev->readChunk(chunkdims,chunkorigin,data);
+}
+int f_io_writestrm(IOFile deviceID,void *data,int *length){
+  IObase *dev=(IObase*)(deviceID);
+  return dev->writeStream(data,*length);
+}
+
+int f_io_readstrm(IOFile deviceID,void *data,int *length){
+  IObase *dev=(IObase*)(deviceID);
+  return dev->readStream(data,*length);
+}
+
+int f_io_pause(Long8 *deviceID){
+  IObase *dev=(IObase*)(deviceID);
+  return dev->pause();
+}
+
+int f_io_resume(Long8 *deviceID){
+  IObase *dev=(IObase*)(deviceID);
+  return dev->resume();
+}
+
+//====================C Interface========================
+int IOclose(IOFile deviceID){
+  IObase *dev=(IObase*)(deviceID);
+  delete dev;
+  deviceID=0;
+  return 1;
+}
+
+int IOisValid(IOFile deviceID){
+  if(!deviceID) return 0;
+  IObase *dev=(IObase*)(deviceID);
+  return dev->isValid();
+}
+
+int IOsizeOf(int type){
+  return IObase::sizeOf(IObase::Int2DataType(type));
+}
+
+int IOnElements(int rank,int *dims){
+  return IObase::nElements(rank,dims);
+}
+
+int IOnBytes(int type,int rank,int *dims){
+  return IObase::nBytes(IObase::Int2DataType(type),rank,dims);
+}
+
+int IOwrite(IOFile deviceID,int typeID,int rank,int *dims,void *data){
+  IObase *dev=(IObase*)(deviceID);
+  return dev->write(dev->Int2DataType(typeID),rank,dims,data);
+}
+
+int IOreadInfo(IOFile deviceID,int *typeID,int *rank,int *dims,int maxdims){
+  IObase *dev=(IObase*)(deviceID);
+  IObase::DataType tid;
+  int r=dev->readInfo(tid,*rank,dims,maxdims);
+  *typeID=tid;
+  return r;
+}
+
+int IOread(IOFile deviceID,void *data){
+  IObase *dev=(IObase*)(deviceID);
+  return dev->read(data);
+}
+
+int IOseek(IOFile deviceID,int dataset_index){
+  IObase *dev=(IObase*)(deviceID);
+  return dev->seek(dataset_index);
+}
+
+int IOnDatasets(IOFile deviceID){
+  IObase *dev=(IObase*)(deviceID);
+  return dev->nDatasets();
+}
+
+int IOwriteAnnotation(IOFile deviceID,char *annotation){
+  IObase *dev=(IObase*)(deviceID);
+  //annotation[size]='\0'; // Yeah! its unsafe I know..!! (but it works!)
+  return dev->writeAnnotation(annotation);
+}
+ 
+int IOreadAnnotation(IOFile deviceID,int index,char *annotation,int maxsize){
+  IObase *dev=(IObase*)(deviceID);
+  return dev->readAnnotation(index,annotation,maxsize);
+}
+
+int IOreadAnnotationInfo(IOFile deviceID,int index,int *size){
+  IObase *dev=(IObase*)(deviceID);
+  return dev->readAnnotationInfo(index,*size);
+}
+
+int IOnAnnotations(IOFile deviceID){
+  IObase *dev=(IObase*)(deviceID);
+  return dev->nAnnotations();
+}
+
+int IOwriteAttribute(IOFile deviceID,char *name,int type,Long length,void *data){
+  IObase *dev=(IObase*)(deviceID);
+  IObase::DataType typeID = IObase::Int2DataType(type);
+  return dev->writeAttribute(name,typeID,length,data);
+}
+
+int IOreadIndexedAttributeInfo(IOFile deviceID,int number,char *name,
+			int *type,Long *nelem,int maxnamelen){
+  IObase *dev=(IObase*)(deviceID);
+  IObase::DataType typeID;
+  int i=dev->readAttributeInfo(number,name,typeID,*nelem,maxnamelen);
+  *type=typeID; // convert from enum
+  return i;
+}
+
+int IOreadAttributeInfo(IOFile deviceID,char *name,int *type,Long *nelem){
+  IObase *dev=(IObase*)(deviceID);
+  IObase::DataType typeID;
+  int i=dev->readAttributeInfo(name,typeID,*nelem);
+  *type=typeID; // convert from enum
+  return i;
+}  
+
+int IOreadAttribute(IOFile deviceID,int number,void *data){
+  IObase *dev=(IObase*)(deviceID);
+  return dev->readAttribute(number,data);
+}  
+
+int IOnAttributes(IOFile deviceID){
+  IObase *dev=(IObase*)(deviceID);
+  return dev->nAttributes();
+}
+
+//=============C Chunking interface---------------
+int IOreserveChunk(IOFile deviceID,int typeID,int rank,int *dims){
+  IObase *dev=(IObase*)(deviceID);
+  return dev->reserveChunk(dev->Int2DataType(typeID),rank,dims);
+}
+
+int IOwriteChunk(IOFile deviceID,int *chunkdims,int *chunkorigin,void *data){
+  IObase *dev=(IObase*)(deviceID);
+  return dev->writeChunk(chunkdims,chunkorigin,data);
+}
+
+int IOreadChunk(IOFile deviceID,int *chunkdims,int *chunkorigin,void *data){
+  IObase *dev=(IObase*)(deviceID);
+  return dev->readChunk(chunkdims,chunkorigin,data);
+}
+
+int IOwriteStream(IOFile deviceID,void *data,int length){
+  IObase *dev=(IObase*)(deviceID);
+  return dev->writeStream(data,length);
+}
+
+int IOreadStream(IOFile deviceID,void *data,int length){
+  IObase *dev=(IObase*)(deviceID);
+  return dev->readStream(data,length);
+}
+int IOpause(IOFile deviceID){
+  IObase *dev=(IObase*)(deviceID);
+  return dev->pause();
+}
+
+int IOresume(IOFile deviceID){
+  IObase *dev=(IObase*)(deviceID);
+  return dev->resume();
+}
diff --git a/src/IEEEIO/IO.h b/src/IEEEIO/IO.h
new file mode 100644
index 0000000..1a17cd3
--- /dev/null
+++ b/src/IEEEIO/IO.h
@@ -0,0 +1,63 @@
+/* IO.h C header */
+#ifndef __IO_H_
+#define __IO_H_
+#ifndef UNITYPE
+#define UNITYPE
+#define BYTE 0
+#define INT8 0
+#define INT16 1
+#define INT32 2
+#define INT64 3
+#define FLOAT32 4
+#define FLOAT64 5
+#define UCHAR 6
+#define UINT8 6
+#define UINT16 7
+#define UINT32 8
+#define UINT64 9
+/* special string types */
+#define CHAR 10
+#define CHAR8 10
+#define STRING 10
+#define UNICODE 11
+#define CHAR16 11
+#endif
+/* error */
+#define TYPEERROR -1
+
+#define IOREAD 0
+#define IOWRITE 1
+#define IOAPPEND 2
+#define IOCREATE 1 /* Same as IOWRITE */
+
+#include "Arch.h"
+
+int IOclose PROTO((IOFile deviceID));
+int IOisValid PROTO((IOFile deviceID));
+int IOsizeOf PROTO((int datatype));
+int IOnElements PROTO((int rank,int *dims));
+int IOnBytes PROTO((int datatype,int rank,int *dims));
+int IOwrite PROTO((IOFile deviceID,int typeID,int rank,int *dims,void *data));
+int IOreadInfo PROTO((IOFile deviceID,int *typeID,int *rank,int *dims,int maxdims));
+int IOread PROTO((IOFile deviceID,void *data));
+int IOseek PROTO((IOFile deviceID,int dataset_index));
+int IOnDatasets PROTO((IOFile deviceID));
+int IOwriteAnnotation PROTO((IOFile deviceID,char *annotation));
+int IOreadAnnotation PROTO((IOFile deviceID,int index,char *annotation,int maxsize));
+int IOreadAnnotationInfo PROTO((IOFile deviceID,int index,int *size));
+int IOnAnnotations PROTO((IOFile deviceID));
+
+int IOwriteAttribute PROTO((IOFile deviceID,char *name,int type,Long length,void *data));
+int IOreadIndexedAttributeInfo PROTO((IOFile deviceID,int number,
+				      char *name,int *type,Long *nelem,int maxnamelen));
+int IOreadAttributeInfo PROTO((IOFile deviceID,char *name,int *type,Long *nelem));
+int IOreadAttribute PROTO((IOFile deviceID,int number,void *data));
+int IOnAttributes PROTO((IOFile deviceID));
+int IOreserveChunk PROTO((IOFile deviceID,int typeID,int rank,int *dims));
+int IOwriteChunk PROTO((IOFile deviceID,int *chunkdims,int *chunkorigin,void *data));
+int IOreadChunk PROTO((IOFile deviceID,int *chunkdims,int *chunkorigin,void *data));
+int IOwriteStream PROTO((IOFile deviceID,void *data,int length));
+int IOreadStream PROTO((IOFile deviceID,void *data,int length));
+int IOpause PROTO((IOFile deviceID));
+int IOresume PROTO((IOFile deviceID));
+#endif
diff --git a/src/IEEEIO/IO.hh b/src/IEEEIO/IO.hh
new file mode 100644
index 0000000..443f2c9
--- /dev/null
+++ b/src/IEEEIO/IO.hh
@@ -0,0 +1,168 @@
+#ifndef __IO_HH_
+#define __IO_HH_
+#include "Arch.h"
+
+class IObase {
+public:
+  //-----------------------Public Enums....................
+  enum AccessMode {Read=0,Write=1,Create=1,Append=2};
+  enum DataType {Byte=0,Int8=0,Int16=1,Int32=2,Int64=3,
+		 Float32=4,Float64=5,
+		 uInt8=6,uChar=6,uInt16=7,uInt32=8,uInt64=9,
+		 Char=10,Char8=10,String=10,Unicode=11,Char16=11, // special string types
+		 Error=-1};
+protected:
+  //-----------------------State variables.................
+  int index,nitems;
+  AccessMode accessmode;
+  char filename[128];
+public:
+  //------------------------core stuff.....................
+  IObase(CONST char *fname,AccessMode access);
+  virtual ~IObase() {}
+  virtual int isValid() { return 0; }
+  // could use overloading to differentiate type here... (but I'm going simple)
+  virtual int write(DataType typeID,int rank,CONST int *dims,void *data)=0;
+  virtual int readInfo(DataType &typeID,int &rank,int *dims,int maxdims=3)=0;
+  virtual int read(void *data)=0;
+  virtual int seek(int dataset_index)=0;
+  virtual int nDatasets()=0;
+  virtual int writeAnnotation(CONST char *annotation)=0;
+  virtual int readAnnotationInfo(int number,int &length)=0; // returns length (-1 if none left)
+  virtual int readAnnotation(int number,char *annotation,int maxsize=128)=0;
+  virtual int nAnnotations()=0;
+  
+  virtual int writeAttribute(CONST char *name,IObase::DataType typeID,Long length,void *data)=0;
+  // returns number
+  virtual int readAttributeInfo(int number,char *name,IObase::DataType &typeID,Long &nelem,int maxnamelen=128)=0;
+  virtual int readAttributeInfo(CONST char *name,IObase::DataType &typeID,Long &nelem)=0; // returns number
+  virtual int readAttribute(int number,void *data)=0;
+  // virtual Long readAttribute(char *name,void *data);
+  virtual int nAttributes()=0;
+
+  //-----------------Chunking Features (for MPI/HPF)................
+  virtual int reserveChunk(IObase::DataType typeID,int rank,CONST int *dims)=0;
+  
+  virtual int writeChunk(CONST int *chunkdims,CONST int *chunkorigin,void *data)=0;
+  // virtual int writeStridedChunk()=0;
+  // virtual int writeHPF(int processorID,int *proclayout,IObase::HPFlayout *arraylayout,void *data)=0;
+  virtual int readChunk(CONST int *chunkdims,CONST int *chunkorigin,void *data)=0;
+  // virtual int readStrided()=0;
+  // virtual int readHPF(int processorID,int *proclayout,IObase::HPFlayout *arraylayout,void *data)=0;
+  //-----------------Streaming interface (for PANDA, Sockets & MPIO).........
+  virtual int reserveStream(IObase::DataType typeID,int rank,CONST int *dims)=0;
+  virtual int writeStream(void *data,int length)=0;
+  virtual int readStream(void *data,int length)=0;
+  //-----------------------Utilities........................
+  // unfortunately you can cast enums to int but not the reverse
+  static DataType Int2DataType(int dt);
+  static int sizeOf(DataType dt);
+  static int nElements(int rank,CONST int *dims); // returns number of elements in dataset
+  static int nBytes(DataType dt,int rank,CONST int *dims); // returns number of bytes (good for malloc)
+  CONST char *name() { return filename; }
+  virtual int pause()   { return 0; }
+  virtual int resume() { return 0; }
+};
+
+//===================Fortran77 Interface 
+#define f_io_close F77NAME(io_close_,io_close,IO_CLOSE)
+#define f_io_isvalid F77NAME(io_isvalid_,io_isvalid,IO_ISVALID)
+#define f_io_sizeof F77NAME(io_sizeof_,io_sizeof,IO_SIZEOF)
+#define f_io_nbytes F77NAME(io_nbytes_,io_nbytes,IO_NBYTES)
+#define f_io_nelements F77NAME(io_nelements_,io_nelements,IO_NELEMENTS)
+#define f_io_write F77NAME(io_write_,io_write,IO_WRITE)
+#define f_io_readinfo F77NAME(io_readinfo_,io_readinfo,IO_READINFO)
+#define f_io_read F77NAME(io_read_,io_read,IO_READ)
+#define f_io_seek F77NAME(io_seek_,io_seek,IO_SEEK)
+#define f_io_numdata F77NAME(io_numdata_,io_numdata,IO_NUMDATA)
+#define f_io_writenote F77NAME(io_writenote_,io_writenote,IO_WRITENOTE)
+#define f_io_readnote F77NAME(io_readnote_,io_readnote,IO_READNOTE)
+#define f_io_numnote F77NAME(io_numnote_,io_numnote,IO_NUMNOTE)
+#define f_io_writeatt F77NAME(io_writeatt_,io_writeatt,IO_WRITEATT)
+#define f_io_attinfo F77NAME(io_attinfo_,io_attinfo,IO_ATTINFO)
+#define f_io_iattinfo F77NAME(io_iattinfo_,io_iattinfo,IO_IATTINFO)
+#define f_io_readatt F77NAME(io_readatt_,io_readatt,IO_READATT)
+#define f_io_numatt F77NAME(io_numatt_,io_numatt,IO_NUMATT)
+#define f_io_reserveck F77NAME(io_reserveck_,io_reserveck,IO_RESERVECK)
+#define f_io_writeck F77NAME(io_writeck_,io_writeck,IO_WRITECK)
+#define f_io_readck F77NAME(io_readck_,io_readck,IO_READCK)
+#define f_io_writestrm F77NAME(io_writestrm_,io_writestrm,IO_WRITESTRM)
+#define f_io_readstrm F77NAME(io_readstrm_,io_readstrm,IO_READSTRM)
+#define f_io_pause F77NAME(io_pause_,io_pause,IO_PAUSE)
+#define f_io_resume F77NAME(io_resume_,io_resume,IO_RESUME)
+extern "C"{ 
+//================Ansi C interface
+#include "IO.h"
+//==== f77 interface
+#ifdef CRAY
+  // Cray fortran uses an _fcd data structure to pass strings.  This is in preparation
+  // for using this on Crays, but since nobody is using them yet, I've left the code for
+  // these versions of the routines blank.  It'll be filled in later if the need arises
+#include <fortran.h>
+  int f_io_close (Long8 *deviceID);
+  int f_io_isvalid (Long8 *deviceID);
+  int f_io_sizeof (int *datatype);
+  int f_io_nelements (int *rank,int *dims);
+  int f_io_nbytes (int *datatype,int *rank,int *dims);
+  int f_io_write (Long8 *deviceID,int *typeID,int *rank,int *dims,void *data);
+  int f_io_readinfo (Long8 *deviceID,int *typeID,int *rank,int *dims,int *maxdims);
+  int f_io_read (Long8 *deviceID,void *data);
+  int f_io_seek (Long8 *deviceID,int *dataset_index);
+  int f_io_numdata (Long8 *deviceID);
+  int f_io_writenote (Long8 *deviceID,_fcd fcannotation);
+  int f_io_noteinfo (Long8 *deviceID,int *number,int *length);
+  int f_io_readnote (Long8 *deviceID,int *number,_fcd fcannotation);
+  int f_io_numnote (Long8 *deviceID);
+
+  int f_io_writeatt (Long8 *deviceID,_fcd fcname,
+		   int *datatype,Long *nelements,void *data);
+  int f_io_attinfo (Long8 *deviceID,_fcd fcname,
+		  int *datatype,Long *nelements);
+  int f_io_iattinfo (Long8 *deviceID,int *index,_fcd fcname,
+		  int *datatype,Long *nelements);
+  int f_io_readatt (Long8 *deviceID,int *number,void *data);
+  int f_io_numatt (Long8 *deviceID);
+  int f_io_reserveck(Long8 *deviceID,int *typeID,int *rank,int *dims);
+  int f_io_writeck(Long8 *deviceID,int *chunkdims,int *chunkorigin,void *data);
+  int f_io_readck(Long8 *deviceID,int *chunkdims,int *chunkorigin,void *data);
+  int f_io_pause(Long8 *deviceID);
+  int f_io_resume(Long8 *deviceID);
+#else
+  int f_io_close (Long8 *deviceID);
+  int f_io_isvalid (Long8 *deviceID);
+  int f_io_sizeof (int *datatype);
+  int f_io_nelements (int *rank,int *dims);
+  int f_io_nbytes (int *datatype,int *rank,int *dims);
+  int f_io_write (Long8 *deviceID,int *typeID,int *rank,int *dims,void *data);
+  int f_io_readinfo (Long8 *deviceID,int *typeID,int *rank,int *dims,int *maxdims);
+  int f_io_read (Long8 *deviceID,void *data);
+  int f_io_seek (Long8 *deviceID,int *dataset_index);
+  int f_io_numdata (Long8 *deviceID);
+
+  int f_io_writenote (Long8 *deviceID,char *annotation,int size);
+  int f_io_noteinfo (Long8 *deviceID,int *number,int length);
+  int f_io_readnote (Long8 *deviceID,int *number,char *annotation,int maxsize);
+  int f_io_numnote (Long8 *deviceID);
+
+  int f_io_writeatt (Long8 *deviceID,char *name,
+		   int *datatype,Long *nelements,void *data,int namesize);
+  int f_io_attinfo (Long8 *deviceID,char *name,
+		  int *datatype,Long *nelements,int namesize);
+  int f_io_iattinfo (Long8 *deviceID,int *index,char *name,
+		  int *datatype,Long *nelements,int namesize);
+  int f_io_readatt (Long8 *deviceID,int *number,void *data);
+  int f_io_numatt (Long8 *deviceID);
+  int f_io_reserveck(Long8 *deviceID,int *typeID,int *rank,int *dims);
+  int f_io_writeck(Long8 *deviceID,int *chunkdims,int *chunkorigin,void *data);
+  int f_io_readck(Long8 *deviceID,int *chunkdims,int *chunkorigin,void *data);
+  int f_io_writestrm(Long8 *deviceID,void *data,int *length);
+  int f_io_readstrm(Long8 *deviceID,void *data,int *length);
+  int f_io_pause(Long8 *deviceID);
+  int f_io_resume(Long8 *deviceID);
+#endif
+}
+
+#endif
+
+
+
diff --git a/src/IEEEIO/make.code.defn b/src/IEEEIO/make.code.defn
new file mode 100644
index 0000000..c945b9e
--- /dev/null
+++ b/src/IEEEIO/make.code.defn
@@ -0,0 +1,84 @@
+# Main make.code.defn file for thorn FlexIO
+# $Header$
+
+# Source files in this directory
+SRCS = IO.cc IEEEIO.cc
+
+# Subdirectories containing source files
+SUBDIRS = 
+
+# Work out some thorn specific compilation flags
+# These are stripped down lines from the original IEEEIO makefile
+
+# The 9000 names of the cygwin tools and T3E...
+TMPUN   := $(shell uname)
+ifeq ($(TMPUN), CYGWIN32_95)
+UNAME   = CYGWIN
+else
+ifeq ($(TMPUN), CYGWIN32_NT)
+UNAME   = CYGWIN
+else
+ifeq ($(TMPUN), CYGWIN_NT-4.0)
+UNAME   = CYGWIN
+else 
+UNAME            := $(shell uname | perl -pe 's/(sn\d\d\d\d|jsimpson)/UNICOS\/mk/')
+endif
+endif
+endif
+
+# 64 Bit Irix
+ifeq ($(UNAME), IRIX64)
+
+CXXFLAGS += -DANSI -ptused 
+
+endif
+
+# 32 Bit Irix
+ifeq ($(UNAME), IRIX)
+
+CXXFLAGS += -DANSI -ptused 
+
+endif
+
+# HP
+ifeq ($(UNAME), HP-UX)
+
+CXXFLAGS += -DANSI -DHP +a1
+
+endif
+
+# Alpha
+ifeq ($(UNAME), OSF1)
+
+CXXFLAGS += -DANSI 
+
+endif
+
+# Linux
+ifeq ($(UNAME), Linux)
+
+CXXFLAGS += -DANSI 
+
+endif
+
+# Macintosh /PowerMach-MachTen
+ifeq ($(UNAME), machten)
+
+CXXFLAGS += -DANSI 
+
+endif
+
+# Cygwin / Win32
+ifeq ($(UNAME), CYGWIN)
+
+CXXFLAGS += -DANSI -DWIN32
+
+endif
+
+# T3E
+ifeq ($(UNAME), UNICOS/mk)
+
+CXXFLAGS += -DANSI -DT3E -hinstantiate=used 
+
+endif
+
diff --git a/src/Output2D.c b/src/Output2D.c
new file mode 100644
index 0000000..23b1298
--- /dev/null
+++ b/src/Output2D.c
@@ -0,0 +1,283 @@
+ /*@@
+   @file      Output2D.c
+   @date      Tue Jan 9 1999
+   @author    Gabrielle Allen
+   @desc 
+   Functions to deal 2D output of GFs
+   @enddesc 
+   @history
+   @hauthor Thomas Radke @hdate 16 Mar 1999
+   @hdesc Converted to Cactus 4.0
+   @hendhistory
+ @@*/
+
+#include <stdio.h>
+#include <string.h>
+
+#include "flesh.h"
+#include "Groups.h"
+#include "declare_parameters.h"
+#include "CactusBase/IOUtil/src/ioGH.h"
+#include "ioFlexGH.h"
+
+
+int IOFlexIO_Output2DVarAs (cGH *GH, const char *var, const char *alias);
+int IOFlexIO_TimeFor2D (cGH *GH, int index);
+
+/*@@
+   @routine    IOFlexIO_Output2DGH
+   @date       Sat March 6 1999
+   @author     Gabrielle Allen
+   @desc 
+   Loops over all variables and outputs them if necessary
+   @enddesc 
+   @calls    CCTK_GetGHExtensionHandle 
+             CCTK_GetNumVars
+	     CCTK_GetImplementationFromVar
+	     CCTK_GetVarName
+             IOFlexIO_TimeFor2D
+	     IOFlexIO_Output2DVarAs  	     
+   @calledby
+   @history 
+ 
+   @endhistory 
+   @var     GH
+   @vdesc   Pointer to CCTK GH
+   @vtype   cGH
+   @vio     in
+   @vcomment 
+   @endvar 
+@@*/
+
+int IOFlexIO_Output2DGH (cGH *GH)
+{
+  DECLARE_PARAMETERS
+  int i;
+  flexioGH *myGH;
+  char *name, *fullname;
+
+  /* Get the GH extension for IOFlexIO */
+  myGH = (flexioGH *) GH->extensions [CCTK_GetGHExtensionHandle ("IOFlexIO")];
+
+  /* Loop over all variables */
+  for (i = 0; i < CCTK_GetNumVars (); i++) {
+
+    if (IOFlexIO_TimeFor2D (GH, i)) {
+
+      name = CCTK_GetVarName (i);
+      fullname = CCTK_GetFullName (i);
+
+      if (IO_verbose) {
+        printf("IOFlexIO Output2DGH : \n");
+        printf("  fullname/name = %s/%s\n", fullname, name);
+      }
+
+      IOFlexIO_Output2DVarAs (GH, fullname, name);
+
+      free (fullname);
+
+      /* Register another 2D output for this GF */
+      myGH->IO_2Dnum[i]++;
+
+      /* Register GF as having 2D output this iteration */
+      myGH->IO_2Dlast[i] = GH->iteration;
+    }
+  }
+
+  return (0);
+}
+
+
+/*@@
+   @routine    IOFlexIO_Output2DVarAs
+   @date       Sat March 6 1999
+   @author     Gabrielle Allen
+   @desc 
+   unconditional output of a variable using the IO 2D output method
+   @enddesc 
+   @calls    CCTK_DecomposeName
+	     CCTK_GetVarIndex
+	     CCTK_GetGHExtensionHandle
+	     IOFlexIO_Write2D  	     
+   @calledby IOFlexIO_Output2DGH
+   @history 
+ 
+   @endhistory 
+   @var     GH
+   @vdesc   Pointer to CCTK GH
+   @vtype   cGH
+   @vio     in
+   @vcomment 
+   @endvar 
+   @var     fullname
+   @vdesc   complete name of variable to output
+   @vtype   const char *
+   @vio     in
+   @vcomment 
+   @endvar 
+   @var     alias
+   @vdesc   alias name of variable to output (used to generate output filename)
+   @vtype   const char *
+   @vio     in
+   @vcomment 
+   @endvar 
+@@*/
+
+int IOFlexIO_Output2DVarAs (cGH *GH, const char *fullname, const char *alias)
+{
+  DECLARE_PARAMETERS
+  int index, first;
+  flexioGH *myGH;
+
+  index = CCTK_GetVarIndex (fullname);
+
+  /* get the handle for IOFlexIO extensions */
+  myGH = (flexioGH *) GH->extensions [CCTK_GetGHExtensionHandle ("IOFlexIO")];
+
+  /* Output with correct file opening behaviour */
+  first = (myGH->IO_2Dnum [index] == 1);
+
+  if (IO_verbose) {
+    printf ("\nIn IOFlexIO Output2DVarAs\n-------------------\n");
+    printf ("  Fullname = -%s-\n", fullname);
+    printf ("  Alias = -%s-\n", alias);
+    printf ("  Index = %d\n", index);
+  }
+
+  /* Do the 2D output */
+  IOFlexIO_Write2D (GH, index, alias, first);
+
+  return (0);
+}
+
+
+ /*@@
+   @routine    IOFlexIO_TimeFor2D
+   @date       Sat March 6 1999
+   @author     Gabrielle Allen
+   @desc 
+   Decides if it is time to output a variable using the IO 2D output
+   method
+   @enddesc 
+   @calls    CCTK_GetGHExtensionHandle 
+             CCTK_GetVarGType
+	     CCTK_WARN
+	     CCTK_QueryGroupStorage_ByIndex
+	     CCTK_GetGroupNameFromVar_ByIndex
+   @calledby IOFlexIO_Output2DGH
+   @history 
+ 
+   @endhistory 
+   @var     GH
+   @vdesc   Pointer to CCTK GH
+   @vtype   cGH
+   @vio     in
+   @vcomment 
+   @endvar 
+   @var     index
+   @vdesc   index of variable
+   @vtype   int
+   @vio     in
+   @vcomment 
+   @endvar 
+@@*/
+
+		  
+int IOFlexIO_TimeFor2D (cGH *GH, int index)
+{
+  ioGH *ioUtilGH;
+  flexioGH *myGH;
+
+  /* Get the GH extension for IOUtil and IOFlexIO */
+  ioUtilGH = (ioGH *) GH->extensions [CCTK_GetGHExtensionHandle ("IO")];
+  myGH = (flexioGH *) GH->extensions [CCTK_GetGHExtensionHandle ("IOFlexIO")];
+
+  /* Check this GF should be output */
+  if (! (CCTK_GetVarGType (index) == GROUP_GF && myGH->IO_2Dnum [index] != 0
+         && GH->iteration % ioUtilGH->IO_2Devery == 0))
+    return (0);
+
+  /* Check GF not already output this iteration */
+  if (myGH->IO_2Dlast [index] == GH->iteration) {
+    CCTK_WARN (2, "Already done 2D output in IOFlexIO");
+    return (0);
+  }
+
+  /* Check GF has storage */
+  if (! CCTK_QueryGroupStorage_ByIndex (GH,
+            CCTK_GetGroupIndexFromVar_ByIndex(index))) {
+    char *fullname = CCTK_GetFullName (index);
+    char *msg = (char *) malloc (80 + strlen (fullname));
+
+    sprintf (msg, "No IOFlexIO 2D output for '%s' (no storage)", fullname);
+    CCTK_WARN (2, msg);
+    free (fullname);
+    free (msg);
+    return (0);
+  }
+
+  return (1);
+}
+
+
+/*@@
+   @routine    IOFlexIO_TriggerOutput2D
+   @date       Sat March 6 1999
+   @author     Gabrielle Allen
+   @desc 
+   Triggers the output a variable using the IO 2D output
+   method
+   @enddesc 
+   @calls    CCTK_GetGHExtensionHandle 
+             CCTK_GetVarName
+	     IOFlexIO_Write2D  	     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+   @var     GH
+   @vdesc   Pointer to CCTK GH
+   @vtype   cGH
+   @vio     in
+   @vcomment 
+   @endvar 
+   @var     index
+   @vdesc   index of variable to output
+   @vtype   int
+   @vio     in
+   @vcomment 
+   @endvar 
+@@*/
+
+int IOFlexIO_TriggerOutput2D (cGH *GH, int index)
+{
+  DECLARE_PARAMETERS
+  int first;
+  flexioGH *myGH;
+  char *varname;
+
+  varname = CCTK_GetVarName (index);
+
+  /* get the handle for IOFlexIO extensions */
+  myGH = (flexioGH *) GH->extensions [CCTK_GetGHExtensionHandle ("IOFlexIO")];
+
+  /* Output with correct file opening behaviour */
+  first =  (myGH->IO_2Dnum [index] == 1);
+
+  if (IO_verbose) {
+    printf("\nIn IOFlexIO TriggerOutput2D\n---------------------\n");
+    printf("  Index = %d\n", index);
+    printf("  Variable = -%s-\n", varname);
+  }
+
+  /* Do the 2D output */
+  IOFlexIO_Write2D (GH, index, varname, first);
+
+  /* Register another 2D output for this GF */
+  myGH->IO_2Dnum [index]++;
+
+  /* Register GF as having 2D output this iteration */
+  myGH->IO_2Dlast [index] = GH->iteration;
+
+  return (0);
+}
diff --git a/src/Output3D.c b/src/Output3D.c
new file mode 100644
index 0000000..64f0b1a
--- /dev/null
+++ b/src/Output3D.c
@@ -0,0 +1,285 @@
+ /*@@
+   @file      Output3D.c
+   @date      Tue Jan 9 1999
+   @author    Gabrielle Allen
+   @desc 
+   Functions to deal 3D output of GFs
+   @enddesc 
+   @history
+   @hauthor Thomas Radke @hdate 16 Mar 1999
+   @hdesc Converted to Cactus 4.0
+   @hendhistory
+ @@*/
+
+#include <stdio.h>
+#include <assert.h>
+
+#include "flesh.h"
+#include "Groups.h"
+#include "declare_parameters.h"
+#include "CactusBase/IOUtil/src/ioGH.h"
+#include "ioFlexGH.h"
+
+int IOFlexIO_Output3DVarAs (cGH *GH, const char *var, const char *alias);
+int IOFlexIO_TimeFor3D (cGH *GH, int index);
+
+/*@@
+   @routine    IOFlexIO_Output3DGH
+   @date       Sat March 6 1999
+   @author     Gabrielle Allen
+   @desc
+   Loops over all variables and outputs them if necessary
+   @enddesc
+   @calls    CCTK_GetGHExtensionHandle
+             CCTK_GetNumVars
+             CCTK_GetImplementationFromVar
+             CCTK_GetVarName
+             IOFlexIO_TimeFor3D
+             IOFlexIO_Output3DVarAs
+   @calledby
+   @history
+
+   @endhistory
+   @var     GH
+   @vdesc   Pointer to CCTK GH
+   @vtype   cGH
+   @vio     in
+   @vcomment
+   @endvar
+@@*/
+
+int IOFlexIO_Output3DGH (cGH *GH)
+{
+  int i;
+  flexioGH *myGH;
+  char *implementation;
+  char *name;
+  char *fullname;
+  DECLARE_PARAMETERS
+
+  /* Get the GH extension for IOFlexIO */
+  myGH = (flexioGH *) GH->extensions [CCTK_GetGHExtensionHandle ("IOFlexIO")];
+
+  /* Loop over all variables */
+  for (i = 0; i < CCTK_GetNumVars (); i++) {
+    if (IOFlexIO_TimeFor3D (GH, i)) {
+      implementation = CCTK_GetImplementationFromVar (i);
+      name           = CCTK_GetVarName (i);
+      fullname       = (char *) malloc (strlen (implementation) +
+                                        strlen (name) + 3);
+      assert (fullname);
+      sprintf (fullname, "%s::%s", implementation, name);
+
+      if (IO_verbose) {
+        printf ("IOFlexIO Output3DGH : \n");
+        printf ("  fullname/name = %s/%s\n", fullname, name);
+      }
+
+      IOFlexIO_Output3DVarAs (GH, fullname, name);
+
+      free (fullname);
+
+      /* Register another 3D output for this GF */
+      myGH->IO_3Dnum [i]++;
+
+      /* Register GF as having 3D output this iteration */
+      myGH->IO_3Dlast [i] = GH->iteration;
+    }
+  }
+
+  return (0);
+}
+
+
+/*@@
+   @routine    IOFlexIO_Output3DVarAs
+   @date       Sat March 6 1999
+   @author     Gabrielle Allen
+   @desc
+   unconditional output of a variable using the IOFlexIO 3D output method
+   @enddesc
+   @calls    CCTK_DecomposeName
+             CCTK_GetVarIndex
+             CCTK_GetGHExtensionHandle
+             IOFlexIO_Write3D
+   @calledby IOFlexIO_Output3DGH
+   @history
+
+   @endhistory
+   @var     GH
+   @vdesc   Pointer to CCTK GH
+   @vtype   cGH
+   @vio     in
+   @vcomment
+   @endvar
+   @var     fullname
+   @vdesc   complete name of variable to output
+   @vtype   const char *
+   @vio     in
+   @vcomment
+   @endvar
+   @var     alias
+   @vdesc   alias name of variable to output (used to generate output filename)
+   @vtype   const char *
+   @vio     in
+   @vcomment
+   @endvar
+@@*/
+
+int IOFlexIO_Output3DVarAs (cGH *GH, const char *fullname, const char *alias)
+{
+  DECLARE_PARAMETERS
+  int index, first;
+  flexioGH *myGH;
+
+  index = CCTK_GetVarIndex(fullname);
+
+  /* Get the GH extension for IOFlexIO */
+  myGH = (flexioGH *) GH->extensions [CCTK_GetGHExtensionHandle ("IOFlexIO")];
+
+  /* Output with correct file opening behaviour */
+  first = (myGH->IO_3Dnum [index] == 1);
+
+  if (IO_verbose) {
+    printf ("\nIn IOFlexIO Output3DVarAs\n-------------------\n");
+    printf ("  Fullname = -%s-\n", fullname);
+    printf ("  Alias = -%s-\n", alias);
+    printf ("  Index = %d\n", index);
+  }
+
+  /* Do the 3D output */
+  IOFlexIO_Write3D (GH, index, alias, first);
+
+  return (0);
+}
+
+
+/*@@
+   @routine    IOFlexIO_TimeFor3D
+   @date       Sat March 6 1999
+   @author     Gabrielle Allen
+   @desc
+   Decides if it is time to output a variable using the IOFlexIO 3D output
+   method
+   @enddesc
+   @calls    CCTK_GetGHExtensionHandle
+             CCTK_GetVarGType
+             CCTK_WARN
+             CCTK_QueryGroupStorage_ByIndex
+             CCTK_GetGroupNameFromVar_ByIndex
+   @calledby IOFlexIO_Output3DGH
+   @history
+
+   @endhistory
+   @var     GH
+   @vdesc   Pointer to CCTK GH
+   @vtype   cGH
+   @vio     in
+   @vcomment
+   @endvar
+   @var     index
+   @vdesc   index of variable
+   @vtype   int
+   @vio     in
+   @vcomment
+   @endvar
+@@*/
+
+int IOFlexIO_TimeFor3D (cGH *GH, int index)
+{
+  ioGH *ioUtilGH;
+  flexioGH *myGH;
+
+  /* Get the GH extension for IOUtil and IOFlexIO */
+  ioUtilGH = (ioGH *) GH->extensions [CCTK_GetGHExtensionHandle ("IO")];
+  myGH = (flexioGH *) GH->extensions [CCTK_GetGHExtensionHandle ("IOFlexIO")];
+
+  /* Check this GF should be output */
+  if (! (myGH->IO_3Dnum [index] != 0 &&
+         GH->iteration % ioUtilGH->IO_3Devery == 0))
+    return (0);
+
+  /* Check GF not already output this iteration */
+  if (myGH->IO_3Dlast [index] == GH->iteration) {
+      CCTK_WARN (2, "Already done 3D output in IOFlexIO");
+      return (0);
+  }
+
+  /* Check GF has storage */
+  if (! CCTK_QueryGroupStorage_ByIndex (GH,
+            CCTK_GetGroupIndexFromVar_ByIndex(index))) {
+    char *fullname = CCTK_GetFullName (index);
+    char *msg = (char *) malloc (80 + strlen (fullname));
+
+    sprintf (msg, "No IOFlexIO 3D output for '%s' (no storage)", fullname);
+    CCTK_WARN (2, msg);
+    free (fullname);
+    free (msg);
+    return (0);
+  }
+
+  return (1);
+}
+
+
+/*@@
+   @routine    IOFlexIO_TriggerOutput3D
+   @date       Sat March 6 1999
+   @author     Gabrielle Allen
+   @desc
+   Triggers the output a variable using the IOFlexIO 3D output
+   method
+   @enddesc
+   @calls    CCTK_GetGHExtensionHandle
+             CCTK_GetVarName
+             IOFlexIO_Write3D
+   @calledby
+   @history
+
+   @endhistory
+   @var     GH
+   @vdesc   Pointer to CCTK GH
+   @vtype   cGH
+   @vio     in
+   @vcomment
+   @endvar
+   @var     index
+   @vdesc   index of variable to output
+   @vtype   int
+   @vio     in
+   @vcomment
+   @endvar
+@@*/
+
+int IOFlexIO_TriggerOutput3D (cGH *GH, int index)
+{
+  DECLARE_PARAMETERS
+  int first;
+  flexioGH *myGH;
+  char *varname;
+
+  varname = CCTK_GetVarName (index);
+
+  /* Get the GH extension for IOFlexIO */
+  myGH = (flexioGH *) GH->extensions [CCTK_GetGHExtensionHandle ("IOFlexIO")];
+
+  /* Output with correct file opening behaviour */
+  first =  (myGH->IO_3Dnum [index] == 1);
+
+  if (IO_verbose) {
+    printf("\nIn IOFlexIO TriggerOutput3D\n---------------------\n");
+    printf("  Index = %d\n", index);
+    printf("  Variable = -%s-\n", varname);
+  }
+
+  /* Do the 3D output */
+  IOFlexIO_Write3D (GH, index, varname, first);
+
+  /* Register another 3D output for this GF */
+  myGH->IO_3Dnum [index]++;
+
+  /* Register GF as having 3D output this iteration */
+  myGH->IO_3Dlast [index] = GH->iteration;
+
+  return (0);
+}
diff --git a/src/RecoverGH.c b/src/RecoverGH.c
new file mode 100644
index 0000000..0391b30
--- /dev/null
+++ b/src/RecoverGH.c
@@ -0,0 +1,343 @@
+ /*@@
+   @file      RecoverGH.c
+   @date      Fri Jun 19 09:14:22 1998
+   @author    Tom Goodale
+   @desc 
+   Contains the routines to do the internal checkpoint recovery.
+
+   Currently can recover from:
+     (1) One file containing recombined data
+     (2) Multiple unrecombined files, where the current
+	     number of processors and outputing processors
+		 match those used to write the data.
+   @enddesc 
+   @history
+   @hauthor Gabrielle Allen @hdate 19 Oct 1998      
+   @hdesc Changed names ready for thorn_IO
+   @endhistory
+   @version $Id$
+ @@*/
+
+static char *rcsid = "$Id$";
+
+#include <stdio.h>
+
+#include "cctk.h"
+#include "flesh.h"
+#include "declare_parameters.h"
+#include "GHExtensions.h"
+#include "WarnLevel.h"
+#include "Comm.h"
+#ifdef CACTUSBASE_PUGH
+#include "CactusBase/pugh/src/include/pugh.h"
+#endif
+#include "CactusBase/IOUtil/src/ioGH.h"
+#include "ioFlexGH.h"
+
+
+/* this one comes from RestoreFile.c */
+int IOFlexIO_RestoreIEEEIOfile (cGH *GH, IOFile ifp,
+                              int IOrecover_ioproc,
+                              int IOrecover_ioproc_every,
+                              int IOrecover_unchunked);
+
+
+ /*@@
+   @routine    IOFlexIO_RecoverGH
+   @date       Fri Jun 19 09:22:52 1998
+   @author     Tom Goodale
+   @desc 
+   Recovers a GH.
+   @enddesc 
+   @calls IOUtil_PrepareFilename IOFlexIO_RestoreIEEEIOfile
+   @calledby   
+   @history 
+   @hauthor Gabrielle Allen
+   @hdate Thu Jul  2 18:17:59 1998
+   @hdesc Restore the physical time and iteration count, pass myproc to
+          IEEEIOparamRestore
+	  Derives the filename from IOUtil_PrepareFilename (in thorn IOUtil)
+   @hauthor Gabrielle Allen @hdate Oct 17 1998
+   @hdesc Added input of (some) GH structure variables
+   @endhistory 
+
+@@*/
+
+int IOFlexIO_RecoverGH (cGH *GH, const char *basename, int called_from)
+{
+#ifdef CACTUSBASE_PUGH
+
+  DECLARE_PARAMETERS
+  IOFile ifp;
+  char ftmp [1024], fname [1024];
+  int proc, nprocs, myproc;
+  int index;
+  int is_IEEEIO_file;
+  Long nels_stored;
+  int nt_stored;
+  int iteration_stored;
+  int file_ioproc, file_ioproc_every;
+  int file_nprocs, file_unchunked;
+  CCTK_INT4 tmpInt;
+  char msg [512];
+  pGH *pughGH;
+  ioGH *ioUtilGH;
+  cTimer total_time, dataset_time, param_time;
+#ifdef MPI
+  CCTK_INT info [3];
+#endif
+
+
+  /* Get the handles for PUGH and IOUtil extensions */
+  pughGH = (pGH *) GH->extensions [CCTK_GetGHExtensionHandle ("PUGH")];
+  ioUtilGH = (ioGH *) GH->extensions [CCTK_GetGHExtensionHandle ("IO")];
+
+  /* identify myself */
+  nprocs = CCTK_GetnProcs (GH);
+  myproc = CCTK_GetMyProc (GH);
+
+  /* initialize timers */
+  CactusResetTimer (&total_time);
+  CactusResetTimer (&dataset_time);
+  CactusResetTimer (&param_time);
+
+  CactusStartTimer (&total_time);
+
+  /* Examine base file to find whether recovering from
+   * one (recombined) file or from multiple files
+   */
+
+  if (myproc == 0) {
+    
+    /* Determine name of base file
+       NOTE: As we don't know whether the file is chunked or not
+             we need to try both file names. */
+    /* at first try with current chunking mode */
+    file_unchunked = ioUtilGH->unchunked;
+    IOUtil_PrepareFilename(GH, basename, fname, called_from, 0, file_unchunked);
+    if (called_from == CP_RECOVER_DATA)
+      strcat (fname, ".chkpt");
+    strcat (fname, ".ieee");
+
+    if (IO_verbose)
+      printf ("Opening file %s\n", fname);
+
+    /* Open file, make sure the file is valid */
+    ifp = IEEEopen (fname, "r");
+    if (IOisValid (ifp))
+      is_IEEEIO_file = 1;
+    else {
+      if (IO_verbose)
+        printf ("Cannot open file '%s'\n", fname);
+
+      /* now try with the other chunking mode */
+      file_unchunked = ! ioUtilGH->unchunked;
+      IOUtil_PrepareFilename (GH, basename, fname, called_from,
+                          0, file_unchunked);
+      if (called_from == CP_RECOVER_DATA)
+        strcat (fname, ".chkpt");
+      strcat (fname, ".ieee");
+      if (IO_verbose)
+        printf ("Trying now file '%s'...\n", fname);
+
+      /* Open file, make sure the file is valid */
+      ifp = IEEEopen (fname, "r");
+      is_IEEEIO_file = IOisValid (ifp);
+    }
+  }
+    
+  if (myproc == 0 && is_IEEEIO_file) {
+    /* Now determine how the data was written */
+
+    /* Read nioprocs used to write data */
+    index = IOreadAttributeInfo (ifp, "GH$ioproc_every", &nt_stored, &nels_stored);
+    if (index >= 0 && nt_stored == FLEXIO_INT4 && nels_stored == 1) {
+      IOreadAttribute (ifp, index, &tmpInt);
+      file_ioproc_every = tmpInt;
+    } else {
+      CCTK_WARN (1, "Unable to restore GH$ioproc_every. "
+                 "Assuming it is nprocs and continuing");
+      file_ioproc_every = nprocs;
+    }
+
+    /* Read nprocs used to write data */
+    index = IOreadAttributeInfo (ifp, "GH$nprocs", &nt_stored, &nels_stored);
+    if (index >= 0 && nt_stored == FLEXIO_INT4 && nels_stored == 1) {
+      IOreadAttribute (ifp, index, &tmpInt);
+      file_nprocs = tmpInt;
+    } else {
+      CCTK_WARN (1, "Unable to restore GH$nprocs. "
+                 "Assuming it is 1 and continuing");
+      file_nprocs = 1;
+    }
+    
+    /* Determine whether data is chunked or unchunked
+       We could derive this from the filename itself but just to be sure ... */
+    index = IOreadAttributeInfo (ifp, "unchunked", &nt_stored, &nels_stored);
+    if (index >= 0 && nt_stored == FLEXIO_INT4 && nels_stored == 1) {
+      IOreadAttribute (ifp, index, &tmpInt);
+      file_unchunked = tmpInt;
+    } else {
+      sprintf (msg, "Unable to restore 'unchunked' attribute. "
+               "Assuming it is %s and continuing",
+               file_unchunked ? "true" : "false");
+      CCTK_WARN (1, msg);
+    }
+
+    /* If we restore from multiple files
+     * the number of processors must match.
+     */
+ 
+    if (file_ioproc_every == nprocs || file_unchunked) {
+      if (IO_verbose)
+        printf ("Recovering from one %s file\n",
+                file_unchunked ? "unchunked" : "chunked");
+    } else {
+      if (file_nprocs != nprocs) {
+        sprintf (msg, "Must restart on %d processors with multiple files "
+                 "or recombine them", file_nprocs);
+        CCTK_WARN (0, msg);
+      }
+      if (IO_verbose)
+        printf ("Recovering from %d chunked files\n",
+                nprocs / file_ioproc_every + (nprocs % file_ioproc_every?1:0));
+    }
+  }
+  
+#ifdef MPI
+  /* Broadcast chunking mode to all processors from processor zero */
+  info [0] = is_IEEEIO_file;
+  info [1] = unchunked;
+  info [2] = file_ioproc_every;
+  CACTUS_MPI_ERROR (MPI_Bcast (info, 3, PUGH_MPI_INT, 0,
+                    pughGH->PUGH_COMM_WORLD));
+  is_IEEEIO_file = info [0];
+  unchunked = info [1];
+  file_ioproc_every = info [2];
+#endif
+
+  /* return here to IOUtil if no valid file could be found */
+  if (! is_IEEEIO_file) {
+    if (myproc == 0) {
+      sprintf (msg, "No valid IEEEIO file '%s' found !", fname);
+      CCTK_WARN (2, msg);
+    }
+    return (-1);
+  }
+
+  /* Determine the IO processors for each node and the corresponding
+     checkpoint file */
+  file_ioproc = myproc - (myproc % file_ioproc_every);
+  IOUtil_PrepareFilename (GH, basename, fname, called_from,
+                      file_ioproc/file_ioproc_every, file_unchunked);
+  if (called_from == CP_RECOVER_DATA)
+    strcat (fname, ".chkpt");
+  strcat (fname, ".ieee");
+
+  /* Open chunked files on other IO processors */
+  if (myproc != 0 && myproc == file_ioproc) {
+      
+    if (IO_verbose)
+      printf ("Opening chunked file '%s' on processor %d.\n",
+              fname, myproc);
+
+    /* Open file, make sure the file is valid */
+    ifp = IEEEopen (fname, "r");
+    if (! IOisValid (ifp)) {
+      sprintf (msg, "Cannot open file '%s' on processor %d",
+               fname, myproc);
+      CCTK_WARN (0, msg);
+    }
+  }
+ 
+  /* Restore the data */
+  if (IO_verbose && myproc == 0)
+    printf ("Recovering %schunked data with ioproc %d, ioproc_every %d.\n",
+            file_unchunked ? "un" : "", file_ioproc, file_ioproc_every);
+
+  CactusStartTimer (&dataset_time);
+  IOFlexIO_RestoreIEEEIOfile (GH, ifp, file_ioproc, file_ioproc_every,
+                            file_unchunked);
+  CactusStopTimer (&dataset_time);
+
+  /* Close the file. */
+  if (myproc == file_ioproc) {
+    if (IO_verbose)
+      printf ("Closing file '%s' after recovery.\n", fname);
+    IOclose (ifp);
+  }
+ 
+  if (called_from == CP_RECOVER_DATA) {
+    /* Must read in parameters and scalars on all processors. */
+    CactusStartTimer (&param_time);
+    for (proc = file_ioproc;
+         proc < file_ioproc+file_ioproc_every && proc < nprocs;
+         proc++) {
+    
+      /* Only have the file open by one proc at any time. */
+      if (proc == myproc) {
+      
+        /* Open file, make sure the file is valid */
+        ifp = IEEEopen (fname, "r");
+        if (! IOisValid (ifp)) {
+          sprintf (msg, "Cannot open checkpoint file '%s' on processor %d",
+                   fname, myproc);
+          CCTK_WARN (0, msg);
+        }
+
+        /* Restore the parameters. */
+        if (IO_verbose)
+          printf ("Recovering parameters on processor %d.\n", myproc);
+/*** FIXME ***/
+#if	0
+        IO_IEEEIOparamRestore (ifp, myproc);
+      
+        /* Restore the structure variables. */
+        if (IO_verbose)
+          printf ("Recovering GH variables.\n");
+        IO_IEEEIOStructRestore (GH, ifp);
+#endif
+
+        /* Restore global variables */
+      
+        /* Get the iteration number. */
+        if (IO_verbose)
+          printf ("Recovering iteration number.\n");
+        index = IOreadAttributeInfo (ifp, "iteration", &nt_stored, &nels_stored);
+	
+        if (index >= 0 && nt_stored == FLEXIO_INT4 && nels_stored == 1) {
+          IOreadAttribute (ifp, index, &iteration_stored);
+          if (IO_verbose)
+            printf ("Iteration number is %d\n", (int) iteration_stored);
+          GH->iteration = iteration_stored;
+        } else
+          printf ("*Warning* Unable to restore iteration number\n");
+      
+        /* Close the file. */
+        if (IO_verbose)
+          printf ("Closing '%s' after recovery.\n", fname);
+        IOclose (ifp);
+      }
+    
+      /* Synchronise all processors */
+      CCTK_Barrier (GH);
+    }
+    CactusStopTimer (&param_time);
+  }
+
+  /* print timing output */
+  if (IO_verbose && called_from == CP_RECOVER_DATA && myproc == 0) {
+    printf (
+      "----------------------------------------------------------------\n");
+/*** FIXME: choose right component of basic[] ***/
+    printf ("Time to restore data:            %10.3lf sec\n",
+            dataset_time.total.basic [0]);
+    printf ("Time to restore parameters:      %10.3lf sec\n",
+            param_time.total.basic [0]);
+    printf ("Time to recover from checkpoint: %10.3lf sec\n",
+            total_time.total.basic [0]);
+  }
+
+  return (0);
+#endif /* CACTUSBASE_PUGH */
+}
diff --git a/src/RestoreFile.c b/src/RestoreFile.c
new file mode 100644
index 0000000..60fb56a
--- /dev/null
+++ b/src/RestoreFile.c
@@ -0,0 +1,520 @@
+ /*@@
+   @file      RestoreFile.c
+   @date      Thu Jun 18 16:34:59 1998
+   @author    Tom Goodale
+   @desc 
+   Routines to restore GFs from a given file
+   @enddesc 
+   @history
+   @hauthor Gabrielle Allen @hdate 19 Oct 1998      
+   @hdesc Changed names ready for thorn_IO
+   @endhistory
+   @version $Id$
+ @@*/
+
+static char *rcsid = "$Id$";
+
+#include <stdio.h>
+#include <stdlib.h>
+
+#include "cctk.h"
+#include "flesh.h"
+#include "Groups.h"
+#include "Misc.h"
+#include "Comm.h"
+#include "GHExtensions.h"
+#include "WarnLevel.h"
+#include "declare_parameters.h"
+#ifdef CACTUSBASE_PUGH
+#include "CactusBase/pugh/src/include/pugh.h"
+#endif
+#include "CactusBase/IOUtil/src/ioGH.h"
+#include "ioFlexGH.h"
+
+
+/* MPI tag base */
+#define STARTUPBASE 1001        
+/* the maximum number of dimensions we can deal with (up to now :-) */
+#define	MAXDIM      3
+
+
+/* local routine getDatasetAttributes() reads in the next dataset's attributes
+   and verifies them:
+
+  * checks if there is a variable with the name given by the name attribute
+  * verifies that this variable still belongs to the same group
+  * checks the group data info:
+    - group type
+    - variable type
+    - ntimelevels
+    - sizes (rank, dimensions) according to chunking mode
+
+ If there is a mismatch a warning (warning level 2) is printed and value of 1
+ is returned to indicate that this dataset should be ignored.
+ If successful, the global variable index, the group type and the timelevel
+ to restore are stored in {*index, *gtype, *timelevel}, and 0 is returned.
+
+ Currently only restoring of CCTK_VARIABLE_REAL variables is supported.
+*/
+
+int GetCommonAttributes (cGH *GH, IOFile ifp, int unchunked, int *index,
+                         int *gtype, int *timelevel)
+{
+  int i, flag;
+  int atype;
+  Long asize;
+  int vtype_stored, rank_stored, dims_stored [MAXDIM];
+  CCTK_INT4 gtype_stored, ntimelevels_stored, timelevel_stored;
+  int vtype, rank, ntimelevels;
+  int result;
+  char fullname [512], groupname [512];
+  char msg [160];
+
+
+  /* read the next dataset's info from the file */
+  result = IOreadInfo (ifp, &vtype_stored, &rank_stored, dims_stored, MAXDIM);
+  CACTUS_IEEEIO_ERROR (result);
+  if (result == 0) {
+    CCTK_WARN (1, "Can't read dataset info");
+    return (1);
+  }
+      
+  /* retrieve the name attribute */
+  i = IOreadAttributeInfo (ifp, "name", &atype, &asize);
+  if (i < 0 || (atype != BYTE && atype != CHAR) || asize >= sizeof (fullname)) {
+    CCTK_WARN (2, "Can't read name attribute");
+    return (1);
+  }
+  CACTUS_IEEEIO_ERROR (IOreadAttribute (ifp, i, fullname));
+
+  /* check if there is a matching variable */
+  *index = CCTK_GetVarIndex (fullname);
+  if (*index < 0) {
+    sprintf (msg, "No matching variable found for '%s'", fullname);
+    CCTK_WARN (2, msg);
+    return (1);
+  }
+
+  /* read and verify the group name */
+  i = IOreadAttributeInfo (ifp, "groupname", &atype, &asize);
+  if (i < 0 || (atype != BYTE && atype != CHAR) || asize >= sizeof (groupname)) {
+    sprintf (msg, "Can't read groupname attribute of '%s'", fullname);
+    CCTK_WARN (2, msg);
+    return (1);
+  }
+  CACTUS_IEEEIO_ERROR (IOreadAttribute (ifp, i, groupname));
+  if (! CCTK_Equals (groupname, CCTK_GetGroupNameFromVar_ByIndex (*index))) {
+    sprintf (msg, "Groupnames don't match for '%s'", fullname);
+    CCTK_WARN (2, msg);
+    return (1);
+  }
+
+  /* read the group type */
+  i = IOreadAttributeInfo (ifp, "grouptype", &atype, &asize);
+  if (i < 0 || atype != FLEXIO_INT4 || asize != 1) {
+    sprintf (msg, "Can't read grouptype attribute for '%s'", fullname);
+    CCTK_WARN (2, msg);
+    return (1);
+  }
+  CACTUS_IEEEIO_ERROR (IOreadAttribute (ifp, i, &gtype_stored));
+
+  /* read the number of timelevels */
+  i = IOreadAttributeInfo (ifp, "ntimelevels", &atype, &asize);
+  if (i < 0 || atype != FLEXIO_INT4 || asize != 1) {
+    sprintf (msg, "Can't read ntimelevels attribute for '%s'", fullname);
+    CCTK_WARN (2, msg);
+    return (1);
+  }
+  CACTUS_IEEEIO_ERROR (IOreadAttribute (ifp, i, &ntimelevels_stored));
+
+  /* read the timelevel to restore */
+  i = IOreadAttributeInfo (ifp, "timelevel", &atype, &asize);
+  if (i < 0 || atype != FLEXIO_INT4 || asize != 1) {
+    sprintf (msg, "Can't read timelevel attribute for '%s'", fullname);
+    CCTK_WARN (2, msg);
+    return (1);
+  }
+  CACTUS_IEEEIO_ERROR (IOreadAttribute (ifp, i, &timelevel_stored));
+  *timelevel = timelevel_stored;
+
+  /* verify group type, variable type, dims, sizes and ntimelevels */
+  CCTK_GetGroupData (CCTK_GetGroupIndex (groupname), gtype, &vtype, &rank,
+                     &i, &ntimelevels);
+  if (*gtype != gtype_stored) {
+    sprintf (msg, "Group types don't match for '%s'", fullname);
+    CCTK_WARN (2, msg);
+    return (1);
+  }
+  if (ntimelevels != ntimelevels_stored) {
+    sprintf (msg, "Number of timelevels don't match for '%s'", fullname);
+    CCTK_WARN (2, msg);
+    return (1);
+  }
+  /* The CCTK variable type defines do not correlate with the IEEEIO defines
+     so compare them explicitely here.
+     Up to now CCTK knows only REALs, but COMPLEX types aren't supported
+     anyway by IOFlexIO. */
+  if ((vtype_stored == FLEXIO_REAL && vtype == CCTK_VARIABLE_REAL) ||
+      (vtype_stored == FLEXIO_INT  && vtype == CCTK_VARIABLE_INTEGER) ||
+      (vtype_stored == FLEXIO_CHAR && vtype == CCTK_VARIABLE_CHAR)) {
+#if 0
+    || (vtype_stored == FLEXIO_COMPLEX && vtype == CCTK_VARIABLE_COMPLEX)) {
+#endif
+    /* everything is okay */
+  } else {
+    sprintf (msg, "Variable types don't match for '%s'", fullname);
+    CCTK_WARN (2, msg);
+    return (1);
+  }
+
+  /* verify the dims and sizes */
+  flag = 0;
+  if (rank != rank_stored)
+    flag = 1;
+  switch (*gtype) {
+    case GROUP_SCALAR:
+      if (dims_stored [0] != 1)
+        flag = 1;
+      break;
+    case GROUP_ARRAY:
+      /*** FIXME: what's the local size of the array ?? ***/
+      /* for (i = 0; i < rank; i++)
+        if (*(int *) CCTK_ArrayGroupSize (GH, groupname, i) != dims_stored [i])
+          flag = 1;
+      */
+      break;
+    case GROUP_GF:
+      if (unchunked) {
+        for (i = 0; i < rank; i++)
+          if (GH->global_shape [i] != dims_stored [i])
+            flag = 1;
+      } else {
+        for (i = 0; i < rank; i++)
+          if (GH->local_shape [i] != dims_stored [i])
+            flag = 1;
+      }
+      break;
+  }
+  if (flag) {
+    sprintf (msg, "Variable sizes don't match for '%s'", fullname);
+    CCTK_WARN (2, msg);
+    return (1);
+  }
+
+  if (! CCTK_QueryGroupStorage_ByIndex (GH, 
+             CCTK_GetGroupIndexFromVar_ByIndex (*index))) {
+    CCTK_WARN (2, "Can't read into variable with no storage");
+    return (1);
+  }
+
+  return (0);
+}
+
+
+/* local routine GetChunkAttributes() reads in the info of the next dataset
+   that should be a chunk.
+   It verifies via the name attribute that this chunk belongs to the
+   current variable given by its index. */
+
+int GetChunkAttributes (cGH *GH, IOFile ifp, int index)
+{
+  int i;
+  int atype;
+  Long asize;
+  int result;
+  int vtype_stored, rank_stored, dims_stored [MAXDIM];
+  char fullname [512];
+  char msg [160];
+  pGH *pughGH;               /* PUGH extension handle */
+
+
+  /* Get the handle for PUGH extensions */
+  pughGH = (pGH *) GH->extensions [CCTK_GetGHExtensionHandle ("PUGH")];
+
+  /* read the next dataset's info from the file */
+  result = IOreadInfo (ifp, &vtype_stored, &rank_stored, dims_stored, MAXDIM);
+  CACTUS_IEEEIO_ERROR (result);
+  if (result == 0) {
+    CCTK_WARN (1, "Can't read dataset info");
+    return (1);
+  }
+      
+  /* retrieve the name attribute */
+  i = IOreadAttributeInfo (ifp, "name", &atype, &asize);
+  if (i < 0 || (atype != BYTE && atype != CHAR) || asize >= sizeof (fullname)) {
+    CCTK_WARN (2, "Can't read name attribute");
+    return (1);
+  }
+  CACTUS_IEEEIO_ERROR (IOreadAttribute (ifp, i, fullname));
+
+  /* check if there is a matching variable */
+  if (index != CCTK_GetVarIndex (fullname)) {
+    sprintf (msg, "No matching variable found for '%s'", fullname);
+    CCTK_WARN (2, msg);
+    return (1);
+  }
+
+  return (0);
+}
+
+ /*@@
+   @routine    IOFlexIO_RestoreIEEEIOfile
+   @date       Fri Jun 19 09:19:48 1998
+   @author     Tom Goodale
+   @desc 
+   Reads in data from an open IEEEIO file.
+   Code was originally in StartupReader.
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+   @hauthor Gabrielle Allen @hdate Oct 17 1998
+   @hdesc Changed logic so that cactus stops if any of the dimensions of the 
+          input file and the current cactus run differ.
+   @hauthor Thomas Radke @hdate May 05 1999
+   @hdesc Added parameter unchunked
+   @endhistory 
+
+@@*/
+int IOFlexIO_RestoreIEEEIOfile (cGH *GH, IOFile ifp, int file_ioproc,
+                              int file_ioproc_every, int file_unchunked)
+{
+#ifdef CACTUSBASE_PUGH
+
+  DECLARE_PARAMETERS
+  int index, gtype;
+  int myproc, nprocs;
+  int nDatasets, currentDataset;
+  int timelevel;             /* current timelevel to be restored */
+  pGH *pughGH;               /* PUGH extension handle */
+  char msg [512];            /* just a message buffer */
+#ifdef MPI
+  int proc;
+  CCTK_INT info [2];         /* communication buffer for MPI */
+#endif
+
+
+  /* Get the handles for PUGH and IO extensions */
+  pughGH = (pGH *) GH->extensions [CCTK_GetGHExtensionHandle ("PUGH")];
+
+  myproc = CCTK_GetMyProc (GH);
+  nprocs = CCTK_GetnProcs (GH);
+
+  /* all IO procs determine the number of datasets in their checkpoint files */
+  if (myproc == file_ioproc) {
+
+    /* Get the number of sets */
+    nDatasets = IOnDatasets (ifp); 
+    if (IO_verbose)
+      printf ("  Input file has %d datasets\n", nDatasets);
+
+  }
+
+  /* In Cactus 3.x we had only datasets containing grid function data.
+     This distributed data was stored as one dataset per processor within
+     the group belonging to one IO processor. So there should be
+     nGF*ioproc_every datasets within each checkpoint file.
+
+     This consistency condition is no longer true
+     because now there might be datasets containing a SCALAR grouptype.
+     These datasets are stored only from the IO processors,
+     and they are just distributed to the non-IO processors again
+     during recovery.
+
+  if (nDatasets % file_ioproc_every != 0)
+    CCTK_WARN (0, "Number of datasets isn't a multiple of nioprocs");
+  */
+
+  /* Now process the datasets.
+     All IO processors read the datasets from their checkpoint file
+     verify their contents and communicate them to the non-IO processors. */
+
+  /* At first the code for the IO processors ... */
+  if (myproc == file_ioproc) {
+
+    /* Seek here once to the beginning of the file, the file pointer
+       is advanced then implicitely by subsequent calls to IOreadInfo() */
+    CACTUS_IEEEIO_ERROR (IOseek (ifp, 0));
+
+    /* Each IO processor loops over all available datasets, checks their
+       consistency and broadcasts them to the non-IO processors. */
+    for (currentDataset = 0; currentDataset < nDatasets; currentDataset++) {
+
+      /* read in the next dataset's attributes and verify them */
+      if (GetCommonAttributes (GH, ifp, file_unchunked, &index, &gtype, &timelevel)) { 
+        sprintf (msg, "Ignoring dataset %d", currentDataset);
+        CCTK_WARN (1, msg);
+      }
+
+      /*** FIXME: process all group types ! ***/
+      else if (gtype != GROUP_GF && gtype != GROUP_SCALAR) {
+        CCTK_WARN (1, "Currently only restoring of GF and SCALAR datasets is supported");
+      }
+
+      /* Read in the data */
+      else {
+
+        if (IO_verbose) {
+          char *varname = CCTK_GetFullName (index);
+
+          printf ("    dataset %d: %s (timelevel %d)\n", currentDataset,
+                  varname, timelevel);
+          free (varname);
+        }
+
+        if (file_ioproc_every == 1) 
+          CACTUS_IEEEIO_ERROR (IOread (ifp, GH->data [index][timelevel]));
+#ifdef MPI
+        else {
+          int npoints;
+          void *buffer;
+          int chunkdims [3], *chunkorigin;
+          int element_size, mpi_type;
+
+          switch (CCTK_GetVarVType (index)) {
+            case CCTK_VARIABLE_CHAR:
+              element_size = sizeof (CCTK_CHAR); mpi_type = PUGH_MPI_CHAR;
+              break;
+            case CCTK_VARIABLE_INTEGER:
+              element_size = sizeof (CCTK_INT); mpi_type = PUGH_MPI_INT;
+              break;
+            case CCTK_VARIABLE_REAL:
+              element_size = sizeof (CCTK_REAL); mpi_type = PUGH_MPI_REAL;
+              break;
+            case CCTK_VARIABLE_COMPLEX:
+              CCTK_WARN (1, "Restoring of complex datatypes not yet supported");
+              continue;
+            default:
+              CCTK_WARN (1, "Unknown datatype in IOFlexIO_RestoreIEEEIOfile");
+              continue;
+          }
+
+          /* read my own data directly into data,
+             read others data into buffer and communicate it */
+          if (! file_unchunked || gtype == GROUP_SCALAR)
+            CACTUS_IEEEIO_ERROR (IOread (ifp, GH->data [index][timelevel]));
+          else {
+            chunkdims [0] = pughGH->rnx [file_ioproc];
+            chunkdims [1] = pughGH->rny [file_ioproc];
+            chunkdims [2] = pughGH->rnz [file_ioproc];
+            chunkorigin   = pughGH->lb [file_ioproc];
+
+            CACTUS_IEEEIO_ERROR (IOreadChunk (ifp, chunkdims, chunkorigin,
+                                 GH->data [index][timelevel]));
+          }
+
+          /* read data for non-IO processors */
+          if (gtype == GROUP_SCALAR) {
+            npoints = 1;
+            buffer = GH->data [index][timelevel];
+          } else {
+            /* allocate memory for the biggest chunk */
+            npoints = pughGH->rnpoints [file_ioproc + 1];
+            for (proc = 2; proc < file_ioproc_every; proc++)
+              if (npoints < pughGH->rnpoints [file_ioproc + proc])
+                npoints = pughGH->rnpoints [file_ioproc + proc];
+            buffer = malloc (npoints * element_size);
+          }
+
+          for (proc = file_ioproc + 1;
+               proc < file_ioproc + file_ioproc_every && proc < nprocs;
+               proc++) {
+
+            if (gtype != GROUP_SCALAR) {
+
+              if (! unchunked) {
+                /* Also increment dataset counter here !!! */
+                currentDataset++;
+                if (GetChunkAttributes (GH, ifp, index)) { 
+                  sprintf (msg, "Ignoring chunk in dataset %d",
+                           currentDataset+1);
+                  CCTK_WARN (1, msg);
+                  continue;
+                }
+
+                CACTUS_IEEEIO_ERROR (IOread (ifp, buffer));
+              } else {
+                chunkdims [0] = pughGH->rnx [proc];
+                chunkdims [1] = pughGH->rny [proc];
+                chunkdims [2] = pughGH->rnz [proc];
+                chunkorigin   = pughGH->lb [proc];
+
+                CACTUS_IEEEIO_ERROR (IOreadChunk (ifp, chunkdims, chunkorigin,
+                                     buffer));
+              }
+
+              npoints = pughGH->rnpoints [proc];
+            }
+
+            /* and finally send the index and the data */
+            info [0] = index; info [1] = timelevel;
+            CACTUS_MPI_ERROR (MPI_Send (info, 2, PUGH_MPI_INT, proc,
+                              STARTUPBASE, pughGH->PUGH_COMM_WORLD));
+            CACTUS_MPI_ERROR (MPI_Send (buffer, npoints, mpi_type, proc,
+                              STARTUPBASE, pughGH->PUGH_COMM_WORLD));
+
+          }
+
+          if (gtype != GROUP_GF)
+            free (buffer);
+        }
+#endif
+      } /* reading data for file_ioproc_every processors */
+
+    } /* end of loop over all datasets */
+
+#ifdef MPI
+    /* Finally an invalid variable index is communicated
+       to indicate completion to the non-IO processors. */
+    info [0] = -1;
+    for (proc = 1; proc < file_ioproc_every; proc++)
+      CACTUS_MPI_ERROR (MPI_Send (info, 2, PUGH_MPI_INT, proc + file_ioproc,
+                        STARTUPBASE, pughGH->PUGH_COMM_WORLD));
+#endif
+
+  } else {
+
+    /* And here the code for non-IO processors: */
+#ifdef MPI
+    int mpi_type;
+    MPI_Status ms;
+
+    /* They don't know how many datasets there are, because the IO processors
+       could skip some on the fly during their consistency checks.
+       The IO Processor sends the index of the variable to be processed next.
+       So, all non-IO processors execute a loop where the termination condition
+       is when an invalid index was received.
+    */
+    while (1) {
+      /* receive the next variable index from the IO processor */
+      CACTUS_MPI_ERROR (MPI_Recv (info, 2, PUGH_MPI_INT, file_ioproc,
+                        STARTUPBASE, pughGH->PUGH_COMM_WORLD, &ms));
+      index = info [0]; timelevel = info [1];
+
+      /* check for termination condition */
+      if (index < 0)
+        break;
+
+      switch (CCTK_GetVarVType (index)) {
+        case CCTK_VARIABLE_CHAR:    mpi_type = PUGH_MPI_CHAR; break;
+        case CCTK_VARIABLE_INTEGER: mpi_type = PUGH_MPI_INT;  break;
+        case CCTK_VARIABLE_REAL:    mpi_type = PUGH_MPI_REAL; break;
+        case CCTK_VARIABLE_COMPLEX:
+          CCTK_WARN (1, "Restoring of complex datatypes not yet supported");
+          continue;
+        default:
+          CCTK_WARN (1, "Unknown datatype in IOFlexIO_RestoreIEEEIOfile");
+          continue;
+      }
+
+      /* receive following data from my IO processor */
+      CACTUS_MPI_ERROR (MPI_Recv (GH->data [index][timelevel], pughGH->npoints,
+                        mpi_type, file_ioproc,
+                        STARTUPBASE, pughGH->PUGH_COMM_WORLD, &ms));
+    }
+#endif
+  }
+
+  return (0);
+#endif /* CACTUSBASE_PUGH */
+}
diff --git a/src/Startup.c b/src/Startup.c
new file mode 100644
index 0000000..a9d5113
--- /dev/null
+++ b/src/Startup.c
@@ -0,0 +1,78 @@
+ /*@@
+   @file      Startup.c
+   @date      Fri May 21 1999
+   @author    Thomas Radke
+   @desc 
+   Startup routines for IOFlexIO.
+   @enddesc 
+   @history
+   @hauthor Thomas Radke @hdate May 21 1999
+   @hdesc Just copied from thorn IOFlexIO.
+   @endhistory
+ @@*/
+
+#include <stdio.h>
+
+#include "flesh.h"
+#include "GHExtensions.h"
+#include "declare_parameters.h"
+
+/* prototypes of functions to be registered */
+int IOFlexIO_Output2DGH (cGH *GH);
+int IOFlexIO_TriggerOutput2D (cGH *GH, int);
+int IOFlexIO_TimeFor2D (cGH *GH, int);
+int IOFlexIO_Output2DVarAs (cGH *GH, const char *var, const char *alias);
+int IOFlexIO_Output3DGH (cGH *GH);
+int IOFlexIO_TriggerOutput3D (cGH *GH, int);
+int IOFlexIO_TimeFor3D (cGH *GH, int);
+int IOFlexIO_Output3DVarAs (cGH *GH, const char *var, const char *alias);
+void *IOFlexIO_SetupGH (tFleshConfig *config, int convergence_level, cGH *GH);
+int IOFlexIO_InitGH (cGH *GH);
+int IOFlexIO_rfrTraverseGH (cGH *GH, int rfrpoint);
+
+
+
+ /*@@
+   @routine   IOFlexIO_Startup
+   @date      Fri May 21 1999
+   @author    Thomas Radke
+   @desc 
+   The startup registration routine for IOFlexIO.
+   Registers the GH extensions needed for IOFlexIO and
+   the registerable routines used for each method of IOFlexIO.
+   IOFlexIO does not overload any functions.
+   @enddesc 
+   @calls     
+   @calledby   
+   @history 
+ 
+   @endhistory 
+
+@@*/
+void IOFlexIO_Startup (void)
+{
+  DECLARE_PARAMETERS
+  int IO_GHExtension;
+  int IOMethod;
+
+  IO_GHExtension = CCTK_RegisterGHExtension ("IOFlexIO");
+  CCTK_RegisterGHExtensionSetupGH (IO_GHExtension, IOFlexIO_SetupGH);
+  CCTK_RegisterGHExtensionInitGH (IO_GHExtension, IOFlexIO_InitGH);
+  CCTK_RegisterGHExtensionrfrTraverseGH (IO_GHExtension, IOFlexIO_rfrTraverseGH);
+
+  /* Register the 2D and 3D IOFlexIO routines as output methods  */
+  if (strlen (output2D) > 0) {
+    IOMethod = CCTK_RegisterIOMethod ("IOFlexIO_2D");
+    CCTK_RegisterIOMethodOutputGH (IOMethod, IOFlexIO_Output2DGH);
+    CCTK_RegisterIOMethodOutputVarAs (IOMethod, IOFlexIO_Output2DVarAs);
+    CCTK_RegisterIOMethodTimeToOutput (IOMethod, IOFlexIO_TimeFor2D);
+    CCTK_RegisterIOMethodTriggerOutput (IOMethod, IOFlexIO_TriggerOutput2D);
+  }
+  if (strlen (output3D) > 0) {
+    IOMethod = CCTK_RegisterIOMethod ("IOFlexIO_3D");
+    CCTK_RegisterIOMethodOutputGH (IOMethod, IOFlexIO_Output3DGH);
+    CCTK_RegisterIOMethodOutputVarAs (IOMethod, IOFlexIO_Output3DVarAs);
+    CCTK_RegisterIOMethodTimeToOutput (IOMethod, IOFlexIO_TimeFor3D);
+    CCTK_RegisterIOMethodTriggerOutput (IOMethod, IOFlexIO_TriggerOutput3D);
+  }
+}
diff --git a/src/Write2D.c b/src/Write2D.c
new file mode 100644
index 0000000..c708440
--- /dev/null
+++ b/src/Write2D.c
@@ -0,0 +1,479 @@
+/*@@
+   @file      Write2D.c
+   @date      Tue Feb  3 19:06:12 1998
+   @author    Paul Walker
+   @desc 
+   Two dimensional slices into IEEEIO files on MPI distributed cactus data.
+   @enddesc 
+   @history
+   @hauthor Gabrielle Allen @hdate 26 Sep 1998
+   @hdesc Changed subroutine name (SliceTwoD -> Write2D), renamed IO parameters
+   @hauthor Gabrielle Allen @hdate Oct 17 1998
+   @hdesc Added the IO_verbose parameter instead of an ifdef
+   @hauthor Gabrielle Allen @hdate 19 Oct 1998      
+   @hdesc Changed names ready for thorn_IO
+   @hauthor Thomas Radke @hdate 16 Mar 1999
+   @hdesc Converted to Cactus 4.0
+   @hauthor Thomas Radke @hdate 17 Mar 1999
+   @hdesc Changed naming: IEEEIO -> IOFlexIO
+   @hendhistory
+   @version $Id$
+ @@*/
+
+static char *rcsid = "$Id$";
+
+#include <stdio.h>
+#include <assert.h>
+#include <stdlib.h>
+
+#include "cctk.h"
+#include "flesh.h"
+#include "Comm.h"
+#include "declare_parameters.h"
+#include "GroupsOnGH.h"
+#ifdef CACTUSBASE_PUGH
+#include "CactusBase/pugh/src/include/pugh.h"
+#endif
+#include "CactusBase/IOUtil/src/ioGH.h"
+#include "ioFlexGH.h"
+
+/* MPI message tag */
+#define BASE 1234
+
+
+/*@@
+   @routine    IOFlexIO_Write2D
+   @date       Sat March 6 1999
+   @author     Gabrielle Allen
+   @desc       Writes the 2D data of a variable into separate IEEEIO files
+   @enddesc
+   @calls    CCTK_GetGHExtensionHandle
+             CCTK_GetMyProc
+             CCTK_GetnProcs
+             CCTK_WARN
+   @calledby IOFlexIO_Output2DVarAs
+             IOFlexIO_TriggerOutput2D
+   @history
+   @endhistory
+   @var     GH
+   @vdesc   Pointer to CCTK GH
+   @vtype   cGH
+   @vio     in
+   @vcomment
+   @endvar
+   @var     index
+   @vdesc   index of variable to output
+   @vtype   int
+   @vio     in
+   @vcomment
+   @endvar
+   @var     alias
+   @vdesc   alias name of variable to output
+   @vtype   const char *
+   @vio     in
+   @vcomment
+   @endvar
+   @var     first
+   @vdesc   flag to indicate whether file should be opened in "w" or "a" mode
+   @vtype   int (boolean)
+   @vio     in
+   @vcomment
+   @endvar
+@@*/
+
+void IOFlexIO_Write2D (cGH *GH, int index, const char *alias, int first)
+{
+  DECLARE_PARAMETERS
+  int i, j, k, l;
+  int myproc, nprocs;
+  pGH *pughGH;
+  ioGH *ioUtilGH;
+  flexioGH *myGH;
+  int dir, ngpoints;
+  CCTK_INT npoints, *nrempoints, bnd [4];
+  void *data;
+  int ni, nj;                        /* "x" and "y" direction ns */
+  /* Attributes */
+  CCTK_REAL origin [2], delta [2];
+  int maxrempt;
+  int dims [2];
+  void *locdat, *alldat, *recd;
+  int hi, lo;
+  int locbnd [4];
+  int timelevel, variable_type, element_size, ioflex_type;
+#ifdef MPI
+  MPI_Status ms;
+  int mpi_type;
+#endif
+
+
+  /* Get the handles for pugh, IOFlexIO and IOUtil extensions */
+  pughGH = (pGH *) GH->extensions [CCTK_GetGHExtensionHandle ("PUGH")];
+  ioUtilGH = (ioGH *) GH->extensions [CCTK_GetGHExtensionHandle ("IO")];
+  myGH = (flexioGH *) GH->extensions [CCTK_GetGHExtensionHandle ("IOFlexIO")];
+
+  /* What processor are we on? */
+  myproc = CCTK_GetMyProc (GH);
+  nprocs = CCTK_GetnProcs (GH);
+
+  if (IO_verbose) {
+    printf ("<><><><><><><><><><><><><><><><><>\n");
+    printf ("2D Slice for [%s]\n", alias);
+#if        0
+    cactus_memtrace ();
+#endif
+  }
+
+  /* Open the files on the first trip through if we are proc. 0 */
+  if (myproc == 0) {
+    if (first) {
+      int len;
+      char *fname;
+      char *msg;
+
+      len = strlen (ioUtilGH->outpfx_2D) + strlen (alias);
+      fname = (char *) malloc (len + 20);
+      msg = (char *) malloc (len + 20 + 40);
+      assert (fname && msg);
+ 
+      /* Create files (open mode "w") */
+      sprintf (fname, "%s/%s_2d_yz.ieee", ioUtilGH->outpfx_2D, alias);
+      myGH->IEEEfile_2D [index][0] = IEEEopen (fname, "w");
+      if (! IOisValid (myGH->IEEEfile_2D [index][0])) {
+        sprintf (msg, "Cannot open output file %s", fname);
+        CCTK_WARN (1, msg);
+        return;
+      }
+      
+      sprintf (fname, "%s/%s_2d_xz.ieee", ioUtilGH->outpfx_2D, alias);
+      myGH->IEEEfile_2D [index][1] = IEEEopen (fname, "w");
+      if (! IOisValid (myGH->IEEEfile_2D [index][1])) {
+        sprintf (msg, "Cannot open output file %s", fname);
+        CCTK_WARN (1, msg);
+        return;
+      }
+      
+      sprintf (fname, "%s/%s_2d_xy.ieee", ioUtilGH->outpfx_2D, alias);
+      myGH->IEEEfile_2D [index][2] = IEEEopen (fname, "w");
+      if (! IOisValid (myGH->IEEEfile_2D [index][2])) {
+        sprintf (msg, "Cannot open output file %s", fname);
+        CCTK_WARN (1, msg);
+        return;
+      }
+
+      free (msg);
+      free (fname);
+    }
+  }
+
+  nrempoints = (CCTK_INT *) malloc (nprocs * sizeof (CCTK_INT));
+  assert (nrempoints);
+
+  timelevel = CCTK_GetNumTimeLevels (index) - 1;
+  if (timelevel > 0)
+    timelevel--;
+  data = CCTK_GetVarDataPtr_ByIndex (GH, timelevel, index);
+
+  variable_type = CCTK_GetVarVType (index);
+  switch (variable_type) {
+    case CCTK_VARIABLE_CHAR:
+      ioflex_type = FLEXIO_CHAR;
+      element_size = sizeof (CCTK_CHAR);
+#ifdef MPI
+      mpi_type = PUGH_MPI_CHAR;
+#endif
+      break;
+
+    case CCTK_VARIABLE_INTEGER:
+      ioflex_type = FLEXIO_INT;
+      element_size = sizeof (CCTK_INT);
+#ifdef MPI
+      mpi_type = PUGH_MPI_INT;
+#endif
+      break;
+
+    case CCTK_VARIABLE_REAL:
+      ioflex_type = FLEXIO_REAL;
+      element_size = ioUtilGH->out_single ?
+                     sizeof (CCTK_REAL4) : sizeof (CCTK_REAL);
+#ifdef MPI
+      mpi_type = ioUtilGH->out_single ? PUGH_MPI_REAL4 : PUGH_MPI_REAL;
+#endif
+      break;
+
+    default:
+      CCTK_WARN (1, "Unsupported variable type in IOFlexIO_Write2D");
+      return;
+  }
+
+  for (dir = 0; dir < 3; dir++) {
+    npoints = 0;
+    /* Figure out the number of points */
+    if ((pughGH->ownership[pughGH->stagger][dir][0] <= pughGH->sp2dxyz[dir]) &&
+        (pughGH->sp2dxyz[dir] < pughGH->ownership[pughGH->stagger][dir][1]))
+      npoints = (pughGH->ownership[pughGH->stagger][(dir+1)%3][1] -
+                 pughGH->ownership[pughGH->stagger][(dir+1)%3][0]) * 
+                (pughGH->ownership[pughGH->stagger][(dir+2)%3][1] -
+                 pughGH->ownership[pughGH->stagger][(dir+2)%3][0]);
+
+    /* Set up the number of global points and the points in the "i"
+       and "j" directions */
+    if (dir == 0) {
+      ngpoints = GH->global_shape[1] * GH->global_shape[2];
+      /* FIXME : origin */
+      ni = GH->global_shape[1]; origin[0] = 0.0; delta[0] = pughGH->dy0;
+      nj = GH->global_shape[2]; origin[1] = 0.0; delta[1] = pughGH->dz0;
+    } else if (dir == 1) {
+      ngpoints = GH->global_shape[0] * GH->global_shape[2];
+      /* FIXME : Origin */
+      ni = GH->global_shape[0]; origin[0] = 0.0; delta[0] = pughGH->dx0;
+      nj = GH->global_shape[2]; origin[1] = 0.0; delta[1] = pughGH->dz0;
+    } else {
+      ngpoints = GH->global_shape[0] * GH->global_shape[1];
+      /* FIXME : Origin */
+      ni = GH->global_shape[0]; origin[0] = 0.0; delta[0] = pughGH->dx0;
+      nj = GH->global_shape[1]; origin[1] = 0.0; delta[1] = pughGH->dy0;
+    }
+
+    if (IO_verbose) {
+      printf ("Npoints in dir %d is %d -> ",dir,npoints);
+      printf ("Global npoints is %d\n",ngpoints);
+    }
+
+#ifdef MPI
+    /* Send the number of points */
+    CACTUS_MPI_ERROR (MPI_Gather (&npoints, 1, PUGH_MPI_INT,
+                      nrempoints, 1, PUGH_MPI_INT, 0, pughGH->PUGH_COMM_WORLD));
+
+    if (IO_verbose && myproc == 0)
+      for (i = 0; i < nprocs; i++)
+        printf ("From proc %d: %d\n", i, (int) nrempoints [i]);
+
+#else
+    nrempoints [0] = npoints;
+#endif
+
+    maxrempt = 0;
+    if (myproc == 0)
+      for (i = 0; i < nprocs; i++)
+        maxrempt = (nrempoints [i] > maxrempt ? nrempoints [i] : maxrempt);
+
+    if (IO_verbose)
+      printf ("Maximum remote points: %d\n", maxrempt);
+
+    /* OK so we want the ordering in bnd to be yz xz or xy but the
+       mod naturally gives is yz zx xy so we need to flip the order.
+       thus put the dir+1 % 3 in at lo, the dir+2 % 3 in at hi, and
+       set lo and hi directionally as follows ... */
+    if (dir == 0) {lo = 0; hi = 2;}
+    if (dir == 1) {lo = 2; hi = 0;}
+    if (dir == 2) {lo = 0; hi = 2;}
+    if (npoints > 0) {
+      /* See above comment ... */
+      bnd [lo]   = GH->lower_bound[(dir+1)%3] + 
+                   pughGH->ownership[pughGH->stagger][(dir+1)%3][0];
+      bnd [lo+1] = GH->lower_bound[(dir+1)%3] + 
+                   pughGH->ownership[pughGH->stagger][(dir+1)%3][1];
+
+      bnd [hi]   = GH->lower_bound[(dir+2)%3] + 
+                   pughGH->ownership[pughGH->stagger][(dir+2)%3][0];
+      bnd [hi+1] = GH->lower_bound[(dir+2)%3] + 
+                   pughGH->ownership[pughGH->stagger][(dir+2)%3][1];
+
+      locbnd [lo]   = pughGH->ownership[pughGH->stagger][(dir+1)%3][0];
+      locbnd [lo+1] = pughGH->ownership[pughGH->stagger][(dir+1)%3][1];
+
+      locbnd [hi]   = pughGH->ownership[pughGH->stagger][(dir+2)%3][0];
+      locbnd [hi+1] = pughGH->ownership[pughGH->stagger][(dir+2)%3][1];
+
+    } else {
+      bnd [lo] = 0;
+      bnd [lo+1] = 0;
+      bnd [hi] = 0;
+      bnd [hi+1] = 0;
+    }
+
+    if (IO_verbose)
+      printf ("Local bounds %d: j -> %d %d  k -> %d %d\n",
+               dir, locbnd [0], locbnd [1], locbnd [2], locbnd [3]);
+
+    /* Subsample the data */
+    if (npoints > 0) {
+      locdat = malloc (npoints * element_size);
+      assert (locdat);
+      l = 0;
+      if (variable_type == CCTK_VARIABLE_CHAR) {
+        CCTK_CHAR *char_locdat = (CCTK_CHAR *) locdat;
+
+        for (k = locbnd [2]; k < locbnd [3]; k++)
+          for (j = locbnd [0]; j < locbnd [1]; j++) {
+            int pt;
+            if (dir==0) pt = DI (pughGH, pughGH->sp2dxyz [dir], j, k);
+            if (dir==1) pt = DI (pughGH, j, pughGH->sp2dxyz [dir], k);
+            if (dir==2) pt = DI (pughGH, j, k, pughGH->sp2dxyz [dir]);
+            assert (l <= npoints && pt <= pughGH->npoints);
+            char_locdat [l++] = ((CCTK_CHAR *) data) [pt];
+          }
+      } else if (variable_type == CCTK_VARIABLE_INTEGER) {
+        CCTK_INT *int_locdat = (CCTK_INT *) locdat;
+
+        for (k = locbnd [2]; k < locbnd [3]; k++)
+          for (j = locbnd [0]; j < locbnd [1]; j++) {
+            int pt;
+            if (dir==0) pt = DI (pughGH, pughGH->sp2dxyz [dir], j, k);
+            if (dir==1) pt = DI (pughGH, j, pughGH->sp2dxyz [dir], k);
+            if (dir==2) pt = DI (pughGH, j, k, pughGH->sp2dxyz [dir]);
+            assert (l <= npoints && pt <= pughGH->npoints);
+            int_locdat [l++] = ((CCTK_INT *) data) [pt];
+          }
+      } else if (variable_type == CCTK_VARIABLE_REAL) {
+        CCTK_REAL4 *real4_locdat = (CCTK_REAL4 *) locdat;
+        CCTK_REAL  *real_locdat = (CCTK_REAL *) locdat;
+
+        for (k = locbnd [2]; k < locbnd [3]; k++)
+          for (j = locbnd [0]; j < locbnd [1]; j++) {
+            int pt;
+            if (dir==0) pt = DI (pughGH, pughGH->sp2dxyz [dir], j, k);
+            if (dir==1) pt = DI (pughGH, j, pughGH->sp2dxyz [dir], k);
+            if (dir==2) pt = DI (pughGH, j, k, pughGH->sp2dxyz [dir]);
+            assert (l <= npoints && pt <= pughGH->npoints);
+            if (ioUtilGH->out_single)
+              real4_locdat [l++] = (CCTK_REAL4) (((CCTK_REAL *) data) [pt]);
+            else
+              real_locdat [l++] = ((CCTK_REAL *) data) [pt];
+          }
+      }
+    } 
+
+    /* Send the data */
+#ifdef MPI
+    /* Only send if we have points and we are not on proc 0 */
+    if (npoints > 0 && myproc != 0) {
+      CACTUS_MPI_ERROR (MPI_Send (bnd, 4, PUGH_MPI_INT, 0, 2*myproc+1+BASE,
+                        pughGH->PUGH_COMM_WORLD));
+      CACTUS_MPI_ERROR (MPI_Send (locdat, (int) npoints, mpi_type, 0,
+                        2*myproc+BASE, pughGH->PUGH_COMM_WORLD));
+    }
+
+    /* Receive the data on processor 0 */
+    if (myproc == 0) {
+      alldat = malloc (ngpoints * element_size);
+      recd   = malloc (maxrempt * element_size);
+      assert (alldat && recd);
+
+      for (i = 0; i < nprocs; i++) {
+        void *ud;
+
+        if (IO_verbose)
+          printf ("Proceesing proc %d\n",i);fflush(stdout);
+
+        if (nrempoints [i] > 0) {
+          if (i == 0)
+            ud = locdat;
+          else {
+
+            if (IO_verbose) {
+              printf ("Receiving %d from %d\n", nrempoints [i], i);
+              fflush (stdout);
+            }
+
+            CACTUS_MPI_ERROR (MPI_Recv (bnd, 4, PUGH_MPI_INT,
+                              i, 2*i+1+BASE, pughGH->PUGH_COMM_WORLD, &ms));
+            CACTUS_MPI_ERROR (MPI_Recv (recd, (int) nrempoints [i], mpi_type,
+                              i, 2*i+BASE, pughGH->PUGH_COMM_WORLD, &ms));
+            ud = recd;
+
+            if (IO_verbose) {
+              printf ("   Received %d points from %d\n",
+                      (int) nrempoints [i], i);
+              printf ("   Data lives at [%d %d] -> [%d %d]\n",
+                      (int) bnd [0], (int) bnd [2], (int) bnd [1], (int)bnd[3]);
+            }
+
+          } /* End proc != 0 */
+
+          if (IO_verbose)
+            printf ("Remote Bounds %d: j -> %d %d  k -> %d %d\n",
+                    dir, (int) bnd [0], (int) bnd[1], (int)bnd[2], (int)bnd[3]);
+
+          l = 0;
+          if (variable_type == CCTK_VARIABLE_CHAR)
+            for (k = bnd [2]; k < bnd [3]; k++)
+              for (j = bnd [0]; j < bnd [1]; j++)
+                ((CCTK_CHAR *) alldat) [j + k*ni] = ((CCTK_CHAR *) ud) [l++];
+          else if (variable_type == CCTK_VARIABLE_INTEGER)
+            for (k = bnd [2]; k < bnd [3]; k++)
+              for (j = bnd [0]; j < bnd [1]; j++)
+                ((CCTK_INT *) alldat) [j + k*ni] = ((CCTK_INT *) ud) [l++];
+          else if (variable_type == CCTK_VARIABLE_REAL && ioUtilGH->out_single)
+            for (k = bnd [2]; k < bnd [3]; k++)
+              for (j = bnd [0]; j < bnd [1]; j++)
+                ((CCTK_REAL4 *) alldat) [j + k*ni] = ((CCTK_REAL4 *) ud) [l++];
+          else if (variable_type == CCTK_VARIABLE_REAL && !ioUtilGH->out_single)
+            for (k = bnd [2]; k < bnd [3]; k++)
+              for (j = bnd [0]; j < bnd [1]; j++)
+                ((CCTK_REAL *) alldat) [j + k*ni] = ((CCTK_REAL *) ud) [l++];
+        } /* End if > 0 */
+      } /* End processor loop */
+      free (recd);
+    }
+#else
+    alldat = locdat;
+#endif
+
+    /* Proc 0 writes */
+    if (myproc == 0) {
+      dims [0] = ni;
+      dims [1] = nj;
+
+      /*** FIXME: since slice center is not yet set up, npoints might be 0 ***/
+      if (npoints <= 0)
+        CCTK_WARN (1, "2D output doesn't work since slice center is not set up");
+      else {
+        /* the data itself */
+        CACTUS_IEEEIO_ERROR (IOwrite (myGH->IEEEfile_2D [index][dir],
+                             ioflex_type, 2, dims, alldat));
+      
+        /* Time and space attributes */
+        CACTUS_IEEEIO_ERROR (IOwriteAttribute (myGH->IEEEfile_2D [index][dir],
+                             "time", FLEXIO_REAL, 1, &GH->time));
+        CACTUS_IEEEIO_ERROR (IOwriteAttribute (myGH->IEEEfile_2D [index][dir],
+                             "origin", FLEXIO_REAL, 2, origin));
+        CACTUS_IEEEIO_ERROR (IOwriteAttribute (myGH->IEEEfile_2D [index][dir],
+                             "delta", FLEXIO_REAL, 2, delta));
+      }
+    }
+
+    /* Free memory for this time through */
+    if (npoints > 0) 
+      free (locdat);
+#ifdef MPI
+    if (myproc == 0)
+      free (alldat);
+#endif /* If not MPI then alldat and locdat are not distinct */
+
+  }
+
+  free (nrempoints);
+
+#if        0
+  if (IO_verbose)
+    cactus_memtrace();
+#endif
+}
+ 
+
+
+#if        0
+void FMODIFIER FORTRAN_NAME(io_write2d_,
+                            IO_WRITE2D,
+                            io_write2d)
+     (CCTK_REAL *x)
+{
+  pGF *GF;
+
+  GF = locateGFbyDataPtr(x);
+  IO_Write2D(GF->parentGH,GF);
+}
+
+#endif
diff --git a/src/Write3D.c b/src/Write3D.c
new file mode 100644
index 0000000..6d72a93
--- /dev/null
+++ b/src/Write3D.c
@@ -0,0 +1,462 @@
+/*@@
+   @file      Write3D.c
+   @date      Fri Feb 21 15:27:03 1997
+   @author    
+   @desc 
+   A horrible blocking write and not-quite-so-horrible
+   slicer for xgraph output. Basically MPI based IO
+   bites the big one, and this is MPI based IO written
+   quickly and in the obvious inefficient fashion.
+   <p>
+   Oh, not it also contains zero-d I/O and a closer also.
+   <p>
+   This routine can write DFSD style HDF (not the default)
+   or Johns IEEEIO 3D files (the default). For 0 and 1D
+   I/O it writes xgraph files.
+   @enddesc 
+   @history
+   @hauthor Gabrielle Allen @hdate 26 Sep 1998
+   @hdesc Changed routine name (WriteGF->Write3D), Lots of tidying up and 
+          renaming of parameters
+   @hauthor Gabrielle Allen @hdate 17 Oct 1998
+   @hdesc Implemented new parameters IO_verbose,IO_parameter,IO_structures,IO_scalars
+   @hauthor Gabrielle Allen @hdate 19 Oct 1998      
+   @hdesc Changed names ready for thorn_IO
+   @hauthor Thomas Radke @hdate 16 Mar 1999
+   @hdesc Converted to Cactus 4.0
+   @hauthor Thomas Radke @hdate 17 Mar 1999
+   @hdesc Changed naming: IEEEIO -> IOFlexIO
+   @hendhistory
+   @version $Id$
+ @@*/
+
+static char *rcsid = "$Id$";
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <assert.h>
+
+#include "cctk.h"
+#include "flesh.h"
+#include "Comm.h"
+#include "Groups.h"
+#include "GHExtensions.h"
+#include "WarnLevel.h"
+#include "declare_parameters.h"
+#ifdef CACTUSBASE_PUGH
+#include "CactusBase/pugh/src/include/pugh.h"
+#endif
+#include "CactusBase/IOUtil/src/ioGH.h"
+#include "ioFlexGH.h"
+
+
+/* local function prototypes */
+void IOFlexIO_Write3D_filename (cGH *GH, char *filename, const char *name);
+void IOFlexIO_Write3D_openFile (cGH *GH, const char *filename, IOFile *iofile, int first);
+void IOFlexIO_Write3D_closeFile (cGH *GH, const char *filename, IOFile *iofile);
+/* external function prototypes */
+void IOFlexIO_DumpVar (cGH *GH, int index, int timelevel, IOFile iof);
+
+
+/*@@
+  @routine IOFlexIO_Write3D
+  @author  Paul Walker 
+  @date    Feb 1997
+  @desc
+  This routine does 3D IO for a grid function which has been
+  distributed with MPI.
+  <p>
+  3D MPI IO! Wow, you may say, that sounds pretty clever. Well it isn't
+  because <b><i>The method I use sucks rocks</i></b>. That is.
+  I spam the data from all processors onto processor 0. Luckily
+  it doesn't just put it all in one memory chunk, but recieves it
+  bit by bit and writes it out chunkwise, using IEEEIO. Alas, 
+  while this happens the other processors <i>sit in a barrier 
+  with their metaphoric thumbs up their asses</i>. 
+  <p>
+  <b>Even Worse</b> it is a blocking send which must come in order
+  from the other processors.
+  <p> 
+  This is inefficient to say the least. I/O costs can go up to 70-80%
+  using this technique on the T3E. So there are also a couple of ways to
+  get around it.
+  <p>
+  The easiest option is one-file-per-processor which is
+  good on some machines, (eg T3E) and bad on others (eg O2K, I think).
+  So you can activate this mode (or de-activate it on the T3E)
+  by setting the parameter onefileperproc to "yes" or "no".
+  This will create a bunch of files in their own subdirectory,
+  and they can be smooshed together using cactus/lib/etc/Recombiner.C
+  <p>
+  Also, you can downsample output. You can also output one file per
+  slice. Various combinations of these work or don't, depending.
+  <p>
+  Finally notice that @seeroutine main guarantees that this is never
+  called on a 1- or 2-D grid.
+  @enddesc
+  @calls    CCTK_GetGHExtensionHandle
+            CCTK_GetMyProc
+            CCTK_GetnProcs
+            CCTK_WARN
+  @calledby IOFlexIO_Output3DVarAs
+            IOFlexIO_TriggerOutput3D
+  @history 
+  @hdate Wed Sep  2 10:11:41 1998 @hauthor Tom Goodale
+  @hdesc Merged in Szu-Wen's Panda calls.
+  @hdate Sep 25 1998 1998 @hauthor Gabrielle Allen
+  @hdesc Split into subroutines
+  @endhistory 
+  @var     GH
+  @vdesc   Pointer to CCTK GH
+  @vtype   cGH
+  @vio     in
+  @vcomment
+  @endvar
+  @var     index
+  @vdesc   index of variable to output
+  @vtype   int
+  @vio     in
+  @vcomment
+  @endvar
+  @var     alias
+  @vdesc   alias name of variable to output
+  @vtype   const char *
+  @vio     in
+  @vcomment
+  @endvar
+  @var     first
+  @vdesc   flag to indicate whether file should be opened in "w" or "a" mode
+  @vtype   int (boolean)
+  @vio     in
+  @vcomment
+  @endvar
+  @@*/
+
+
+void IOFlexIO_Write3D (cGH *GH, int index, const char *alias, int first)
+{
+  DECLARE_PARAMETERS
+  int myproc, nprocs;
+  int timelevel;
+  ioGH *ioUtilGH;
+  flexioGH *myGH;
+  IOFile *iofile;
+  char *filename;
+
+
+  /* Get the handle for IOUtil and IOFlexIO extensions */
+  ioUtilGH = (ioGH *) GH->extensions [CCTK_GetGHExtensionHandle ("IO")];
+  myGH = (flexioGH *) GH->extensions [CCTK_GetGHExtensionHandle ("IOFlexIO")];
+
+  /* Get filename and reference to file handle */
+  filename = myGH->IEEEfname_3D [index];
+  iofile = &myGH->IEEEfile_3D [index];
+
+  if (IO_verbose) {
+    printf ("-------------------------------------------------------\n");
+    printf ("3D I/O on Grid Function %s\n", alias);
+    printf ("  downsample (x,y,z): (%d,%d,%d)  ioproc : %d\n", 
+            ioUtilGH->downsample_x, ioUtilGH->downsample_y,
+            ioUtilGH->downsample_z, ioUtilGH->ioproc);
+    fflush (stdout);
+  }
+
+  /* What processor are we on? */
+  myproc = CCTK_GetMyProc (GH);
+  nprocs = CCTK_GetnProcs (GH);
+
+
+  /* build the filename for output */
+  /* This code is semi-repeated in DumpGH.c */
+  /* The output directory is also created. */
+  if (first || onefileperslice)
+    IOFlexIO_Write3D_filename (GH, filename, alias);
+
+  /* Close the file if it is already open (if IEEEfile_3D is set). 
+   *  This shouldn't be needed, since we do it below, 
+   *  but better safe than sorry.
+   *  Also I have some vague recollection of needing it in some
+   *  very strange case ... So leave it here for now (PW 11.5.98)
+   */
+  if (*iofile && onefileperslice) {
+    if (IO_verbose)
+      printf ("Closing IEEEfile from previous iteration.\n");
+    CACTUS_IEEEIO_ERROR (IOclose (*iofile));
+    *iofile = NULL; /* prevent attempts to close it twice! */
+  }
+
+
+  /* open the output file (only if we are an output processor) */
+  if (myproc == ioUtilGH->ioproc)
+    IOFlexIO_Write3D_openFile (GH, filename, iofile, first);
+  
+  /* get the current timelevel */
+  timelevel = CCTK_GetNumTimeLevels (index) - 1;
+  if (timelevel > 0)
+    timelevel--;
+
+  /* output the data */
+  IOFlexIO_DumpVar (GH, index, timelevel, *iofile);
+
+  /* output the scalars and parameters */
+  /* Don't do this if the output is being byte compared, as the
+   * thorn lists might be different.
+   */
+  if (myproc == ioUtilGH->ioproc) {
+    /* output parameters necessary for filereader datafiles */
+    if (first)
+      IOFlexIO_IEEEIOStructDump (GH, *iofile);
+#if        0
+    if (ioUtilGH->parameters)
+      IOFlexIO_IEEEIOparamDump (*iofile);
+#endif
+    /* close the file */
+    IOFlexIO_Write3D_closeFile (GH, filename, iofile);
+  }
+
+  if (IO_verbose)
+    printf ("-------------------\n");
+}
+
+
+/*@@
+  @routine IOFlexIO_Write3D_closeFile
+  @author  Paul Walker
+  @date    Feb 1997
+  @desc
+  Closes or pauses the IEEEIO output file.
+  @enddesc
+  @history 
+  @hauthor Gabrielle Allen
+  @hdate Sep 1998
+  @hdesc Split into subroutine
+  @endhistory 
+  @@*/
+
+
+void IOFlexIO_Write3D_closeFile (cGH *GH, const char *filename, IOFile *iofile)
+{
+  DECLARE_PARAMETERS
+  ioGH *ioUtilGH;
+  flexioGH *myGH;
+
+  ioUtilGH = (ioGH *) GH->extensions [CCTK_GetGHExtensionHandle ("IO")];
+  myGH = (flexioGH *) GH->extensions [CCTK_GetGHExtensionHandle ("IOFlexIO")];
+
+  if (*iofile == NULL)
+    return;
+
+  if (onefileperslice) {    
+    if (IO_verbose)
+      printf ("Closing IEEEfile from this iteration\n");
+    CACTUS_IEEEIO_ERROR (IOclose (*iofile));
+    *iofile = NULL;
+    return;
+  }
+  
+  if (myGH->reuse_fh) {
+    if (IO_verbose)
+      printf ("Pausing file %s\n", filename);
+    IEEEbufferOff (*iofile); 
+    CACTUS_IEEEIO_ERROR (IOpause (*iofile));
+  } else {
+    if (IO_verbose)
+      printf ("Closing file %s\n", filename);
+    CACTUS_IEEEIO_ERROR (IOclose (*iofile));
+    *iofile = NULL;
+  }
+}
+
+
+/*@@
+  @routine IOFlexIO_Write3D_openFile
+  @author  Paul Walker
+  @date    Feb 1997
+  @desc
+  Opens or resumes the IEEEIO output file
+  @enddesc
+  @history 
+  @hauthor Gabrielle Allen
+  @hdate Sep 1998
+  @hdesc Split into subroutine
+  @endhistory 
+  @@*/
+
+void IOFlexIO_Write3D_openFile (cGH *GH, const char *filename, IOFile *iofile, int first)
+{
+  DECLARE_PARAMETERS
+  ioGH *ioUtilGH;
+  flexioGH *myGH;
+
+  ioUtilGH = (ioGH *) GH->extensions [CCTK_GetGHExtensionHandle ("IO")];
+  myGH = (flexioGH *) GH->extensions [CCTK_GetGHExtensionHandle ("IOFlexIO")];
+
+  /* If we are the first time through, or writing one file per slice, we need
+     to open a file in write mode. Or, if we are reusing file handles we will
+     need to resume the file */
+  
+  if (first || onefileperslice) {
+    
+    /* First time through or one file per slice; open anew */
+    if (IO_verbose)
+      printf ("Opening file %s\n", filename);
+    
+    *iofile = IEEEopen (filename, "w");
+    if (! IOisValid (*iofile)) {
+      char *msg = (char *) malloc (strlen (filename) + 80);
+
+      sprintf (msg, "The file %s could not be opened for writing.", filename);
+      CCTK_WARN (1, msg);
+      free (msg);
+    }
+  } 
+  
+  else if (myGH->reuse_fh) {
+    
+    /* resume the file (reopen descriptor)
+       This allows descriptor reuse without requiring expensive destruction
+       and reallocation of the associated datastructures */
+    
+    if (*iofile) {
+      
+      if (IO_verbose)
+        printf ("Resuming file %s\n", filename);
+
+      CACTUS_IEEEIO_ERROR (IOresume (*iofile));
+    }
+  } else {
+    /* File is closed, not onefileperslice and not first time.
+       Therefore we must be conserving file handles and failed
+       to use the pause()/resume methods(), So append.
+    */
+    if (IO_verbose)
+      printf ("Re-opening file %s in append mode\n", filename);
+    
+    *iofile = IEEEopen (filename, "a");
+    if (! IOisValid (*iofile)) {
+      char *msg = (char *) malloc (strlen (filename) + 80);
+
+      sprintf (msg, "The file %s could not be opened for appending.\n",
+               filename);
+      CCTK_WARN (1, msg);
+      free (msg);
+    }
+  }
+  
+  /* Turn on buffer cache (with default size chosen by IEEEIO lib) */
+  IEEEbufferOn (*iofile, 0);
+}
+
+
+/*@@
+  @routine IOFlexIO_Write3D_filename
+  @author  Paul Walker
+  @date    Feb 1997
+  @desc
+  Generates the filename for output
+  @enddesc
+  @history 
+  @hauthor Gabrielle Allen
+  @hdate Sep 1998
+  @hdesc Split into subroutine
+  @endhistory 
+  @@*/
+
+void IOFlexIO_Write3D_filename (cGH *GH, char *filename, const char *name)
+{
+  DECLARE_PARAMETERS
+  ioGH *ioUtilGH;       /* handle for IOUtil extensions */
+  pGH *pughGH;          /* handle for PUGH extensions */
+  int nprocs;
+  int myproc;
+  char extra [256];        /* Extra stuff in fname based on mode */
+  char extradir [256];        /* Extra stuff for an output dir */
+  char createdir [256];        /* Text for system call to create output directory */
+
+  ioUtilGH = (ioGH *) GH->extensions [CCTK_GetGHExtensionHandle ("IO")];
+  pughGH = (pGH *) GH->extensions [CCTK_GetGHExtensionHandle ("PUGH")];
+
+  extra [0] = '\0';
+  extradir [0] = '\0';
+  
+  nprocs = CCTK_GetnProcs (GH);
+  myproc = CCTK_GetMyProc (GH);
+
+#if        0
+  /* Say if we have a chunked data in the output file(s) */
+  if (! ioUtilGH->unchunked && nprocs > 1)
+    sprintf (extra, "%s_chunked", extra);
+#endif
+  
+  if (onefileperslice)
+    sprintf (extra, "%s.time_%7.3f", extra, GH->time);
+  
+ 
+  /* OUTPUT ONE FILE FOR EACH N PROCESSORS 
+   *  ------------------------------------- 
+   *  
+   *  If only one output file, the single file for each GF is written 
+   *  in the normal output dir (that is extradir = ""). Otherwise
+   *  a directory is created for each output GF to hold the multiple
+   *  file 
+   */
+
+  if (ioUtilGH->ioproc_every < nprocs) {
+    
+    /* Add the output processor number to the extra string */
+    sprintf (extra, "%s.file_%d", extra, myproc / ioUtilGH->ioproc_every);
+    
+    /* If necessary create the output directory */
+    if (myproc == 0) {
+      
+      sprintf (createdir, "mkdir -p %s/%s_3d", ioUtilGH->outpfx_3D, name);
+      
+      if (IO_verbose) {
+        printf ("Creating output directory with command\n %s\n", createdir);
+        fflush (stdout);
+      }
+      
+      if (system (createdir) != 0)
+        printf ("Could not create directory\n");
+    }
+    
+#ifdef MPI
+    /* Wait for all processors to catch up */
+    CACTUS_MPI_ERROR (MPI_Barrier (pughGH->PUGH_COMM_WORLD));
+#endif
+    
+    /* extradir is the relative output directory */
+    sprintf (extradir, "%s_3d/", name);
+    
+  }
+  
+  
+  /* CREATE THE COMPLETE OUTPUT FILENAME 
+     ----------------------------------- */
+  
+  sprintf (filename, "%s/%s%s_3d%s.ieee", ioUtilGH->outpfx_3D,
+           extradir, name, extra);
+
+  /* And be sure to replace any spaces in the filename with an _ */
+  do
+    if (*filename == ' ')
+      *filename = '_';
+  while (*++filename);
+
+}
+
+
+#if        0
+void FMODIFIER FORTRAN_NAME(io_write3d_,
+                            IO_WRITE3D,
+                            io_write3d)
+     (Double *x)
+{
+  pGF *GF;
+
+  GF = locateGFbyDataPtr(x);
+  IO_Write3D(GF->parentGH,GF);
+
+}
+
+#endif
diff --git a/src/ioFlexGH.h b/src/ioFlexGH.h
new file mode 100644
index 0000000..573acbf
--- /dev/null
+++ b/src/ioFlexGH.h
@@ -0,0 +1,85 @@
+ /*@@
+   @header    IOFlexIOGH.h
+   @date      Tue 9th Jan 1999
+   @author    Gabrielle Allen
+   @desc 
+   The extensions to the GH structure from IO.
+   @history
+   @hauthor Thomas Radke @hdate 16 Mar 1999
+   @hdesc Added parameters for 2D and 3D output
+   @hauthor Thomas Radke @hdate 17 Mar 1999
+   @hdesc Changed naming: IEEEIO -> IOFlexIO
+   @hauthor Thomas Radke @hdate 30 Mar 1999
+   @hdesc Undefined DI macro
+   @endhistory
+   @version $Header$
+ @@*/
+
+#include <string.h>
+
+#include "IEEEIO/IO.h"
+#include "IEEEIO/IEEEIO.h"
+
+
+/* define the IOFlexIO datatypes according to CCTK_??? datatypes */
+#define FLEXIO_REAL4   FLOAT32
+
+#ifdef  CCTK_REAL_PRECISION_16
+#error  "128-bit precision floating point types are not supported by IOFlexIO !"
+#elif   CCTK_REAL_PRECISION_8
+#define FLEXIO_REAL    FLOAT64
+#elif   CCTK_REAL_PRECISION_4
+#define FLEXIO_REAL    FLOAT32
+#endif
+
+#ifdef  CCTK_INTEGER_PRECISION_8
+#define FLEXIO_INT     INT64
+#elif   CCTK_INTEGER_PRECISION_4
+#define FLEXIO_INT     INT32
+#elif   CCTK_INTEGER_PRECISION_2
+#define FLEXIO_INT     INT16
+#endif
+
+#define FLEXIO_INT4    INT32
+#define FLEXIO_CHAR    CHAR
+
+
+/* Check error flags from IOFlexIO */
+#define CACTUS_IEEEIO_ERROR(fn_call)                                          \
+    do {                                                                      \
+                                                                              \
+      int error_code = fn_call;                                               \
+                                                                              \
+      if (error_code < 0) {                                                   \
+        char *msg = (char*) malloc(strlen(__FILE__) + strlen(#fn_call) + 160);\
+                                                                              \
+        sprintf (msg,                                                         \
+                 "IEEEIO call %s returned error code %d, file %s, line %d\n", \
+                 #fn_call, error_code, __FILE__, __LINE__);                   \
+        CCTK_WARN (1, msg);                                                   \
+        free (msg);                                                           \
+      }                                                                       \
+    } while (0)
+
+
+typedef struct IOFlexIOGH {
+
+  /* The number of times output */
+  int   *IO_2Dnum;
+  int   *IO_3Dnum;
+
+  /* The last iteration output */
+  int   *IO_2Dlast;
+  int   *IO_3Dlast;
+
+  /* IEEEIO file pointers for 2D output (3 for each variable) */
+  IOFile **IEEEfile_2D;
+
+  /* IEEEIO file name and pointer for 3D output */
+  char  **IEEEfname_3D;
+  IOFile *IEEEfile_3D;
+
+  /* reuse filehandles during IEEEIO */
+  int reuse_fh;
+
+} flexioGH;
diff --git a/src/make.code.defn b/src/make.code.defn
new file mode 100644
index 0000000..2bbe405
--- /dev/null
+++ b/src/make.code.defn
@@ -0,0 +1,3 @@
+SRCS = 	Startup.c GHExtension.c Output2D.c Write2D.c Output3D.c Write3D.c DumpVar.c DumpGH.c RecoverGH.c RestoreFile.c
+
+SUBDIRS = IEEEIO