Added support for higher-order interpolation in space and time.

darcs-hash:20010324213842-f6438-3ccfdd7797b28055f08d28b77e33205c69c60e27.gz

16 files changed, 606 insertions, 687 deletions

diff --git a/Carpet/CarpetLib/src/copy_3d_real8.F77 b/Carpet/CarpetLib/src/copy_3d_real8.F77
index ae5d1c34f..e00dd8d54 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.2 2001/03/22 18:42:05 eschnett Exp $
+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 $
       
 #include "cctk.h"
       
@@ -43,6 +43,12 @@ c     bbox(:,3) is stride
      $        .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).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 extents")
+         end if
       end do
       
       if (srciext.ne.(srcbbox(1,2)-srcbbox(1,1))/srcbbox(1,3)+1
diff --git a/Carpet/CarpetLib/src/data.cc b/Carpet/CarpetLib/src/data.cc
index 23cd28b57..b75c770ff 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.8 2001/03/22 18:42:05 eschnett Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/data.cc,v 1.9 2001/03/24 22:38:48 eschnett Exp $
 
  ***************************************************************************/
 
@@ -198,109 +198,64 @@ void data<T,D>
 
 template<class T, int D>
 void data<T,D>
-::interpolate_from_innerloop (const generic_data<D>* gsrc, const ibbox& box)
+::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 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 (has_storage());
+  assert (all(box.lower()>=extent().lower()));
+  assert (all(box.upper()<=extent().upper()));
   assert (all(box.stride()==extent().stride()));
   assert (all((box.lower()-extent().lower())%box.stride() == 0));
-  
-  assert (proc() == src->proc());
+  vector<const data*> srcs(gsrcs.size());
+  for (int t=0; t<(int)srcs.size(); ++t) srcs[t] = (const data*)gsrcs[t];
+  assert (srcs.size() == tls.size() && srcs.size()>0);
+  for (int t=0; t<(int)srcs.size(); ++t) {
+    assert (srcs[t]->has_storage());
+    assert (all(box.lower()>=srcs[t]->extent().lower()));
+    assert (all(box.upper()<=srcs[t]->extent().upper()));
+    assert (proc() == srcs[t]->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;
-    
-    // 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
-  
-}
-
-
-
-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 (order_space == 1);
   
-  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));
+  vector<T> src_fac(srcs.size());
+  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]);
+      }
+    }
+  }
   
   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(src1->extent());
+      = ibbox(pos,pos,extent().stride()).expanded_for(srcs[0]->extent());
     
     // interpolate from box to position
-    T src1_sum = 0;
-    T src2_sum = 0;
+    T 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];
+	const ivect str = srcs[0]->extent().stride();
+	const T f = prod(vect<T,D>(str - abs(pos - frompos)) / vect<T,D>(str));
+	for (int t=0; t<(int)src_fac.size(); ++t) {
+	  sum += f * src_fac[t] * (*srcs[t])[frompos];
+	}
       }
-    (*this)[pos] = src1_fac * src1_sum + src2_fac * src2_sum;
+    (*this)[pos] = sum;
     
   } // for pos
   
@@ -381,6 +336,18 @@ void data<CCTK_REAL8,3>
 
 
 extern "C" {
+  
+  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]);
+  
+  
+  
   void CCTK_FCALL CCTK_FNAME(prolongate_3d_real8)
     (const CCTK_REAL8* src,
      const int& srciext, const int& srcjext, const int& srckext,
@@ -389,7 +356,7 @@ extern "C" {
      const int srcbbox[3][3],
      const int dstbbox[3][3],
      const int regbbox[3][3]);
-  void CCTK_FCALL CCTK_FNAME(restrict_3d_real8)
+  void CCTK_FCALL CCTK_FNAME(prolongate_3d_real8_o3)
     (const CCTK_REAL8* src,
      const int& srciext, const int& srcjext, const int& srckext,
      CCTK_REAL8* dst,
@@ -397,77 +364,7 @@ extern "C" {
      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;
-  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());
-  
-  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];
-    
-    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 {
-    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,
@@ -477,40 +374,67 @@ extern "C" {
      const int srcbbox[3][3],
      const int dstbbox[3][3],
      const int regbbox[3][3]);
+  void CCTK_FCALL CCTK_FNAME(prolongate_3d_real8_2tl_o3)
+    (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]);
+  
+  void CCTK_FCALL CCTK_FNAME(prolongate_3d_real8_3tl)
+    (const CCTK_REAL8* src1, const int& t1,
+     const CCTK_REAL8* src2, const int& t2,
+     const CCTK_REAL8* src3, const int& t3,
+     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]);
+  void CCTK_FCALL CCTK_FNAME(prolongate_3d_real8_3tl_o3)
+    (const CCTK_REAL8* src1, const int& t1,
+     const CCTK_REAL8* src2, const int& t2,
+     const CCTK_REAL8* src3, const int& t3,
+     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)
+::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 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 (has_storage());
+  assert (all(box.lower()>=extent().lower()));
+  assert (all(box.upper()<=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 (all((box.lower()-extent().lower())%box.stride() == 0));
+  vector<const data*> srcs(gsrcs.size());
+  for (int t=0; t<(int)srcs.size(); ++t) srcs[t] = (const data*)gsrcs[t];
+  assert (srcs.size() == tls.size() && srcs.size()>0);
+  for (int t=0; t<(int)srcs.size(); ++t) {
+    assert (srcs[t]->has_storage());
+    assert (all(box.lower()>=srcs[t]->extent().lower()));
+    assert (all(box.upper()<=srcs[t]->extent().upper()));
+  }
   
-  assert (proc() == src1->proc());
+  assert (proc() == srcs[0]->proc());
   
   int rank;
   MPI_Comm_rank (dist::comm, &rank);
   assert (rank == proc());
   
-  assert (src1->extent() == src2->extent());
-  
-  const ibbox& sext = src1->extent();
+  const ibbox& sext = srcs[0]->extent();
   const ibbox& dext = extent();
   
   int srcshp[3], dstshp[3];
@@ -534,16 +458,102 @@ void data<CCTK_REAL8,3>
   }
   
   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,
+  if (all(sext.stride() < dext.stride())) {
+    
+    assert (srcs.size() == 1);
+    CCTK_FNAME(restrict_3d_real8)
+      ((const CCTK_REAL8*)srcs[0]->storage(),
        srcshp[0], srcshp[1], srcshp[2],
-       (CCTK_REAL8*)storage(), t,
+       (CCTK_REAL8*)storage(),
        dstshp[0], dstshp[1], dstshp[2],
-       srcbbox,
-       dstbbox,
-       regbbox);
+       srcbbox, dstbbox, regbbox);
+    
+  } else if (all(sext.stride() > dext.stride())) {
+    
+    switch (srcs.size()) {
+      
+    case 1:
+      // One timelevel
+      switch (order_space) {
+      case 1:
+	CCTK_FNAME(prolongate_3d_real8)
+	  ((const CCTK_REAL8*)srcs[0]->storage(),
+	   srcshp[0], srcshp[1], srcshp[2],
+	   (CCTK_REAL8*)storage(),
+	   dstshp[0], dstshp[1], dstshp[2],
+	   srcbbox, dstbbox, regbbox);
+	break;
+      case 3:
+	CCTK_FNAME(prolongate_3d_real8_o3)
+	  ((const CCTK_REAL8*)srcs[0]->storage(),
+	   srcshp[0], srcshp[1], srcshp[2],
+	   (CCTK_REAL8*)storage(),
+	   dstshp[0], dstshp[1], dstshp[2],
+	   srcbbox, dstbbox, regbbox);
+	break;
+      default:
+	abort();
+      }
+      break;
+      
+    case 2:
+      // Two timelevels
+      switch (order_space) {
+      case 1:
+	CCTK_FNAME(prolongate_3d_real8_2tl)
+	  ((const CCTK_REAL8*)srcs[0]->storage(), tls[0],
+	   (const CCTK_REAL8*)srcs[1]->storage(), tls[1],
+	   srcshp[0], srcshp[1], srcshp[2],
+	   (CCTK_REAL8*)storage(), tl,
+	   dstshp[0], dstshp[1], dstshp[2],
+	   srcbbox, dstbbox, regbbox);
+	break;
+      case 3:
+	CCTK_FNAME(prolongate_3d_real8_2tl_o3)
+	  ((const CCTK_REAL8*)srcs[0]->storage(), tls[0],
+	   (const CCTK_REAL8*)srcs[1]->storage(), tls[1],
+	   srcshp[0], srcshp[1], srcshp[2],
+	   (CCTK_REAL8*)storage(), tl,
+	   dstshp[0], dstshp[1], dstshp[2],
+	   srcbbox, dstbbox, regbbox);
+	break;
+      default:
+	abort();
+      }
+      break;
+      
+    case 3:
+      // Three timelevels
+      switch (order_space) {
+      case 1:
+	CCTK_FNAME(prolongate_3d_real8_3tl)
+	  ((const CCTK_REAL8*)srcs[0]->storage(), tls[0],
+	   (const CCTK_REAL8*)srcs[1]->storage(), tls[1],
+	   (const CCTK_REAL8*)srcs[2]->storage(), tls[2],
+	   srcshp[0], srcshp[1], srcshp[2],
+	   (CCTK_REAL8*)storage(), tl,
+	   dstshp[0], dstshp[1], dstshp[2],
+	   srcbbox, dstbbox, regbbox);
+	break;
+      case 3:
+	CCTK_FNAME(prolongate_3d_real8_3tl_o3)
+	  ((const CCTK_REAL8*)srcs[0]->storage(), tls[0],
+	   (const CCTK_REAL8*)srcs[1]->storage(), tls[1],
+	   (const CCTK_REAL8*)srcs[2]->storage(), tls[2],
+	   srcshp[0], srcshp[1], srcshp[2],
+	   (CCTK_REAL8*)storage(), tl,
+	   dstshp[0], dstshp[1], dstshp[2],
+	   srcbbox, dstbbox, regbbox);
+	break;
+      default:
+	abort();
+      }
+      break;
+      
+    default:
+      abort();
+    }
+    
   } else {
     abort();
   }
@@ -564,7 +574,7 @@ void data<T,D>::write_ascii_output_element (ofstream& file, const ivect& index)
 // Output
 template<class T,int D>
 ostream& data<T,D>::out (ostream& os) const {
-  os << "data<" STR(T) "," << D << ">:"
+  os << "data<T," << D << ">:"
      << "extent=" << extent() << ","
      << "stride=" << stride() << ",size=" << size();
   return os;
diff --git a/Carpet/CarpetLib/src/data.hh b/Carpet/CarpetLib/src/data.hh
index c4594c56d..5182d0379 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.4 2001/03/22 18:42:05 eschnett Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/data.hh,v 1.5 2001/03/24 22:38:48 eschnett Exp $
 
  ***************************************************************************/
 
@@ -99,13 +99,12 @@ public:
   }
   
   // Data manipulators
-  void copy_from_innerloop (const generic_data<D>* src,
+  void copy_from_innerloop (const generic_data<D>* gsrc,
 			    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 interpolate_from_innerloop (const vector<const generic_data<D>*> gsrcs,
+				   const vector<int> tls,
+				   const ibbox& box, const int tl,
+				   const int order_space);
   
   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 bf483c62f..aa25b17b8 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.3 2001/03/22 18:42:05 eschnett Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/defs.cc,v 1.4 2001/03/24 22:38:48 eschnett Exp $
 
  ***************************************************************************/
 
@@ -79,6 +79,7 @@ ostream& operator<< (ostream& os, const vector<T>& v) {
 template ostream& operator<< (ostream& os, const list<bbox<int,3> >& l);
 template ostream& operator<< (ostream& os, const set<bbox<int,3> >& s);
 template ostream& operator<< (ostream& os, const set<bboxset<int,3> >& s);
+template ostream& operator<< (ostream& os, const vector<int>& v);
 template ostream& operator<< (ostream& os, const vector<list<bbox<int,3> > >& v);
 template ostream& operator<< (ostream& os, const vector<vector<bbox<int,3> > >& v);
 template ostream& operator<< (ostream& os, const vector<vector<vector<bbox<int,3> > > >& v);
diff --git a/Carpet/CarpetLib/src/gdata.cc b/Carpet/CarpetLib/src/gdata.cc
index 5418db23a..6f2a2ecd3 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.9 2001/03/22 18:42:05 eschnett Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gdata.cc,v 1.10 2001/03/24 22:38:48 eschnett Exp $
 
  ***************************************************************************/
 
@@ -88,78 +88,38 @@ void generic_data<D>::copy_from (const generic_data* src, const ibbox& box)
 
 template<int D>
 void generic_data<D>
-::interpolate_from (const generic_data* src, const ibbox& box)
+::interpolate_from (const vector<const generic_data*> srcs,
+		    const vector<int> tls,
+		    const ibbox& box, const int tl,
+		    const int order_space)
 {
-  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 (has_storage());
+  assert (all(box.lower()>=extent().lower()));
+  assert (all(box.upper()<=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;
-    
+  assert (srcs.size() == tls.size() && srcs.size()>0);
+  for (int t=0; t<(int)srcs.size(); ++t) {
+    assert (srcs[t]->has_storage());
+    assert (all(box.lower()>=srcs[t]->extent().lower()));
+    assert (all(box.upper()<=srcs[t]->extent().upper()));
   }
-}
-
-
-
-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()) {
+  if (proc() == srcs[0]->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);
+      interpolate_from_innerloop (srcs, tls, box, tl, order_space);
     }
     
   } 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->allocate (box, srcs[0]->proc());
+    tmp->interpolate_from (srcs, tls, box, tl, order_space);
     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 1a16f92c2..ddb1c0603 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.5 2001/03/22 18:42:05 eschnett Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gdata.hh,v 1.6 2001/03/24 22:38:48 eschnett Exp $
 
  ***************************************************************************/
 
@@ -135,19 +135,18 @@ public:
 
   // Data manipulators
   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);
+  void interpolate_from (const vector<const generic_data*> srcs,
+			 const vector<int> tls,
+			 const ibbox& box, const int tl,
+			 const int order_space);
 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;
+  interpolate_from_innerloop (const vector<const generic_data*> srcs,
+			      const vector<int> tls,
+			      const ibbox& box, const int tl,
+			      const int order_space) = 0;
   
 public:
   
diff --git a/Carpet/CarpetLib/src/ggf.cc b/Carpet/CarpetLib/src/ggf.cc
index 09966ac69..2fbb0eae3 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.6 2001/03/22 18:42:05 eschnett Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/ggf.cc,v 1.7 2001/03/24 22:38:48 eschnett Exp $
 
  ***************************************************************************/
 
@@ -20,6 +20,7 @@
  ***************************************************************************/
 
 #include <assert.h>
+#include <stdlib.h>
 
 #include <iostream>
 #include <string>
@@ -170,7 +171,7 @@ void generic_gf<D>::cycle (int rl, int c, int ml) {
 
 // Operations
 
-// Copy region for a component (between time levels)
+// Copy a region
 template<int D>
 void generic_gf<D>::copycat (int tl1, int rl1, int c1, int ml1,
 			     const ibbox dh<D>::dboxes::* recv_list,
@@ -194,7 +195,7 @@ void generic_gf<D>::copycat (int tl1, int rl1, int c1, int ml1,
     (storage[tl2-tmin][rl2][c2][ml2], recv);
 }
 
-// Copy regions for a component (between multigrid levels)
+// Copy regions
 template<int D>
 void generic_gf<D>::copycat (int tl1, int rl1, int c1, int ml1,
 			     const iblist dh<D>::dboxes::* recv_list,
@@ -217,12 +218,12 @@ void generic_gf<D>::copycat (int tl1, int rl1, int c1, int ml1,
        r!=recv.end(); ++r, ++s) {
     // (use the send boxes for communication)
     // copy the content
-    storage[tl1-tmin][rl1][c1][ml1]->interpolate_from
+    storage[tl1-tmin][rl1][c1][ml1]->copy_from
       (storage[tl2-tmin][rl2][c2][ml2], *r);
   }
 }
 
-// Copy regions for a level (between refinement levels)
+// Copy regions
 template<int D>
 void generic_gf<D>::copycat (int tl1, int rl1, int c1, int ml1,
 			     const iblistvect dh<D>::dboxes::* recv_listvect,
@@ -246,54 +247,75 @@ void generic_gf<D>::copycat (int tl1, int rl1, int c1, int ml1,
       	 r!=recv.end(); ++r, ++s) {
       // (use the send boxes for communication)
       // copy the content
-      storage[tl1-tmin][rl1][c1][ml1]->interpolate_from
+      storage[tl1-tmin][rl1][c1][ml1]->copy_from
       	(storage[tl2-tmin][rl2][c2][ml2], *r);
     }
   }
 }
 
-// Interpolate a component (between time levels)
+// Interpolate a region
 template<int D>
 void generic_gf<D>::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)
+			      const vector<int> tl2s, int rl2, int ml2,
+			      const ibbox dh<D>::dboxes::* send_list,
+			      int order_space)
 {
   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 (tl2a>=tmin && tl2a<=tmax);
-  assert (tl2b>=tmin && tl2b<=tmax);
+  for (int i=0; i<(int)tl2s.size(); ++i) {
+    assert (tl2s[i]>=tmin && tl2s[i]<=tmax);
+  }
   assert (rl2>=0 && rl2<h.reflevels());
   const int c2=c1;
   assert (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) {
+    gsrcs[i] = storage[tl2s[i]-tmin][rl2][c2][ml2];
+    tls[i] = tl2s[i] * t.get_delta(rl2,ml2);
+  }
+  const int tl = tl1 * t.get_delta(rl1,ml1);
+  
   const ibbox recv = d.boxes[rl1][c1][ml1].*recv_list;
   const ibbox send = d.boxes[rl2][c2][ml2].*send_list;
   assert (recv.size()==send.size());
   // interpolate the content
   assert (recv==send);
   storage[tl1-tmin][rl1][c1][ml1]->interpolate_from
-    (storage[tl2a-tmin][rl2][c2][ml2], tl2a,
-     storage[tl2b-tmin][rl2][c2][ml2], tl2b, recv, tl1);
+    (gsrcs, tls, recv, tl, order_space);
 }
 
-// Interpolate a component (between multigrid levels)
+// Interpolate regions
 template<int D>
 void generic_gf<D>::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)
+			      const vector<int> tl2s, int rl2, int ml2,
+			      const iblist dh<D>::dboxes::* send_list,
+			      int order_space)
 {
   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 (tl2a>=tmin && tl2a<=tmax);
-  assert (tl2b>=tmin && tl2b<=tmax);
+  for (int i=0; i<(int)tl2s.size(); ++i) {
+    assert (tl2s[i]>=tmin && tl2s[i]<=tmax);
+  }
   assert (rl2>=0 && rl2<h.reflevels());
   const int c2=c1;
   assert (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) {
+    gsrcs[i] = storage[tl2s[i]-tmin][rl2][c2][ml2];
+    tls[i] = tl2s[i] * t.get_delta(rl2,ml2);
+  }
+  const int tl = tl1 * t.get_delta(rl1,ml1);
+  
   const iblist recv = d.boxes[rl1][c1][ml1].*recv_list;
   const iblist send = d.boxes[rl2][c2][ml2].*send_list;
   assert (recv.size()==send.size());
@@ -303,28 +325,38 @@ 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[tl2a-tmin][rl2][c2][ml2], tl2a,
-       storage[tl2b-tmin][rl2][c2][ml2], tl2b, *r, tl1);
+      (gsrcs, tls, *r, tl, order_space);
   }
 }
 
-// Interpolate a level (between refinement levels)
+// Interpolate regions
 template<int D>
 void generic_gf<D>::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)
+			      const vector<int> tl2s, int rl2, int ml2,
+			      const iblistvect dh<D>::dboxes::* send_listvect,
+			      int order_space)
 {
   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 (tl2a>=tmin && tl2a<=tmax);
-  assert (tl2b>=tmin && tl2b<=tmax);
+  for (int i=0; i<(int)tl2s.size(); ++i) {
+    assert (tl2s[i]>=tmin && tl2s[i]<=tmax);
+  }
   assert (rl2>=0 && rl2<h.reflevels());
   // walk all components
   for (int c2=0; c2<h.components(rl2); ++c2) {
     assert (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) {
+      gsrcs[i] = storage[tl2s[i]-tmin][rl2][c2][ml2];
+      tls[i] = tl2s[i] * t.get_delta(rl2,ml2);
+    }
+    const int tl = tl1 * t.get_delta(rl1,ml1);
+    
     const iblist recv = (d.boxes[rl1][c1][ml1].*recv_listvect)[c2];
     const iblist send = (d.boxes[rl2][c2][ml2].*send_listvect)[c1];
     assert (recv.size()==send.size());
@@ -334,8 +366,7 @@ 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[tl2a-tmin][rl2][c2][ml2], tl2a,
-      	 storage[tl2b-tmin][rl2][c2][ml2], tl2b, *r, tl1);
+      	(gsrcs, tls, *r, tl, order_space);
     }
   }
 }
@@ -360,59 +391,75 @@ 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) {
+void generic_gf<D>::ref_bnd_prolongate (int tl, int rl, int c, int ml,
+					int order_space, int order_time) {
   // Interpolate
   assert (rl>=1);
-  double time =
-    (t.time(tl,rl,ml) - t.get_time(rl-1,ml)) / (double)t.get_delta(rl-1, ml);
-  const int tl2 = (int)floor(time + 1e-6);
-  assert (tl2>=tmin && tl2<=tmax);
-  if (time==tl2) {
-    copycat (tl ,rl  ,c,ml, &dh<D>::dboxes::recv_ref_bnd_coarse,
-      	     tl2,rl-1,  ml, &dh<D>::dboxes::send_ref_bnd_fine);
+  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 {
-    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,
-      	      tl2a, tl2b, rl-1,ml, &dh<D>::dboxes::send_ref_bnd_fine);
+    // 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;
   }
+  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);
 }
 
 // Restrict a multigrid level
 template<int D>
-void generic_gf<D>::mg_restrict (int tl, int rl, int c, int ml) {
+void generic_gf<D>::mg_restrict (int tl, int rl, int c, int ml,
+				 int order_space) {
   // Require same times
   assert (t.get_time(rl,ml) == t.get_time(rl,ml-1));
-  copycat (tl,rl,c,ml,   &dh<D>::dboxes::recv_mg_coarse,
-      	   tl,rl,  ml-1, &dh<D>::dboxes::send_mg_fine);
+  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);
 }
 
 // Prolongate a multigrid level
 template<int D>
-void generic_gf<D>::mg_prolongate (int tl, int rl, int c, int ml) {
+void generic_gf<D>::mg_prolongate (int tl, int rl, int c, int ml,
+				   int order_space) {
   // Require same times
   assert (t.get_time(rl,ml) == t.get_time(rl,ml+1));
-  copycat (tl,rl,c,ml,   &dh<D>::dboxes::recv_mg_coarse,
-      	   tl,rl,  ml+1, &dh<D>::dboxes::send_mg_fine);
+  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);
 }
 
 // Restrict a refinement level
 template<int D>
-void generic_gf<D>::ref_restrict (int tl, int rl, int c, int ml) {
+void generic_gf<D>::ref_restrict (int tl, int rl, int c, int ml,
+				  int order_space) {
   // Require same times
   assert (t.get_time(rl,ml) == t.get_time(rl+1,ml));
-  copycat (tl,rl  ,c,ml, &dh<D>::dboxes::recv_ref_fine,
-      	   tl,rl+1,  ml, &dh<D>::dboxes::send_ref_coarse);
+  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);
 }
 
 // Prolongate a refinement level
 template<int D>
-void generic_gf<D>::ref_prolongate (int tl, int rl, int c, int ml) {
+void generic_gf<D>::ref_prolongate (int tl, int rl, int c, int ml,
+				    int order_space) {
   // Require same times
   assert (t.get_time(rl,ml) == t.get_time(rl-1,ml));
-  copycat (tl,rl  ,c,ml, &dh<D>::dboxes::recv_ref_coarse,
-      	   tl,rl-1,  ml, &dh<D>::dboxes::send_ref_fine);
+  const vector<int> tl2s(1,tl);
+  intercat (tl  ,rl  ,c,ml, &dh<D>::dboxes::recv_ref_coarse,
+	    tl2s,rl-1,  ml, &dh<D>::dboxes::send_ref_fine,
+	    order_space);
 }
 
 
diff --git a/Carpet/CarpetLib/src/ggf.hh b/Carpet/CarpetLib/src/ggf.hh
index 6e998c145..6cacc0724 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.3 2001/03/22 18:42:06 eschnett Exp $
+    $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/ggf.hh,v 1.4 2001/03/24 22:38:48 eschnett Exp $
 
  ***************************************************************************/
 
@@ -26,6 +26,7 @@
 
 #include <iostream>
 #include <string>
+#include <vector>
 
 #include "defs.hh"
 #include "dh.hh"
@@ -110,41 +111,44 @@ protected:
   
 protected:
   
-  // Copy region for a component (between time levels)
+  // Copy a region
   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)
+  // Copy regions
   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)
+  // Copy regions
   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)
+  // Interpolate a region
   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);
+		 const vector<int> tl2s, int rl2, int ml2,
+		 const ibbox dh<D>::dboxes::* send_list,
+		 int order_space);
 
-  // Interpolate a component (between multigrid levels)
+  // Interpolate regions
   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);
+		 const vector<int> tl2s, int rl2, int ml2,
+		 const iblist dh<D>::dboxes::* send_list,
+		 int order_space);
 
-  // Interpolate a level (between refinement levels)
+  // Interpolate regions
   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);
+		 const vector<int> tl2s, int rl2, int ml2,
+		 const iblistvect dh<D>::dboxes::* send_listvect,
+		 int order_space);
 
 
 
@@ -163,19 +167,20 @@ 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);
+  void ref_bnd_prolongate (int tl, int rl, int c, int ml,
+			   int order_space=1, int order_time=1);
 
   // Restrict a multigrid level
-  void mg_restrict (int tl, int rl, int c, int ml);
+  void mg_restrict (int tl, int rl, int c, int ml, int order_space=1);
 
   // Prolongate a multigrid level
-  void mg_prolongate (int tl, int rl, int c, int ml);
+  void mg_prolongate (int tl, int rl, int c, int ml, int order_space=1);
 
   // Restrict a refinement level
-  void ref_restrict (int tl, int rl, int c, int ml);
+  void ref_restrict (int tl, int rl, int c, int ml, int order_space=1);
 
   // Prolongate a refinement level
-  void ref_prolongate (int tl, int rl, int c, int ml);
+  void ref_prolongate (int tl, int rl, int c, int ml, int order_space=1);
   
   
   
diff --git a/Carpet/CarpetLib/src/make.code.defn b/Carpet/CarpetLib/src/make.code.defn
index dd1756780..b072b8c62 100644
--- a/Carpet/CarpetLib/src/make.code.defn
+++ b/Carpet/CarpetLib/src/make.code.defn
@@ -1,9 +1,16 @@
 # Main make.code.defn file for thorn CarpetLib   -*-Makefile-*-
-# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/make.code.defn,v 1.3 2001/03/22 18:42:06 eschnett Exp $
+# $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/make.code.defn,v 1.4 2001/03/24 22:38:48 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 \
-	copy_3d_real8.F77 prolongate_3d_real8.F77 prolongate_3d_real8_2tl.F77 restrict_3d_real8.F77
+	copy_3d_real8.F77			\
+	prolongate_3d_real8.F77			\
+	prolongate_3d_real8_o3.F77		\
+	prolongate_3d_real8_2tl.F77		\
+	prolongate_3d_real8_2tl_o3.F77		\
+	prolongate_3d_real8_3tl.F77		\
+	prolongate_3d_real8_3tl_o3.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 04053c524..531a3a6a7 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.2 2001/03/22 18:42:06 eschnett Exp $
+c     $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8.F77,v 1.3 2001/03/24 22:38:48 eschnett Exp $
       
 #include "cctk.h"
       
@@ -53,6 +53,12 @@ c     bbox(:,3) is stride
      $        .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).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 extents")
+         end if
       end do
       
       if (srciext.ne.(srcbbox(1,2)-srcbbox(1,1))/srcbbox(1,3)+1
diff --git a/Carpet/CarpetLib/src/prolongate_3d_real8_2tl.F77 b/Carpet/CarpetLib/src/prolongate_3d_real8_2tl.F77
index f9b380c5d..0cb837b3d 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.2 2001/03/22 18:42:06 eschnett Exp $
+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 $
       
 #include "cctk.h"
       
@@ -62,6 +62,12 @@ c     bbox(:,3) is stride
      $        .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).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 extents")
+         end if
       end do
       
       if (srciext.ne.(srcbbox(1,2)-srcbbox(1,1))/srcbbox(1,3)+1
@@ -93,6 +99,11 @@ c     bbox(:,3) is stride
       
       
       
+c     Linear (first order) interpolation
+      if (t1.eq.t2) then
+         call CCTK_WARN (0, "Internal error: arrays have same time")
+      end if
+      
       s1fac = (t - t2) * one / (t1 - t2)
       s2fac = (t - t1) * one / (t2 - t1)
       
diff --git a/Carpet/CarpetLib/src/prolongate_3d_real8_2tl_o3.F77 b/Carpet/CarpetLib/src/prolongate_3d_real8_2tl_o3.F77
index 0f0509cff..554c4fc38 100644
--- a/Carpet/CarpetLib/src/prolongate_3d_real8_2tl_o3.F77
+++ b/Carpet/CarpetLib/src/prolongate_3d_real8_2tl_o3.F77
@@ -1,20 +1,7 @@
 c     -*-Fortran-*-
-c     $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_2tl_o3.F77,v 1.13 2004/03/11 12:03:09 schnetter Exp $
+c     $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_2tl_o3.F77,v 1.1 2001/03/24 22:38:48 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_o3 (
      $     src1, t1, src2, t2, srciext, srcjext, srckext,
@@ -26,24 +13,19 @@ 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
       integer srcbbox(3,3), dstbbox(3,3), regbbox(3,3)
       
-      integer offsetlo, offsethi
-      
       integer regiext, regjext, regkext
       
       integer dstifac, dstjfac, dstkfac
@@ -53,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(4), jfac(4), kfac(4)
       integer ii, jj, kk
-      integer fac
+      CCTK_REAL8 fi, fj, fk
+      CCTK_REAL8 fac
       CCTK_REAL8 res
       integer d
       
-      character msg*1000
-      
       
       
       do d=1,3
@@ -80,41 +58,15 @@ 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
-         regkext = (regbbox(d,2) - regbbox(d,1)) / regbbox(d,3) + 1
-         dstkfac = srcbbox(d,3) / dstbbox(d,3)
-         srckoff = (regbbox(d,1) - srcbbox(d,1)) / dstbbox(d,3)
-         offsetlo = regbbox(d,3)
-         if (mod(srckoff + 0, dstkfac).eq.0) then
-            offsetlo = 0
-            if (regkext.gt.1) then
-               offsetlo = regbbox(d,3)
-            end if
-         end if
-         offsethi = regbbox(d,3)
-         if (mod(srckoff + regkext-1, dstkfac).eq.0) then
-            offsethi = 0
-            if (regkext.gt.1) then
-               offsethi = regbbox(d,3)
-            end if
-         end if
-         if (regbbox(d,1)-offsetlo.lt.srcbbox(d,1)
-     $        .or. regbbox(d,2)+offsethi.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 extent")
+         if (regbbox(d,1)-1.lt.srcbbox(d,1)
+     $        .or. regbbox(d,1)-1.lt.dstbbox(d,1)
+     $        .or. regbbox(d,2)+1.gt.srcbbox(d,2)
+     $        .or. regbbox(d,2)+1.gt.dstbbox(d,2)) then
+            call CCTK_WARN (0, "Internal error: region extent is not contained in array extents")
          end if
       end do
       
@@ -151,65 +103,74 @@ 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 / (6*dstifac**3 * 6*dstjfac**3 * 6*dstkfac**3)
-      
       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
-               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
+               
+               i0 = (srcioff + i) / dstifac + 1
+               j0 = (srcjoff + j) / dstjfac + 1
+               k0 = (srckoff + k) / dstkfac + 1
+               
+               fi = mod(srcioff + i, dstifac) * one / dstifac
+               fj = mod(srcjoff + j, dstjfac) * one / dstjfac
+               fk = mod(srckoff + k, dstkfac) * one / dstkfac
                
                res = 0
                
-               do kk=1,4
-                  do jj=1,4
-                     do ii=1,4
+               do kk=-1,2
+                  do jj=-1,2
+                     do ii=-1,2
                         
-                        fac = ifac(ii) * jfac(jj) * kfac(kk)
+                        fac = 1
+                        
+                        if (ii.eq.-1) then
+                           fac = fac * (fi  ) * (fi-1) * (fi-2) / (-6)
+                        else if (ii.eq.0) then
+                           fac = fac * (fi+1) * (fi-1) * (fi-2) / 2
+                        else if (ii.eq.1) then
+                           fac = fac * (fi+1) * (fi  ) * (fi-2) / (-2)
+                        else
+                           fac = fac * (fi+1) * (fi  ) * (fi-1) / 6
+                        end if
+                        
+                        if (jj.eq.-1) then
+                           fac = fac * (fj  ) * (fj-1) * (fj-2) / (-6)
+                        else if (jj.eq.0) then
+                           fac = fac * (fj+1) * (fj-1) * (fj-2) / 2
+                        else if (jj.eq.1) then
+                           fac = fac * (fj+1) * (fj  ) * (fj-2) / (-2)
+                        else
+                           fac = fac * (fj+1) * (fj  ) * (fj-1) / 6
+                        end if
+                        
+                        if (kk.eq.-1) then
+                           fac = fac * (fk  ) * (fk-1) * (fk-2) / (-6)
+                        else if (kk.eq.0) then
+                           fac = fac * (fk+1) * (fk-1) * (fk-2) / 2
+                        else if (kk.eq.1) then
+                           fac = fac * (fk+1) * (fk  ) * (fk-2) / (-2)
+                        else
+                           fac = fac * (fk+1) * (fk  ) * (fk-1) / 6
+                        end if
                         
                         if (fac.ne.0) then
-                           CHKIDX (i0+ii-1, j0+jj-1, k0+kk-1, \
-                                   srciext,srcjext,srckext, "source")
                            res = res
-     $                          + fac * s1fac * src1(i0+ii-1, j0+jj-1, k0+kk-1)
-     $                          + fac * s2fac * src2(i0+ii-1, j0+jj-1, k0+kk-1)
+     $                          + 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
                
             end do
          end do
diff --git a/Carpet/CarpetLib/src/prolongate_3d_real8_3tl.F77 b/Carpet/CarpetLib/src/prolongate_3d_real8_3tl.F77
index 9f3e3d944..17834ddd1 100644
--- a/Carpet/CarpetLib/src/prolongate_3d_real8_3tl.F77
+++ b/Carpet/CarpetLib/src/prolongate_3d_real8_3tl.F77
@@ -1,20 +1,7 @@
 c     -*-Fortran-*-
-c     $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_3tl.F77,v 1.9 2004/03/11 12:03:09 schnetter Exp $
+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 $
       
 #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_3tl (
      $     src1, t1, src2, t2, src3, t3, srciext, srcjext, srckext,
@@ -26,19 +13,16 @@ 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
       CCTK_REAL8 src3(srciext,srcjext,srckext)
-      CCTK_REAL8 t3
+      integer t3
       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
@@ -53,18 +37,14 @@ c     bbox(:,3) is stride
       
       CCTK_REAL8 s1fac, s2fac, s3fac
       
-      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
@@ -80,24 +60,15 @@ 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")
+            call CCTK_WARN (0, "Internal error: region extent is not contained in array extents")
          end if
       end do
       
@@ -134,61 +105,63 @@ c     Quadratic (second order) interpolation
       if (t1.eq.t2 .or. t1.eq.t3 .or. t2.eq.t3) then
          call CCTK_WARN (0, "Internal error: arrays have same time")
       end if
-      if (t.lt.min(t1,t2,t3)-eps .or. t.gt.max(t1,t2,t3)+eps) then
-         call CCTK_WARN (0, "Internal error: extrapolation in time")
-      end if
       
-      s1fac = (t - t2) * (t - t3) / ((t1 - t2) * (t1 - t3))
-      s2fac = (t - t1) * (t - t3) / ((t2 - t1) * (t2 - t3))
-      s3fac = (t - t1) * (t - t2) / ((t3 - t1) * (t3 - t2))
+      s1fac = (t - t2) * (t - t3) * one / ((t1 - t2) * (t1 - t3))
+      s2fac = (t - t1) * (t - t3) * one / ((t2 - t1) * (t2 - t3))
+      s3fac = (t - t1) * (t - t2) * one / ((t3 - t1) * (t3 - t2))
       
       
       
 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 = 1
                         
-                        fac = ifac(ii) * jfac(jj) * kfac(kk)
+                        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 * s3fac * src3(i0+ii, j0+jj, k0+kk)
+     $                          + fac * s1fac * src1(i0+ii,j0+jj,k0+kk)
+     $                          + fac * s2fac * src2(i0+ii,j0+jj,k0+kk)
+     $                          + fac * s3fac * src3(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/prolongate_3d_real8_3tl_o3.F77 b/Carpet/CarpetLib/src/prolongate_3d_real8_3tl_o3.F77
index c44e1119f..ad360647c 100644
--- a/Carpet/CarpetLib/src/prolongate_3d_real8_3tl_o3.F77
+++ b/Carpet/CarpetLib/src/prolongate_3d_real8_3tl_o3.F77
@@ -1,20 +1,7 @@
 c     -*-Fortran-*-
-c     $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_3tl_o3.F77,v 1.13 2004/03/11 12:03:09 schnetter Exp $
+c     $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_3tl_o3.F77,v 1.1 2001/03/24 22:38:48 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_3tl_o3 (
      $     src1, t1, src2, t2, src3, t3, srciext, srcjext, srckext,
@@ -26,26 +13,21 @@ 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
       CCTK_REAL8 src3(srciext,srcjext,srckext)
-      CCTK_REAL8 t3
+      integer t3
       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
       integer srcbbox(3,3), dstbbox(3,3), regbbox(3,3)
       
-      integer offsetlo, offsethi
-      
       integer regiext, regjext, regkext
       
       integer dstifac, dstjfac, dstkfac
@@ -55,18 +37,14 @@ c     bbox(:,3) is stride
       
       CCTK_REAL8 s1fac, s2fac, s3fac
       
-      CCTK_REAL8 dstdiv
       integer i, j, k
       integer i0, j0, k0
-      integer fi, fj, fk
-      integer ifac(4), jfac(4), kfac(4)
       integer ii, jj, kk
-      integer fac
+      CCTK_REAL8 fi, fj, fk
+      CCTK_REAL8 fac
       CCTK_REAL8 res
       integer d
       
-      character msg*1000
-      
       
       
       do d=1,3
@@ -82,41 +60,15 @@ 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
-         regkext = (regbbox(d,2) - regbbox(d,1)) / regbbox(d,3) + 1
-         dstkfac = srcbbox(d,3) / dstbbox(d,3)
-         srckoff = (regbbox(d,1) - srcbbox(d,1)) / dstbbox(d,3)
-         offsetlo = regbbox(d,3)
-         if (mod(srckoff + 0, dstkfac).eq.0) then
-            offsetlo = 0
-            if (regkext.gt.1) then
-               offsetlo = regbbox(d,3)
-            end if
-         end if
-         offsethi = regbbox(d,3)
-         if (mod(srckoff + regkext-1, dstkfac).eq.0) then
-            offsethi = 0
-            if (regkext.gt.1) then
-               offsethi = regbbox(d,3)
-            end if
-         end if
-         if (regbbox(d,1)-offsetlo.lt.srcbbox(d,1)
-     $        .or. regbbox(d,2)+offsethi.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 extent")
+         if (regbbox(d,1)-1.lt.srcbbox(d,1)
+     $        .or. regbbox(d,1)-1.lt.dstbbox(d,1)
+     $        .or. regbbox(d,2)+1.gt.srcbbox(d,2)
+     $        .or. regbbox(d,2)+1.gt.dstbbox(d,2)) then
+            call CCTK_WARN (0, "Internal error: region extent is not contained in array extents")
          end if
       end do
       
@@ -153,67 +105,76 @@ c     Quadratic (second order) interpolation
       if (t1.eq.t2 .or. t1.eq.t3 .or. t2.eq.t3) then
          call CCTK_WARN (0, "Internal error: arrays have same time")
       end if
-      if (t.lt.min(t1,t2,t3)-eps .or. t.gt.max(t1,t2,t3)+eps) then
-         call CCTK_WARN (0, "Internal error: extrapolation in time")
-      end if
       
-      s1fac = (t - t2) * (t - t3) / ((t1 - t2) * (t1 - t3))
-      s2fac = (t - t1) * (t - t3) / ((t2 - t1) * (t2 - t3))
-      s3fac = (t - t1) * (t - t2) / ((t3 - t1) * (t3 - t2))
+      s1fac = (t - t2) * (t - t3) * one / ((t1 - t2) * (t1 - t3))
+      s2fac = (t - t1) * (t - t3) * one / ((t2 - t1) * (t2 - t3))
+      s3fac = (t - t1) * (t - t2) * one / ((t3 - t1) * (t3 - t2))
       
       
       
 c     Loop over fine region
-      dstdiv = one / (6*dstifac**3 * 6*dstjfac**3 * 6*dstkfac**3)
-      
       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
-               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
+               
+               i0 = (srcioff + i) / dstifac + 1
+               j0 = (srcjoff + j) / dstjfac + 1
+               k0 = (srckoff + k) / dstkfac + 1
+               
+               fi = mod(srcioff + i, dstifac) * one / dstifac
+               fj = mod(srcjoff + j, dstjfac) * one / dstjfac
+               fk = mod(srckoff + k, dstkfac) * one / dstkfac
                
                res = 0
                
-               do kk=1,4
-                  do jj=1,4
-                     do ii=1,4
+               do kk=-1,2
+                  do jj=-1,2
+                     do ii=-1,2
                         
-                        fac = ifac(ii) * jfac(jj) * kfac(kk)
+                        fac = 1
+                        
+                        if (ii.eq.-1) then
+                           fac = fac * (fi  ) * (fi-1) * (fi-2) / (-6)
+                        else if (ii.eq.0) then
+                           fac = fac * (fi+1) * (fi-1) * (fi-2) / 2
+                        else if (ii.eq.1) then
+                           fac = fac * (fi+1) * (fi  ) * (fi-2) / (-2)
+                        else
+                           fac = fac * (fi+1) * (fi  ) * (fi-1) / 6
+                        end if
+                        
+                        if (jj.eq.-1) then
+                           fac = fac * (fj  ) * (fj-1) * (fj-2) / (-6)
+                        else if (jj.eq.0) then
+                           fac = fac * (fj+1) * (fj-1) * (fj-2) / 2
+                        else if (jj.eq.1) then
+                           fac = fac * (fj+1) * (fj  ) * (fj-2) / (-2)
+                        else
+                           fac = fac * (fj+1) * (fj  ) * (fj-1) / 6
+                        end if
+                        
+                        if (kk.eq.-1) then
+                           fac = fac * (fk  ) * (fk-1) * (fk-2) / (-6)
+                        else if (kk.eq.0) then
+                           fac = fac * (fk+1) * (fk-1) * (fk-2) / 2
+                        else if (kk.eq.1) then
+                           fac = fac * (fk+1) * (fk  ) * (fk-2) / (-2)
+                        else
+                           fac = fac * (fk+1) * (fk  ) * (fk-1) / 6
+                        end if
                         
                         if (fac.ne.0) then
-                           CHKIDX (i0+ii-1, j0+jj-1, k0+kk-1, \
-                                   srciext,srcjext,srckext, "source")
                            res = res
-     $                          + fac * s1fac * src1(i0+ii-1, j0+jj-1, k0+kk-1)
-     $                          + fac * s2fac * src2(i0+ii-1, j0+jj-1, k0+kk-1)
-     $                          + fac * s3fac * src3(i0+ii-1, j0+jj-1, k0+kk-1)
+     $                          + fac * s1fac * src1(i0+ii,j0+jj,k0+kk)
+     $                          + fac * s2fac * src2(i0+ii,j0+jj,k0+kk)
+     $                          + fac * s3fac * src3(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
                
             end do
          end do
diff --git a/Carpet/CarpetLib/src/prolongate_3d_real8_o3.F77 b/Carpet/CarpetLib/src/prolongate_3d_real8_o3.F77
index 21d82a733..c539ee7de 100644
--- a/Carpet/CarpetLib/src/prolongate_3d_real8_o3.F77
+++ b/Carpet/CarpetLib/src/prolongate_3d_real8_o3.F77
@@ -1,20 +1,7 @@
 c     -*-Fortran-*-
-c     $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_o3.F77,v 1.9 2004/03/11 12:03:09 schnetter Exp $
+c     $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_o3.F77,v 1.1 2001/03/24 22:38:48 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_o3 (
      $     src, srciext, srcjext, srckext,
@@ -35,8 +22,6 @@ c     bbox(:,2) is upper boundary (inclusive)
 c     bbox(:,3) is stride
       integer srcbbox(3,3), dstbbox(3,3), regbbox(3,3)
       
-      integer offsetlo, offsethi
-      
       integer regiext, regjext, regkext
       
       integer dstifac, dstjfac, dstkfac
@@ -44,18 +29,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(4), jfac(4), kfac(4)
       integer ii, jj, kk
-      integer fac
+      CCTK_REAL8 fi, fj, fk
+      CCTK_REAL8 fac
       CCTK_REAL8 res
       integer d
       
-      character msg*1000
-      
       
       
       do d=1,3
@@ -71,41 +52,15 @@ 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
-         regkext = (regbbox(d,2) - regbbox(d,1)) / regbbox(d,3) + 1
-         dstkfac = srcbbox(d,3) / dstbbox(d,3)
-         srckoff = (regbbox(d,1) - srcbbox(d,1)) / dstbbox(d,3)
-         offsetlo = regbbox(d,3)
-         if (mod(srckoff + 0, dstkfac).eq.0) then
-            offsetlo = 0
-            if (regkext.gt.1) then
-               offsetlo = regbbox(d,3)
-            end if
-         end if
-         offsethi = regbbox(d,3)
-         if (mod(srckoff + regkext-1, dstkfac).eq.0) then
-            offsethi = 0
-            if (regkext.gt.1) then
-               offsethi = regbbox(d,3)
-            end if
-         end if
-         if (regbbox(d,1)-offsetlo.lt.srcbbox(d,1)
-     $        .or. regbbox(d,2)+offsethi.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 extent")
+         if (regbbox(d,1)-1.lt.srcbbox(d,1)
+     $        .or. regbbox(d,1)-1.lt.dstbbox(d,1)
+     $        .or. regbbox(d,2)+1.gt.srcbbox(d,2)
+     $        .or. regbbox(d,2)+1.gt.dstbbox(d,2)) then
+            call CCTK_WARN (0, "Internal error: region extent is not contained in array extents")
          end if
       end do
       
@@ -139,53 +94,65 @@ c     This could be handled, but is likely to point to an error elsewhere
       
       
 c     Loop over fine region
-      dstdiv = one / (6*dstifac**3 * 6*dstjfac**3 * 6*dstkfac**3)
-      
       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
-               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
+               
+               i0 = (srcioff + i) / dstifac + 1
+               j0 = (srcjoff + j) / dstjfac + 1
+               k0 = (srckoff + k) / dstkfac + 1
+               
+               fi = mod(srcioff + i, dstifac) * one / dstifac
+               fj = mod(srcjoff + j, dstjfac) * one / dstjfac
+               fk = mod(srckoff + k, dstkfac) * one / dstkfac
                
                res = 0
                
-               do kk=1,4
-                  do jj=1,4
-                     do ii=1,4
+               do kk=-1,2
+                  do jj=-1,2
+                     do ii=-1,2
+                        
+                        fac = 1
+                        
+                        if (ii.eq.-1) then
+                           fac = fac * (fi  ) * (fi-1) * (fi-2) / (-6)
+                        else if (ii.eq.0) then
+                           fac = fac * (fi+1) * (fi-1) * (fi-2) / 2
+                        else if (ii.eq.1) then
+                           fac = fac * (fi+1) * (fi  ) * (fi-2) / (-2)
+                        else
+                           fac = fac * (fi+1) * (fi  ) * (fi-1) / 6
+                        end if
+                        
+                        if (jj.eq.-1) then
+                           fac = fac * (fj  ) * (fj-1) * (fj-2) / (-6)
+                        else if (jj.eq.0) then
+                           fac = fac * (fj+1) * (fj-1) * (fj-2) / 2
+                        else if (jj.eq.1) then
+                           fac = fac * (fj+1) * (fj  ) * (fj-2) / (-2)
+                        else
+                           fac = fac * (fj+1) * (fj  ) * (fj-1) / 6
+                        end if
                         
-                        fac = ifac(ii) * jfac(jj) * kfac(kk)
+                        if (kk.eq.-1) then
+                           fac = fac * (fk  ) * (fk-1) * (fk-2) / (-6)
+                        else if (kk.eq.0) then
+                           fac = fac * (fk+1) * (fk-1) * (fk-2) / 2
+                        else if (kk.eq.1) then
+                           fac = fac * (fk+1) * (fk  ) * (fk-2) / (-2)
+                        else
+                           fac = fac * (fk+1) * (fk  ) * (fk-1) / 6
+                        end if
                         
                         if (fac.ne.0) then
-                           CHKIDX (i0+ii-1, j0+jj-1, k0+kk-1, \
-                                   srciext,srcjext,srckext, "source")
-                           res = res + fac * src(i0+ii-1, j0+jj-1, k0+kk-1)
+                           res = res + fac * src(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
                
             end do
          end do
diff --git a/Carpet/CarpetLib/src/restrict_3d_real8.F77 b/Carpet/CarpetLib/src/restrict_3d_real8.F77
index 68ea98b11..d580f423c 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.2 2001/03/22 18:42:06 eschnett Exp $
+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 $
       
 #include "cctk.h"
       
@@ -48,6 +48,12 @@ c     bbox(:,3) is stride
      $        .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).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 extents")
+         end if
       end do
       
       if (srciext.ne.(srcbbox(1,2)-srcbbox(1,1))/srcbbox(1,3)+1