Brought in latest differences from the SGI version. This is work

Brought in latest differences from the SGI version.  This is work
towards a code that compiles on both architectures.

darcs-hash:20010322174200-f6438-23ab5f26cf84d2666312791c6bdb5a0fc1d0390a.gz

42 files changed, 1064 insertions, 610 deletions

diff --git a/Carpet/Carpet/options/carpet-harpo-gcc b/Carpet/Carpet/options/carpet-harpo-gcc
index ba494a449..bdb26744a 100644
--- a/Carpet/Carpet/options/carpet-harpo-gcc
+++ b/Carpet/Carpet/options/carpet-harpo-gcc
@@ -1,4 +1,4 @@
-# $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/options/carpet-harpo-gcc,v 1.4 2002/09/08 20:38:35 schnetter Exp $
+# $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/options/carpet-harpo-gcc,v 1.2 2001/03/22 18:42:00 eschnett Exp $
 
 CC			= gcc
 CXX			= g++
@@ -13,14 +13,14 @@ LD			= g++
 #LD			= /u1/eschnett/bin/g++
 
 CFLAGS			=
-CXXFLAGS		= -ftemplate-depth-30 -DMPIPP_H
+CXXFLAGS		= -DTMPL_EXPLICIT -ftemplate-depth-30 -DMPIPP_H
 FFLAGS			=
 F90FLAGS		=
 
 C_OPTIMISE_FLAGS	= -O3 -funroll-loops
 CXX_OPTIMISE_FLAGS	= -O3 -funroll-loops
-F90_OPTIMISE_FLAGS	= -O3 -funroll-loops -malign-double
-F77_OPTIMISE_FLAGS	= -O3 -funroll-loops -malign-double
+F90_OPTIMISE_FLAGS	= -O3 -funroll-loops
+F77_OPTIMISE_FLAGS	= -O3 -funroll-loops
 
 C_WARN_FLAGS		= -Wall -Wpointer-arith -Wbad-function-cast -Wcast-qual -Wcast-align -Wstrict-prototypes -Wmissing-declarations -Winline
 CXX_WARN_FLAGS		= -Wall -Wpointer-arith -Wbad-function-cast -Wcast-qual -Wcast-align -Wstrict-prototypes -Wmissing-declarations -Winline
diff --git a/Carpet/Carpet/options/carpet-harpo-sgi b/Carpet/Carpet/options/carpet-harpo-sgi
index 64fcb11b8..a7e740646 100644
--- a/Carpet/Carpet/options/carpet-harpo-sgi
+++ b/Carpet/Carpet/options/carpet-harpo-sgi
@@ -1,30 +1,26 @@
-# $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/options/carpet-harpo-sgi,v 1.5 2002/06/07 20:44:03 schnetter Exp $
+# $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/options/carpet-harpo-sgi,v 1.2 2001/03/22 18:42:00 eschnett Exp $
 
-CPP			= /lib/cpp -DMPIPP_H
+CPP			= /lib/cpp -DTMPL_EXPLICIT -DMPIPP_H
 
-CFLAGS			= -mips4 -r12000 -64 -LANG:restrict=ON
-CXXFLAGS		= -mips4 -r12000 -64 -LANG:restrict=ON -LANG:std -LANG:vla=ON -no_auto_include -ptused -DMPIPP_H
-F90FLAGS		= -mips4 -r12000 -64
-F77FLAGS		= -mips4 -r12000 -64
+CXXFLAGS		= -LANG:std -LANG:vla=ON -64 -mips4 -r12000 -no_auto_include -ptused -DTMPL_EXPLICIT -DMPIPP_H
 LDFLAGS			= -LANG:std -LANG:vla=ON -64 -mips4 -r12000 -Wl,"-woff 84","-woff 85"
 
-C_OPTIMISE_FLAGS	= -O3 -OPT:Olimit=0
-CXX_OPTIMISE_FLAGS	= -O3 -OPT:Olimit=0
-F77_OPTIMISE_FLAGS	= -O3 -OPT:Olimit=0
-F90_OPTIMISE_FLAGS	= -O3 -OPT:Olimit=0
+C_OPTIMISE_FLAGS	= -O3 -INLINE -LNO -OPT:Olimit=100000
+CXX_OPTIMISE_FLAGS	= -O3 -INLINE -LNO -OPT:Olimit=100000
+F77_OPTIMISE_FLAGS	= -O3 -INLINE -LNO -OPT:Olimit=100000
+F90_OPTIMISE_FLAGS	= -O3 -INLINE -LNO -OPT:Olimit=100000
 
 C_DEBUG_FLAGS		= -g3
 CXX_DEBUG_FLAGS		= -g3
 F77_DEBUG_FLAGS		= -g3
 F90_DEBUG_FLAGS		= -g3
 
-SYS_INC_DIRS		= /usr/local/hdf/include /u1/eschnett/FlexIO
-LIBDIRS			= /usr/local/hdf/lib /u1/eschnett/FlexIO
-LIBS			= AMR h5io hdfio hlio ieeeio   ffio fpe fortran ftn ftn90   hdf5 mfhdf df jpeg z
+SYS_INC_DIRS		= /usr/local/hdf/include /u1/eschnett/scratch/FlexIO
+LIBDIRS			= /usr/local/hdf/lib /u1/eschnett/scratch/FlexIO
+LIBS			= AMR h5io hdfio hlio ieeeio  hdf5 mfhdf df jpeg z
 
 HDF5			= yes
 HDF5_DIR		= /usr/local/hdf5_64
-
 MPI			= NATIVE
 
 OPTIMISE		= yes
diff --git a/Carpet/Carpet/options/carpet-lilypond b/Carpet/Carpet/options/carpet-lilypond
index 0b20c83bf..e01ba7751 100644
--- a/Carpet/Carpet/options/carpet-lilypond
+++ b/Carpet/Carpet/options/carpet-lilypond
@@ -1,9 +1,9 @@
-# $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/options/carpet-lilypond,v 1.3 2001/03/19 21:29:45 eschnett Exp $
+# $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/options/carpet-lilypond,v 1.4 2001/03/22 18:42:00 eschnett Exp $
 
+CPP			= gcc -E -DTMPL_EXPLICIT -DCARPET_REAL
 F77			= g77
 F90			= g77
 
-CPPFLAGS		= -DTMPL_EXPLICIT -DCARPET_REAL
 CXXFLAGS		= -DTMPL_EXPLICIT -DCARPET_REAL -ftemplate-depth-30
 
 C_OPTIMISE_FLAGS	= -O3 -funroll-loops -march=pentiumpro
diff --git a/Carpet/Carpet/src/carpet.cc b/Carpet/Carpet/src/carpet.cc
index a317a7dfb..010e783d2 100644
--- a/Carpet/Carpet/src/carpet.cc
+++ b/Carpet/Carpet/src/carpet.cc
@@ -1,5 +1,3 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Attic/carpet.cc,v 1.16 2001/03/19 21:29:59 eschnett Exp $
-
 // It is assumed that the number of components of all arrays is equal
 // to the number of components of the grid functions, and that their
 // distribution onto the processors is the same, and that all
@@ -8,14 +6,15 @@
 // Scalar variables currently exist in one single incarnation for all
 // refinement levels and all components.
 
+#include <assert.h>
+#include <math.h>
+#include <stdarg.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+
 #include <algorithm>
-#include <cassert>
-#include <cmath>
 #include <complex>
-#include <cstdarg>
-#include <cstdio>
-#include <cstdlib>
-#include <cstring>
 
 #include <mpi.h>
 
@@ -36,12 +35,14 @@
 
 #include "carpet.hh"
 
-static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Attic/carpet.cc,v 1.16 2001/03/19 21:29:59 eschnett Exp $";
+static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/Attic/carpet.cc,v 1.17 2001/03/22 18:42:05 eschnett Exp $";
 
 
 
 namespace Carpet {
   
+  using namespace std;
+  
   static void Recompose (cGH* cgh);
   static void CycleTimeLevels (cGH* cgh);
   static void Restrict (cGH* cgh);
@@ -138,7 +139,8 @@ namespace Carpet {
     
     // Refinement information
     maxreflevels = max_refinement_levels;
-    maxreflevelfact = (int)floor(pow(refinement_factor, maxreflevels-1) + 0.5);
+    maxreflevelfact
+      = (int)floor(pow((double)refinement_factor, maxreflevels-1) + 0.5);
     
     // Ghost zones
     vect<int,dim> lghosts, ughosts;
@@ -212,9 +214,7 @@ namespace Carpet {
 	   multigrid_factor, vertex_centered,
 	   arrext);
 	
-	arrdata[group].tt = new th<dim>
-	  (*hh,
-	   (int)floor(pow(refinement_factor, max_refinement_levels-1) + 0.5));
+	arrdata[group].tt = new th<dim> (*hh, maxreflevelfact);
 	
 	vect<int,dim> alghosts, aughosts;
 	for (int d=0; d<dim; ++d) {
@@ -642,7 +642,11 @@ namespace Carpet {
 #include "typecase"
 #undef TYPECASE
 	  default:
-	    abort();
+	    CCTK_VWarn
+	      (0, __LINE__, __FILE__, CCTK_THORNSTRING,
+	       "Carpet does not support the type of the variable \"%s\".\n"
+	       "Either enable support for this type, "
+	       "or change the type of this variable.", CCTK_FullName(n));
 	  }
 	}
       }
@@ -666,7 +670,11 @@ namespace Carpet {
 #include "typecase"
 #undef TYPECASE
 	default:
-	  abort();
+	  CCTK_VWarn
+	    (0, __LINE__, __FILE__, CCTK_THORNSTRING,
+	     "Carpet does not support the type of the variable \"%s\".\n"
+	     "Either enable support for this type, "
+	     "or change the type of this variable.", CCTK_FullName(n));
 	}
       }
       break;
@@ -688,7 +696,11 @@ namespace Carpet {
 #include "typecase"
 #undef TYPECASE
 	default:
-	  abort();
+	  CCTK_VWarn
+	    (0, __LINE__, __FILE__, CCTK_THORNSTRING,
+	     "Carpet does not support the type of the variable \"%s\".\n"
+	     "Either enable support for this type, "
+	     "or change the type of this variable.", CCTK_FullName(n));
 	}
       }
       break;
@@ -732,7 +744,11 @@ namespace Carpet {
 #include "typecase"
 #undef TYPECASE
 	  default:
-	    abort();
+	    CCTK_VWarn
+	      (0, __LINE__, __FILE__, CCTK_THORNSTRING,
+	       "Carpet does not support the type of the variable \"%s\".\n"
+	       "Either enable support for this type, "
+	       "or change the type of this variable.", CCTK_FullName(n));
 	  }
 	  scdata[group][var][ti] = 0;
 	}
@@ -754,7 +770,11 @@ namespace Carpet {
 #include "typecase"
 #undef TYPECASE
 	default:
-	  abort();
+	  CCTK_VWarn
+	    (0, __LINE__, __FILE__, CCTK_THORNSTRING,
+	     "Carpet does not support the type of the variable \"%s\".\n"
+	     "Either enable support for this type, "
+	     "or change the type of this variable.", CCTK_FullName(n));
 	}
 	arrdata[group].data[var] = 0;
       }
@@ -776,7 +796,11 @@ namespace Carpet {
 #include "typecase"
 #undef TYPECASE
 	default:
-	  abort();
+	  CCTK_VWarn
+	    (0, __LINE__, __FILE__, CCTK_THORNSTRING,
+	     "Carpet does not support the type of the variable \"%s\".\n"
+	     "Either enable support for this type, "
+	     "or change the type of this variable.", CCTK_FullName(n));
 	}
 	gfdata[group].data[var] = 0;
       }
@@ -889,6 +913,7 @@ namespace Carpet {
 	  case CCTK_SCALAR: {
 	    assert (group<(int)scdata.size());
 	    assert (var<(int)scdata[group].size());
+	    const int num_tl = CCTK_NumTimeLevelsFromVarI(n);
 	    void* tmpdata = scdata[group][var][0];
 	    for (int ti=0; ti<num_tl-1; ++ti) {
 	      // TODO: Which refinement level to use?
@@ -902,11 +927,7 @@ namespace Carpet {
 	    assert (group<(int)arrdata.size());
 	    assert (var<(int)arrdata[group].data.size());
 	    for (int c=0; c<arrdata[group].hh->components(reflevel); ++c) {
-	      for (int ti=0; ti<num_tl-1; ++ti) {
-		const int tmin = min(0, 2 - num_tl);
-		const int tl = tmin + ti;
- 		arrdata[group].data[var]->copy(tl, reflevel, c, mglevel);
-	      }
+	      arrdata[group].data[var]->cycle (reflevel, c, mglevel);
 	    }
 	    break;
 	  }
@@ -914,11 +935,7 @@ namespace Carpet {
 	    assert (group<(int)gfdata.size());
 	    assert (var<(int)gfdata[group].data.size());
 	    for (int c=0; c<hh->components(reflevel); ++c) {
-	      for (int ti=0; ti<CCTK_NumTimeLevelsFromVarI(n)-1; ++ti) {
-		const int tmin = min(0, 2 - CCTK_NumTimeLevelsFromVarI(n));
-		const int tl   = tmin + ti;
-		gfdata[group].data[var]->copy(tl, reflevel, c, mglevel);
-	      }
+	      gfdata[group].data[var]->cycle (reflevel, c, mglevel);
 	    }
 	    break;
 	  }
@@ -1143,7 +1160,7 @@ namespace Carpet {
     // Change
     reflevel = rl;
     const bbox<int,dim>& base = hh->baseextent;
-    reflevelfact = (int)floor(pow(hh->reffact, reflevel)+0.5);
+    reflevelfact = (int)floor(pow((double)hh->reffact, reflevel)+0.5);
     cgh->cctk_delta_time = base_delta_time / reflevelfact;
     for (int d=0; d<dim; ++d) {
       cgh->cctk_gsh[d]
diff --git a/Carpet/Carpet/src/typecase b/Carpet/Carpet/src/typecase
index f964d3300..dbe0d8823 100644
--- a/Carpet/Carpet/src/typecase
+++ b/Carpet/Carpet/src/typecase
@@ -1,7 +1,7 @@
 // Instantiate type cases for all available types    -*-C++-*-
 // (C) 2001 Erik Schnetter <schnetter@uni-tuebingen.de>
 
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/typecase,v 1.2 2001/03/19 21:30:00 eschnett Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/Carpet/src/typecase,v 1.3 2001/03/22 18:42:05 eschnett Exp $
 
 // Usage:
 // Define the macro TYPECASE(N,T) to be a typecase for the type T with name N,
@@ -51,7 +51,9 @@
 #endif
 
 #if !defined(CARPET_BYTE) && !defined(CARPET_INT2) && !defined(CARPET_INT4) && !defined(CARPET_INT8) && !defined(CARPET_REAL4) && !defined(CARPET_REAL8) && !defined(CARPET_REAL16) && !defined(CARPET_COMPLEX8) && !defined(CARPET_COMPLEX16) && !defined(CARPET_COMPLEX32)
-// Assume the user just wants REAL
+// Assume the user just wants INT and REAL
+#  undef CARPET_INT
+#  define CARPET_INT
 #  undef CARPET_REAL
 #  define CARPET_REAL
 #endif
@@ -123,7 +125,7 @@ TYPECASE(CCTK_VARIABLE_INT2, CCTK_INT2)
 #ifdef CARPET_INT4
 TYPECASE(CCTK_VARIABLE_INT4, CCTK_INT4)
 #endif
-#ifdef CARPET_INT4
+#ifdef CARPET_INT8
 #  ifdef CCTK_INT8
 TYPECASE(CCTK_VARIABLE_INT8, CCTK_INT8)
 #  endif
diff --git a/Carpet/CarpetIOASCII/schedule.ccl b/Carpet/CarpetIOASCII/schedule.ccl
index da776c850..aff6826a2 100644
--- a/Carpet/CarpetIOASCII/schedule.ccl
+++ b/Carpet/CarpetIOASCII/schedule.ccl
@@ -1,5 +1,5 @@
 # Schedule definitions for thorn CarpetIOASCII
-# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetIOASCII/schedule.ccl,v 1.2 2001/03/18 05:20:24 eschnett Exp $
+# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetIOASCII/schedule.ccl,v 1.3 2001/03/22 18:42:05 eschnett Exp $
 
 schedule CarpetIOASCIIStartup at STARTUP after IOUtil_Startup
 {
diff --git a/Carpet/CarpetIOASCII/src/ioascii.cc b/Carpet/CarpetIOASCII/src/ioascii.cc
index f76f5e2e5..99aa1491d 100644
--- a/Carpet/CarpetIOASCII/src/ioascii.cc
+++ b/Carpet/CarpetIOASCII/src/ioascii.cc
@@ -1,11 +1,13 @@
-#include <cassert>
-#include <climits>
-#include <cstdio>
-#include <cstdlib>
-#include <cstring>
-#include <fstream>
+#include <alloca.h>
+#include <assert.h>
+#include <limits.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
 #include <sys/stat.h>
 #include <sys/types.h>
+
+#include <fstream>
 #include <vector>
 
 #include "cctk.h"
@@ -22,10 +24,11 @@
 
 #include "ioascii.hh"
 
-static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetIOASCII/src/ioascii.cc,v 1.10 2001/03/18 22:37:04 eschnett Exp $";
+static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetIOASCII/src/ioascii.cc,v 1.11 2001/03/22 18:42:05 eschnett Exp $";
 
 
 
+using namespace std;
 using namespace Carpet;
 
 
@@ -218,13 +221,16 @@ int CarpetIOASCII<outdim>
 	if (desired) {
 	  
 	  // Invent a file name
-	  char filename[strlen(myoutdir)+strlen(alias)+100];
+	  char* const filename
+	    = (char*)alloca(strlen(myoutdir)+strlen(alias)+100);
 	  sprintf (filename, "%s/%s.", myoutdir, alias);
 	  for (int d=0; d<outdim; ++d) {
 	    assert (dirs[d]>=0 && dirs[d]<3);
-	    sprintf (filename, "%s%c", filename, "xyz"[dirs[d]]);
+	    const char* const coords = "xyz";
+	    sprintf (filename, "%s%c", filename, coords[dirs[d]]);
 	  }
-	  sprintf (filename, "%s%c", filename, "lpv"[outdim-1]);
+	  const char* const suffixes = "lpv";
+	  sprintf (filename, "%s%c", filename, suffixes[outdim-1]);
 	  
 	  // If this is the first time, then write a nice header on
 	  // the root processor
@@ -575,7 +581,7 @@ bool CheckForVariable (cGH* const cgh,
   const int numvars = CCTK_NumVars();
   assert (vindex>=0 && vindex<numvars);
   
-  bool flags[numvars];
+  bool* const flags = (bool*)alloca(numvars * sizeof(bool));
   
   for (int i=0; i<numvars; ++i) {
     flags[i] = false;
@@ -594,6 +600,6 @@ void SetFlag (int index, const char* optstring, void* arg)
 
 
 // Explicit instantiation for all output dimensions
-template CarpetIOASCII<1>;
-template CarpetIOASCII<2>;
-template CarpetIOASCII<3>;
+template class CarpetIOASCII<1>;
+template class CarpetIOASCII<2>;
+template class CarpetIOASCII<3>;
diff --git a/Carpet/CarpetIOASCII/src/ioascii.hh b/Carpet/CarpetIOASCII/src/ioascii.hh
index 061e4c18c..a0a311c41 100644
--- a/Carpet/CarpetIOASCII/src/ioascii.hh
+++ b/Carpet/CarpetIOASCII/src/ioascii.hh
@@ -1,9 +1,11 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetIOASCII/src/ioascii.hh,v 1.4 2001/03/16 21:32:17 eschnett Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetIOASCII/src/ioascii.hh,v 1.5 2001/03/22 18:42:05 eschnett Exp $
 
 #include <vector>
 
 #include "cctk.h"
 
+using namespace std;
+
 
 
 // scheduled functions
diff --git a/Carpet/CarpetLib/src/bbox.cc b/Carpet/CarpetLib/src/bbox.cc
index de281f582..f90c14d4e 100644
--- a/Carpet/CarpetLib/src/bbox.cc
+++ b/Carpet/CarpetLib/src/bbox.cc
@@ -5,7 +5,7 @@
     copyright            : (C) 2000 by Erik Schnetter
     email                : schnetter@astro.psu.edu
 
-    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/bbox.cc,v 1.4 2001/03/12 16:54:25 eschnett Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/bbox.cc,v 1.5 2001/03/22 18:42:05 eschnett Exp $
 
  ***************************************************************************/
 
@@ -18,7 +18,8 @@
  *                                                                         *
  ***************************************************************************/
 
-#include <cassert>
+#include <assert.h>
+
 #include <iostream>
 
 #include "defs.hh"
@@ -28,6 +29,8 @@
 #  include "bbox.hh"
 #endif
 
+using namespace std;
+
 
 
 // Constructors
@@ -223,19 +226,37 @@ bbox<T,D>::iterator bbox<T,D>::end () const {
 
 
 // Output
+#ifndef SGI
+// This doesn't work on SGIs.  Is this legal C++?
 template<class T,int D>
 ostream& operator<< (ostream& os, const bbox<T,D>& b) {
   os << "(" << b.lower() << ":" << b.upper() << ":" << b.stride() << ")";
   return os;
 }
+#else
+ostream& operator<< (ostream& os, const bbox<int,1>& b) {
+  os << "(" << b.lower() << ":" << b.upper() << ":" << b.stride() << ")";
+  return os;
+}
+ostream& operator<< (ostream& os, const bbox<int,2>& b) {
+  os << "(" << b.lower() << ":" << b.upper() << ":" << b.stride() << ")";
+  return os;
+}
+ostream& operator<< (ostream& os, const bbox<int,3>& b) {
+  os << "(" << b.lower() << ":" << b.upper() << ":" << b.stride() << ")";
+  return os;
+}
+#endif
 
 
 
 #if defined(TMPL_EXPLICIT)
 template class bbox<int,1>;
 template ostream& operator<< (ostream& os, const bbox<int,1>& b);
+
 template class bbox<int,2>;
 template ostream& operator<< (ostream& os, const bbox<int,2>& b);
+
 template class bbox<int,3>;
 template ostream& operator<< (ostream& os, const bbox<int,3>& b);
 #endif
diff --git a/Carpet/CarpetLib/src/bbox.hh b/Carpet/CarpetLib/src/bbox.hh
index db43519b4..396a9c5ca 100644
--- a/Carpet/CarpetLib/src/bbox.hh
+++ b/Carpet/CarpetLib/src/bbox.hh
@@ -5,7 +5,7 @@
     copyright            : (C) 2000 by Erik Schnetter
     email                : schnetter@astro.psu.edu
 
-    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/bbox.hh,v 1.6 2001/03/14 11:00:26 eschnett Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/bbox.hh,v 1.7 2001/03/22 18:42:05 eschnett Exp $
 
  ***************************************************************************/
 
@@ -26,13 +26,15 @@
 #include "defs.hh"
 #include "vect.hh"
 
+using namespace std;
+
 
 
 // Forward definition
 template<class T, int D> class bbox;
 
 // Output
-template<class T,int D>
+template<class T, int D>
 ostream& operator<< (ostream& os, const bbox<T,D>& b);
 
 
@@ -102,7 +104,7 @@ public:
   // Iterators
   class iterator {
   protected:
-    bbox box;
+    const bbox& box;
     vect<T,D> pos;
   public:
     iterator (const bbox& box, const vect<T,D>& pos);
@@ -115,7 +117,7 @@ public:
   iterator end () const;
   
   // Output
-  friend ostream& operator<< <>(ostream& os, const bbox& b);
+  friend ostream& operator<< <>(ostream& os, const bbox& s);
 };
 
 
diff --git a/Carpet/CarpetLib/src/bboxset.cc b/Carpet/CarpetLib/src/bboxset.cc
index 489856298..f18385c1b 100644
--- a/Carpet/CarpetLib/src/bboxset.cc
+++ b/Carpet/CarpetLib/src/bboxset.cc
@@ -5,7 +5,7 @@
     copyright            : (C) 2000 by Erik Schnetter
     email                : schnetter@astro.psu.edu
 
-    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/bboxset.cc,v 1.4 2001/03/19 21:30:19 eschnett Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/bboxset.cc,v 1.5 2001/03/22 18:42:05 eschnett Exp $
 
  ***************************************************************************/
 
@@ -18,7 +18,8 @@
  *                                                                         *
  ***************************************************************************/
 
-#include <cassert>
+#include <assert.h>
+
 #include <iostream>
 #include <set>
 
@@ -28,6 +29,8 @@
 #  include "bboxset.hh"
 #endif
 
+using namespace std;
+
 
 
 // Constructors
@@ -210,6 +213,8 @@ bboxset<T,D>& bboxset<T,D>::operator&= (const bboxset& s) {
 
 
 // Difference
+#ifndef SGI
+// This doesn't work on SGIs.  Is this legal C++?
 template<class T, int D>
 bboxset<T,D> operator- (const bbox<T,D>& b1, const bbox<T,D>& b2) {
   assert (b1.is_aligned_with(b2));
@@ -243,6 +248,104 @@ bboxset<T,D> operator- (const bbox<T,D>& b1, const bbox<T,D>& b2) {
   assert (r.invariant());
   return r;
 }
+#else
+bboxset<int,1> operator- (const bbox<int,1>& b1, const bbox<int,1>& b2) {
+  assert (b1.is_aligned_with(b2));
+  if (b1.empty()) return bboxset<int,1>();
+  if (b2.empty()) return bboxset<int,1>(b1);
+  const vect<int,1> str = b1.stride();
+  bboxset<int,1> r;
+  for (int d=0; d<1; ++d) {
+    // make resulting bboxes as large as possible in x-direction (for
+    // better consumption by Fortranly ordered arrays)
+    vect<int,1> lb, ub;
+    bbox<int,1> b;
+    for (int dd=0; dd<1; ++dd) {
+      if (dd<d) {
+	lb[dd] = b2.lower()[dd];
+	ub[dd] = b2.upper()[dd];
+      } else if (dd>d) {
+	lb[dd] = b1.lower()[dd];
+	ub[dd] = b1.upper()[dd];
+      }
+    }
+    lb[d] = b1.lower()[d];
+    ub[d] = b2.lower()[d] - str[d];
+    b = bbox<int,1>(lb,ub,str) & b1;
+    r += b;
+    lb[d] = b2.upper()[d] + str[d];
+    ub[d] = b1.upper()[d];
+    b = bbox<int,1>(lb,ub,str) & b1;
+    r += b;
+  }
+  assert (r.invariant());
+  return r;
+}
+bboxset<int,2> operator- (const bbox<int,2>& b1, const bbox<int,2>& b2) {
+  assert (b1.is_aligned_with(b2));
+  if (b1.empty()) return bboxset<int,2>();
+  if (b2.empty()) return bboxset<int,2>(b1);
+  const vect<int,2> str = b1.stride();
+  bboxset<int,2> r;
+  for (int d=0; d<2; ++d) {
+    // make resulting bboxes as large as possible in x-direction (for
+    // better consumption by Fortranly ordered arrays)
+    vect<int,2> lb, ub;
+    bbox<int,2> b;
+    for (int dd=0; dd<2; ++dd) {
+      if (dd<d) {
+	lb[dd] = b2.lower()[dd];
+	ub[dd] = b2.upper()[dd];
+      } else if (dd>d) {
+	lb[dd] = b1.lower()[dd];
+	ub[dd] = b1.upper()[dd];
+      }
+    }
+    lb[d] = b1.lower()[d];
+    ub[d] = b2.lower()[d] - str[d];
+    b = bbox<int,2>(lb,ub,str) & b1;
+    r += b;
+    lb[d] = b2.upper()[d] + str[d];
+    ub[d] = b1.upper()[d];
+    b = bbox<int,2>(lb,ub,str) & b1;
+    r += b;
+  }
+  assert (r.invariant());
+  return r;
+}
+bboxset<int,3> operator- (const bbox<int,3>& b1, const bbox<int,3>& b2) {
+  assert (b1.is_aligned_with(b2));
+  if (b1.empty()) return bboxset<int,3>();
+  if (b2.empty()) return bboxset<int,3>(b1);
+  const vect<int,3> str = b1.stride();
+  bboxset<int,3> r;
+  for (int d=0; d<3; ++d) {
+    // make resulting bboxes as large as possible in x-direction (for
+    // better consumption by Fortranly ordered arrays)
+    vect<int,3> lb, ub;
+    bbox<int,3> b;
+    for (int dd=0; dd<3; ++dd) {
+      if (dd<d) {
+	lb[dd] = b2.lower()[dd];
+	ub[dd] = b2.upper()[dd];
+      } else if (dd>d) {
+	lb[dd] = b1.lower()[dd];
+	ub[dd] = b1.upper()[dd];
+      }
+    }
+    lb[d] = b1.lower()[d];
+    ub[d] = b2.lower()[d] - str[d];
+    b = bbox<int,3>(lb,ub,str) & b1;
+    r += b;
+    lb[d] = b2.upper()[d] + str[d];
+    ub[d] = b1.upper()[d];
+    b = bbox<int,3>(lb,ub,str) & b1;
+    r += b;
+  }
+  assert (r.invariant());
+  return r;
+}
+#endif
 
 template<class T, int D>
 bboxset<T,D> bboxset<T,D>::operator- (const box& b) const {
@@ -281,21 +384,38 @@ bboxset<T,D> bboxset<T,D>::operator- (const bboxset& s) const {
   return r;
 }
 
+#ifndef SGI
+// This doesn't work on SGIs.  Is this legal C++?
 template<class T, int D>
 bboxset<T,D> operator- (const bbox<T,D>& b, const bboxset<T,D>& s) {
   bboxset<T,D> r = bboxset<T,D>(b) - s;
   assert (r.invariant());
   return r;
 }
+#else
+bboxset<int,1> operator- (const bbox<int,1>& b, const bboxset<int,1>& s) {
+  bboxset<int,1> r = bboxset<int,1>(b) - s;
+  assert (r.invariant());
+  return r;
+}
+bboxset<int,2> operator- (const bbox<int,2>& b, const bboxset<int,2>& s) {
+  bboxset<int,2> r = bboxset<int,2>(b) - s;
+  assert (r.invariant());
+  return r;
+}
+bboxset<int,3> operator- (const bbox<int,3>& b, const bboxset<int,3>& s) {
+  bboxset<int,3> r = bboxset<int,3>(b) - s;
+  assert (r.invariant());
+  return r;
+}
+#endif
 
 
 
 // Output
 template<class T,int D>
 ostream& operator<< (ostream& os, const bboxset<T,D>& s) {
-//   os << "bboxset<" STR(T) "," << D << ">:size=" << s.size() << ","
-//      << "set=" << s.bs;
-  os << s.bs;
+  os << "bboxset<T," << D << ">:size=" << s.size() << "," << "set=" << s.bs;
   return os;
 }
 
@@ -303,7 +423,11 @@ ostream& operator<< (ostream& os, const bboxset<T,D>& s) {
 
 #if defined(TMPL_EXPLICIT)
 template class bboxset<int,3>;
-template bboxset<int,3>	operator- (const bbox<int,3>& b1, const bbox<int,3>& b3);
-template bboxset<int,3>	operator- (const bbox<int,3>& b, const bboxset<int,3>& s);
-template ostream& operator<< (ostream& os, const bboxset<int,3>& b);
+
+template
+bboxset<int,3> operator- (const bbox<int,3>& b1, const bbox<int,3>& b3);
+template
+bboxset<int,3> operator- (const bbox<int,3>& b, const bboxset<int,3>& s);
+template
+ostream& operator<< (ostream& os, const bboxset<int,3>& b);
 #endif
diff --git a/Carpet/CarpetLib/src/bboxset.hh b/Carpet/CarpetLib/src/bboxset.hh
index dae3d3f3f..4883cc6e9 100644
--- a/Carpet/CarpetLib/src/bboxset.hh
+++ b/Carpet/CarpetLib/src/bboxset.hh
@@ -5,7 +5,7 @@
     copyright            : (C) 2000 by Erik Schnetter
     email                : schnetter@astro.psu.edu
 
-    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/bboxset.hh,v 1.3 2001/03/12 16:54:25 eschnett Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/bboxset.hh,v 1.4 2001/03/22 18:42:05 eschnett Exp $
 
  ***************************************************************************/
 
@@ -21,7 +21,8 @@
 #ifndef BBOXSET_HH
 #define BBOXSET_HH
 
-#include <cassert>
+#include <assert.h>
+
 #include <iostream>
 #include <set>
 
@@ -29,6 +30,8 @@
 #include "defs.hh"
 #include "vect.hh"
 
+using namespace std;
+
 
 
 // Forward definition
@@ -105,8 +108,8 @@ public:
   // friend bboxset operator- <T,D>(const box& b, const bboxset& s);
   
   // Iterators
-  typedef bset::const_iterator const_iterator;
-  typedef bset::iterator       iterator;
+  typedef typename bset::const_iterator const_iterator;
+  typedef typename bset::iterator       iterator;
   
   iterator begin () const { return bs.begin(); }
   iterator end () const   { return bs.end(); }
diff --git a/Carpet/CarpetLib/src/copy_3d_real8.F77 b/Carpet/CarpetLib/src/copy_3d_real8.F77
index 716ffd135..ae5d1c34f 100644
--- a/Carpet/CarpetLib/src/copy_3d_real8.F77
+++ b/Carpet/CarpetLib/src/copy_3d_real8.F77
@@ -1,10 +1,7 @@
 c     -*-Fortran-*-
-c     $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/copy_3d_real8.F77,v 1.8 2004/03/11 12:03:09 schnetter Exp $
+c     $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/copy_3d_real8.F77,v 1.2 2001/03/22 18:42:05 eschnett Exp $
       
 #include "cctk.h"
-#include "cctk_Parameters.h"
-
-
       
       subroutine copy_3d_real8 (
      $     src, srciext, srcjext, srckext,
@@ -13,8 +10,6 @@ c     $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/copy_3d_real8
       
       implicit none
       
-      DECLARE_CCTK_PARAMETERS
-      
       integer srciext, srcjext, srckext
       CCTK_REAL8 src(srciext,srcjext,srckext)
       integer dstiext, dstjext, dstkext
@@ -32,8 +27,6 @@ c     bbox(:,3) is stride
       integer i, j, k
       integer d
       
-      character msg*1000
-      
       
       
       do d=1,3
@@ -46,25 +39,10 @@ c     bbox(:,3) is stride
             call CCTK_WARN (0, "Internal error: strides disagree")
          end if
          if (mod(srcbbox(d,1), srcbbox(d,3)).ne.0
-     $        .or. mod(dstbbox(d,1), dstbbox(d,3)).ne.0
+     $     .or. mod(dstbbox(d,1), dstbbox(d,3)).ne.0
      $        .or. mod(regbbox(d,1), regbbox(d,3)).ne.0) then
             call CCTK_WARN (0, "Internal error: array origins are not integer multiples of the strides")
          end if
-         if (regbbox(d,1).gt.regbbox(d,2)) then
-c     This could be handled, but is likely to point to an error elsewhere
-            call CCTK_WARN (0, "Internal error: region extent is empty")
-         end if
-         if (mod(srcbbox(d,2) - srcbbox(d,1), srcbbox(d,3)).ne.0
-     $        .or. mod(dstbbox(d,2) - dstbbox(d,1), dstbbox(d,3)).ne.0
-     $        .or. mod(regbbox(d,2) - regbbox(d,1), regbbox(d,3)).ne.0) then
-            call CCTK_WARN (0, "Internal error: array extents are not integer multiples of the strides")
-         end if
-         if (regbbox(d,1).lt.srcbbox(d,1)
-     $        .or. regbbox(d,1).lt.dstbbox(d,1)
-     $        .or. regbbox(d,2).gt.srcbbox(d,2)
-     $        .or. regbbox(d,2).gt.dstbbox(d,2)) then
-            call CCTK_WARN (0, "Internal error: region extent is not contained in array extent")
-         end if
       end do
       
       if (srciext.ne.(srcbbox(1,2)-srcbbox(1,1))/srcbbox(1,3)+1
@@ -93,19 +71,12 @@ c     This could be handled, but is likely to point to an error elsewhere
       
       
 c     Loop over region
-      do k = 1, regkext
-         do j = 1, regjext
-            do i = 1, regiext
-               
-               if (check_array_accesses.ne.0) then
-                  call checkindex (srcioff+i, srcjoff+j+1, srckoff+k+1, 1,1,1,
-     $                 "source")
-                  call checkindex (dstioff+i, dstjoff+j+1, dstkoff+k+1, 1,1,1,
-     $                 "destination")
-               end if
+      do k = 0, regkext-1
+         do j = 0, regjext-1
+            do i = 0, regiext-1
                
-               dst (dstioff+i, dstjoff+j, dstkoff+k)
-     $              = src (srcioff+i, srcjoff+j, srckoff+k)
+               dst (dstioff+i+1, dstjoff+j+1, dstkoff+k+1)
+     $              = src (srcioff+i+1, srcjoff+j+1, srckoff+k+1)
                
             end do
          end do
diff --git a/Carpet/CarpetLib/src/data.cc b/Carpet/CarpetLib/src/data.cc
index 932b692e6..23cd28b57 100644
--- a/Carpet/CarpetLib/src/data.cc
+++ b/Carpet/CarpetLib/src/data.cc
@@ -5,7 +5,7 @@
     copyright            : (C) 2000 by Erik Schnetter
     email                : schnetter@astro.psu.edu
 
-    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/data.cc,v 1.7 2001/03/19 21:30:19 eschnett Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/data.cc,v 1.8 2001/03/22 18:42:05 eschnett Exp $
 
  ***************************************************************************/
 
@@ -18,12 +18,15 @@
  *                                                                         *
  ***************************************************************************/
 
-#include <cassert>
+#include <assert.h>
+
 #include <fstream>
 #include <string>
 
 #include <mpi.h>
 
+#include "cctk.h"
+
 #include "bbox.hh"
 #include "defs.hh"
 #include "dist.hh"
@@ -33,6 +36,8 @@
 #  include "data.hh"
 #endif
 
+using namespace std;
+
 
 
 // Constructors
@@ -65,7 +70,7 @@ data<T,D>* data<T,D>::make_typed () const {
 // Storage management
 template<class T, int D>
 void data<T,D>::allocate (const ibbox& extent, const int proc,
-			  void* const mem=0) {
+			  void* const mem) {
   assert (!_has_storage);
   _has_storage = true;
   // data
@@ -110,7 +115,7 @@ void data<T,D>::transfer_from (generic_data<D>* gsrc) {
 
 // Processor management
 template<class T, int D>
-void data<T,D>::change_processor (const int newproc, void* const mem=0) {
+void data<T,D>::change_processor (const int newproc, void* const mem) {
   if (newproc == _proc) {
     assert (!mem);
     return;
@@ -162,7 +167,9 @@ void data<T,D>::change_processor (const int newproc, void* const mem=0) {
 
 // Data manipulators
 template<class T, int D>
-void data<T,D>::copy_from (const generic_data<D>* gsrc, const ibbox& box) {
+void data<T,D>
+::copy_from_innerloop (const generic_data<D>* gsrc, const ibbox& box)
+{
   const data* src = (const data*)gsrc;
   assert (has_storage() && src->has_storage());
   assert (all(box.lower()>=extent().lower()
@@ -174,164 +181,371 @@ void data<T,D>::copy_from (const generic_data<D>* gsrc, const ibbox& box) {
   assert (all((box.lower()-extent().lower())%box.stride() == 0
 	      && (box.lower()-src->extent().lower())%box.stride() == 0));
   
-  if (proc() == src->proc()) {
-    // copy on same processor
+  assert (proc() == src->proc());
+  
+  int rank;
+  MPI_Comm_rank (dist::comm, &rank);
+  assert (rank == proc());
+  
+  for (ibbox::iterator it=box.begin(); it!=box.end(); ++it) {
+    const ivect index = *it;
+    (*this)[index] = (*src)[index];
+  }
+  
+}
+
+
+
+template<class T, int D>
+void data<T,D>
+::interpolate_from_innerloop (const generic_data<D>* gsrc, const ibbox& box)
+{
+  const data* src = (const data*)gsrc;
+  assert (has_storage() && src->has_storage());
+  assert (all(box.lower()>=extent().lower()
+	      && box.upper()<=extent().upper()));
+  assert (all(box.lower()>=extent().lower()
+	      && box.lower()>=src->extent().lower()));
+  assert (all(box.upper()<=extent().upper()
+	      && box.upper()<=src->extent().upper()));
+  assert (all(box.stride()==extent().stride()));
+  assert (all((box.lower()-extent().lower())%box.stride() == 0));
+  
+  assert (proc() == src->proc());
+  
+  int rank;
+  MPI_Comm_rank (dist::comm, &rank);
+  assert (rank == proc());
+  
+  for (ibbox::iterator posi=box.begin(); posi!=box.end(); ++posi) {
+    const ivect& pos = *posi;
     
-    int rank;
-    MPI_Comm_rank (dist::comm, &rank);
-    if (rank == proc()) {
-      
-      for (ibbox::iterator it=box.begin(); it!=box.end(); ++it) {
-        const ivect index = *it;
-        (*this)[index] = (*src)[index];
+    // get box around current position
+    const ibbox frombox
+      = ibbox(pos,pos,extent().stride()).expanded_for(src->extent());
+    
+    // interpolate from box to position
+    T sum = 0;
+    for (ibbox::iterator fromposi=frombox.begin();
+	 fromposi!=frombox.end(); ++fromposi)
+      {
+	const ivect& frompos = *fromposi;
+	
+	// interpolation weight
+	const ivect str = src->extent().stride();
+	const T f = prod(vect<T,D>(str - abs(pos - frompos))
+			 / vect<T,D>(str));
+	sum += f * (*src)[frompos];
       }
-      
-    }
+    (*this)[pos] = sum;
     
-  } else {
+  } // for pos
+  
+}
+
+
+
+template<class T, int D>
+void data<T,D>
+::interpolate_from_innerloop (const generic_data<D>* gsrc1, const int t1,
+			      const generic_data<D>* gsrc2, const int t2,
+			      const ibbox& box, const int t)
+{
+  const data* src1 = (const data*)gsrc1;
+  const data* src2 = (const data*)gsrc2;
+  assert (has_storage() && src1->has_storage() && src2->has_storage());
+  assert (all(box.lower()>=extent().lower()
+	      && box.upper()<=extent().upper()));
+  assert (all(box.lower()>=extent().lower()
+	      && box.lower()>=src1->extent().lower()
+	      && box.lower()>=src2->extent().lower()));
+  assert (all(box.upper()<=extent().upper()
+	      && box.upper()<=src1->extent().upper()
+	      && box.upper()<=src2->extent().upper()));
+  assert (all(box.stride()==extent().stride()));
+  assert (all((box.lower()-extent().lower())%box.stride() == 0
+	      && (box.lower()-src1->extent().lower())%box.stride() == 0
+	      && (box.lower()-src2->extent().lower())%box.stride() == 0));
+  assert (src1->proc() == src2->proc());
+  
+  assert (proc() == src1->proc());
+  
+  int rank;
+  MPI_Comm_rank (dist::comm, &rank);
+  assert (rank == proc());
+  
+  const T one = 1;
+  const T src1_fac = (T)((t - t2) * one / (t1 - t2));
+  const T src2_fac = (T)((t - t1) * one / (t2 - t1));
+  
+  for (ibbox::iterator posi=box.begin(); posi!=box.end(); ++posi) {
+    const ivect& pos = *posi;
     
-    // copy to different processor
-    data* const tmp = new data(box, src->proc());
-    tmp->copy_from (src, box);
-    tmp->change_processor (proc());
-    copy_from (tmp, box);
-    delete tmp;
+    // get box around current position
+    const ibbox frombox
+      = ibbox(pos,pos,extent().stride()).expanded_for(src1->extent());
     
-  }
+    // interpolate from box to position
+    T src1_sum = 0;
+    T src2_sum = 0;
+    for (ibbox::iterator fromposi=frombox.begin();
+	 fromposi!=frombox.end(); ++fromposi)
+      {
+	const ivect& frompos = *fromposi;
+	
+	// interpolation weight
+	const ivect str = src1->extent().stride();
+	const T f = prod(vect<T,D>(str - abs(pos - frompos))
+			 / vect<T,D>(str));
+	src1_sum += f * (*src1)[frompos];
+	src2_sum += f * (*src2)[frompos];
+      }
+    (*this)[pos] = src1_fac * src1_sum + src2_fac * src2_sum;
+    
+  } // for pos
+  
 }
 
-template<class T, int D>
-void data<T,D>::interpolate_from (const generic_data<D>* gsrc,
-                                  const ibbox& box)
+
+
+extern "C" {
+  void CCTK_FCALL CCTK_FNAME(copy_3d_real8)
+    (const CCTK_REAL8* src,
+     const int& srciext, const int& srcjext, const int& srckext,
+     CCTK_REAL8* dst,
+     const int& dstiext, const int& dstjext, const int& dstkext,
+     const int srcbbox[3][3],
+     const int dstbbox[3][3],
+     const int regbbox[3][3]);
+}
+
+template<>
+void data<CCTK_REAL8,3>
+::copy_from_innerloop (const generic_data<3>* gsrc, const ibbox& box)
 {
   const data* src = (const data*)gsrc;
   assert (has_storage() && src->has_storage());
   assert (all(box.lower()>=extent().lower()
-	      && box.upper()<=extent().upper()));
-  assert (all(box.lower()>=extent().lower()
 	      && box.lower()>=src->extent().lower()));
   assert (all(box.upper()<=extent().upper()
 	      && box.upper()<=src->extent().upper()));
   assert (all(box.stride()==extent().stride()
-	      /* && box.stride()<=src->extent().stride() */ ));
+	      && box.stride()==src->extent().stride()));
   assert (all((box.lower()-extent().lower())%box.stride() == 0
-	      /* && (box.lower()-src->extent().lower())%box.stride() == 0 */ ));
+	      && (box.lower()-src->extent().lower())%box.stride() == 0));
+  
+  assert (proc() == src->proc());
   
-  if (proc() == src->proc()) {
-    // interpolate on same processor
+  int rank;
+  MPI_Comm_rank (dist::comm, &rank);
+  assert (rank == proc());
+  
+  const ibbox& sext = src->extent();
+  const ibbox& dext = extent();
+  
+  int srcshp[3], dstshp[3];
+  int srcbbox[3][3], dstbbox[3][3], regbbox[3][3];
+  
+  for (int d=0; d<3; ++d) {
+    srcshp[d] = (sext.shape() / sext.stride())[d];
+    dstshp[d] = (dext.shape() / dext.stride())[d];
     
-    int rank;
-    MPI_Comm_rank (dist::comm, &rank);
-    if (rank == proc()) {
-      
-      for (ibbox::iterator posi=box.begin(); posi!=box.end(); ++posi) {
-        const ivect& pos = *posi;
-	
-        // get box around current position
-        const ibbox frombox
-	  = ibbox(pos,pos,extent().stride()).expanded_for(src->extent());
-	
-        // interpolate from box to position
-        T sum = 0;
-        for (ibbox::iterator fromposi=frombox.begin();
-	     fromposi!=frombox.end(); ++fromposi)
-     	  {
-	    const ivect& frompos = *fromposi;
-	    
-	    // interpolation weight
-	    const ivect str = src->extent().stride();
-	    const T f = prod(vect<T,D>(str - abs(pos - frompos))
-			     / vect<T,D>(str));
-	    sum += f * (*src)[frompos];
-	  }
-        (*this)[pos] = sum;
-	
-      } // for pos
-      
-    }
+    srcbbox[0][d] = sext.lower()[d];
+    srcbbox[1][d] = sext.upper()[d];
+    srcbbox[2][d] = sext.stride()[d];
     
-  } else {
-    // interpolate from other processor
+    dstbbox[0][d] = dext.lower()[d];
+    dstbbox[1][d] = dext.upper()[d];
+    dstbbox[2][d] = dext.stride()[d];
     
-    data* const tmp = new data(box, src->proc());
-    tmp->interpolate_from (src, box);
-    tmp->change_processor (proc());
-    copy_from (tmp, box);
-    delete tmp;
+    regbbox[0][d] = box.lower()[d];
+    regbbox[1][d] = box.upper()[d];
+    regbbox[2][d] = box.stride()[d];
+  }
+  
+  assert (all(dext.stride() == box.stride()));
+  if (all(sext.stride() == dext.stride())) {
+    CCTK_FNAME(copy_3d_real8) ((const CCTK_REAL8*)src->storage(),
+			       srcshp[0], srcshp[1], srcshp[2],
+			       (CCTK_REAL8*)storage(),
+			       dstshp[0], dstshp[1], dstshp[2],
+			       srcbbox,
+			       dstbbox,
+			       regbbox);
     
+  } else {
+    abort();
   }
 }
 
-template<class T, int D>
-void data<T,D>::interpolate_from (const generic_data<D>* gsrc,
-                                  const double sfact,
-				  const generic_data<D>* gtrc,
-				  const double tfact,
-				  const ibbox& box)
+
+
+extern "C" {
+  void CCTK_FCALL CCTK_FNAME(prolongate_3d_real8)
+    (const CCTK_REAL8* src,
+     const int& srciext, const int& srcjext, const int& srckext,
+     CCTK_REAL8* dst,
+     const int& dstiext, const int& dstjext, const int& dstkext,
+     const int srcbbox[3][3],
+     const int dstbbox[3][3],
+     const int regbbox[3][3]);
+  void CCTK_FCALL CCTK_FNAME(restrict_3d_real8)
+    (const CCTK_REAL8* src,
+     const int& srciext, const int& srcjext, const int& srckext,
+     CCTK_REAL8* dst,
+     const int& dstiext, const int& dstjext, const int& dstkext,
+     const int srcbbox[3][3],
+     const int dstbbox[3][3],
+     const int regbbox[3][3]);
+}
+
+template<>
+void data<CCTK_REAL8,3>
+::interpolate_from_innerloop (const generic_data<3>* gsrc, const ibbox& box)
 {
   const data* src = (const data*)gsrc;
-  const data* trc = (const data*)gtrc;
-  assert (has_storage() && src->has_storage() && trc->has_storage());
+  assert (has_storage() && src->has_storage());
   assert (all(box.lower()>=extent().lower()
 	      && box.upper()<=extent().upper()));
   assert (all(box.lower()>=extent().lower()
-	      && box.lower()>=src->extent().lower()
-	      && box.lower()>=trc->extent().lower()));
+	      && box.lower()>=src->extent().lower()));
   assert (all(box.upper()<=extent().upper()
-	      && box.upper()<=src->extent().upper()
-	      && box.upper()<=trc->extent().upper()));
-  assert (all(box.stride()==extent().stride()
-	      /* && box.stride()<=src->extent().stride() */
-	      /* && trc->extent().stride()==src->extent().stride() */ ));
-  assert (all((box.lower()-extent().lower())%box.stride() == 0
-	      && (box.lower()-src->extent().lower())%box.stride() == 0
-	      && (box.lower()-trc->extent().lower())%box.stride() == 0));
-  assert (src->proc() == trc->proc());
+	      && box.upper()<=src->extent().upper()));
+  assert (all(box.stride()==extent().stride()));
+  assert (all((box.lower()-extent().lower())%box.stride() == 0));
   
-  if (proc() == src->proc()) {
-    // interpolate on same processor
+  assert (proc() == src->proc());
+  
+  int rank;
+  MPI_Comm_rank (dist::comm, &rank);
+  assert (rank == proc());
+  
+  const ibbox& sext = src->extent();
+  const ibbox& dext = extent();
+  
+  int srcshp[3], dstshp[3];
+  int srcbbox[3][3], dstbbox[3][3], regbbox[3][3];
+  
+  for (int d=0; d<3; ++d) {
+    srcshp[d] = (sext.shape() / sext.stride())[d];
+    dstshp[d] = (dext.shape() / dext.stride())[d];
     
-    int rank;
-    MPI_Comm_rank (dist::comm, &rank);
-    if (rank == proc()) {
-      
-      for (ibbox::iterator posi=box.begin(); posi!=box.end(); ++posi) {
-        const ivect& pos = *posi;
-	
-        // get box around current position
-        const ibbox frombox
-	  = ibbox(pos,pos,extent().stride()).expanded_for(src->extent());
-	
-        // interpolate from box to position
-        T src_sum = 0;
-        T trc_sum = 0;
-        for (ibbox::iterator fromposi=frombox.begin();
-	     fromposi!=frombox.end(); ++fromposi)
-  	    {
-	      const ivect& frompos = *fromposi;
-	      
-	      // interpolation weight
-	      const ivect str = src->extent().stride();
-	      const T f = prod(vect<T,D>(str - abs(pos - frompos))
-			       / vect<T,D>(str));
-	      src_sum += f * (*src)[frompos];
-	      trc_sum += f * (*trc)[frompos];
-	    }
-        (*this)[pos] = (T)sfact * src_sum + (T)tfact * trc_sum;
-	
-      } // for pos
-      
-    }
+    srcbbox[0][d] = sext.lower()[d];
+    srcbbox[1][d] = sext.upper()[d];
+    srcbbox[2][d] = sext.stride()[d];
+    
+    dstbbox[0][d] = dext.lower()[d];
+    dstbbox[1][d] = dext.upper()[d];
+    dstbbox[2][d] = dext.stride()[d];
     
+    regbbox[0][d] = box.lower()[d];
+    regbbox[1][d] = box.upper()[d];
+    regbbox[2][d] = box.stride()[d];
+  }
+  
+  assert (all(dext.stride() == box.stride()));
+  if (all(sext.stride() > dext.stride())) {
+    CCTK_FNAME(prolongate_3d_real8) ((const CCTK_REAL8*)src->storage(),
+				     srcshp[0], srcshp[1], srcshp[2],
+				     (CCTK_REAL8*)storage(),
+				     dstshp[0], dstshp[1], dstshp[2],
+				     srcbbox,
+				     dstbbox,
+				     regbbox);
+  } else if (all(sext.stride() < dext.stride())) {
+    CCTK_FNAME(restrict_3d_real8) ((const CCTK_REAL8*)src->storage(),
+				   srcshp[0], srcshp[1], srcshp[2],
+				   (CCTK_REAL8*)storage(),
+				   dstshp[0], dstshp[1], dstshp[2],
+				   srcbbox,
+				   dstbbox,
+				   regbbox);
   } else {
-    // interpolate from other processor
+    abort();
+  }
+}
+
+
+
+extern "C" {
+  void CCTK_FCALL CCTK_FNAME(prolongate_3d_real8_2tl)
+    (const CCTK_REAL8* src1, const int& t1,
+     const CCTK_REAL8* src2, const int& t2,
+     const int& srciext, const int& srcjext, const int& srckext,
+     CCTK_REAL8* dst, const int& t,
+     const int& dstiext, const int& dstjext, const int& dstkext,
+     const int srcbbox[3][3],
+     const int dstbbox[3][3],
+     const int regbbox[3][3]);
+}
+
+template<>
+void data<CCTK_REAL8,3>
+::interpolate_from_innerloop (const generic_data<3>* gsrc1, const int t1,
+			      const generic_data<3>* gsrc2, const int t2,
+			      const ibbox& box, const int t)
+{
+  const data* src1 = (const data*)gsrc1;
+  const data* src2 = (const data*)gsrc2;
+  assert (has_storage() && src1->has_storage() && src2->has_storage());
+  assert (all(box.lower()>=extent().lower()
+	      && box.upper()<=extent().upper()));
+  assert (all(box.lower()>=extent().lower()
+	      && box.lower()>=src1->extent().lower()
+	      && box.lower()>=src2->extent().lower()));
+  assert (all(box.upper()<=extent().upper()
+	      && box.upper()<=src1->extent().upper()
+	      && box.upper()<=src2->extent().upper()));
+  assert (all(box.stride()==extent().stride()));
+  assert (all((box.lower()-extent().lower())%box.stride() == 0
+	      && (box.lower()-src1->extent().lower())%box.stride() == 0
+	      && (box.lower()-src2->extent().lower())%box.stride() == 0));
+  assert (src1->proc() == src2->proc());
+  
+  assert (proc() == src1->proc());
+  
+  int rank;
+  MPI_Comm_rank (dist::comm, &rank);
+  assert (rank == proc());
+  
+  assert (src1->extent() == src2->extent());
+  
+  const ibbox& sext = src1->extent();
+  const ibbox& dext = extent();
+  
+  int srcshp[3], dstshp[3];
+  int srcbbox[3][3], dstbbox[3][3], regbbox[3][3];
+  
+  for (int d=0; d<3; ++d) {
+    srcshp[d] = (sext.shape() / sext.stride())[d];
+    dstshp[d] = (dext.shape() / dext.stride())[d];
     
-    data* const tmp = new data(box, src->proc());
-    tmp->interpolate_from (src, sfact, trc, tfact, box);
-    tmp->change_processor (proc());
-    copy_from (tmp, box);
-    delete tmp;
+    srcbbox[0][d] = sext.lower()[d];
+    srcbbox[1][d] = sext.upper()[d];
+    srcbbox[2][d] = sext.stride()[d];
     
+    dstbbox[0][d] = dext.lower()[d];
+    dstbbox[1][d] = dext.upper()[d];
+    dstbbox[2][d] = dext.stride()[d];
+    
+    regbbox[0][d] = box.lower()[d];
+    regbbox[1][d] = box.upper()[d];
+    regbbox[2][d] = box.stride()[d];
+  }
+  
+  assert (all(dext.stride() == box.stride()));
+  if (all(sext.stride() > dext.stride())) {
+    CCTK_FNAME(prolongate_3d_real8_2tl)
+      ((const CCTK_REAL8*)src1->storage(), t1,
+       (const CCTK_REAL8*)src2->storage(), t2,
+       srcshp[0], srcshp[1], srcshp[2],
+       (CCTK_REAL8*)storage(), t,
+       dstshp[0], dstshp[1], dstshp[2],
+       srcbbox,
+       dstbbox,
+       regbbox);
+  } else {
+    abort();
   }
 }
 
@@ -361,7 +575,7 @@ ostream& data<T,D>::out (ostream& os) const {
 #if defined(TMPL_EXPLICIT)
 
 #define INSTANTIATE(T)				\
-template data<T,3>;
+template class data<T,3>;
 
 #include "instantiate"
 
diff --git a/Carpet/CarpetLib/src/data.hh b/Carpet/CarpetLib/src/data.hh
index 84e4c92b7..c4594c56d 100644
--- a/Carpet/CarpetLib/src/data.hh
+++ b/Carpet/CarpetLib/src/data.hh
@@ -5,7 +5,7 @@
     copyright            : (C) 2000 by Erik Schnetter
     email                : schnetter@astro.psu.edu
 
-    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/data.hh,v 1.3 2001/03/10 20:55:06 eschnett Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/data.hh,v 1.4 2001/03/22 18:42:05 eschnett Exp $
 
  ***************************************************************************/
 
@@ -21,7 +21,8 @@
 #ifndef DATA_HH
 #define DATA_HH
 
-#include <cassert>
+#include <assert.h>
+
 #include <string>
 
 #include "defs.hh"
@@ -30,6 +31,8 @@
 #include "gdata.hh"
 #include "vect.hh"
 
+using namespace std;
+
 
 
 // Forward definition
@@ -50,7 +53,7 @@ class data: public generic_data<D> {
   typedef bbox<int,D> ibbox;
 
   // Fields
-  T* restrict _storage;		// the data (if located on this processor)
+  T* _storage;			// the data (if located on this processor)
 
 public:
   
@@ -74,12 +77,12 @@ public:
   virtual void change_processor (const int newproc, void* const mem=0);
 
   // Accessors
-  virtual const T* storage () const {
+  virtual const void* storage () const {
     assert (_has_storage);
     return _storage;
   }
 
-  virtual T* storage () {
+  virtual void* storage () {
     assert (_has_storage);
     return _storage;
   }
@@ -96,14 +99,13 @@ public:
   }
   
   // Data manipulators
-  virtual void copy_from (const generic_data<D>* gsrc, const ibbox& b);
-  virtual void interpolate_from (const generic_data<D>* gsrc,
-				 const ibbox& box);
-  virtual void interpolate_from (const generic_data<D>* gsrc,
-                                 const double sfact,
-				 const generic_data<D>* gtrc,
-				 const double tfact,
-				 const ibbox& box);
+  void copy_from_innerloop (const generic_data<D>* src,
+			    const ibbox& box);
+  void interpolate_from_innerloop (const generic_data<D>* src,
+				   const ibbox& box);
+  void interpolate_from_innerloop (const generic_data<D>* src1, const int t1,
+				   const generic_data<D>* src2, const int t2,
+				   const ibbox& box, const int t);
   
   void write_ascii_output_element (ofstream& file, const ivect& index) const;
 //   void write_ieee (const string name, const int time,
diff --git a/Carpet/CarpetLib/src/defs.cc b/Carpet/CarpetLib/src/defs.cc
index bbfdf990a..bf483c62f 100644
--- a/Carpet/CarpetLib/src/defs.cc
+++ b/Carpet/CarpetLib/src/defs.cc
@@ -5,7 +5,7 @@
     copyright            : (C) 2000 by Erik Schnetter
     email                : schnetter@astro.psu.edu
 
-    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/defs.cc,v 1.2 2001/03/19 21:30:19 eschnett Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/defs.cc,v 1.3 2001/03/22 18:42:05 eschnett Exp $
 
  ***************************************************************************/
 
@@ -18,7 +18,8 @@
  *                                                                         *
  ***************************************************************************/
 
-#include <cassert>
+#include <assert.h>
+
 #include <iostream>
 #include <list>
 #include <set>
@@ -28,6 +29,8 @@
 #  include "defs.hh"
 #endif
 
+using namespace std;
+
 
 
 // List output
diff --git a/Carpet/CarpetLib/src/defs.hh b/Carpet/CarpetLib/src/defs.hh
index 9071bad2b..ec59237f8 100644
--- a/Carpet/CarpetLib/src/defs.hh
+++ b/Carpet/CarpetLib/src/defs.hh
@@ -5,7 +5,7 @@
     copyright            : (C) 2000 by Erik Schnetter
     email                : schnetter@astro.psu.edu
 
-    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/defs.hh,v 1.3 2001/03/18 05:20:24 eschnett Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/defs.hh,v 1.4 2001/03/22 18:42:05 eschnett Exp $
 
  ***************************************************************************/
 
@@ -31,6 +31,8 @@
 #include <set>
 #include <vector>
 
+using namespace std;
+
 // Stringification
 #define STR(s) #s
 
diff --git a/Carpet/CarpetLib/src/dh.cc b/Carpet/CarpetLib/src/dh.cc
index 09962b96d..c9e03753d 100644
--- a/Carpet/CarpetLib/src/dh.cc
+++ b/Carpet/CarpetLib/src/dh.cc
@@ -6,7 +6,7 @@
     copyright            : (C) 2000 by Erik Schnetter
     email                : schnetter@astro.psu.edu
 
-    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/dh.cc,v 1.9 2001/03/19 21:30:19 eschnett Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/dh.cc,v 1.10 2001/03/22 18:42:05 eschnett Exp $
 
  ***************************************************************************/
 
@@ -19,7 +19,7 @@
  *                                                                         *
  ***************************************************************************/
 
-#include <cassert>
+#include <assert.h>
 
 #include "defs.hh"
 #include "dist.hh"
@@ -32,6 +32,8 @@
 
 #undef DEBUG_OUTPUT
 
+using namespace std;
+
 
 
 // Constructors
diff --git a/Carpet/CarpetLib/src/dh.hh b/Carpet/CarpetLib/src/dh.hh
index 91634d1ca..aa9e331da 100644
--- a/Carpet/CarpetLib/src/dh.hh
+++ b/Carpet/CarpetLib/src/dh.hh
@@ -6,7 +6,7 @@
     copyright            : (C) 2000 by Erik Schnetter
     email                : schnetter@astro.psu.edu
 
-    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/dh.hh,v 1.3 2001/03/10 20:55:06 eschnett Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/dh.hh,v 1.4 2001/03/22 18:42:05 eschnett Exp $
 
  ***************************************************************************/
 
@@ -22,7 +22,8 @@
 #ifndef DH_HH
 #define DH_HH
 
-#include <cassert>
+#include <assert.h>
+
 #include <iostream>
 #include <list>
 #include <string>
@@ -34,6 +35,8 @@
 #include "gh.hh"
 #include "vect.hh"
 
+using namespace std;
+
 
 
 // Forward declaration
@@ -57,6 +60,7 @@ class dh {
   typedef list<ibbox>    iblist;
   typedef vector<iblist> iblistvect; // vector of lists
   
+public:
   struct dboxes {
     ibbox exterior;		// whole region (including boundaries)
     
@@ -77,6 +81,7 @@ class dh {
     iblistvect recv_ref_bnd_coarse; // received from coarser grids
     ibset sync_not;		// not received while syncing (outer boundary)
   };
+private:
   
   struct dbases {
     ibbox exterior;		// whole region (including boundaries)
diff --git a/Carpet/CarpetLib/src/dist.cc b/Carpet/CarpetLib/src/dist.cc
index 9aa613764..99cc2e83a 100644
--- a/Carpet/CarpetLib/src/dist.cc
+++ b/Carpet/CarpetLib/src/dist.cc
@@ -5,7 +5,7 @@
     copyright            : (C) 2000 by Erik Schnetter
     email                : schnetter@astro.psu.edu
 
-    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/dist.cc,v 1.3 2001/03/07 13:00:57 eschnett Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/dist.cc,v 1.4 2001/03/22 18:42:05 eschnett Exp $
 
  ***************************************************************************/
 
@@ -18,7 +18,7 @@
  *                                                                         *
  ***************************************************************************/
 
-#include <cassert>
+#include <assert.h>
 
 #include <mpi.h>
 
@@ -31,6 +31,8 @@
 #  include "dist.hh"
 #endif
 
+using namespace std;
+
 
 
 namespace dist {
diff --git a/Carpet/CarpetLib/src/dist.hh b/Carpet/CarpetLib/src/dist.hh
index dac83076f..9c488e0fd 100644
--- a/Carpet/CarpetLib/src/dist.hh
+++ b/Carpet/CarpetLib/src/dist.hh
@@ -5,7 +5,7 @@
     copyright            : (C) 2000 by Erik Schnetter
     email                : schnetter@astro.psu.edu
 
-    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/dist.hh,v 1.3 2001/03/07 13:00:57 eschnett Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/dist.hh,v 1.4 2001/03/22 18:42:05 eschnett Exp $
 
  ***************************************************************************/
 
@@ -21,9 +21,9 @@
 #ifndef DIST_HH
 #define DIST_HH
 
-#include <cassert>
-#include <cstdio>
-#include <cstdlib>
+#include <assert.h>
+#include <stdio.h>
+#include <stdlib.h>
 
 #include <complex>
 
@@ -31,6 +31,8 @@
 
 #include "defs.hh"
 
+using namespace std;
+
 
 
 namespace dist {
diff --git a/Carpet/CarpetLib/src/gdata.cc b/Carpet/CarpetLib/src/gdata.cc
index eb19701ce..5418db23a 100644
--- a/Carpet/CarpetLib/src/gdata.cc
+++ b/Carpet/CarpetLib/src/gdata.cc
@@ -5,7 +5,7 @@
     copyright            : (C) 2000 by Erik Schnetter
     email                : schnetter@astro.psu.edu
 
-    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gdata.cc,v 1.8 2001/03/19 21:30:19 eschnett Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gdata.cc,v 1.9 2001/03/22 18:42:05 eschnett Exp $
 
  ***************************************************************************/
 
@@ -18,7 +18,8 @@
  *                                                                         *
  ***************************************************************************/
 
-#include <cassert>
+#include <assert.h>
+
 #include <fstream>
 #include <iomanip>
 
@@ -31,6 +32,8 @@
 #  include "gdata.hh"
 #endif
 
+using namespace std;
+
 
 
 // Constructors
@@ -45,6 +48,127 @@ generic_data<D>::~generic_data () { }
 
 
 
+// Data manipulators
+template<int D>
+void generic_data<D>::copy_from (const generic_data* src, const ibbox& box)
+{
+  assert (has_storage() && src->has_storage());
+  assert (all(box.lower()>=extent().lower()
+	      && box.lower()>=src->extent().lower()));
+  assert (all(box.upper()<=extent().upper()
+	      && box.upper()<=src->extent().upper()));
+  assert (all(box.stride()==extent().stride()
+	      && box.stride()==src->extent().stride()));
+  assert (all((box.lower()-extent().lower())%box.stride() == 0
+	      && (box.lower()-src->extent().lower())%box.stride() == 0));
+  
+  if (proc() == src->proc()) {
+    // copy on same processor
+    
+    int rank;
+    MPI_Comm_rank (dist::comm, &rank);
+    if (rank == proc()) {
+      copy_from_innerloop (src, box);
+    }
+    
+  } else {
+    
+    // copy to different processor
+    generic_data* const tmp = make_typed();
+    tmp->allocate (box, src->proc());
+    tmp->copy_from (src, box);
+    tmp->change_processor (proc());
+    copy_from (tmp, box);
+    delete tmp;
+    
+  }
+}
+
+
+
+template<int D>
+void generic_data<D>
+::interpolate_from (const generic_data* src, const ibbox& box)
+{
+  assert (has_storage() && src->has_storage());
+  assert (all(box.lower()>=extent().lower()
+	      && box.upper()<=extent().upper()));
+  assert (all(box.lower()>=extent().lower()
+	      && box.lower()>=src->extent().lower()));
+  assert (all(box.upper()<=extent().upper()
+	      && box.upper()<=src->extent().upper()));
+  assert (all(box.stride()==extent().stride()));
+  assert (all((box.lower()-extent().lower())%box.stride() == 0));
+  
+  if (proc() == src->proc()) {
+    // interpolate on same processor
+    
+    int rank;
+    MPI_Comm_rank (dist::comm, &rank);
+    if (rank == proc()) {
+      interpolate_from_innerloop (src, box);
+    }
+    
+  } else {
+    // interpolate from other processor
+    
+    generic_data* const tmp = make_typed();
+    tmp->allocate (box, src->proc());
+    tmp->interpolate_from (src, box);
+    tmp->change_processor (proc());
+    copy_from (tmp, box);
+    delete tmp;
+    
+  }
+}
+
+
+
+template<int D>
+void generic_data<D>
+::interpolate_from (const generic_data* src1, const int t1,
+		    const generic_data* src2, const int t2,
+		    const ibbox& box, const int t)
+{
+  assert (has_storage() && src1->has_storage() && src2->has_storage());
+  assert (all(box.lower()>=extent().lower()
+	      && box.upper()<=extent().upper()));
+  assert (all(box.lower()>=extent().lower()
+	      && box.lower()>=src1->extent().lower()
+	      && box.lower()>=src2->extent().lower()));
+  assert (all(box.upper()<=extent().upper()
+	      && box.upper()<=src1->extent().upper()
+	      && box.upper()<=src2->extent().upper()));
+  assert (all(box.stride()==extent().stride()));
+  assert (all((box.lower()-extent().lower())%box.stride() == 0
+	      && (box.lower()-src1->extent().lower())%box.stride() == 0
+	      && (box.lower()-src2->extent().lower())%box.stride() == 0));
+  assert (src1->proc() == src2->proc());
+  
+  if (proc() == src1->proc()) {
+    // interpolate on same processor
+    
+    int rank;
+    MPI_Comm_rank (dist::comm, &rank);
+    if (rank == proc()) {
+      interpolate_from_innerloop (src1, t1, src2, t2, box, t);
+    }
+    
+  } else {
+    // interpolate from other processor
+    
+    generic_data* const tmp = make_typed();
+    tmp->allocate (box, src1->proc());
+    tmp->interpolate_from (src1, t1, src2, t2, box, t);
+    tmp->change_processor (proc());
+    copy_from (tmp, box);
+    delete tmp;
+    
+  }
+}
+
+
+
 // Output
 template<int D>
 template<int DD>
@@ -75,7 +199,8 @@ void generic_data<D>::write_ascii (const string name, const int time,
 	   << "   component " << c << "   multigrid level " << ml << endl
 	   << "# column format: it tl rl c ml";
       assert (D>=1 && D<=3);
-      for (int d=0; d<D; ++d) file << " " << "xyz"[d];
+      const char* const coords = "xyz";
+      for (int d=0; d<D; ++d) file << " " << coords[d];
       file << " data" << endl;
       
       const vect<int,DD> lo = extent().lower()[dirs];
diff --git a/Carpet/CarpetLib/src/gdata.hh b/Carpet/CarpetLib/src/gdata.hh
index e765cc7c6..1a16f92c2 100644
--- a/Carpet/CarpetLib/src/gdata.hh
+++ b/Carpet/CarpetLib/src/gdata.hh
@@ -5,7 +5,7 @@
     copyright            : (C) 2000 by Erik Schnetter
     email                : schnetter@astro.psu.edu
 
-    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gdata.hh,v 1.4 2001/03/14 11:00:26 eschnett Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gdata.hh,v 1.5 2001/03/22 18:42:05 eschnett Exp $
 
  ***************************************************************************/
 
@@ -21,8 +21,9 @@
 #ifndef GDATA_HH
 #define GDATA_HH
 
-#include <cassert>
-#include <cstdlib>
+#include <assert.h>
+#include <stdlib.h>
+
 #include <fstream>
 #include <iostream>
 #include <string>
@@ -32,6 +33,8 @@
 #include "bbox.hh"
 #include "vect.hh"
 
+using namespace std;
+
 
 
 // Forward declaration
@@ -131,13 +134,22 @@ public:
   }
 
   // Data manipulators
-  virtual void copy_from (const generic_data* src,
-                          const ibbox& b) = 0;
-  virtual void interpolate_from (const generic_data* src,
-				 const ibbox& box) = 0;
-  virtual void interpolate_from (const generic_data* src, const double sfact,
-				 const generic_data* trc, const double tfact,
-				 const ibbox& box) = 0;
+  void copy_from (const generic_data* src, const ibbox& box);
+  void interpolate_from (const generic_data* src, const ibbox& box);
+  void interpolate_from (const generic_data* src1, const int t1,
+			 const generic_data* src2, const int t2,
+			 const ibbox& box, const int t);
+protected:
+  virtual void
+  copy_from_innerloop (const generic_data* src, const ibbox& box) = 0;
+  virtual void
+  interpolate_from_innerloop (const generic_data* src, const ibbox& box) =0;
+  virtual void
+  interpolate_from_innerloop (const generic_data* src1, const int t1,
+			      const generic_data* src2, const int t2,
+			      const ibbox& box, const int t) = 0;
+  
+public:
   
   // Output
   template<int DD>
diff --git a/Carpet/CarpetLib/src/gf.cc b/Carpet/CarpetLib/src/gf.cc
index 86a2d518d..6847d8506 100644
--- a/Carpet/CarpetLib/src/gf.cc
+++ b/Carpet/CarpetLib/src/gf.cc
@@ -6,7 +6,7 @@
     copyright            : (C) 2000 by Erik Schnetter
     email                : schnetter@astro.psu.edu
 
-    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gf.cc,v 1.2 2001/03/19 21:30:19 eschnett Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gf.cc,v 1.3 2001/03/22 18:42:05 eschnett Exp $
 
  ***************************************************************************/
 
@@ -19,7 +19,7 @@
  *                                                                         *
  ***************************************************************************/
 
-#include <cassert>
+#include <assert.h>
 
 #include "defs.hh"
 
@@ -27,6 +27,8 @@
 #  include "gf.hh"
 #endif
 
+using namespace std;
+
 
 
 // Constructors
@@ -81,6 +83,7 @@ ostream& gf<T,D>::out (ostream& os) const {
 template class gf<T,3>;
 
 #include "instantiate"
+
 #undef INSTANTIATE
 
 #endif
diff --git a/Carpet/CarpetLib/src/gf.hh b/Carpet/CarpetLib/src/gf.hh
index 298c26589..bf916bb82 100644
--- a/Carpet/CarpetLib/src/gf.hh
+++ b/Carpet/CarpetLib/src/gf.hh
@@ -6,7 +6,7 @@
     copyright            : (C) 2000 by Erik Schnetter
     email                : schnetter@astro.psu.edu
 
-    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gf.hh,v 1.1 2001/03/01 13:40:10 eschnett Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gf.hh,v 1.2 2001/03/22 18:42:05 eschnett Exp $
 
  ***************************************************************************/
 
@@ -22,8 +22,9 @@
 #ifndef GF_HH
 #define GF_HH
 
-#include <cassert>
-#include <cmath>
+#include <assert.h>
+#include <math.h>
+
 #include <iostream>
 #include <string>
 
@@ -36,6 +37,8 @@
 #include "th.hh"
 #include "vect.hh"
 
+using namespace std;
+
 
 
 // A real grid function
diff --git a/Carpet/CarpetLib/src/ggf.cc b/Carpet/CarpetLib/src/ggf.cc
index 53aa15179..09966ac69 100644
--- a/Carpet/CarpetLib/src/ggf.cc
+++ b/Carpet/CarpetLib/src/ggf.cc
@@ -6,7 +6,7 @@
     copyright            : (C) 2000 by Erik Schnetter
     email                : schnetter@astro.psu.edu
 
-    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/ggf.cc,v 1.5 2001/03/19 21:30:19 eschnett Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/ggf.cc,v 1.6 2001/03/22 18:42:05 eschnett Exp $
 
  ***************************************************************************/
 
@@ -19,7 +19,8 @@
  *                                                                         *
  ***************************************************************************/
 
-#include <cassert>
+#include <assert.h>
+
 #include <iostream>
 #include <string>
 
@@ -31,6 +32,8 @@
 #  include "ggf.hh"
 #endif
 
+using namespace std;
+
 
 
 // Constructors
@@ -150,6 +153,19 @@ void generic_gf<D>::recompose () {
   } // for tl
 }
 
+// Cycle the time levels by rotating the data sets
+template<int D>
+void generic_gf<D>::cycle (int rl, int c, int ml) {
+  assert (rl>=0 && rl<h.reflevels());
+  assert (c>=0 && c<h.components(rl));
+  assert (ml>=0 && ml<h.mglevels(rl,c));
+  generic_data<D>* tmpdata = storage[tmin-tmin][rl][c][ml];
+  for (int tl=tmin; tl<=tmax-1; ++tl) {
+    storage[tl-tmin][rl][c][ml] = storage[tl+1-tmin][rl][c][ml];
+  }
+  storage[tmax-tmin][rl][c][ml] = tmpdata;
+}
+
 
 
 // Operations
@@ -240,17 +256,15 @@ void generic_gf<D>::copycat (int tl1, int rl1, int c1, int ml1,
 template<int D>
 void generic_gf<D>::intercat (int tl1, int rl1, int c1, int ml1,
 			      const ibbox dh<D>::dboxes::* recv_list,
-			      int tl2, const double fact2,
-			      int tl3, const double fact3,
-			      int rl2, int ml2,
+			      int tl2a, int tl2b, int rl2, int ml2,
 			      const ibbox dh<D>::dboxes::* send_list)
 {
   assert (tl1>=tmin && tl1<=tmax);
   assert (rl1>=0 && rl1<h.reflevels());
   assert (c1>=0 && c1<h.components(rl1));
   assert (ml1>=0 && ml1<h.mglevels(rl1,c1));
-  assert (tl2>=tmin && tl2<=tmax);
-  assert (tl3>=tmin && tl3<=tmax);
+  assert (tl2a>=tmin && tl2a<=tmax);
+  assert (tl2b>=tmin && tl2b<=tmax);
   assert (rl2>=0 && rl2<h.reflevels());
   const int c2=c1;
   assert (ml2<h.mglevels(rl2,c2));
@@ -260,25 +274,23 @@ void generic_gf<D>::intercat (int tl1, int rl1, int c1, int ml1,
   // interpolate the content
   assert (recv==send);
   storage[tl1-tmin][rl1][c1][ml1]->interpolate_from
-    (storage[tl2-tmin][rl2][c2][ml2], fact2,
-     storage[tl3-tmin][rl2][c2][ml2], fact3, recv);
+    (storage[tl2a-tmin][rl2][c2][ml2], tl2a,
+     storage[tl2b-tmin][rl2][c2][ml2], tl2b, recv, tl1);
 }
 
 // Interpolate a component (between multigrid levels)
 template<int D>
 void generic_gf<D>::intercat (int tl1, int rl1, int c1, int ml1,
 			      const iblist dh<D>::dboxes::* recv_list,
-			      int tl2, const double fact2,
-			      int tl3, const double fact3,
-			      int rl2, int ml2,
+			      int tl2a, int tl2b, int rl2, int ml2,
 			      const iblist dh<D>::dboxes::* send_list)
 {
   assert (tl1>=tmin && tl1<=tmax);
   assert (rl1>=0 && rl1<h.reflevels());
   assert (c1>=0 && c1<h.components(rl1));
   assert (ml1>=0 && ml1<h.mglevels(rl1,c1));
-  assert (tl2>=tmin && tl2<=tmax);
-  assert (tl3>=tmin && tl3<=tmax);
+  assert (tl2a>=tmin && tl2a<=tmax);
+  assert (tl2b>=tmin && tl2b<=tmax);
   assert (rl2>=0 && rl2<h.reflevels());
   const int c2=c1;
   assert (ml2<h.mglevels(rl2,c2));
@@ -291,8 +303,8 @@ void generic_gf<D>::intercat (int tl1, int rl1, int c1, int ml1,
     // (use the send boxes for communication)
     // interpolate the content
     storage[tl1-tmin][rl1][c1][ml1]->interpolate_from
-      (storage[tl2-tmin][rl2][c2][ml2], fact2,
-       storage[tl3-tmin][rl2][c2][ml2], fact3, *r);
+      (storage[tl2a-tmin][rl2][c2][ml2], tl2a,
+       storage[tl2b-tmin][rl2][c2][ml2], tl2b, *r, tl1);
   }
 }
 
@@ -300,18 +312,16 @@ void generic_gf<D>::intercat (int tl1, int rl1, int c1, int ml1,
 template<int D>
 void generic_gf<D>::intercat (int tl1, int rl1, int c1, int ml1,
 			      const iblistvect dh<D>::dboxes::* recv_listvect,
-			      int tl2, const double fact2,
-			      int tl3, const double fact3,
-			      int rl2, int ml2,
+			      int tl2a, int tl2b, int rl2, int ml2,
 			      const iblistvect dh<D>::dboxes::* send_listvect)
 {
   assert (tl1>=tmin && tl1<=tmax);
   assert (rl1>=0 && rl1<h.reflevels());
   assert (c1>=0 && c1<h.components(rl1));
   assert (ml1>=0 && ml1<h.mglevels(rl1,c1));
-  assert (tl2>=tmin && tl2<=tmax);
+  assert (tl2a>=tmin && tl2a<=tmax);
+  assert (tl2b>=tmin && tl2b<=tmax);
   assert (rl2>=0 && rl2<h.reflevels());
-  assert (tl3>=tmin && tl3<=tmax);
   // walk all components
   for (int c2=0; c2<h.components(rl2); ++c2) {
     assert (ml2<h.mglevels(rl2,c2));
@@ -324,8 +334,8 @@ void generic_gf<D>::intercat (int tl1, int rl1, int c1, int ml1,
       // (use the send boxes for communication)
       // interpolate the content
       storage[tl1-tmin][rl1][c1][ml1]->interpolate_from
-      	(storage[tl2-tmin][rl2][c2][ml2], fact2,
-      	 storage[tl3-tmin][rl2][c2][ml2], fact3, *r);
+      	(storage[tl2a-tmin][rl2][c2][ml2], tl2a,
+      	 storage[tl2b-tmin][rl2][c2][ml2], tl2b, *r, tl1);
     }
   }
 }
@@ -361,13 +371,11 @@ void generic_gf<D>::ref_bnd_prolongate (int tl, int rl, int c, int ml) {
     copycat (tl ,rl  ,c,ml, &dh<D>::dboxes::recv_ref_bnd_coarse,
       	     tl2,rl-1,  ml, &dh<D>::dboxes::send_ref_bnd_fine);
   } else {
-    int tl3=tl2+1;
-    assert (tl3>=tmin && tl3<=tmax);
-    const double fact2 = 1 - (time - tl2);
-    const double fact3 = 1 - fact2;
+    const int tl2a=tl2;
+    const int tl2b=tl2+1;
+    assert (tl2b>=tmin && tl2b<=tmax);
     intercat (tl,rl,c,ml, &dh<D>::dboxes::recv_ref_bnd_coarse,
-      	      tl2,fact2, tl3,fact3,
-      	      rl-1,ml, &dh<D>::dboxes::send_ref_bnd_fine);
+      	      tl2a, tl2b, rl-1,ml, &dh<D>::dboxes::send_ref_bnd_fine);
   }
 }
 
diff --git a/Carpet/CarpetLib/src/ggf.hh b/Carpet/CarpetLib/src/ggf.hh
index f59d4e415..6e998c145 100644
--- a/Carpet/CarpetLib/src/ggf.hh
+++ b/Carpet/CarpetLib/src/ggf.hh
@@ -6,7 +6,7 @@
     copyright            : (C) 2000 by Erik Schnetter
     email                : schnetter@astro.psu.edu
 
-    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/ggf.hh,v 1.2 2001/03/15 09:59:43 eschnett Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/ggf.hh,v 1.3 2001/03/22 18:42:06 eschnett Exp $
 
  ***************************************************************************/
 
@@ -22,7 +22,8 @@
 #ifndef GGF_HH
 #define GGF_HH
 
-#include <cassert>
+#include <assert.h>
+
 #include <iostream>
 #include <string>
 
@@ -32,6 +33,8 @@
 #include "gh.hh"
 #include "th.hh"
 
+using namespace std;
+
 
 
 // Forward declaration
@@ -83,12 +86,15 @@ public:
   virtual ~generic_gf ();
 
   // Comparison
-  virtual bool operator== (const generic_gf<D>& f) const;
+  bool operator== (const generic_gf<D>& f) const;
 
 
 
   // Modifiers
-  virtual void recompose ();
+  void recompose ();
+
+  // Cycle the time levels by rotating the data sets
+  void cycle (int rl, int c, int ml);
   
   
   
@@ -105,46 +111,40 @@ protected:
 protected:
   
   // Copy region for a component (between time levels)
-  virtual void copycat (int tl1, int rl1, int c1, int ml1,
-			const ibbox dh<D>::dboxes::* recv_list,
-			int tl2, int rl2, int ml2,
-			const ibbox dh<D>::dboxes::* send_list);
+  void copycat (int tl1, int rl1, int c1, int ml1,
+		const ibbox dh<D>::dboxes::* recv_list,
+		int tl2, int rl2, int ml2,
+		const ibbox dh<D>::dboxes::* send_list);
 
   // Copy regions for a component (between multigrid levels)
-  virtual void copycat (int tl1, int rl1, int c1, int ml1,
-			const iblist dh<D>::dboxes::* recv_list,
-			int tl2, int rl2, int ml2,
-			const iblist dh<D>::dboxes::* send_list);
+  void copycat (int tl1, int rl1, int c1, int ml1,
+		const iblist dh<D>::dboxes::* recv_list,
+		int tl2, int rl2, int ml2,
+		const iblist dh<D>::dboxes::* send_list);
 
   // Copy regions for a level (between refinement levels)
-  virtual void copycat (int tl1, int rl1, int c1, int ml1,
-			const iblistvect dh<D>::dboxes::* recv_listvect,
-			int tl2, int rl2, int ml2,
-			const iblistvect dh<D>::dboxes::* send_listvect);
+  void copycat (int tl1, int rl1, int c1, int ml1,
+		const iblistvect dh<D>::dboxes::* recv_listvect,
+		int tl2, int rl2, int ml2,
+		const iblistvect dh<D>::dboxes::* send_listvect);
   
   // Interpolate a component (between time levels)
-  virtual void intercat (int tl1, int rl1, int c1, int ml1,
-			 const ibbox dh<D>::dboxes::* recv_list,
-			 int tl2, const double fact2,
-			 int tl3, const double fact3,
-			 int rl2, int ml2,
-			 const ibbox dh<D>::dboxes::* send_list);
+  void intercat (int tl1, int rl1, int c1, int ml1,
+		 const ibbox dh<D>::dboxes::* recv_list,
+		 int tl2a, int tl2b, int rl2, int ml2,
+		 const ibbox dh<D>::dboxes::* send_list);
 
   // Interpolate a component (between multigrid levels)
-  virtual void intercat (int tl1, int rl1, int c1, int ml1,
-			 const iblist dh<D>::dboxes::* recv_list,
-			 int tl2, const double fact2,
-			 int tl3, const double fact3,
-			 int rl2, int ml2,
-			 const iblist dh<D>::dboxes::* send_list);
+  void intercat (int tl1, int rl1, int c1, int ml1,
+		 const iblist dh<D>::dboxes::* recv_list,
+		 int tl2a, int tl2b, int rl2, int ml2,
+		 const iblist dh<D>::dboxes::* send_list);
 
   // Interpolate a level (between refinement levels)
-  virtual void intercat (int tl1, int rl1, int c1, int ml1,
-			 const iblistvect dh<D>::dboxes::* recv_listvect,
-			 int tl2, const double fact2,
-			 int tl3, const double fact3,
-			 int rl2, int ml2,
-			 const iblistvect dh<D>::dboxes::* send_listvect);
+  void intercat (int tl1, int rl1, int c1, int ml1,
+		 const iblistvect dh<D>::dboxes::* recv_listvect,
+		 int tl2a, int tl2b, int rl2, int ml2,
+		 const iblistvect dh<D>::dboxes::* send_listvect);
 
 
 
@@ -157,25 +157,25 @@ public:
   // synchronised.  They don't need to be prolongated.
 
   // Copy a component from the next time level
-  virtual void copy (int tl, int rl, int c, int ml);
+  void copy (int tl, int rl, int c, int ml);
 
   // Synchronise the boundaries of a component
-  virtual void sync (int tl, int rl, int c, int ml);
+  void sync (int tl, int rl, int c, int ml);
 
   // Prolongate the boundaries of a component
-  virtual void ref_bnd_prolongate (int tl, int rl, int c, int ml);
+  void ref_bnd_prolongate (int tl, int rl, int c, int ml);
 
   // Restrict a multigrid level
-  virtual void mg_restrict (int tl, int rl, int c, int ml);
+  void mg_restrict (int tl, int rl, int c, int ml);
 
   // Prolongate a multigrid level
-  virtual void mg_prolongate (int tl, int rl, int c, int ml);
+  void mg_prolongate (int tl, int rl, int c, int ml);
 
   // Restrict a refinement level
-  virtual void ref_restrict (int tl, int rl, int c, int ml);
+  void ref_restrict (int tl, int rl, int c, int ml);
 
   // Prolongate a refinement level
-  virtual void ref_prolongate (int tl, int rl, int c, int ml);
+  void ref_prolongate (int tl, int rl, int c, int ml);
   
   
   
diff --git a/Carpet/CarpetLib/src/gh.cc b/Carpet/CarpetLib/src/gh.cc
index b1608b6b7..ee10251b1 100644
--- a/Carpet/CarpetLib/src/gh.cc
+++ b/Carpet/CarpetLib/src/gh.cc
@@ -7,7 +7,7 @@
     copyright            : (C) 2000 by Erik Schnetter
     email                : schnetter@astro.psu.edu
 
-    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gh.cc,v 1.4 2001/03/19 21:30:19 eschnett Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gh.cc,v 1.5 2001/03/22 18:42:06 eschnett Exp $
 
  ***************************************************************************/
 
@@ -20,8 +20,8 @@
  *                                                                         *
  ***************************************************************************/
 
-#include <cassert>
-#include <cstdlib>
+#include <assert.h>
+#include <stdlib.h>
 #include <iostream>
 
 #include "defs.hh"
@@ -32,6 +32,8 @@
 #  include "gh.hh"
 #endif
 
+using namespace std;
+
 
 
   // Constructors
diff --git a/Carpet/CarpetLib/src/gh.hh b/Carpet/CarpetLib/src/gh.hh
index ff90a9476..841abbdfb 100644
--- a/Carpet/CarpetLib/src/gh.hh
+++ b/Carpet/CarpetLib/src/gh.hh
@@ -7,7 +7,7 @@
     copyright            : (C) 2000 by Erik Schnetter
     email                : schnetter@astro.psu.edu
 
-    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gh.hh,v 1.2 2001/03/07 13:00:57 eschnett Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gh.hh,v 1.3 2001/03/22 18:42:06 eschnett Exp $
 
  ***************************************************************************/
 
@@ -23,7 +23,8 @@
 #ifndef GH_HH
 #define GH_HH
 
-#include <cassert>
+#include <assert.h>
+
 #include <iostream>
 #include <list>
 #include <vector>
@@ -33,6 +34,8 @@
 #include "dist.hh"
 #include "vect.hh"
 
+using namespace std;
+
 
 
 // Forward declaration
diff --git a/Carpet/CarpetLib/src/instantiate b/Carpet/CarpetLib/src/instantiate
index aafc645d2..f85d27c58 100644
--- a/Carpet/CarpetLib/src/instantiate
+++ b/Carpet/CarpetLib/src/instantiate
@@ -1,7 +1,7 @@
 // Instantiate templates for all available types    -*-C++-*-
 // (C) 2001 Erik Schnetter <schnetter@uni-tuebingen.de>
 
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/instantiate,v 1.2 2001/03/19 21:30:19 eschnett Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/instantiate,v 1.3 2001/03/22 18:42:06 eschnett Exp $
 
 // Usage:
 // Define the macro INSTANTIATE(T) to instantiate for the type T,
@@ -55,7 +55,9 @@
 #endif
 
 #if !defined(CARPET_BYTE) && !defined(CARPET_INT2) && !defined(CARPET_INT4) && !defined(CARPET_INT8) && !defined(CARPET_REAL4) && !defined(CARPET_REAL8) && !defined(CARPET_REAL16) && !defined(CARPET_COMPLEX8) && !defined(CARPET_COMPLEX16) && !defined(CARPET_COMPLEX32)
-// Assume the user just wants REAL
+// Assume the user just wants INT and REAL
+#  undef CARPET_INT
+#  define CARPET_INT
 #  undef CARPET_REAL
 #  define CARPET_REAL
 #endif
diff --git a/Carpet/CarpetLib/src/make.code.defn b/Carpet/CarpetLib/src/make.code.defn
index 121d29247..dd1756780 100644
--- a/Carpet/CarpetLib/src/make.code.defn
+++ b/Carpet/CarpetLib/src/make.code.defn
@@ -1,8 +1,9 @@
 # Main make.code.defn file for thorn CarpetLib   -*-Makefile-*-
-# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/make.code.defn,v 1.2 2001/03/17 22:37:56 eschnett Exp $
+# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/make.code.defn,v 1.3 2001/03/22 18:42:06 eschnett Exp $
 
 # Source files in this directory
-SRCS = bbox.cc bboxset.cc data.cc defs.cc dh.cc dist.cc gdata.cc gf.cc ggf.cc gh.cc th.cc vect.cc
+SRCS = bbox.cc bboxset.cc data.cc defs.cc dh.cc dist.cc gdata.cc gf.cc ggf.cc gh.cc th.cc vect.cc \
+	copy_3d_real8.F77 prolongate_3d_real8.F77 prolongate_3d_real8_2tl.F77 restrict_3d_real8.F77
 # Note: skipping io.cc
 
 # Subdirectories containing source files
diff --git a/Carpet/CarpetLib/src/prolongate_3d_real8.F77 b/Carpet/CarpetLib/src/prolongate_3d_real8.F77
index 8eae50332..04053c524 100644
--- a/Carpet/CarpetLib/src/prolongate_3d_real8.F77
+++ b/Carpet/CarpetLib/src/prolongate_3d_real8.F77
@@ -1,21 +1,8 @@
 c     -*-Fortran-*-
-c     $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8.F77,v 1.11 2004/03/11 12:03:09 schnetter Exp $
+c     $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8.F77,v 1.2 2001/03/22 18:42:06 eschnett Exp $
       
 #include "cctk.h"
-
-
-
-#define CHKIDX(i,j,k, imax,jmax,kmax, where)                                  \
-      if ((i).lt.1 .or. (i).gt.(imax)                                         \
-           .or. (j).lt.1 .or. (j).gt.(jmax)                                   \
-           .or. (k).lt.1 .or. (k).gt.(kmax)) then                           &&\
-         write (msg, '(a, " array index out of bounds: shape is (",i4,",",i4,",",i4,"), index is (",i4,",",i4,",",i4,")")') \
-              (where), (imax), (jmax), (kmax), (i), (j), (k)                &&\
-         call CCTK_WARN (0, msg(1:len_trim(msg)))                           &&\
-      end if
-
-
-
+      
       subroutine prolongate_3d_real8 (
      $     src, srciext, srcjext, srckext,
      $     dst, dstiext, dstjext, dstkext,
@@ -23,9 +10,6 @@ c     $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d
       
       implicit none
       
-      CCTK_REAL8 one
-      parameter (one = 1)
-      
       integer srciext, srcjext, srckext
       CCTK_REAL8 src(srciext,srcjext,srckext)
       integer dstiext, dstjext, dstkext
@@ -42,18 +26,14 @@ c     bbox(:,3) is stride
       integer srcioff, srcjoff, srckoff
       integer dstioff, dstjoff, dstkoff
       
-      CCTK_REAL8 dstdiv
       integer i, j, k
       integer i0, j0, k0
       integer fi, fj, fk
-      integer ifac(2), jfac(2), kfac(2)
       integer ii, jj, kk
       integer fac
       CCTK_REAL8 res
       integer d
       
-      character msg*1000
-      
       
       
       do d=1,3
@@ -69,25 +49,10 @@ c     bbox(:,3) is stride
             call CCTK_WARN (0, "Internal error: destination strides are not integer multiples of the source strides")
          end if
          if (mod(srcbbox(d,1), srcbbox(d,3)).ne.0
-     $        .or. mod(dstbbox(d,1), dstbbox(d,3)).ne.0
+     $     .or. mod(dstbbox(d,1), dstbbox(d,3)).ne.0
      $        .or. mod(regbbox(d,1), regbbox(d,3)).ne.0) then
             call CCTK_WARN (0, "Internal error: array origins are not integer multiples of the strides")
          end if
-         if (regbbox(d,1).gt.regbbox(d,2)) then
-c     This could be handled, but is likely to point to an error elsewhere
-            call CCTK_WARN (0, "Internal error: region extent is empty")
-         end if
-         if (mod(srcbbox(d,2) - srcbbox(d,1), srcbbox(d,3)).ne.0
-     $        .or. mod(dstbbox(d,2) - dstbbox(d,1), dstbbox(d,3)).ne.0
-     $        .or. mod(regbbox(d,2) - regbbox(d,1), regbbox(d,3)).ne.0) then
-            call CCTK_WARN (0, "Internal error: array extents are not integer multiples of the strides")
-         end if
-         if (regbbox(d,1).lt.srcbbox(d,1)
-     $        .or. regbbox(d,1).lt.dstbbox(d,1)
-     $        .or. regbbox(d,2).gt.srcbbox(d,2)
-     $        .or. regbbox(d,2).gt.dstbbox(d,2)) then
-            call CCTK_WARN (0, "Internal error: region extent is not contained in array extent")
-         end if
       end do
       
       if (srciext.ne.(srcbbox(1,2)-srcbbox(1,1))/srcbbox(1,3)+1
@@ -120,71 +85,52 @@ c     This could be handled, but is likely to point to an error elsewhere
       
       
 c     Loop over fine region
-      dstdiv = one / (dstifac * dstjfac * dstkfac)
-      
       do k = 0, regkext-1
-         k0 = (srckoff + k) / dstkfac
-         fk = mod(srckoff + k, dstkfac)
-         kfac(1) = (fk-dstkfac) * (-1)
-         kfac(2) = (fk        ) * 1
-         
          do j = 0, regjext-1
-            j0 = (srcjoff + j) / dstjfac
-            fj = mod(srcjoff + j, dstjfac)
-            jfac(1) = (fj-dstjfac) * (-1)
-            jfac(2) = (fj        ) * 1
-            
             do i = 0, regiext-1
-               i0 = (srcioff + i) / dstifac
+               
+               i0 = (srcioff + i) / dstifac + 1
+               j0 = (srcjoff + j) / dstjfac + 1
+               k0 = (srckoff + k) / dstkfac + 1
+               
                fi = mod(srcioff + i, dstifac)
-               ifac(1) = (fi-dstifac) * (-1)
-               ifac(2) = (fi        ) * 1
+               fj = mod(srcjoff + j, dstjfac)
+               fk = mod(srckoff + k, dstkfac)
                
                res = 0
                
-               do kk=1,2
-                  do jj=1,2
-                     do ii=1,2
+               do kk=0,1
+                  do jj=0,1
+                     do ii=0,1
                         
-                        fac = ifac(ii) * jfac(jj) * kfac(kk)
+                        fac = 1
+                        
+                        if (ii.eq.0) then
+                           fac = fac * (dstifac - fi)
+                        else
+                           fac = fac *            fi
+                        end if
+                        if (jj.eq.0) then
+                           fac = fac * (dstjfac - fj)
+                        else
+                           fac = fac *            fj
+                        end if
+                        if (kk.eq.0) then
+                           fac = fac * (dstkfac - fk)
+                        else
+                           fac = fac *            fk
+                        end if
                         
                         if (fac.ne.0) then
-                           CHKIDX (i0+ii, j0+jj, k0+kk, \
-                                   srciext,srcjext,srckext, "source")
-                           res = res + fac * src(i0+ii, j0+jj, k0+kk)
+                           res = res + fac * src(i0+ii,j0+jj,k0+kk)
                         end if
                         
                      end do
                   end do
                end do
                
-c$$$               fac = ifac(1) * jfac(1) * kfac(1)
-c$$$               if (fac.ne.0) res = res + fac * src(i0+1, j0+1, k0+1)
-c$$$               
-c$$$               fac = ifac(2) * jfac(1) * kfac(1)
-c$$$               if (fac.ne.0) res = res + fac * src(i0+2, j0+1, k0+1)
-c$$$               
-c$$$               fac = ifac(1) * jfac(2) * kfac(1)
-c$$$               if (fac.ne.0) res = res + fac * src(i0+1, j0+2, k0+1)
-c$$$               
-c$$$               fac = ifac(2) * jfac(2) * kfac(1)
-c$$$               if (fac.ne.0) res = res + fac * src(i0+2, j0+2, k0+1)
-c$$$               
-c$$$               fac = ifac(1) * jfac(1) * kfac(2)
-c$$$               if (fac.ne.0) res = res + fac * src(i0+1, j0+1, k0+2)
-c$$$               
-c$$$               fac = ifac(2) * jfac(1) * kfac(2)
-c$$$               if (fac.ne.0) res = res + fac * src(i0+2, j0+1, k0+2)
-c$$$               
-c$$$               fac = ifac(1) * jfac(2) * kfac(2)
-c$$$               if (fac.ne.0) res = res + fac * src(i0+1, j0+2, k0+2)
-c$$$               
-c$$$               fac = ifac(2) * jfac(2) * kfac(2)
-c$$$               if (fac.ne.0) res = res + fac * src(i0+2, j0+2, k0+2)
-               
-               CHKIDX (dstioff+i+1, dstjoff+j+1, dstkoff+k+1, \
-                       dstiext,dstjext,dstkext, "destination")
-               dst (dstioff+i+1, dstjoff+j+1, dstkoff+k+1) = dstdiv * res
+               dst (dstioff+i+1, dstjoff+j+1, dstkoff+k+1)
+     $              = res / (dstifac * dstjfac * dstkfac)
                
             end do
          end do
diff --git a/Carpet/CarpetLib/src/prolongate_3d_real8_2tl.F77 b/Carpet/CarpetLib/src/prolongate_3d_real8_2tl.F77
index 8fd69178f..f9b380c5d 100644
--- a/Carpet/CarpetLib/src/prolongate_3d_real8_2tl.F77
+++ b/Carpet/CarpetLib/src/prolongate_3d_real8_2tl.F77
@@ -1,20 +1,7 @@
 c     -*-Fortran-*-
-c     $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_2tl.F77,v 1.12 2004/03/11 12:03:09 schnetter Exp $
+c     $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_2tl.F77,v 1.2 2001/03/22 18:42:06 eschnett Exp $
       
 #include "cctk.h"
-
-
-
-#define CHKIDX(i,j,k, imax,jmax,kmax, where)                                  \
-      if ((i).lt.1 .or. (i).gt.(imax)                                         \
-           .or. (j).lt.1 .or. (j).gt.(jmax)                                   \
-           .or. (k).lt.1 .or. (k).gt.(kmax)) then                           &&\
-         write (msg, '(a, " array index out of bounds: shape is (",i4,",",i4,",",i4,"), index is (",i4,",",i4,",",i4,")")') \
-              (where), (imax), (jmax), (kmax), (i), (j), (k)                &&\
-         call CCTK_WARN (0, msg(1:len_trim(msg)))                           &&\
-      end if
-
-
       
       subroutine prolongate_3d_real8_2tl (
      $     src1, t1, src2, t2, srciext, srcjext, srckext,
@@ -26,17 +13,14 @@ c     $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d
       CCTK_REAL8 one
       parameter (one = 1)
       
-      CCTK_REAL8 eps
-      parameter (eps = 1.0d-10)
-      
       integer srciext, srcjext, srckext
       CCTK_REAL8 src1(srciext,srcjext,srckext)
-      CCTK_REAL8 t1
+      integer t1
       CCTK_REAL8 src2(srciext,srcjext,srckext)
-      CCTK_REAL8 t2
+      integer t2
       integer dstiext, dstjext, dstkext
       CCTK_REAL8 dst(dstiext,dstjext,dstkext)
-      CCTK_REAL8 t
+      integer t
 c     bbox(:,1) is lower boundary (inclusive)
 c     bbox(:,2) is upper boundary (inclusive)
 c     bbox(:,3) is stride
@@ -51,18 +35,14 @@ c     bbox(:,3) is stride
       
       CCTK_REAL8 s1fac, s2fac
       
-      CCTK_REAL8 dstdiv
       integer i, j, k
       integer i0, j0, k0
       integer fi, fj, fk
-      integer ifac(2), jfac(2), kfac(2)
       integer ii, jj, kk
       integer fac
       CCTK_REAL8 res
       integer d
       
-      character msg*1000
-      
       
       
       do d=1,3
@@ -78,25 +58,10 @@ c     bbox(:,3) is stride
             call CCTK_WARN (0, "Internal error: destination strides are not integer multiples of the source strides")
          end if
          if (mod(srcbbox(d,1), srcbbox(d,3)).ne.0
-     $        .or. mod(dstbbox(d,1), dstbbox(d,3)).ne.0
+     $     .or. mod(dstbbox(d,1), dstbbox(d,3)).ne.0
      $        .or. mod(regbbox(d,1), regbbox(d,3)).ne.0) then
             call CCTK_WARN (0, "Internal error: array origins are not integer multiples of the strides")
          end if
-         if (regbbox(d,1).gt.regbbox(d,2)) then
-c     This could be handled, but is likely to point to an error elsewhere
-            call CCTK_WARN (0, "Internal error: region extent is empty")
-         end if
-         if (mod(srcbbox(d,2) - srcbbox(d,1), srcbbox(d,3)).ne.0
-     $        .or. mod(dstbbox(d,2) - dstbbox(d,1), dstbbox(d,3)).ne.0
-     $        .or. mod(regbbox(d,2) - regbbox(d,1), regbbox(d,3)).ne.0) then
-            call CCTK_WARN (0, "Internal error: array extents are not integer multiples of the strides")
-         end if
-         if (regbbox(d,1).lt.srcbbox(d,1)
-     $        .or. regbbox(d,1).lt.dstbbox(d,1)
-     $        .or. regbbox(d,2).gt.srcbbox(d,2)
-     $        .or. regbbox(d,2).gt.dstbbox(d,2)) then
-            call CCTK_WARN (0, "Internal error: region extent is not contained in array extent")
-         end if
       end do
       
       if (srciext.ne.(srcbbox(1,2)-srcbbox(1,1))/srcbbox(1,3)+1
@@ -128,63 +93,60 @@ c     This could be handled, but is likely to point to an error elsewhere
       
       
       
-c     Linear (first order) interpolation
-      if (t1.eq.t2) then
-         call CCTK_WARN (0, "Internal error: arrays have same time")
-      end if
-      if (t.lt.min(t1,t2)-eps .or. t.gt.max(t1,t2)+eps) then
-         call CCTK_WARN (0, "Internal error: extrapolation in time")
-      end if
-      
-      s1fac = (t - t2) / (t1 - t2)
-      s2fac = (t - t1) / (t2 - t1)
+      s1fac = (t - t2) * one / (t1 - t2)
+      s2fac = (t - t1) * one / (t2 - t1)
       
       
       
 c     Loop over fine region
-      dstdiv = one / (dstifac * dstjfac * dstkfac)
-      
       do k = 0, regkext-1
-         k0 = (srckoff + k) / dstkfac
-         fk = mod(srckoff + k, dstkfac)
-         kfac(1) = (fk-dstkfac) * (-1)
-         kfac(2) = (fk        ) * 1
-         
          do j = 0, regjext-1
-            j0 = (srcjoff + j) / dstjfac
-            fj = mod(srcjoff + j, dstjfac)
-            jfac(1) = (fj-dstjfac) * (-1)
-            jfac(2) = (fj        ) * 1
-            
             do i = 0, regiext-1
-               i0 = (srcioff + i) / dstifac
+               
+               i0 = (srcioff + i) / dstifac + 1
+               j0 = (srcjoff + j) / dstjfac + 1
+               k0 = (srckoff + k) / dstkfac + 1
+               
                fi = mod(srcioff + i, dstifac)
-               ifac(1) = (fi-dstifac) * (-1)
-               ifac(2) = (fi        ) * 1
+               fj = mod(srcjoff + j, dstjfac)
+               fk = mod(srckoff + k, dstkfac)
                
                res = 0
                
-               do kk=1,2
-                  do jj=1,2
-                     do ii=1,2
+               do kk=0,1
+                  do jj=0,1
+                     do ii=0,1
                         
-                        fac = ifac(ii) * jfac(jj) * kfac(kk)
+                        fac = 1
+                        
+                        if (ii.eq.0) then
+                           fac = fac * (dstifac - fi)
+                        else
+                           fac = fac *            fi
+                        end if
+                        if (jj.eq.0) then
+                           fac = fac * (dstjfac - fj)
+                        else
+                           fac = fac *            fj
+                        end if
+                        if (kk.eq.0) then
+                           fac = fac * (dstkfac - fk)
+                        else
+                           fac = fac *            fk
+                        end if
                         
                         if (fac.ne.0) then
-                           CHKIDX (i0+ii, j0+jj, k0+kk, \
-                                   srciext,srcjext,srckext, "source")
                            res = res
-     $                          + fac * s1fac * src1(i0+ii, j0+jj, k0+kk)
-     $                          + fac * s2fac * src2(i0+ii, j0+jj, k0+kk)
+     $                          + fac * s1fac * src1(i0+ii,j0+jj,k0+kk)
+     $                          + fac * s2fac * src2(i0+ii,j0+jj,k0+kk)
                         end if
                         
                      end do
                   end do
                end do
                
-               CHKIDX (dstioff+i+1, dstjoff+j+1, dstkoff+k+1, \
-                       dstiext,dstjext,dstkext, "destination")
-               dst (dstioff+i+1, dstjoff+j+1, dstkoff+k+1) = dstdiv * res
+               dst (dstioff+i+1, dstjoff+j+1, dstkoff+k+1)
+     $              = res / (dstifac * dstjfac * dstkfac)
                
             end do
          end do
diff --git a/Carpet/CarpetLib/src/restrict_3d_real8.F77 b/Carpet/CarpetLib/src/restrict_3d_real8.F77
index 81f4cfd0a..68ea98b11 100644
--- a/Carpet/CarpetLib/src/restrict_3d_real8.F77
+++ b/Carpet/CarpetLib/src/restrict_3d_real8.F77
@@ -1,20 +1,7 @@
 c     -*-Fortran-*-
-c     $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/restrict_3d_real8.F77,v 1.7 2004/03/11 12:03:09 schnetter Exp $
+c     $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/restrict_3d_real8.F77,v 1.2 2001/03/22 18:42:06 eschnett Exp $
       
 #include "cctk.h"
-
-
-
-#define CHKIDX(i,j,k, imax,jmax,kmax, where)                                  \
-      if ((i).lt.1 .or. (i).gt.(imax)                                         \
-           .or. (j).lt.1 .or. (j).gt.(jmax)                                   \
-           .or. (k).lt.1 .or. (k).gt.(kmax)) then                           &&\
-         write (msg, '(a, " array index out of bounds: shape is (",i4,",",i4,",",i4,"), index is (",i4,",",i4,",",i4,")")') \
-              (where), (imax), (jmax), (kmax), (i), (j), (k)                &&\
-         call CCTK_WARN (0, msg(1:len_trim(msg)))                           &&\
-      end if
-
-
       
       subroutine restrict_3d_real8 (
      $     src, srciext, srcjext, srckext,
@@ -42,8 +29,6 @@ c     bbox(:,3) is stride
       integer i, j, k
       integer d
       
-      character msg*1000
-      
       
       
       do d=1,3
@@ -59,25 +44,10 @@ c     bbox(:,3) is stride
             call CCTK_WARN (0, "Internal error: source strides are not integer multiples of the destination strides")
          end if
          if (mod(srcbbox(d,1), srcbbox(d,3)).ne.0
-     $        .or. mod(dstbbox(d,1), dstbbox(d,3)).ne.0
+     $     .or. mod(dstbbox(d,1), dstbbox(d,3)).ne.0
      $        .or. mod(regbbox(d,1), regbbox(d,3)).ne.0) then
             call CCTK_WARN (0, "Internal error: array origins are not integer multiples of the strides")
          end if
-         if (regbbox(d,1).gt.regbbox(d,2)) then
-c     This could be handled, but is likely to point to an error elsewhere
-            call CCTK_WARN (0, "Internal error: region extent is empty")
-         end if
-         if (mod(srcbbox(d,2) - srcbbox(d,1), srcbbox(d,3)).ne.0
-     $        .or. mod(dstbbox(d,2) - dstbbox(d,1), dstbbox(d,3)).ne.0
-     $        .or. mod(regbbox(d,2) - regbbox(d,1), regbbox(d,3)).ne.0) then
-            call CCTK_WARN (0, "Internal error: array extents are not integer multiples of the strides")
-         end if
-         if (regbbox(d,1).lt.srcbbox(d,1)
-     $        .or. regbbox(d,1).lt.dstbbox(d,1)
-     $        .or. regbbox(d,2).gt.srcbbox(d,2)
-     $        .or. regbbox(d,2).gt.dstbbox(d,2)) then
-            call CCTK_WARN (0, "Internal error: region extent is not contained in array extent")
-         end if
       end do
       
       if (srciext.ne.(srcbbox(1,2)-srcbbox(1,1))/srcbbox(1,3)+1
@@ -114,10 +84,6 @@ c     Loop over coarse region
          do j = 0, regjext-1
             do i = 0, regiext-1
                
-               CHKIDX (srcioff+srcifac*i+1, srcjoff+srcjfac*j+1, srckoff+srckfac*k+1, \
-                       srciext,srcjext,srckext, "source")
-               CHKIDX (dstioff+i+1, dstjoff+j+1, dstkoff+k+1, \
-                       dstiext,dstjext,dstkext, "destination")
                dst (dstioff+i+1, dstjoff+j+1, dstkoff+k+1)
      $              = src (srcioff+srcifac*i+1, srcjoff+srcjfac*j+1, srckoff+srckfac*k+1)
                
diff --git a/Carpet/CarpetLib/src/th.cc b/Carpet/CarpetLib/src/th.cc
index aab64bc78..17aaff924 100644
--- a/Carpet/CarpetLib/src/th.cc
+++ b/Carpet/CarpetLib/src/th.cc
@@ -6,7 +6,7 @@
     copyright            : (C) 2000 by Erik Schnetter
     email                : schnetter@astro.psu.edu
 
-    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/th.cc,v 1.2 2001/03/19 21:30:19 eschnett Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/th.cc,v 1.3 2001/03/22 18:42:06 eschnett Exp $
 
  ***************************************************************************/
 
@@ -19,7 +19,8 @@
  *                                                                         *
  ***************************************************************************/
 
-#include <cassert>
+#include <assert.h>
+
 #include <iostream>
 
 #include "defs.hh"
@@ -29,6 +30,8 @@
 #  include "th.hh"
 #endif
 
+using namespace std;
+
 
 
 // Constructors
diff --git a/Carpet/CarpetLib/src/th.hh b/Carpet/CarpetLib/src/th.hh
index 6089e37cd..7c39a7d18 100644
--- a/Carpet/CarpetLib/src/th.hh
+++ b/Carpet/CarpetLib/src/th.hh
@@ -6,7 +6,7 @@
     copyright            : (C) 2000 by Erik Schnetter
     email                : schnetter@astro.psu.edu
 
-    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/th.hh,v 1.1 2001/03/01 13:40:10 eschnett Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/th.hh,v 1.2 2001/03/22 18:42:06 eschnett Exp $
 
  ***************************************************************************/
 
@@ -22,13 +22,16 @@
 #ifndef TH_HH
 #define TH_HH
 
-#include <cassert>
+#include <assert.h>
+
 #include <iostream>
 #include <vector>
 
 #include "defs.hh"
 #include "gh.hh"
 
+using namespace std;
+
 
 
 // Forward declaration
diff --git a/Carpet/CarpetLib/src/vect.cc b/Carpet/CarpetLib/src/vect.cc
index 6c8d012a2..5a622fd1c 100644
--- a/Carpet/CarpetLib/src/vect.cc
+++ b/Carpet/CarpetLib/src/vect.cc
@@ -5,7 +5,7 @@
     copyright            : (C) 2000 by Erik Schnetter
     email                : schnetter@astro.psu.edu
 
-    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/vect.cc,v 1.2 2001/03/19 21:30:19 eschnett Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/vect.cc,v 1.3 2001/03/22 18:42:06 eschnett Exp $
 
  ***************************************************************************/
 
@@ -18,7 +18,8 @@
  *                                                                         *
  ***************************************************************************/
 
-#include <cassert>
+#include <assert.h>
+
 #include <iostream>
 
 #include "defs.hh"
@@ -27,9 +28,13 @@
 #  include "vect.hh"
 #endif
 
+using namespace std;
+
 
 
 // Output
+#ifndef SGI
+// This doesn't work on SGIs.  Is this legal C++?
 template<class T,int D>
 ostream& operator<< (ostream& os, const vect<T,D>& a) {
   os << "[";
@@ -40,6 +45,35 @@ ostream& operator<< (ostream& os, const vect<T,D>& a) {
   os << "]";
   return os;
 }
+#else
+ostream& operator<< (ostream& os, const vect<int,1>& a) {
+  os << "[";
+  for (int d=0; d<1; ++d) {
+    if (d>0) os << ",";
+    os << a[d];
+  }
+  os << "]";
+  return os;
+}
+ostream& operator<< (ostream& os, const vect<int,2>& a) {
+  os << "[";
+  for (int d=0; d<2; ++d) {
+    if (d>0) os << ",";
+    os << a[d];
+  }
+  os << "]";
+  return os;
+}
+ostream& operator<< (ostream& os, const vect<int,3>& a) {
+  os << "[";
+  for (int d=0; d<3; ++d) {
+    if (d>0) os << ",";
+    os << a[d];
+  }
+  os << "]";
+  return os;
+}
+#endif
 
 
 
diff --git a/Carpet/CarpetLib/src/vect.hh b/Carpet/CarpetLib/src/vect.hh
index b4b2a2bba..d5ad5cbab 100644
--- a/Carpet/CarpetLib/src/vect.hh
+++ b/Carpet/CarpetLib/src/vect.hh
@@ -5,7 +5,7 @@
     copyright            : (C) 2000 by Erik Schnetter
     email                : schnetter@astro.psu.edu
 
-    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/vect.hh,v 1.2 2001/03/17 22:37:56 eschnett Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/vect.hh,v 1.3 2001/03/22 18:42:06 eschnett Exp $
 
  ***************************************************************************/
 
@@ -21,10 +21,13 @@
 #ifndef VECT_HH
 #define VECT_HH
 
-#include <cassert>
-#include <cmath>
+#include <assert.h>
+#include <math.h>
+
 #include <iostream>
 
+using namespace std;
+
 
 
 // Forward definition
@@ -201,25 +204,25 @@ public:
   
   // Non-modifying operators
   vect operator+ () const {
-    vect r(*this);
+    vect r;
     for (int d=0; d<D; ++d) r[d]=+r[d];
     return r;
   }
   
   vect operator- () const {
-    vect r(*this);
+    vect r;
     for (int d=0; d<D; ++d) r[d]=-r[d];
     return r;
   }
   
   vect<bool,D> operator! () const {
-    vect r(*this);
+    vect<bool,D> r;
     for (int d=0; d<D; ++d) r[d]=!r[d];
     return r;
   }
   
   vect operator~ () const {
-    vect r(*this);
+    vect r;
     for (int d=0; d<D; ++d) r[d]=~r[d];
     return r;
   }
@@ -273,14 +276,14 @@ public:
   }
   
   vect<bool,D> operator&& (const T x) const {
-    vect r(*this);
-    for (int d=0; d<D; ++d) r[d]=r[d]&&x;
+    vect<bool,D> r;
+    for (int d=0; d<D; ++d) r[d]=(*this)[d]&&x;
     return r;
   }
   
   vect<bool,D> operator|| (const T x) const {
-    vect r(*this);
-    for (int d=0; d<D; ++d) r[d]=r[d]||x;
+    vect<bool,D> r;
+    for (int d=0; d<D; ++d) r[d]=(*this)[d]||x;
     return r;
   }
   
@@ -402,8 +405,6 @@ public:
   };
 #endif
   
-  // Output
-  friend ostream& operator<< <>(ostream& os, const vect& a);
 };
 
 
diff --git a/Carpet/CarpetSlab/src/carpetslab.cc b/Carpet/CarpetSlab/src/carpetslab.cc
index 30948ea12..4d8141e0d 100644
--- a/Carpet/CarpetSlab/src/carpetslab.cc
+++ b/Carpet/CarpetSlab/src/carpetslab.cc
@@ -1,7 +1,8 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetSlab/src/Attic/carpetslab.cc,v 1.5 2001/03/17 16:14:52 eschnett Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetSlab/src/Attic/carpetslab.cc,v 1.6 2001/03/22 18:42:06 eschnett Exp $
 
-#include <cassert>
-#include <cstdlib>
+#include <alloca.h>
+#include <assert.h>
+#include <stdlib.h>
 
 #include "cctk.h"
 
@@ -17,7 +18,7 @@
 
 #include "carpetslab.hh"
 
-static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetSlab/src/Attic/carpetslab.cc,v 1.5 2001/03/17 16:14:52 eschnett Exp $";
+static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetSlab/src/Attic/carpetslab.cc,v 1.6 2001/03/22 18:42:06 eschnett Exp $";
 
 
 
@@ -267,7 +268,7 @@ namespace CarpetSlab {
     assert (hsize);
     
     // Calculate more convenient representation of the direction
-    int dirs[hdim];
+    int* const dirs = (int*)alloca(hdim * sizeof(int));
     // The following if statement is written according to the
     // definition of "dir".
     if (hdim==1) {
diff --git a/CarpetAttic/CarpetIOFlexIO/schedule.ccl b/CarpetAttic/CarpetIOFlexIO/schedule.ccl
index ad6119623..259abe7a0 100644
--- a/CarpetAttic/CarpetIOFlexIO/schedule.ccl
+++ b/CarpetAttic/CarpetIOFlexIO/schedule.ccl
@@ -1,5 +1,5 @@
 # Schedule definitions for thorn CarpetIOFlexIO
-# $Header: /home/eschnett/C/carpet/Carpet/CarpetAttic/CarpetIOFlexIO/schedule.ccl,v 1.3 2001/03/18 05:20:24 eschnett Exp $
+# $Header: /home/eschnett/C/carpet/Carpet/CarpetAttic/CarpetIOFlexIO/schedule.ccl,v 1.4 2001/03/22 18:42:05 eschnett Exp $
 
 schedule CarpetIOFlexIOStartup at STARTUP after IOUtil_Startup
 {
diff --git a/CarpetAttic/CarpetIOFlexIO/src/ioflexio.cc b/CarpetAttic/CarpetIOFlexIO/src/ioflexio.cc
index 2628e236e..e1a6282ca 100644
--- a/CarpetAttic/CarpetIOFlexIO/src/ioflexio.cc
+++ b/CarpetAttic/CarpetIOFlexIO/src/ioflexio.cc
@@ -1,11 +1,13 @@
-#include <cassert>
-#include <climits>
-#include <cstdio>
-#include <cstdlib>
-#include <cstring>
-#include <fstream>
+#include <alloca.h>
+#include <assert.h>
+#include <limits.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
 #include <sys/stat.h>
 #include <sys/types.h>
+
+#include <fstream>
 #include <vector>
 
 #include <AMRwriter.hh>
@@ -31,12 +33,13 @@
 
 #include "ioflexio.hh"
 
-static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/CarpetAttic/CarpetIOFlexIO/src/ioflexio.cc,v 1.6 2001/03/19 21:30:09 eschnett Exp $";
+static const char* rcsid = "$Header: /home/eschnett/C/carpet/Carpet/CarpetAttic/CarpetIOFlexIO/src/ioflexio.cc,v 1.7 2001/03/22 18:42:05 eschnett Exp $";
 
 
 
 namespace CarpetIOFlexIO {
   
+  using namespace std;
   using namespace Carpet;
   
   
@@ -165,7 +168,7 @@ namespace CarpetIOFlexIO {
       }
       extension = GetStringParameter ("out3D_extension", extension);
       
-      char filename[strlen(myoutdir)+strlen(alias)+strlen(extension)+100];
+      char* const filename = (char*)alloca(strlen(myoutdir)+strlen(alias)+strlen(extension)+100);
       sprintf (filename, "%s/%s%s", myoutdir, alias, extension);
       
       IObase* writer = 0;
@@ -446,7 +449,7 @@ namespace CarpetIOFlexIO {
       }
       extension = GetStringParameter ("in3D_extension", extension);
       
-      char filename[strlen(myindir)+strlen(alias)+strlen(extension)+100];
+      char* const filename = (char*)alloca(strlen(myindir)+strlen(alias)+strlen(extension)+100);
       sprintf (filename, "%s/%s.%s", myindir, alias, extension);
       
       IObase* reader = 0;
@@ -652,7 +655,7 @@ namespace CarpetIOFlexIO {
     const int numvars = CCTK_NumVars();
     assert (vindex>=0 && vindex<numvars);
     
-    bool flags[numvars];
+    bool* const flags = (bool*)alloca(numvars * sizeof(bool));
     
     for (int i=0; i<numvars; ++i) {
       flags[i] = false;
diff --git a/CarpetExtra/SpaceTimeToy/src/SpaceTimeToy.F77 b/CarpetExtra/SpaceTimeToy/src/SpaceTimeToy.F77
index 387059dc2..b9fb12aaf 100644
--- a/CarpetExtra/SpaceTimeToy/src/SpaceTimeToy.F77
+++ b/CarpetExtra/SpaceTimeToy/src/SpaceTimeToy.F77
@@ -1,5 +1,5 @@
 c     -*-Fortran-*-
-c     $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/SpaceTimeToy/src/SpaceTimeToy.F77,v 1.5 2001/03/21 22:57:40 eschnett Exp $
+c     $Header: /home/eschnett/C/carpet/Carpet/CarpetExtra/SpaceTimeToy/src/SpaceTimeToy.F77,v 1.6 2001/03/22 18:42:06 eschnett Exp $
       
 #include "cctk.h" 
 #include "cctk_Parameters.h"
@@ -182,7 +182,7 @@ c     Apply boundary condition
       
       call Cart3dSymGN (ierr, cctkGH, "spacetimetoy::spacetimeevolve")
       if (ierr .lt. 0)  then
-         call CCTK_WARN (0, "Error while applying boundary condition")
+         call CCTK_WARN (0, "Error while applying symmetry condition")
       end if
       
       end