Changed handling of interpolation orders; they are now stored in the

Changed handling of interpolation orders; they are now stored in the
grid functions and don't have to be passed around.

Speeded up some prolongation operators.  They are still slow.

darcs-hash:20011209154308-07bb3-caa74a89a87c290852f2b59160ed26d9361f3cf1.gz

23 files changed, 252 insertions, 223 deletions

diff --git a/Carpet/CarpetLib/src/copy_3d_real8.F77 b/Carpet/CarpetLib/src/copy_3d_real8.F77
index e00dd8d54..8a07051dc 100644
--- a/Carpet/CarpetLib/src/copy_3d_real8.F77
+++ b/Carpet/CarpetLib/src/copy_3d_real8.F77
@@ -1,5 +1,5 @@
 c     -*-Fortran-*-
-c     $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/copy_3d_real8.F77,v 1.3 2001/03/24 22:38:48 eschnett Exp $
+c     $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/copy_3d_real8.F77,v 1.4 2001/12/09 16:43:08 schnetter Exp $
       
 #include "cctk.h"
       
@@ -39,7 +39,7 @@ 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
@@ -47,7 +47,7 @@ c     bbox(:,3) is stride
      $        .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 extents")
+            call CCTK_WARN (0, "Internal error: region extent is not contained in array extent")
          end if
       end do
       
diff --git a/Carpet/CarpetLib/src/data.cc b/Carpet/CarpetLib/src/data.cc
index d21eb0bdd..efca6c699 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.13 2001/07/04 12:29:51 schnetter Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/data.cc,v 1.14 2001/12/09 16:43:09 schnetter Exp $
 
  ***************************************************************************/
 
@@ -202,7 +202,8 @@ void data<T,D>
 ::interpolate_from_innerloop (const vector<const generic_data<D>*> gsrcs,
 			      const vector<int> tls,
 			      const ibbox& box, const int tl,
-			      const int order_space)
+			      const int order_space,
+			      const int order_time)
 {
   assert (has_storage());
   assert (all(box.lower()>=extent().lower()));
@@ -223,14 +224,15 @@ void data<T,D>
   MPI_Comm_rank (dist::comm, &rank);
   assert (rank == proc());
   
-  assert (order_space == 1);
+  assert (order_space <= 1);
   
-  vector<T> src_fac(srcs.size());
+  assert (srcs.size() > order_time);
+  vector<T> src_fac(order_time+1);
   for (int t=0; t<(int)src_fac.size(); ++t) {
     src_fac[t] = 1;
     for (int tt=0; tt<(int)src_fac.size(); ++tt) {
       if (tt!=t) {
-	src_fac[t] *= (T)(t - tls[tt]) / (T)(tls[t] - tls[tt]);
+	src_fac[t] *= (T)(tl - tls[tt]) / (T)(tls[t] - tls[tt]);
       }
     }
   }
@@ -413,7 +415,8 @@ void data<CCTK_REAL8,3>
 ::interpolate_from_innerloop (const vector<const generic_data<3>*> gsrcs,
 			      const vector<int> tls,
 			      const ibbox& box, const int tl,
-			      const int order_space)
+			      const int order_space,
+			      const int order_time)
 {
   assert (has_storage());
   assert (all(box.lower()>=extent().lower()));
@@ -471,11 +474,12 @@ void data<CCTK_REAL8,3>
     
   } else if (all(sext.stride() > dext.stride())) {
     
-    switch (srcs.size()) {
+    switch (order_time) {
       
-    case 1:
-      // One timelevel
+    case 0:
+      assert (srcs.size()>=1);
       switch (order_space) {
+      case 0:
       case 1:
 	CCTK_FNAME(prolongate_3d_real8)
 	  ((const CCTK_REAL8*)srcs[0]->storage(),
@@ -484,6 +488,7 @@ void data<CCTK_REAL8,3>
 	   dstshp[0], dstshp[1], dstshp[2],
 	   srcbbox, dstbbox, regbbox);
 	break;
+      case 2:
       case 3:
 	CCTK_FNAME(prolongate_3d_real8_o3)
 	  ((const CCTK_REAL8*)srcs[0]->storage(),
@@ -497,9 +502,10 @@ void data<CCTK_REAL8,3>
       }
       break;
       
-    case 2:
-      // Two timelevels
+    case 1:
+      assert (srcs.size()>=2);
       switch (order_space) {
+      case 0:
       case 1:
 	CCTK_FNAME(prolongate_3d_real8_2tl)
 	  ((const CCTK_REAL8*)srcs[0]->storage(), tls[0],
@@ -509,6 +515,7 @@ void data<CCTK_REAL8,3>
 	   dstshp[0], dstshp[1], dstshp[2],
 	   srcbbox, dstbbox, regbbox);
 	break;
+      case 2:
       case 3:
 	CCTK_FNAME(prolongate_3d_real8_2tl_o3)
 	  ((const CCTK_REAL8*)srcs[0]->storage(), tls[0],
@@ -523,9 +530,10 @@ void data<CCTK_REAL8,3>
       }
       break;
       
-    case 3:
-      // Three timelevels
+    case 2:
+      assert (srcs.size()>=3);
       switch (order_space) {
+      case 0:
       case 1:
 	CCTK_FNAME(prolongate_3d_real8_3tl)
 	  ((const CCTK_REAL8*)srcs[0]->storage(), tls[0],
@@ -536,6 +544,7 @@ void data<CCTK_REAL8,3>
 	   dstshp[0], dstshp[1], dstshp[2],
 	   srcbbox, dstbbox, regbbox);
 	break;
+      case 2:
       case 3:
 	CCTK_FNAME(prolongate_3d_real8_3tl_o3)
 	  ((const CCTK_REAL8*)srcs[0]->storage(), tls[0],
diff --git a/Carpet/CarpetLib/src/data.hh b/Carpet/CarpetLib/src/data.hh
index 8d49998a7..00e239e22 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.7 2001/03/30 00:50:21 eschnett Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/data.hh,v 1.8 2001/12/09 16:43:09 schnetter Exp $
 
  ***************************************************************************/
 
@@ -96,7 +96,8 @@ public:
   void interpolate_from_innerloop (const vector<const generic_data<D>*> gsrcs,
 				   const vector<int> tls,
 				   const ibbox& box, const int tl,
-				   const int order_space);
+				   const int order_space,
+				   const int order_time);
   
   void write_ascii_output_element (ostream& os, const ivect& index) const;
 //   void write_ieee (const string name, const int time,
diff --git a/Carpet/CarpetLib/src/dgdh.cc b/Carpet/CarpetLib/src/dgdh.cc
index 188ca8885..4287621c9 100644
--- a/Carpet/CarpetLib/src/dgdh.cc
+++ b/Carpet/CarpetLib/src/dgdh.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/Attic/dgdh.cc,v 1.1 2001/06/12 14:56:57 schnetter Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/Attic/dgdh.cc,v 1.2 2001/12/09 16:43:09 schnetter Exp $
 
  ***************************************************************************/
 
@@ -36,10 +36,10 @@ using namespace std;
 
 
 // Constructors
-dimgeneric_dh::dimgeneric_dh (int prolongation_order)
-  : prolongation_order(prolongation_order)
+dimgeneric_dh::dimgeneric_dh (int prolongation_order_space)
+  : prolongation_order_space(prolongation_order_space)
 {
-  assert (prolongation_order>=0);
+  assert (prolongation_order_space>=0);
 }
 
 // Destructors
@@ -49,8 +49,8 @@ dimgeneric_dh::~dimgeneric_dh ()
 
 // Helpers
 int dimgeneric_dh::prolongation_stencil_size () const {
-  assert (prolongation_order>=0);
-  return prolongation_order/2;
+  assert (prolongation_order_space>=0);
+  return prolongation_order_space/2;
 }
 
 
diff --git a/Carpet/CarpetLib/src/dgdh.hh b/Carpet/CarpetLib/src/dgdh.hh
index ff7ed8b1e..222e2cfab 100644
--- a/Carpet/CarpetLib/src/dgdh.hh
+++ b/Carpet/CarpetLib/src/dgdh.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/Attic/dgdh.hh,v 1.1 2001/06/12 14:56:58 schnetter Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/Attic/dgdh.hh,v 1.2 2001/12/09 16:43:09 schnetter Exp $
 
  ***************************************************************************/
 
@@ -48,12 +48,12 @@ class dimgeneric_dh {
 public:				// should be readonly
   
   // Fields
-  int prolongation_order;	// order of spatial prolongation operator
+  int prolongation_order_space;	// order of spatial prolongation operator
   
 public:
   
   // Constructors
-  dimgeneric_dh (int prolongation_order);
+  dimgeneric_dh (int prolongation_order_space);
   
   // Destructors
   virtual ~dimgeneric_dh ();
diff --git a/Carpet/CarpetLib/src/dggf.cc b/Carpet/CarpetLib/src/dggf.cc
index 43c72a1df..998c5b90e 100644
--- a/Carpet/CarpetLib/src/dggf.cc
+++ b/Carpet/CarpetLib/src/dggf.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/Attic/dggf.cc,v 1.1 2001/06/12 14:56:58 schnetter Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/Attic/dggf.cc,v 1.2 2001/12/09 16:43:09 schnetter Exp $
 
  ***************************************************************************/
 
@@ -39,8 +39,10 @@ using namespace std;
 
 // Constructors
 dimgeneric_gf::dimgeneric_gf (const string name, th& t,
-			      const int tmin, const int tmax)
-  : name(name), t(t), tmin(tmin), tmax(tmax)
+			      const int tmin, const int tmax,
+			      const int prolongation_order_time)
+  : name(name), t(t),
+    tmin(tmin), tmax(tmax), prolongation_order_time(prolongation_order_time)
 {
   assert (tmin<=tmax+1);
 }
diff --git a/Carpet/CarpetLib/src/dggf.hh b/Carpet/CarpetLib/src/dggf.hh
index 76727afee..28004463e 100644
--- a/Carpet/CarpetLib/src/dggf.hh
+++ b/Carpet/CarpetLib/src/dggf.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/Attic/dggf.hh,v 1.2 2001/07/02 13:22:13 schnetter Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/Attic/dggf.hh,v 1.3 2001/12/09 16:43:09 schnetter Exp $
 
  ***************************************************************************/
 
@@ -53,12 +53,14 @@ public:				// should be readonly
 
   th &t;			// time hierarchy
   int tmin, tmax;		// timelevels
+  int prolongation_order_time;	// order of temporal prolongation operator
 
 public:
 
   // Constructors
   dimgeneric_gf (const string name, th& t, 
-		 const int tmin, const int tmax);
+		 const int tmin, const int tmax,
+		 const int prolongation_order_time);
   
   // Destructors
   virtual ~dimgeneric_gf ();
@@ -99,24 +101,19 @@ public:
   virtual void sync (int tl, int rl, int c, int ml) = 0;
 
   // Prolongate the boundaries of a component
-  virtual void ref_bnd_prolongate (int tl, int rl, int c, int ml,
-				   int order_space=1, int order_time=1) = 0;
+  virtual void ref_bnd_prolongate (int tl, int rl, int c, int ml) = 0;
   
   // Restrict a multigrid level
-  virtual void mg_restrict (int tl, int rl, int c, int ml, int order_space=1)
-    = 0;
+  virtual void mg_restrict (int tl, int rl, int c, int ml) = 0;
 
   // Prolongate a multigrid level
-  virtual void mg_prolongate (int tl, int rl, int c, int ml, int order_space=1)
-    = 0;
+  virtual void mg_prolongate (int tl, int rl, int c, int ml) = 0;
 
   // Restrict a refinement level
-  virtual void ref_restrict (int tl, int rl, int c, int ml, int order_space=1)
-    = 0;
+  virtual void ref_restrict (int tl, int rl, int c, int ml) = 0;
 
   // Prolongate a refinement level
-  virtual void ref_prolongate (int tl, int rl, int c, int ml,
-			       int order_space=1) = 0;
+  virtual void ref_prolongate (int tl, int rl, int c, int ml) = 0;
   
   
   
diff --git a/Carpet/CarpetLib/src/dh.cc b/Carpet/CarpetLib/src/dh.cc
index f1512be42..2f3863dac 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.16 2001/07/04 12:29:51 schnetter Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/dh.cc,v 1.17 2001/12/09 16:43:09 schnetter Exp $
 
  ***************************************************************************/
 
@@ -39,8 +39,8 @@ using namespace std;
 // Constructors
 template<int D>
 dh<D>::dh (gh<D>& h, const ivect& lghosts, const ivect& ughosts,
-	   int prolongation_order)
-  : dimgeneric_dh(prolongation_order),
+	   int prolongation_order_space)
+  : dimgeneric_dh(prolongation_order_space),
     h(h),
     lghosts(lghosts), ughosts(ughosts)
 {
diff --git a/Carpet/CarpetLib/src/dh.hh b/Carpet/CarpetLib/src/dh.hh
index 7ba680691..73c2a54d6 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.8 2001/06/12 14:56:59 schnetter Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/dh.hh,v 1.9 2001/12/09 16:43:09 schnetter Exp $
 
  ***************************************************************************/
 
@@ -109,7 +109,7 @@ public:
   
   // Constructors
   dh (gh<D>& h, const ivect& lghosts, const ivect& ughosts,
-      int prolongation_order);
+      int prolongation_order_space);
   
   // Destructors
   virtual ~dh ();
diff --git a/Carpet/CarpetLib/src/gdata.cc b/Carpet/CarpetLib/src/gdata.cc
index 347bfb52c..41e265238 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.14 2001/07/04 12:29:52 schnetter Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gdata.cc,v 1.15 2001/12/09 16:43:10 schnetter Exp $
 
  ***************************************************************************/
 
@@ -90,7 +90,8 @@ void generic_data<D>
 ::interpolate_from (const vector<const generic_data*> srcs,
 		    const vector<int> tls,
 		    const ibbox& box, const int tl,
-		    const int order_space)
+		    const int order_space,
+		    const int order_time)
 {
   assert (has_storage());
   assert (all(box.lower()>=extent().lower()));
@@ -110,7 +111,7 @@ void generic_data<D>
     int rank;
     MPI_Comm_rank (dist::comm, &rank);
     if (rank == proc()) {
-      interpolate_from_innerloop (srcs, tls, box, tl, order_space);
+      interpolate_from_innerloop (srcs, tls, box, tl, order_space, order_time);
     }
     
   } else {
@@ -118,7 +119,7 @@ void generic_data<D>
     
     generic_data* const tmp = make_typed();
     tmp->allocate (box, srcs[0]->proc());
-    tmp->interpolate_from (srcs, tls, box, tl, order_space);
+    tmp->interpolate_from (srcs, tls, box, tl, order_space, order_time);
     tmp->change_processor (proc());
     copy_from (tmp, box);
     delete tmp;
diff --git a/Carpet/CarpetLib/src/gdata.hh b/Carpet/CarpetLib/src/gdata.hh
index 8e561c073..e827737f8 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.9 2001/06/12 14:56:59 schnetter Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gdata.hh,v 1.10 2001/12/09 16:43:10 schnetter Exp $
 
  ***************************************************************************/
 
@@ -102,7 +102,8 @@ public:
   void interpolate_from (const vector<const generic_data*> srcs,
 			 const vector<int> tls,
 			 const ibbox& box, const int tl,
-			 const int order_space);
+			 const int order_space,
+			 const int order_time);
 protected:
   virtual void
   copy_from_innerloop (const generic_data* src, const ibbox& box) = 0;
@@ -110,7 +111,8 @@ protected:
   interpolate_from_innerloop (const vector<const generic_data*> srcs,
 			      const vector<int> tls,
 			      const ibbox& box, const int tl,
-			      const int order_space) = 0;
+			      const int order_space,
+			      const int order_time) = 0;
   
 public:
   
diff --git a/Carpet/CarpetLib/src/gf.cc b/Carpet/CarpetLib/src/gf.cc
index 30eb4a5b4..4a026a574 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.7 2001/07/04 12:29:52 schnetter Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gf.cc,v 1.8 2001/12/09 16:43:10 schnetter Exp $
 
  ***************************************************************************/
 
@@ -34,8 +34,8 @@ using namespace std;
 // Constructors
 template<class T,int D>
 gf<T,D>::gf (const string name, th& t, dh<D>& d,
-	     const int tmin, const int tmax)
-  : generic_gf<D>(name, t, d, tmin, tmax)
+	     const int tmin, const int tmax, const int prolongation_order_time)
+  : generic_gf<D>(name, t, d, tmin, tmax, prolongation_order_time)
 {
   recompose();
 }
diff --git a/Carpet/CarpetLib/src/gf.hh b/Carpet/CarpetLib/src/gf.hh
index 4a51bc98d..df500e074 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.4 2001/06/12 14:56:59 schnetter Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gf.hh,v 1.5 2001/12/09 16:43:10 schnetter Exp $
 
  ***************************************************************************/
 
@@ -61,7 +61,8 @@ class gf: public generic_gf<D> {
 public:
   
   // Constructors
-  gf (const string name, th& t, dh<D>& d, const int tmin, const int tmax);
+  gf (const string name, th& t, dh<D>& d,
+      const int tmin, const int tmax, const int prolongation_order_time);
   
   // Destructors
   virtual ~gf ();
diff --git a/Carpet/CarpetLib/src/ggf.cc b/Carpet/CarpetLib/src/ggf.cc
index a81bb9bb2..0408417d4 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.12 2001/07/04 12:29:52 schnetter Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/ggf.cc,v 1.13 2001/12/09 16:43:10 schnetter Exp $
 
  ***************************************************************************/
 
@@ -40,8 +40,9 @@ using namespace std;
 // Constructors
 template<int D>
 generic_gf<D>::generic_gf (const string name, th& t, dh<D>& d,
-                           const int tmin, const int tmax)
-  : dimgeneric_gf(name, t, tmin, tmax),
+                           const int tmin, const int tmax,
+			   const int prolongation_order_time)
+  : dimgeneric_gf(name, t, tmin, tmax, prolongation_order_time),
     h(d.h), d(d),
     storage(tmax-tmin+1)
 {
@@ -70,7 +71,9 @@ bool generic_gf<D>::operator== (const generic_gf<D>& f) const {
 template<int D>
 void generic_gf<D>::recompose () {
   // Retain storage that might be needed
-  fdata oldstorage(tmax-tmin+1);
+  fdata oldstorage;
+  
+  oldstorage.resize(tmax-tmin+1);
   for (int tl=tmin; tl<=tmax; ++tl) {
     oldstorage[tl-tmin].resize
       (min(h.reflevels(), (int)storage[tl-tmin].size()));
@@ -88,9 +91,7 @@ void generic_gf<D>::recompose () {
     } // for rl
   } // for tl
   
-  // Delete storage
-  storage.clear();
-  
+  // Resize structure
   storage.resize(tmax-tmin+1);
   for (int tl=tmin; tl<=tmax; ++tl) {
     storage[tl-tmin].resize(h.reflevels());
@@ -98,6 +99,13 @@ void generic_gf<D>::recompose () {
       storage[tl-tmin][rl].resize(h.components(rl));
       for (int c=0; c<h.components(rl); ++c) {
       	storage[tl-tmin][rl][c].resize(h.mglevels(rl,c));
+      } // for c
+    } // for rl
+  } // for tl
+  
+  for (int rl=0; rl<h.reflevels(); ++rl) {
+    for (int tl=tmin; tl<=tmax; ++tl) {
+      for (int c=0; c<h.components(rl); ++c) {
       	for (int ml=0; ml<h.mglevels(rl,c); ++ml) {
 	  
 	  storage[tl-tmin][rl][c][ml] = typed_data();
@@ -258,8 +266,7 @@ template<int D>
 void generic_gf<D>::intercat (int tl1, int rl1, int c1, int ml1,
 			      const ibbox dh<D>::dboxes::* recv_list,
 			      const vector<int> tl2s, int rl2, int ml2,
-			      const ibbox dh<D>::dboxes::* send_list,
-			      int order_space)
+			      const ibbox dh<D>::dboxes::* send_list)
 {
   assert (tl1>=tmin && tl1<=tmax);
   assert (rl1>=0 && rl1<h.reflevels());
@@ -270,11 +277,14 @@ void generic_gf<D>::intercat (int tl1, int rl1, int c1, int ml1,
   }
   assert (rl2>=0 && rl2<h.reflevels());
   const int c2=c1;
-  assert (ml2<h.mglevels(rl2,c2));
+  assert (ml2>=0 && ml2<h.mglevels(rl2,c2));
   
   vector<const generic_data<D>*> gsrcs(tl2s.size());
   vector<int> tls(tl2s.size());
   for (int i=0; i<(int)gsrcs.size(); ++i) {
+    assert (rl2<(int)storage[tl2s[i]-tmin].size());
+    assert (c2<(int)storage[tl2s[i]-tmin][rl2].size());
+    assert (ml2<(int)storage[tl2s[i]-tmin][rl2][ml2].size());
     gsrcs[i] = storage[tl2s[i]-tmin][rl2][c2][ml2];
     tls[i] = t.time(tl2s[i],rl2,ml2);
   }
@@ -286,7 +296,8 @@ 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
-    (gsrcs, tls, recv, tl, order_space);
+    (gsrcs, tls, recv, tl,
+     d.prolongation_order_space, prolongation_order_time);
 }
 
 // Interpolate regions
@@ -294,8 +305,7 @@ template<int D>
 void generic_gf<D>::intercat (int tl1, int rl1, int c1, int ml1,
 			      const iblist dh<D>::dboxes::* recv_list,
 			      const vector<int> tl2s, int rl2, int ml2,
-			      const iblist dh<D>::dboxes::* send_list,
-			      int order_space)
+			      const iblist dh<D>::dboxes::* send_list)
 {
   assert (tl1>=tmin && tl1<=tmax);
   assert (rl1>=0 && rl1<h.reflevels());
@@ -306,11 +316,14 @@ void generic_gf<D>::intercat (int tl1, int rl1, int c1, int ml1,
   }
   assert (rl2>=0 && rl2<h.reflevels());
   const int c2=c1;
-  assert (ml2<h.mglevels(rl2,c2));
+  assert (ml2>=0 && ml2<h.mglevels(rl2,c2));
   
   vector<const generic_data<D>*> gsrcs(tl2s.size());
   vector<int> tls(tl2s.size());
   for (int i=0; i<(int)gsrcs.size(); ++i) {
+    assert (rl2<(int)storage[tl2s[i]-tmin].size());
+    assert (c2<(int)storage[tl2s[i]-tmin][rl2].size());
+    assert (ml2<(int)storage[tl2s[i]-tmin][rl2][ml2].size());
     gsrcs[i] = storage[tl2s[i]-tmin][rl2][c2][ml2];
     tls[i] = t.time(tl2s[i],rl2,ml2);
   }
@@ -325,7 +338,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
-      (gsrcs, tls, *r, tl, order_space);
+      (gsrcs, tls, *r, tl,
+       d.prolongation_order_space, prolongation_order_time);
   }
 }
 
@@ -334,8 +348,7 @@ template<int D>
 void generic_gf<D>::intercat (int tl1, int rl1, int c1, int ml1,
 			      const iblistvect dh<D>::dboxes::* recv_listvect,
 			      const vector<int> tl2s, int rl2, int ml2,
-			      const iblistvect dh<D>::dboxes::* send_listvect,
-			      int order_space)
+			      const iblistvect dh<D>::dboxes::* send_listvect)
 {
   assert (tl1>=tmin && tl1<=tmax);
   assert (rl1>=0 && rl1<h.reflevels());
@@ -347,11 +360,14 @@ void generic_gf<D>::intercat (int tl1, int rl1, int c1, int ml1,
   assert (rl2>=0 && rl2<h.reflevels());
   // walk all components
   for (int c2=0; c2<h.components(rl2); ++c2) {
-    assert (ml2<h.mglevels(rl2,c2));
+    assert (ml2>=0 && ml2<h.mglevels(rl2,c2));
     
     vector<const generic_data<D>*> gsrcs(tl2s.size());
     vector<int> tls(tl2s.size());
     for (int i=0; i<(int)gsrcs.size(); ++i) {
+      assert (rl2<(int)storage[tl2s[i]-tmin].size());
+      assert (c2<(int)storage[tl2s[i]-tmin][rl2].size());
+      assert (ml2<(int)storage[tl2s[i]-tmin][rl2][ml2].size());
       gsrcs[i] = storage[tl2s[i]-tmin][rl2][c2][ml2];
       tls[i] = t.time(tl2s[i],rl2,ml2);
     }
@@ -366,7 +382,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
-      	(gsrcs, tls, *r, tl, order_space);
+      	(gsrcs, tls, *r, tl,
+	 d.prolongation_order_space, prolongation_order_time);
     }
   }
 }
@@ -391,75 +408,83 @@ void generic_gf<D>::sync (int tl, int rl, int c, int ml) {
 
 // Prolongate the boundaries of a component
 template<int D>
-void generic_gf<D>::ref_bnd_prolongate (int tl, int rl, int c, int ml,
-					int order_space, int order_time) {
+void generic_gf<D>::ref_bnd_prolongate (int tl, int rl, int c, int ml) {
   // Interpolate
   assert (rl>=1);
-  const int tmod
-    = ((t.time(tl,rl,ml) - t.get_time(rl-1,ml)) % t.get_delta(rl-1, ml)
-       + t.get_delta(rl-1, ml)) % t.get_delta(rl-1, ml);
+//   const int tmod
+//     = ((t.time(tl,rl,ml) - t.get_time(rl-1,ml)) % t.get_delta(rl-1,ml)
+//        + t.get_delta(rl-1,ml)) % t.get_delta(rl-1,ml);
   vector<int> tl2s;
-  if (tmod == 0) {
-    // No interpolation in time
-    tl2s.resize(1);
-    tl2s[0] = tl;
-  } else {
+//   if (tmod == 0) {
+//     // No interpolation in time
+//     tl2s.resize(1);
+//     tl2s[0] = tl;
+//   } else {
     // Interpolation in time
-    assert (tmax-tmin+1 >= order_time+1);
-    tl2s.resize(order_time+1);
-    for (int i=0; i<=order_time; ++i) tl2s[i] = tmax - i;
-  }
+    assert (tmax-tmin+1 >= prolongation_order_time+1);
+    tl2s.resize(prolongation_order_time+1);
+    for (int i=0; i<=prolongation_order_time; ++i) tl2s[i] = tmax - i;
+//   }
   intercat (tl  ,rl  ,c,ml, &dh<D>::dboxes::recv_ref_bnd_coarse,
-	    tl2s,rl-1,  ml, &dh<D>::dboxes::send_ref_bnd_fine,
-	    order_space);
+	    tl2s,rl-1,  ml, &dh<D>::dboxes::send_ref_bnd_fine);
 }
 
 // Restrict a multigrid level
 template<int D>
-void generic_gf<D>::mg_restrict (int tl, int rl, int c, int ml,
-				 int order_space) {
+void generic_gf<D>::mg_restrict (int tl, int rl, int c, int ml)
+{
   // Require same times
   assert (t.get_time(rl,ml) == t.get_time(rl,ml-1));
   const vector<int> tl2s(1,tl);
   intercat (tl  ,rl,c,ml,   &dh<D>::dboxes::recv_mg_coarse,
-	    tl2s,rl,  ml-1, &dh<D>::dboxes::send_mg_fine,
-	    order_space);
+	    tl2s,rl,  ml-1, &dh<D>::dboxes::send_mg_fine);
 }
 
 // Prolongate a multigrid level
 template<int D>
-void generic_gf<D>::mg_prolongate (int tl, int rl, int c, int ml,
-				   int order_space) {
+void generic_gf<D>::mg_prolongate (int tl, int rl, int c, int ml)
+{
   // Require same times
   assert (t.get_time(rl,ml) == t.get_time(rl,ml+1));
   const vector<int> tl2s(1,tl);
   intercat (tl  ,rl,c,ml,   &dh<D>::dboxes::recv_mg_coarse,
-	    tl2s,rl,  ml+1, &dh<D>::dboxes::send_mg_fine,
-	    order_space);
+	    tl2s,rl,  ml+1, &dh<D>::dboxes::send_mg_fine);
 }
 
 // Restrict a refinement level
 template<int D>
-void generic_gf<D>::ref_restrict (int tl, int rl, int c, int ml,
-				  int order_space) {
+void generic_gf<D>::ref_restrict (int tl, int rl, int c, int ml)
+{
   // Require same times
   assert (t.get_time(rl,ml) == t.get_time(rl+1,ml));
   const vector<int> tl2s(1,tl);
   intercat (tl  ,rl  ,c,ml, &dh<D>::dboxes::recv_ref_fine,
-	    tl2s,rl+1,  ml, &dh<D>::dboxes::send_ref_coarse,
-	    order_space);
+	    tl2s,rl+1,  ml, &dh<D>::dboxes::send_ref_coarse);
 }
 
 // Prolongate a refinement level
 template<int D>
-void generic_gf<D>::ref_prolongate (int tl, int rl, int c, int ml,
-				    int order_space) {
+void generic_gf<D>::ref_prolongate (int tl, int rl, int c, int ml)
+{
   // Require same times
   assert (t.get_time(rl,ml) == t.get_time(rl-1,ml));
-  const vector<int> tl2s(1,tl);
+  assert (rl>=1);
+//   const int tmod
+//     = ((t.time(tl,rl,ml) - t.get_time(rl-1,ml)) % t.get_delta(rl-1,ml)
+//        + t.get_delta(rl-1,ml)) % t.get_delta(rl-1,ml);
+  vector<int> tl2s;
+//   if (tmod == 0) {
+//     // No interpolation in time
+//     tl2s.resize(1);
+//     tl2s[0] = tl;
+//   } else {
+    // Interpolation in time
+    assert (tmax-tmin+1 >= prolongation_order_time+1);
+    tl2s.resize(prolongation_order_time+1);
+    for (int i=0; i<=prolongation_order_time; ++i) tl2s[i] = tmax - i;
+//   }
   intercat (tl  ,rl  ,c,ml, &dh<D>::dboxes::recv_ref_coarse,
-	    tl2s,rl-1,  ml, &dh<D>::dboxes::send_ref_fine,
-	    order_space);
+	    tl2s,rl-1,  ml, &dh<D>::dboxes::send_ref_fine);
 }
 
 
diff --git a/Carpet/CarpetLib/src/ggf.hh b/Carpet/CarpetLib/src/ggf.hh
index 15d817b1b..f269d1772 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.7 2001/07/02 13:22:14 schnetter Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/ggf.hh,v 1.8 2001/12/09 16:43:10 schnetter Exp $
 
  ***************************************************************************/
 
@@ -80,7 +80,8 @@ public:
 
   // Constructors
   generic_gf (const string name, th& t, dh<D>& d,
-              const int tmin, const int tmax);
+              const int tmin, const int tmax,
+	      const int prolongation_order_time);
 
   // Destructors
   virtual ~generic_gf ();
@@ -132,22 +133,19 @@ protected:
   void intercat (int tl1, int rl1, int c1, int ml1,
 		 const ibbox dh<D>::dboxes::* recv_list,
 		 const vector<int> tl2s, int rl2, int ml2,
-		 const ibbox dh<D>::dboxes::* send_list,
-		 int order_space);
+		 const ibbox dh<D>::dboxes::* send_list);
 
   // Interpolate regions
   void intercat (int tl1, int rl1, int c1, int ml1,
 		 const iblist dh<D>::dboxes::* recv_list,
 		 const vector<int> tl2s, int rl2, int ml2,
-		 const iblist dh<D>::dboxes::* send_list,
-		 int order_space);
+		 const iblist dh<D>::dboxes::* send_list);
 
   // Interpolate regions
   void intercat (int tl1, int rl1, int c1, int ml1,
 		 const iblistvect dh<D>::dboxes::* recv_listvect,
 		 const vector<int> tl2s, int rl2, int ml2,
-		 const iblistvect dh<D>::dboxes::* send_listvect,
-		 int order_space);
+		 const iblistvect dh<D>::dboxes::* send_listvect);
 
 
 
@@ -166,20 +164,19 @@ public:
   void sync (int tl, int rl, int c, int ml);
 
   // Prolongate the boundaries of a component
-  void ref_bnd_prolongate (int tl, int rl, int c, int ml,
-			   int order_space=1, int order_time=1);
+  void ref_bnd_prolongate (int tl, int rl, int c, int ml);
 
   // Restrict a multigrid level
-  void mg_restrict (int tl, int rl, int c, int ml, int order_space=1);
+  void mg_restrict (int tl, int rl, int c, int ml);
 
   // Prolongate a multigrid level
-  void mg_prolongate (int tl, int rl, int c, int ml, int order_space=1);
+  void mg_prolongate (int tl, int rl, int c, int ml);
 
   // Restrict a refinement level
-  void ref_restrict (int tl, int rl, int c, int ml, int order_space=1);
+  void ref_restrict (int tl, int rl, int c, int ml);
 
   // Prolongate a refinement level
-  void ref_prolongate (int tl, int rl, int c, int ml, int order_space=1);
+  void ref_prolongate (int tl, int rl, int c, int ml);
   
   
   
diff --git a/Carpet/CarpetLib/src/prolongate_3d_real8.F77 b/Carpet/CarpetLib/src/prolongate_3d_real8.F77
index 81939404a..cff34ddfb 100644
--- a/Carpet/CarpetLib/src/prolongate_3d_real8.F77
+++ b/Carpet/CarpetLib/src/prolongate_3d_real8.F77
@@ -1,5 +1,5 @@
 c     -*-Fortran-*-
-c     $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8.F77,v 1.4 2001/03/26 02:27:46 eschnett Exp $
+c     $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8.F77,v 1.5 2001/12/09 16:43:10 schnetter Exp $
       
 #include "cctk.h"
       
@@ -29,12 +29,13 @@ c     bbox(:,3) is stride
       integer srcioff, srcjoff, srckoff
       integer dstioff, dstjoff, dstkoff
       
+      CCTK_REAL8 dstdiv
       integer i, j, k
       integer i0, j0, k0
-      CCTK_REAL8 fi, fj, fk
-      CCTK_REAL8 ifac(2), jfac(2), kfac(2)
+      integer fi, fj, fk
+      integer ifac(2), jfac(2), kfac(2)
       integer ii, jj, kk
-      CCTK_REAL8 fac
+      integer fac
       CCTK_REAL8 res
       integer d
       
@@ -53,7 +54,7 @@ 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
@@ -61,7 +62,7 @@ c     bbox(:,3) is stride
      $        .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 extents")
+            call CCTK_WARN (0, "Internal error: region extent is not contained in array extent")
          end if
       end do
       
@@ -95,26 +96,25 @@ c     bbox(:,3) is stride
       
       
 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
-               j0 = (srcjoff + j) / dstjfac
-               k0 = (srckoff + k) / dstkfac
-               
-               fi = mod(srcioff + i, dstifac) * one / dstifac
-               fj = mod(srcjoff + j, dstjfac) * one / dstjfac
-               fk = mod(srckoff + k, dstkfac) * one / dstkfac
-               
-               ifac(1) = (fi-1) / (-1)
-               ifac(2) = (fi  ) / 1
-               
-               jfac(1) = (fj-1) / (-1)
-               jfac(2) = (fj  ) / 1
-               
-               kfac(1) = (fk-1) / (-1)
-               kfac(2) = (fk  ) / 1
+               fi = mod(srcioff + i, dstifac)
+               ifac(1) = (fi-dstifac) * (-1)
+               ifac(2) = (fi        ) * 1
                
                res = 0
                
@@ -132,7 +132,7 @@ c     Loop over fine region
                   end do
                end do
                
-               dst (dstioff+i+1, dstjoff+j+1, dstkoff+k+1) = res
+               dst (dstioff+i+1, dstjoff+j+1, dstkoff+k+1) = dstdiv * res
                
             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 0cb837b3d..7b6e43da4 100644
--- a/Carpet/CarpetLib/src/prolongate_3d_real8_2tl.F77
+++ b/Carpet/CarpetLib/src/prolongate_3d_real8_2tl.F77
@@ -1,5 +1,5 @@
 c     -*-Fortran-*-
-c     $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_2tl.F77,v 1.3 2001/03/24 22:38:48 eschnett Exp $
+c     $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_2tl.F77,v 1.4 2001/12/09 16:43:11 schnetter Exp $
       
 #include "cctk.h"
       
@@ -35,9 +35,11 @@ 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
@@ -58,7 +60,7 @@ 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
@@ -66,7 +68,7 @@ c     bbox(:,3) is stride
      $        .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 extents")
+            call CCTK_WARN (0, "Internal error: region extent is not contained in array extent")
          end if
       end do
       
@@ -110,41 +112,33 @@ c     Linear (first order) interpolation
       
       
 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 + 1
-               j0 = (srcjoff + j) / dstjfac + 1
-               k0 = (srckoff + k) / dstkfac + 1
-               
+               i0 = (srcioff + i) / dstifac
                fi = mod(srcioff + i, dstifac)
-               fj = mod(srcjoff + j, dstjfac)
-               fk = mod(srckoff + k, dstkfac)
+               ifac(1) = (fi-dstifac) * (-1)
+               ifac(2) = (fi        ) * 1
                
                res = 0
                
-               do kk=0,1
-                  do jj=0,1
-                     do ii=0,1
-                        
-                        fac = 1
+               do kk=1,2
+                  do jj=1,2
+                     do ii=1,2
                         
-                        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
+                        fac = ifac(ii) * jfac(jj) * kfac(kk)
                         
                         if (fac.ne.0) then
                            res = res
@@ -156,8 +150,7 @@ c     Loop over fine region
                   end do
                end do
                
-               dst (dstioff+i+1, dstjoff+j+1, dstkoff+k+1)
-     $              = res / (dstifac * dstjfac * dstkfac)
+               dst (dstioff+i+1, dstjoff+j+1, dstkoff+k+1) = dstdiv * res
                
             end do
          end do
diff --git a/Carpet/CarpetLib/src/prolongate_3d_real8_2tl_o3.F77 b/Carpet/CarpetLib/src/prolongate_3d_real8_2tl_o3.F77
index 4f17e67c0..081150e87 100644
--- a/Carpet/CarpetLib/src/prolongate_3d_real8_2tl_o3.F77
+++ b/Carpet/CarpetLib/src/prolongate_3d_real8_2tl_o3.F77
@@ -1,5 +1,5 @@
 c     -*-Fortran-*-
-c     $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_2tl_o3.F77,v 1.3 2001/03/28 18:56:09 eschnett Exp $
+c     $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_2tl_o3.F77,v 1.4 2001/12/09 16:43:11 schnetter Exp $
       
 #include "cctk.h"
       
@@ -58,7 +58,7 @@ 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
@@ -66,7 +66,7 @@ c     bbox(:,3) is stride
      $        .or. regbbox(d,2)+regbbox(d,3).gt.srcbbox(d,2)
      $        .or. regbbox(d,1).lt.dstbbox(d,1)
      $        .or. regbbox(d,2).gt.dstbbox(d,2)) then
-            call CCTK_WARN (0, "Internal error: region extent is not contained in array extents")
+            call CCTK_WARN (0, "Internal error: region extent is not contained in array extent")
          end if
       end do
       
diff --git a/Carpet/CarpetLib/src/prolongate_3d_real8_3tl.F77 b/Carpet/CarpetLib/src/prolongate_3d_real8_3tl.F77
index 17834ddd1..0ec7eb6ce 100644
--- a/Carpet/CarpetLib/src/prolongate_3d_real8_3tl.F77
+++ b/Carpet/CarpetLib/src/prolongate_3d_real8_3tl.F77
@@ -1,5 +1,5 @@
 c     -*-Fortran-*-
-c     $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_3tl.F77,v 1.1 2001/03/24 22:38:48 eschnett Exp $
+c     $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_3tl.F77,v 1.2 2001/12/09 16:43:11 schnetter Exp $
       
 #include "cctk.h"
       
@@ -60,7 +60,7 @@ 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
@@ -68,7 +68,7 @@ c     bbox(:,3) is stride
      $        .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 extents")
+            call CCTK_WARN (0, "Internal error: region extent is not contained in array extent")
          end if
       end do
       
diff --git a/Carpet/CarpetLib/src/prolongate_3d_real8_3tl_o3.F77 b/Carpet/CarpetLib/src/prolongate_3d_real8_3tl_o3.F77
index 7c16e8d35..eb4a006ec 100644
--- a/Carpet/CarpetLib/src/prolongate_3d_real8_3tl_o3.F77
+++ b/Carpet/CarpetLib/src/prolongate_3d_real8_3tl_o3.F77
@@ -1,5 +1,5 @@
 c     -*-Fortran-*-
-c     $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_3tl_o3.F77,v 1.3 2001/03/28 18:56:09 eschnett Exp $
+c     $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_3tl_o3.F77,v 1.4 2001/12/09 16:43:11 schnetter Exp $
       
 #include "cctk.h"
       
@@ -60,7 +60,7 @@ 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
@@ -68,7 +68,7 @@ c     bbox(:,3) is stride
      $        .or. regbbox(d,2)+regbbox(d,3).gt.srcbbox(d,2)
      $        .or. regbbox(d,1).lt.dstbbox(d,1)
      $        .or. regbbox(d,2).gt.dstbbox(d,2)) then
-            call CCTK_WARN (0, "Internal error: region extent is not contained in array extents")
+            call CCTK_WARN (0, "Internal error: region extent is not contained in array extent")
          end if
       end do
       
diff --git a/Carpet/CarpetLib/src/prolongate_3d_real8_o3.F77 b/Carpet/CarpetLib/src/prolongate_3d_real8_o3.F77
index 5ee417885..906b67dd6 100644
--- a/Carpet/CarpetLib/src/prolongate_3d_real8_o3.F77
+++ b/Carpet/CarpetLib/src/prolongate_3d_real8_o3.F77
@@ -1,5 +1,5 @@
 c     -*-Fortran-*-
-c     $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_o3.F77,v 1.3 2001/03/28 18:56:09 eschnett Exp $
+c     $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_o3.F77,v 1.4 2001/12/09 16:43:11 schnetter Exp $
       
 #include "cctk.h"
       
@@ -29,12 +29,13 @@ c     bbox(:,3) is stride
       integer srcioff, srcjoff, srckoff
       integer dstioff, dstjoff, dstkoff
       
+      CCTK_REAL8 dstdiv
       integer i, j, k
       integer i0, j0, k0
-      CCTK_REAL8 fi, fj, fk
-      CCTK_REAL8 ifac(4), jfac(4), kfac(4)
+      integer fi, fj, fk
+      integer ifac(4), jfac(4), kfac(4)
       integer ii, jj, kk
-      CCTK_REAL8 fac
+      integer fac
       CCTK_REAL8 res
       integer d
       
@@ -53,7 +54,7 @@ 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
@@ -61,7 +62,7 @@ c     bbox(:,3) is stride
      $        .or. regbbox(d,2)+regbbox(d,3).gt.srcbbox(d,2)
      $        .or. regbbox(d,1).lt.dstbbox(d,1)
      $        .or. regbbox(d,2).gt.dstbbox(d,2)) then
-            call CCTK_WARN (0, "Internal error: region extent is not contained in array extents")
+            call CCTK_WARN (0, "Internal error: region extent is not contained in array extent")
          end if
       end do
       
@@ -95,32 +96,31 @@ c     bbox(:,3) is stride
       
       
 c     Loop over fine region
+      dstdiv = one / (6*dstifac * 6*dstjfac * 6*dstkfac)
+      
       do k = 0, regkext-1
+         k0 = (srckoff + k) / dstkfac
+         fk = mod(srckoff + k, dstkfac)
+         kfac(1) = (fk        ) * (fk-dstkfac) * (fk-2*dstkfac) * (-1)
+         kfac(2) = (fk+dstkfac) * (fk-dstkfac) * (fk-2*dstkfac) * 3
+         kfac(3) = (fk+dstkfac) * (fk        ) * (fk-2*dstkfac) * (-3)
+         kfac(4) = (fk+dstkfac) * (fk        ) * (fk-  dstkfac) * 1
+         
          do j = 0, regjext-1
+            j0 = (srcjoff + j) / dstjfac
+            fj = mod(srcjoff + j, dstjfac)
+            jfac(1) = (fj        ) * (fj-dstjfac) * (fj-2*dstjfac) * (-1)
+            jfac(2) = (fj+dstjfac) * (fj-dstjfac) * (fj-2*dstjfac) * 3
+            jfac(3) = (fj+dstjfac) * (fj        ) * (fj-2*dstjfac) * (-3)
+            jfac(4) = (fj+dstjfac) * (fj        ) * (fj-  dstjfac) * 1
+            
             do i = 0, regiext-1
-               
                i0 = (srcioff + i) / dstifac
-               j0 = (srcjoff + j) / dstjfac
-               k0 = (srckoff + k) / dstkfac
-               
-               fi = mod(srcioff + i, dstifac) * one / dstifac
-               fj = mod(srcjoff + j, dstjfac) * one / dstjfac
-               fk = mod(srckoff + k, dstkfac) * one / dstkfac
-               
-               ifac(1) = (fi  ) * (fi-1) * (fi-2) / (-6)
-               ifac(2) = (fi+1) * (fi-1) * (fi-2) / 2
-               ifac(3) = (fi+1) * (fi  ) * (fi-2) / (-2)
-               ifac(4) = (fi+1) * (fi  ) * (fi-1) / 6
-               
-               jfac(1) = (fj  ) * (fj-1) * (fj-2) / (-6)
-               jfac(2) = (fj+1) * (fj-1) * (fj-2) / 2
-               jfac(3) = (fj+1) * (fj  ) * (fj-2) / (-2)
-               jfac(4) = (fj+1) * (fj  ) * (fj-1) / 6
-               
-               kfac(1) = (fk  ) * (fk-1) * (fk-2) / (-6)
-               kfac(2) = (fk+1) * (fk-1) * (fk-2) / 2
-               kfac(3) = (fk+1) * (fk  ) * (fk-2) / (-2)
-               kfac(4) = (fk+1) * (fk  ) * (fk-1) / 6
+               fi = mod(srcioff + i, dstifac)
+               ifac(1) = (fi        ) * (fi-dstifac) * (fi-2*dstifac) * (-1)
+               ifac(2) = (fi+dstifac) * (fi-dstifac) * (fi-2*dstifac) * 3
+               ifac(3) = (fi+dstifac) * (fi        ) * (fi-2*dstifac) * (-3)
+               ifac(4) = (fi+dstifac) * (fi        ) * (fi-  dstifac) * 1
                
                res = 0
                
@@ -138,7 +138,7 @@ c     Loop over fine region
                   end do
                end do
                
-               dst (dstioff+i+1, dstjoff+j+1, dstkoff+k+1) = res
+               dst (dstioff+i+1, dstjoff+j+1, dstkoff+k+1) = dstdiv * res
                
             end do
          end do
diff --git a/Carpet/CarpetLib/src/restrict_3d_real8.F77 b/Carpet/CarpetLib/src/restrict_3d_real8.F77
index d580f423c..a8ba6c64c 100644
--- a/Carpet/CarpetLib/src/restrict_3d_real8.F77
+++ b/Carpet/CarpetLib/src/restrict_3d_real8.F77
@@ -1,5 +1,5 @@
 c     -*-Fortran-*-
-c     $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/restrict_3d_real8.F77,v 1.3 2001/03/24 22:38:48 eschnett Exp $
+c     $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/restrict_3d_real8.F77,v 1.4 2001/12/09 16:43:11 schnetter Exp $
       
 #include "cctk.h"
       
@@ -44,7 +44,7 @@ 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
@@ -52,7 +52,7 @@ c     bbox(:,3) is stride
      $        .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 extents")
+            call CCTK_WARN (0, "Internal error: region extent is not contained in array extent")
          end if
       end do
       
diff --git a/Carpet/CarpetLib/src/th.cc b/Carpet/CarpetLib/src/th.cc
index 18e26de1e..1e01d4787 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.5 2001/06/12 14:57:00 schnetter Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/th.cc,v 1.6 2001/12/09 16:43:11 schnetter Exp $
 
  ***************************************************************************/
 
@@ -35,7 +35,8 @@ using namespace std;
 
 
 // Constructors
-th::th (dimgeneric_gh* h, const int basedelta) : h(h), delta(basedelta) {
+th::th (dimgeneric_gh* h, const int basedelta)
+  : h(h), delta(basedelta) {
   h->add(this);
 }