Handle empty bboxes.

Handle empty bboxes.
*.F77: Better error checking whether the active region is contained in
the source and destination arrays.

*.F77: Temporarily activated per-gridpoint checking of array accesses.

bbox.cc bboxset.cc: Handle empty bboxes better -- either handle them
correctly, or abort.

gdata.cc: Recognise empty regions.

ggf.cc: Remove line that was commented out for a long time.

dh.cc: Choose send and recv regions for syncing and prolongation so
that they don't overlap.

dh.cc: Handle empty bboxes.

darcs-hash:20030225215700-07bb3-a7296dd92353c003bc0bd3ff435e4939f8041eae.gz

13 files changed, 484 insertions, 161 deletions

diff --git a/Carpet/CarpetLib/src/bbox.cc b/Carpet/CarpetLib/src/bbox.cc
index 0369aec9e..dcae93c02 100644
--- a/Carpet/CarpetLib/src/bbox.cc
+++ b/Carpet/CarpetLib/src/bbox.cc
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/bbox.cc,v 1.12 2003/02/24 17:11:29 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/bbox.cc,v 1.13 2003/02/25 22:57:00 schnetter Exp $
 
 #include <assert.h>
 
@@ -117,6 +117,7 @@ bbox<T,D> bbox<T,D>::operator& (const bbox& b) const {
 // Containment
 template<class T, int D>
 bool bbox<T,D>::is_contained_in (const bbox& b) const {
+  if (empty()) return true;
   // no alignment check
   return all(lower()>=b.lower() && upper()<=b.upper());
 }
@@ -140,7 +141,7 @@ bbox<T,D> bbox<T,D>::expand (const vect<T,D>& lo, const vect<T,D>& hi) const {
 // Find the smallest b-compatible box around *this
 template<class T, int D>
 bbox<T,D> bbox<T,D>::expanded_for (const bbox& b) const {
-  assert (! empty());
+  if (empty()) return bbox(b.lower(), b.lower()-b.stride(), b.stride());
   const vect<T,D> str = b.stride();
   const vect<T,D> loff = ((lower() - b.lower()) % str + str) % str;
   const vect<T,D> uoff = ((upper() - b.lower()) % str + str) % str;
@@ -152,7 +153,7 @@ bbox<T,D> bbox<T,D>::expanded_for (const bbox& b) const {
 // Find the largest b-compatible box inside *this
 template<class T, int D>
 bbox<T,D> bbox<T,D>::contracted_for (const bbox& b) const {
-  assert (! empty());
+  if (empty()) return bbox(b.lower(), b.lower()-b.stride(), b.stride());
   const vect<T,D> str = b.stride();
   const vect<T,D> loff = ((lower() - b.lower()) % str + str) % str;
   const vect<T,D> uoff = ((upper() - b.lower()) % str + str) % str;
diff --git a/Carpet/CarpetLib/src/bboxset.cc b/Carpet/CarpetLib/src/bboxset.cc
index 8de16a015..b440bb59d 100644
--- a/Carpet/CarpetLib/src/bboxset.cc
+++ b/Carpet/CarpetLib/src/bboxset.cc
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/bboxset.cc,v 1.11 2003/01/03 15:49:36 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/bboxset.cc,v 1.12 2003/02/25 22:57:00 schnetter Exp $
 
 #include <assert.h>
 
@@ -318,7 +318,9 @@ bool bboxset<T,D>::operator!= (const bboxset<T,D>& s) const {
 // Output
 template<class T,int D>
 void bboxset<T,D>::output (ostream& os) const {
-  os << "bboxset<T," << D << ">:size=" << size() << "," << "set=" << bs;
+  T Tdummy;
+  os << "bboxset<" << typestring(Tdummy) << "," << D << ">:"
+     << "size=" << size() << "," << "set=" << bs;
 }
 
 
diff --git a/Carpet/CarpetLib/src/copy_3d_real8.F77 b/Carpet/CarpetLib/src/copy_3d_real8.F77
index 230814c04..d115418b8 100644
--- a/Carpet/CarpetLib/src/copy_3d_real8.F77
+++ b/Carpet/CarpetLib/src/copy_3d_real8.F77
@@ -1,7 +1,20 @@
 c     -*-Fortran-*-
-c     $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/copy_3d_real8.F77,v 1.5 2003/02/24 17:43:10 schnetter Exp $
+c     $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/copy_3d_real8.F77,v 1.6 2003/02/25 22:57:00 schnetter 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 copy_3d_real8 (
      $     src, srciext, srcjext, srckext,
@@ -27,6 +40,8 @@ c     bbox(:,3) is stride
       integer i, j, k
       integer d
       
+      character msg*1000
+      
       
       
       do d=1,3
@@ -43,16 +58,16 @@ 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).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 (regbbox(d,1).lt.srcbbox(d,1)
      $        .or. regbbox(d,1).lt.dstbbox(d,1)
      $        .or. regbbox(d,2).gt.srcbbox(d,2)
      $        .or. regbbox(d,2).gt.dstbbox(d,2)) then
             call CCTK_WARN (0, "Internal error: region extent is not contained in array extent")
          end if
-         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
       end do
       
       if (srciext.ne.(srcbbox(1,2)-srcbbox(1,1))/srcbbox(1,3)+1
@@ -85,6 +100,10 @@ c     Loop over region
          do j = 0, regjext-1
             do i = 0, regiext-1
                
+               CHKIDX (srcioff+i+1, srcjoff+j+1, srckoff+k+1, \
+                       srciext,srcjext,srckext, "source")
+               CHKIDX (dstioff+i+1, dstjoff+j+1, dstkoff+k+1, \
+                       dstiext,dstjext,dstkext, "destination")
                dst (dstioff+i+1, dstjoff+j+1, dstkoff+k+1)
      $              = src (srcioff+i+1, srcjoff+j+1, srckoff+k+1)
                
diff --git a/Carpet/CarpetLib/src/dh.cc b/Carpet/CarpetLib/src/dh.cc
index 68f4c18c4..7b856565a 100644
--- a/Carpet/CarpetLib/src/dh.cc
+++ b/Carpet/CarpetLib/src/dh.cc
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/dh.cc,v 1.25 2003/01/10 23:07:16 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/dh.cc,v 1.26 2003/02/25 22:57:00 schnetter Exp $
 
 #include <assert.h>
 
@@ -76,6 +76,7 @@ void dh<D>::recompose () {
 	  if (h.outer_boundaries[rl][c][d][0]) ldist[d] = 0;
 	  if (h.outer_boundaries[rl][c][d][1]) udist[d] = 0;
 	}
+        assert (! intr.empty());
        	boxes[rl][c][ml].exterior = intr.expand(ldist, udist);
 	
        	// Boundaries (ghost zones only)
@@ -116,10 +117,7 @@ void dh<D>::recompose () {
   for (int rl=0; rl<h.reflevels(); ++rl) {
     for (int c=0; c<h.components(rl); ++c) {
       for (int ml=0; ml<h.mglevels(rl,c); ++ml) {
-      	const ibbox intr = boxes[rl][c][ml].interior;
-      	const ibbox extr = boxes[rl][c][ml].exterior;
-
-      	const ibset bnds = boxes[rl][c][ml].boundaries;
+      	const ibset& bnds = boxes[rl][c][ml].boundaries;
 	
       	// Sync boxes
       	for (int cc=0; cc<h.components(rl); ++cc) {
@@ -133,6 +131,16 @@ void dh<D>::recompose () {
 	  }
       	} // for cc
       	
+      } // for ml
+    } // for c
+  } // for rl
+  
+  for (int rl=0; rl<h.reflevels(); ++rl) {
+    for (int c=0; c<h.components(rl); ++c) {
+      for (int ml=0; ml<h.mglevels(rl,c); ++ml) {
+      	const ibbox& intr = boxes[rl][c][ml].interior;
+      	const ibbox& extr = boxes[rl][c][ml].exterior;
+	
       	// Multigrid boxes
       	if (ml>0) {
       	  const ibbox intrf = boxes[rl][c][ml-1].interior;
@@ -142,7 +150,9 @@ void dh<D>::recompose () {
       	    // (the restriction must fill all of the interior of the
       	    // coarse grid, and may use the exterior of the fine grid)
       	    const ibbox recv = intr;
+            assert (! recv.empty());
       	    const ibbox send = recv.expanded_for(extrf);
+            assert (! send.empty());
       	    assert (send.is_contained_in(extrf));
       	    boxes[rl][c][ml-1].send_mg_coarse.push_back(send);
       	    boxes[rl][c][ml  ].recv_mg_fine  .push_back(recv);
@@ -153,30 +163,57 @@ void dh<D>::recompose () {
       	    // grid, and may fill only the interior of the fine grid,
       	    // and the bbox must be as large as possible)
       	    const ibbox recv = extr.contracted_for(intrf) & intrf;
+            assert (! recv.empty());
       	    const ibbox send = recv.expanded_for(extr);
+            assert (! send.empty());
       	    boxes[rl][c][ml-1].recv_mg_coarse.push_back(recv);
       	    boxes[rl][c][ml  ].send_mg_fine  .push_back(send);
       	  }
       	} // if not finest multigrid level
-      	
+        
+      } // for ml
+    } // for c
+  } // for rl
+  
+  for (int rl=0; rl<h.reflevels(); ++rl) {
+    for (int c=0; c<h.components(rl); ++c) {
+      for (int ml=0; ml<h.mglevels(rl,c); ++ml) {
+      	const ibbox& intr = boxes[rl][c][ml].interior;
+      	const ibbox& extr = boxes[rl][c][ml].exterior;
+	
       	// Refinement boxes
       	if (rl<h.reflevels()-1) {
       	  for (int cc=0; cc<h.components(rl+1); ++cc) {
       	    const ibbox intrf = boxes[rl+1][cc][ml].interior;
-// 	    const ibbox extrf = boxes[rl+1][cc][ml].exterior;
       	    // Prolongation (interior)
       	    {
       	      // (the prolongation may use the exterior of the coarse
       	      // grid, and must fill all of the interior of the fine
       	      // grid)
-      	      const ibbox recv = extr.contracted_for(intrf) & intrf;
-      	      const ibbox send = recv.expanded_for(extr);
-      	      assert (send.is_contained_in(extr));
-      	      boxes[rl+1][cc][ml].recv_ref_coarse[c ].push_back(recv);
-      	      boxes[rl  ][c ][ml].send_ref_fine  [cc].push_back(send);
+              const int pss = prolongation_stencil_size();
+              ibset recvs = (extr.expand(-pss,-pss).contracted_for(intrf)
+                             & intrf);
+              const iblistvect& rrc = boxes[rl+1][cc][ml].recv_ref_coarse;
+              for (typename iblistvect::const_iterator lvi=rrc.begin();
+                   lvi!=rrc.end(); ++lvi) {
+                for (typename iblist::const_iterator li=lvi->begin();
+                     li!=lvi->end(); ++li) {
+                  recvs -= *li;
+                }
+              }
+              assert (recvs.setsize() <= 1);
+              if (recvs.setsize() == 1) {
+                const ibbox recv = *recvs.begin();
+                const ibbox send = recv.expanded_for(extr);
+                assert (! send.empty());
+                assert (send.is_contained_in(extr));
+                boxes[rl+1][cc][ml].recv_ref_coarse[c ].push_back(recv);
+                boxes[rl  ][c ][ml].send_ref_fine  [cc].push_back(send);
+              }
       	    }
       	    // Prolongation (boundaries)
       	    {
+              const int pss = prolongation_stencil_size();
       	      ibset bndsf = boxes[rl+1][cc][ml].boundaries;
       	      // coarsify boundaries of fine component
       	      for (typename ibset::const_iterator bi=bndsf.begin();
@@ -185,13 +222,34 @@ void dh<D>::recompose () {
 		// (the prolongation may use the exterior of the
 		// coarse grid, and must fill all of the boundary of
 		// the fine grid)
-		const int pss = prolongation_stencil_size();
-		const ibbox recv = (extr.expand(-pss,-pss).contracted_for(bndf)
-				    & bndf);
-		const ibbox send = recv.expanded_for(extr);
-		assert (send.is_contained_in(extr));
-		boxes[rl+1][cc][ml].recv_ref_bnd_coarse[c ].push_back(recv);
-		boxes[rl  ][c ][ml].send_ref_bnd_fine  [cc].push_back(send);
+                ibset recvs = (extr.expand(-pss,-pss).contracted_for(bndf)
+                               & bndf);
+                const iblistvect& rrbc
+                  = boxes[rl+1][cc][ml].recv_ref_bnd_coarse;
+                for (typename iblistvect::const_iterator lvi=rrbc.begin();
+                     lvi!=rrbc.end(); ++lvi) {
+                  for (typename iblist::const_iterator li=lvi->begin();
+                       li!=lvi->end(); ++li) {
+                    recvs -= *li;
+                  }
+                }
+                const iblistvect& rs = boxes[rl+1][cc][ml].recv_sync;
+                for (typename iblistvect::const_iterator lvi=rs.begin();
+                     lvi!=rs.end(); ++lvi) {
+                  for (typename iblist::const_iterator li=lvi->begin();
+                       li!=lvi->end(); ++li) {
+                    recvs -= *li;
+                  }
+                }
+                assert (recvs.setsize() <= 1);
+                if (recvs.setsize() == 1) {
+                  const ibbox recv = *recvs.begin();
+                  const ibbox send = recv.expanded_for(extr);
+                  assert (! send.empty());
+                  assert (send.is_contained_in(extr));
+                  boxes[rl+1][cc][ml].recv_ref_bnd_coarse[c ].push_back(recv);
+                  boxes[rl  ][c ][ml].send_ref_bnd_fine  [cc].push_back(send);
+                }
       	      }
       	    }
       	    // Restriction (interior)
@@ -201,23 +259,25 @@ void dh<D>::recompose () {
       	      // grid, and the bbox must be as large as possible)
 	      // (the restriction must not fill points that are used
 	      // to prolongate the boundaries)
-      	      const ibbox recv = intrf.contracted_for(intr) & intr;
-	      
-	      const iblist& sendlist = boxes[rl][c][ml].send_ref_bnd_fine[cc];
-	      ibset recv2s (recv);
-	      for (typename iblist::const_iterator sli = sendlist.begin();
-		   sli != sendlist.end();
-		   ++sli) {
-		recv2s -= *sli;
-	      }
-	      assert (recv2s.setsize() == 1);
-	      const ibbox recv2 = *recv2s.begin();
-//       	      const ibbox send = recv.expanded_for(intrf);
-      	      const ibbox send = recv2.expanded_for(intrf);
-      	      boxes[rl+1][cc][ml].send_ref_coarse[c ].push_back(send);
-      	      boxes[rl  ][c ][ml].recv_ref_fine  [cc].push_back(recv);
+              ibset recvs = intrf.contracted_for(intr) & intr;
+              const iblist& sendlist
+                = boxes[rl][c][ml].send_ref_bnd_fine[cc];
+              for (typename iblist::const_iterator sli = sendlist.begin();
+                   sli != sendlist.end();
+                   ++sli) {
+                recvs -= *sli;
+              }
+              assert (recvs.setsize() <= 1);
+              if (recvs.setsize() == 1) {
+                const ibbox recv = *recvs.begin();
+                assert (! recv.empty());
+                const ibbox send = recv.expanded_for(intrf);
+                assert (! send.empty());
+                boxes[rl+1][cc][ml].send_ref_coarse[c ].push_back(send);
+                boxes[rl  ][c ][ml].recv_ref_fine  [cc].push_back(recv);
+              }
       	    }
-      	
+            
       	  } // for cc
       	} // if not finest refinement level
 	
@@ -231,44 +291,41 @@ void dh<D>::recompose () {
 	
 	// Boundaries that are not synced, or are neither synced nor
 	// prolonged to from coarser grids (outer boundaries)
-	ibset& sync_not = boxes[rl][c][ml].sync_not;
-	ibset& recv_not = boxes[rl][c][ml].recv_not;
-	
-	// The whole boundary
-	sync_not = boxes[rl][c][ml].boundaries;
-	recv_not = boxes[rl][c][ml].boundaries;
-	
-	// Subtract boxes received during synchronisation
-	const iblistvect& recv_sync = boxes[rl][c][ml].recv_sync;
-	for (typename iblistvect::const_iterator lvi=recv_sync.begin();
-	     lvi!=recv_sync.end(); ++lvi) {
-	  for (typename iblist::const_iterator li=lvi->begin();
-	       li!=lvi->end(); ++li) {
-	    sync_not -= *li;
-	    recv_not -= *li;
-	  }
-	}
-	
-	// Subtract all boxes received
-	const iblistvect& recv_ref_bnd_coarse
-	  = boxes[rl][c][ml].recv_ref_bnd_coarse;
-	for (typename iblistvect::const_iterator
-	       lvi=recv_ref_bnd_coarse.begin();
-	     lvi!=recv_ref_bnd_coarse.end(); ++lvi) {
-	  for (typename iblist::const_iterator li=lvi->begin();
-	       li!=lvi->end(); ++li) {
-	    recv_not -= *li;
-	  }
-	}
-	
-	// Check that no boundaries are left over
-	if (rl==0) assert (sync_not.empty());
-	assert (recv_not.empty());
+        ibset& sync_not = boxes[rl][c][ml].sync_not;
+        ibset& recv_not = boxes[rl][c][ml].recv_not;
+        
+        // The whole boundary
+        sync_not = boxes[rl][c][ml].boundaries;
+        recv_not = boxes[rl][c][ml].boundaries;
+        
+        // Subtract boxes received during synchronisation
+        const iblistvect& recv_sync = boxes[rl][c][ml].recv_sync;
+        for (typename iblistvect::const_iterator lvi=recv_sync.begin();
+             lvi!=recv_sync.end(); ++lvi) {
+          for (typename iblist::const_iterator li=lvi->begin();
+               li!=lvi->end(); ++li) {
+            sync_not -= *li;
+            recv_not -= *li;
+          }
+        }
+        
+        // Subtract boxes received during prolongation
+        const iblistvect& recv_ref_bnd_coarse
+          = boxes[rl][c][ml].recv_ref_bnd_coarse;
+        for (typename iblistvect::const_iterator
+               lvi=recv_ref_bnd_coarse.begin();
+             lvi!=recv_ref_bnd_coarse.end(); ++lvi) {
+          for (typename iblist::const_iterator li=lvi->begin();
+               li!=lvi->end(); ++li) {
+            recv_not -= *li;
+          }
+        }
 	
       } // for ml
     } // for c
   } // for rl
   
+  // Calculate bases
   bases.resize(h.reflevels());
   for (int rl=0; rl<h.reflevels(); ++rl) {
     if (h.components(rl)==0) {
@@ -334,6 +391,84 @@ void dh<D>::recompose () {
     }
   } // if output_bboxes
   
+  for (int rl=0; rl<h.reflevels(); ++rl) {
+    for (int c=0; c<h.components(rl); ++c) {
+      for (int ml=0; ml<h.mglevels(rl,c); ++ml) {
+	
+	// Assert that all boundaries are synced or received
+        {
+          const ibset& sync_not = boxes[rl][c][ml].sync_not;
+          const ibset& recv_not = boxes[rl][c][ml].recv_not;
+          
+          // Check that no boundaries are left over
+          if (rl==0) assert (sync_not.empty());
+          assert (recv_not.empty());
+        }
+	
+        // Assert that the interior is received exactly once during
+        // prolongation, and that nothing else is received
+        {
+          if (rl==0) {
+            const iblistvect& recv_ref_coarse
+              = boxes[rl][c][ml].recv_ref_coarse;
+            assert (recv_ref_coarse.empty());
+          } else {              // rl!=0
+            const iblistvect& recv_ref_coarse
+              = boxes[rl][c][ml].recv_ref_coarse;
+            ibset intr = boxes[rl][c][ml].interior;
+            for (typename iblistvect::const_iterator
+                   lvi=recv_ref_coarse.begin();
+                 lvi!=recv_ref_coarse.end(); ++lvi) {
+              for (typename iblist::const_iterator li=lvi->begin();
+                   li!=lvi->end(); ++li) {
+                const int old_sz = intr.size();
+                const int this_sz = li->size();
+                intr -= *li;
+                const int new_sz = intr.size();
+                assert (new_sz + this_sz == old_sz);
+              }
+            }
+            assert (intr.empty());
+          }
+        }
+        
+        // Assert that the boundaries are received at most once during
+        // prolongation and synchronisation, and that nothing else is
+        // received
+        {
+          const iblistvect& recv_sync = boxes[rl][c][ml].recv_sync;
+          const iblistvect& recv_ref_bnd_coarse
+            = boxes[rl][c][ml].recv_ref_bnd_coarse;
+          ibset bnds = boxes[rl][c][ml].boundaries;
+          for (typename iblistvect::const_iterator lvi=recv_sync.begin();
+               lvi!=recv_sync.end(); ++lvi) {
+            for (typename iblist::const_iterator li=lvi->begin();
+                 li!=lvi->end(); ++li) {
+              const int old_sz = bnds.size();
+              const int this_sz = li->size();
+              bnds -= *li;
+              const int new_sz = bnds.size();
+              assert (new_sz + this_sz == old_sz);
+            }
+          }
+          for (typename iblistvect::const_iterator
+                 lvi=recv_ref_bnd_coarse.begin();
+               lvi!=recv_ref_bnd_coarse.end(); ++lvi) {
+            for (typename iblist::const_iterator li=lvi->begin();
+                 li!=lvi->end(); ++li) {
+              const int old_sz = bnds.size();
+              const int this_sz = li->size();
+              bnds -= *li;
+              const int new_sz = bnds.size();
+              assert (new_sz + this_sz == old_sz);
+            }
+          }
+        }
+        
+      } // for ml
+    } // for c
+  } // for rl
+  
   for (typename list<ggf<D>*>::iterator f=gfs.begin();
        f!=gfs.end(); ++f) {
     (*f)->recompose();
diff --git a/Carpet/CarpetLib/src/gdata.cc b/Carpet/CarpetLib/src/gdata.cc
index 86f5f89c3..f17534838 100644
--- a/Carpet/CarpetLib/src/gdata.cc
+++ b/Carpet/CarpetLib/src/gdata.cc
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gdata.cc,v 1.21 2003/01/03 15:49:36 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/gdata.cc,v 1.22 2003/02/25 22:57:00 schnetter Exp $
 
 #include <assert.h>
 
@@ -43,6 +43,8 @@ void gdata<D>::copy_from (const gdata* src, const ibbox& box)
   assert (all((box.lower()-extent().lower())%box.stride() == 0
 	      && (box.lower()-src->extent().lower())%box.stride() == 0));
   
+  if (box.empty()) return;
+  
   if (proc() == src->proc()) {
     // copy on same processor
     
@@ -87,6 +89,9 @@ void gdata<D>
     assert (all(box.upper()<=srcs[t]->extent().upper()));
   }
   
+  assert (! box.empty());
+  if (box.empty()) return;
+  
   if (proc() == srcs[0]->proc()) {
     // interpolate on same processor
     
diff --git a/Carpet/CarpetLib/src/ggf.cc b/Carpet/CarpetLib/src/ggf.cc
index fcc9e546b..79a7b105f 100644
--- a/Carpet/CarpetLib/src/ggf.cc
+++ b/Carpet/CarpetLib/src/ggf.cc
@@ -1,4 +1,4 @@
-// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/ggf.cc,v 1.20 2003/01/03 15:49:36 schnetter Exp $
+// $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/ggf.cc,v 1.21 2003/02/25 22:57:00 schnetter Exp $
 
 #include <assert.h>
 #include <stdlib.h>
@@ -32,8 +32,6 @@ ggf<D>::ggf (const string name, th<D>& t, dh<D>& d,
   assert (t.h == &d.h);
 
   d.add(this);
-
-//   recompose();
 }
 
 // Destructors
diff --git a/Carpet/CarpetLib/src/prolongate_3d_real8.F77 b/Carpet/CarpetLib/src/prolongate_3d_real8.F77
index 5ac6d8ce9..2dd5768e6 100644
--- a/Carpet/CarpetLib/src/prolongate_3d_real8.F77
+++ b/Carpet/CarpetLib/src/prolongate_3d_real8.F77
@@ -1,8 +1,21 @@
 c     -*-Fortran-*-
-c     $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8.F77,v 1.9 2003/02/24 17:43:10 schnetter Exp $
+c     $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8.F77,v 1.10 2003/02/25 22:57:00 schnetter Exp $
       
 #include "cctk.h"
-      
+
+
+
+#define CHKIDX(i,j,k, imax,jmax,kmax, where)                                  \
+      if ((i).lt.1 .or. (i).gt.(imax)                                         \
+           .or. (j).lt.1 .or. (j).gt.(jmax)                                   \
+           .or. (k).lt.1 .or. (k).gt.(kmax)) then                           &&\
+         write (msg, '(a, " array index out of bounds: shape is (",i4,",",i4,",",i4,"), index is (",i4,",",i4,",",i4,")")') \
+              (where), (imax), (jmax), (kmax), (i), (j), (k)                &&\
+         call CCTK_WARN (0, msg(1:len_trim(msg)))                           &&\
+      end if
+
+
+
       subroutine prolongate_3d_real8 (
      $     src, srciext, srcjext, srckext,
      $     dst, dstiext, dstjext, dstkext,
@@ -34,11 +47,13 @@ c     bbox(:,3) is stride
       integer i0, j0, k0
       integer fi, fj, fk
       integer ifac(2), jfac(2), kfac(2)
-c$$$      integer ii, jj, kk
+      integer ii, jj, kk
       integer fac
       CCTK_REAL8 res
       integer d
       
+      character msg*1000
+      
       
       
       do d=1,3
@@ -58,16 +73,16 @@ c$$$      integer ii, jj, kk
      $        .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 (regbbox(d,1).lt.srcbbox(d,1)
      $        .or. regbbox(d,1).lt.dstbbox(d,1)
      $        .or. regbbox(d,2).gt.srcbbox(d,2)
      $        .or. regbbox(d,2).gt.dstbbox(d,2)) then
             call CCTK_WARN (0, "Internal error: region extent is not contained in array extent")
          end if
-         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
       end do
       
       if (srciext.ne.(srcbbox(1,2)-srcbbox(1,1))/srcbbox(1,3)+1
@@ -122,44 +137,48 @@ c     Loop over fine region
                
                res = 0
                
-c$$$               do kk=1,2
-c$$$                  do jj=1,2
-c$$$                     do ii=1,2
-c$$$                        
-c$$$                        fac = ifac(ii) * jfac(jj) * kfac(kk)
-c$$$                        
-c$$$                        if (fac.ne.0) then
-c$$$                           res = res + fac * src(i0+ii, j0+jj, k0+kk)
-c$$$                        end if
-c$$$                        
-c$$$                     end do
-c$$$                  end do
-c$$$               end do
-               
-               fac = ifac(1) * jfac(1) * kfac(1)
-               if (fac.ne.0) res = res + fac * src(i0+1, j0+1, k0+1)
-               
-               fac = ifac(2) * jfac(1) * kfac(1)
-               if (fac.ne.0) res = res + fac * src(i0+2, j0+1, k0+1)
-               
-               fac = ifac(1) * jfac(2) * kfac(1)
-               if (fac.ne.0) res = res + fac * src(i0+1, j0+2, k0+1)
-               
-               fac = ifac(2) * jfac(2) * kfac(1)
-               if (fac.ne.0) res = res + fac * src(i0+2, j0+2, k0+1)
-               
-               fac = ifac(1) * jfac(1) * kfac(2)
-               if (fac.ne.0) res = res + fac * src(i0+1, j0+1, k0+2)
-               
-               fac = ifac(2) * jfac(1) * kfac(2)
-               if (fac.ne.0) res = res + fac * src(i0+2, j0+1, k0+2)
-               
-               fac = ifac(1) * jfac(2) * kfac(2)
-               if (fac.ne.0) res = res + fac * src(i0+1, j0+2, k0+2)
+               do kk=1,2
+                  do jj=1,2
+                     do ii=1,2
+                        
+                        fac = ifac(ii) * jfac(jj) * kfac(kk)
+                        
+                        if (fac.ne.0) then
+                           CHKIDX (i0+ii, j0+jj, k0+kk, \
+                                   srciext,srcjext,srckext, "source")
+                           res = res + fac * src(i0+ii, j0+jj, k0+kk)
+                        end if
+                        
+                     end do
+                  end do
+               end do
                
-               fac = ifac(2) * jfac(2) * kfac(2)
-               if (fac.ne.0) res = res + fac * src(i0+2, j0+2, k0+2)
+c$$$               fac = ifac(1) * jfac(1) * kfac(1)
+c$$$               if (fac.ne.0) res = res + fac * src(i0+1, j0+1, k0+1)
+c$$$               
+c$$$               fac = ifac(2) * jfac(1) * kfac(1)
+c$$$               if (fac.ne.0) res = res + fac * src(i0+2, j0+1, k0+1)
+c$$$               
+c$$$               fac = ifac(1) * jfac(2) * kfac(1)
+c$$$               if (fac.ne.0) res = res + fac * src(i0+1, j0+2, k0+1)
+c$$$               
+c$$$               fac = ifac(2) * jfac(2) * kfac(1)
+c$$$               if (fac.ne.0) res = res + fac * src(i0+2, j0+2, k0+1)
+c$$$               
+c$$$               fac = ifac(1) * jfac(1) * kfac(2)
+c$$$               if (fac.ne.0) res = res + fac * src(i0+1, j0+1, k0+2)
+c$$$               
+c$$$               fac = ifac(2) * jfac(1) * kfac(2)
+c$$$               if (fac.ne.0) res = res + fac * src(i0+2, j0+1, k0+2)
+c$$$               
+c$$$               fac = ifac(1) * jfac(2) * kfac(2)
+c$$$               if (fac.ne.0) res = res + fac * src(i0+1, j0+2, k0+2)
+c$$$               
+c$$$               fac = ifac(2) * jfac(2) * kfac(2)
+c$$$               if (fac.ne.0) res = res + fac * src(i0+2, j0+2, k0+2)
                
+               CHKIDX (dstioff+i+1, dstjoff+j+1, dstkoff+k+1, \
+                       dstiext,dstjext,dstkext, "destination")
                dst (dstioff+i+1, dstjoff+j+1, dstkoff+k+1) = dstdiv * res
                
             end do
diff --git a/Carpet/CarpetLib/src/prolongate_3d_real8_2tl.F77 b/Carpet/CarpetLib/src/prolongate_3d_real8_2tl.F77
index 8fab444ae..a89ccfbba 100644
--- a/Carpet/CarpetLib/src/prolongate_3d_real8_2tl.F77
+++ b/Carpet/CarpetLib/src/prolongate_3d_real8_2tl.F77
@@ -1,7 +1,20 @@
 c     -*-Fortran-*-
-c     $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_2tl.F77,v 1.8 2003/02/24 17:43:10 schnetter Exp $
+c     $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_2tl.F77,v 1.9 2003/02/25 22:57:00 schnetter Exp $
       
 #include "cctk.h"
+
+
+
+#define CHKIDX(i,j,k, imax,jmax,kmax, where)                                  \
+      if ((i).lt.1 .or. (i).gt.(imax)                                         \
+           .or. (j).lt.1 .or. (j).gt.(jmax)                                   \
+           .or. (k).lt.1 .or. (k).gt.(kmax)) then                           &&\
+         write (msg, '(a, " array index out of bounds: shape is (",i4,",",i4,",",i4,"), index is (",i4,",",i4,",",i4,")")') \
+              (where), (imax), (jmax), (kmax), (i), (j), (k)                &&\
+         call CCTK_WARN (0, msg(1:len_trim(msg)))                           &&\
+      end if
+
+
       
       subroutine prolongate_3d_real8_2tl (
      $     src1, t1, src2, t2, srciext, srcjext, srckext,
@@ -48,6 +61,8 @@ c     bbox(:,3) is stride
       CCTK_REAL8 res
       integer d
       
+      character msg*1000
+      
       
       
       do d=1,3
@@ -67,16 +82,16 @@ 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).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 (regbbox(d,1).lt.srcbbox(d,1)
      $        .or. regbbox(d,1).lt.dstbbox(d,1)
      $        .or. regbbox(d,2).gt.srcbbox(d,2)
      $        .or. regbbox(d,2).gt.dstbbox(d,2)) then
             call CCTK_WARN (0, "Internal error: region extent is not contained in array extent")
          end if
-         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
       end do
       
       if (srciext.ne.(srcbbox(1,2)-srcbbox(1,1))/srcbbox(1,3)+1
@@ -151,6 +166,8 @@ c     Loop over fine region
                         fac = ifac(ii) * jfac(jj) * kfac(kk)
                         
                         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)
@@ -160,6 +177,8 @@ c     Loop over fine region
                   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
                
             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 f48345562..4d9f1f899 100644
--- a/Carpet/CarpetLib/src/prolongate_3d_real8_2tl_o3.F77
+++ b/Carpet/CarpetLib/src/prolongate_3d_real8_2tl_o3.F77
@@ -1,7 +1,20 @@
 c     -*-Fortran-*-
-c     $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_2tl_o3.F77,v 1.10 2003/02/24 17:43:10 schnetter Exp $
+c     $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_2tl_o3.F77,v 1.11 2003/02/25 22:57:00 schnetter 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,
@@ -50,6 +63,8 @@ c     bbox(:,3) is stride
       CCTK_REAL8 res
       integer d
       
+      character msg*1000
+      
       
       
       do d=1,3
@@ -69,23 +84,33 @@ 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).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
          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) offsetlo = 0
+         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) offsethi = 0
+         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")
          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
       end do
       
       if (srciext.ne.(srcbbox(1,2)-srcbbox(1,1))/srcbbox(1,3)+1
@@ -166,6 +191,8 @@ c     Loop over fine region
                         fac = ifac(ii) * jfac(jj) * kfac(kk)
                         
                         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)
@@ -175,6 +202,8 @@ c     Loop over fine region
                   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
                
             end do
diff --git a/Carpet/CarpetLib/src/prolongate_3d_real8_3tl.F77 b/Carpet/CarpetLib/src/prolongate_3d_real8_3tl.F77
index 2275e4b7f..820dbe89c 100644
--- a/Carpet/CarpetLib/src/prolongate_3d_real8_3tl.F77
+++ b/Carpet/CarpetLib/src/prolongate_3d_real8_3tl.F77
@@ -1,7 +1,20 @@
 c     -*-Fortran-*-
-c     $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_3tl.F77,v 1.6 2003/02/24 17:43:10 schnetter Exp $
+c     $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_3tl.F77,v 1.7 2003/02/25 22:57:00 schnetter 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,
@@ -50,6 +63,8 @@ c     bbox(:,3) is stride
       CCTK_REAL8 res
       integer d
       
+      character msg*1000
+      
       
       
       do d=1,3
@@ -69,16 +84,16 @@ 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).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 (regbbox(d,1).lt.srcbbox(d,1)
      $        .or. regbbox(d,1).lt.dstbbox(d,1)
      $        .or. regbbox(d,2).gt.srcbbox(d,2)
      $        .or. regbbox(d,2).gt.dstbbox(d,2)) then
             call CCTK_WARN (0, "Internal error: region extent is not contained in array extent")
          end if
-         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
       end do
       
       if (srciext.ne.(srcbbox(1,2)-srcbbox(1,1))/srcbbox(1,3)+1
@@ -154,6 +169,8 @@ c     Loop over fine region
                         fac = ifac(ii) * jfac(jj) * kfac(kk)
                         
                         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)
@@ -164,6 +181,8 @@ c     Loop over fine region
                   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
                
             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 b5bb94384..053328536 100644
--- a/Carpet/CarpetLib/src/prolongate_3d_real8_3tl_o3.F77
+++ b/Carpet/CarpetLib/src/prolongate_3d_real8_3tl_o3.F77
@@ -1,7 +1,20 @@
 c     -*-Fortran-*-
-c     $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_3tl_o3.F77,v 1.10 2003/02/24 17:43:10 schnetter Exp $
+c     $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_3tl_o3.F77,v 1.11 2003/02/25 22:57:00 schnetter 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,
@@ -52,6 +65,8 @@ c     bbox(:,3) is stride
       CCTK_REAL8 res
       integer d
       
+      character msg*1000
+      
       
       
       do d=1,3
@@ -71,23 +86,33 @@ 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).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
          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) offsetlo = 0
+         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) offsethi = 0
+         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")
          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
       end do
       
       if (srciext.ne.(srcbbox(1,2)-srcbbox(1,1))/srcbbox(1,3)+1
@@ -169,6 +194,8 @@ c     Loop over fine region
                         fac = ifac(ii) * jfac(jj) * kfac(kk)
                         
                         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)
@@ -179,6 +206,8 @@ c     Loop over fine region
                   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
                
             end do
diff --git a/Carpet/CarpetLib/src/prolongate_3d_real8_o3.F77 b/Carpet/CarpetLib/src/prolongate_3d_real8_o3.F77
index b728980b7..65d6f1f75 100644
--- a/Carpet/CarpetLib/src/prolongate_3d_real8_o3.F77
+++ b/Carpet/CarpetLib/src/prolongate_3d_real8_o3.F77
@@ -1,7 +1,20 @@
 c     -*-Fortran-*-
-c     $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_o3.F77,v 1.7 2003/02/24 17:43:10 schnetter Exp $
+c     $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/prolongate_3d_real8_o3.F77,v 1.8 2003/02/25 22:57:00 schnetter 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,
@@ -41,6 +54,8 @@ c     bbox(:,3) is stride
       CCTK_REAL8 res
       integer d
       
+      character msg*1000
+      
       
       
       do d=1,3
@@ -60,23 +75,33 @@ 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).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
          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) offsetlo = 0
+         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) offsethi = 0
+         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")
          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
       end do
       
       if (srciext.ne.(srcbbox(1,2)-srcbbox(1,1))/srcbbox(1,3)+1
@@ -144,6 +169,8 @@ c     Loop over fine region
                         fac = ifac(ii) * jfac(jj) * kfac(kk)
                         
                         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)
                         end if
                         
@@ -151,6 +178,8 @@ c     Loop over fine region
                   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
                
             end do
diff --git a/Carpet/CarpetLib/src/restrict_3d_real8.F77 b/Carpet/CarpetLib/src/restrict_3d_real8.F77
index b846b18a6..782594a84 100644
--- a/Carpet/CarpetLib/src/restrict_3d_real8.F77
+++ b/Carpet/CarpetLib/src/restrict_3d_real8.F77
@@ -1,7 +1,20 @@
 c     -*-Fortran-*-
-c     $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/restrict_3d_real8.F77,v 1.5 2003/02/24 17:43:10 schnetter Exp $
+c     $Header: /home/eschnett/C/carpet/Carpet/Carpet/CarpetLib/src/restrict_3d_real8.F77,v 1.6 2003/02/25 22:57:00 schnetter Exp $
       
 #include "cctk.h"
+
+
+
+#define CHKIDX(i,j,k, imax,jmax,kmax, where)                                  \
+      if ((i).lt.1 .or. (i).gt.(imax)                                         \
+           .or. (j).lt.1 .or. (j).gt.(jmax)                                   \
+           .or. (k).lt.1 .or. (k).gt.(kmax)) then                           &&\
+         write (msg, '(a, " array index out of bounds: shape is (",i4,",",i4,",",i4,"), index is (",i4,",",i4,",",i4,")")') \
+              (where), (imax), (jmax), (kmax), (i), (j), (k)                &&\
+         call CCTK_WARN (0, msg(1:len_trim(msg)))                           &&\
+      end if
+
+
       
       subroutine restrict_3d_real8 (
      $     src, srciext, srcjext, srckext,
@@ -29,6 +42,8 @@ c     bbox(:,3) is stride
       integer i, j, k
       integer d
       
+      character msg*1000
+      
       
       
       do d=1,3
@@ -48,16 +63,16 @@ 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).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 (regbbox(d,1).lt.srcbbox(d,1)
      $        .or. regbbox(d,1).lt.dstbbox(d,1)
      $        .or. regbbox(d,2).gt.srcbbox(d,2)
      $        .or. regbbox(d,2).gt.dstbbox(d,2)) then
             call CCTK_WARN (0, "Internal error: region extent is not contained in array extent")
          end if
-         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
       end do
       
       if (srciext.ne.(srcbbox(1,2)-srcbbox(1,1))/srcbbox(1,3)+1
@@ -94,6 +109,10 @@ c     Loop over coarse region
          do j = 0, regjext-1
             do i = 0, regiext-1
                
+               CHKIDX (srcioff+srcifac*i+1, srcjoff+srcjfac*j+1, srckoff+srckfac*k+1, \
+                       srciext,srcjext,srckext, "source")
+               CHKIDX (dstioff+i+1, dstjoff+j+1, dstkoff+k+1, \
+                       dstiext,dstjext,dstkext, "destination")
                dst (dstioff+i+1, dstjoff+j+1, dstkoff+k+1)
      $              = src (srcioff+srcifac*i+1, srcjoff+srcjfac*j+1, srckoff+srckfac*k+1)