global: Change the way in which the grid hierarchy is stored

Change the way in which the grid hierarchy is stored.  The new hierarchy is
        map
        mglevel
        reflevel
        component
        timelevel
i.e., mglevel moved from the bottom to almost the top.  This is
because mglevel used to be a true multigrid level, but is now meant to
be a convergence level.

Do not allocate all storage all the time.  Allow storage to be
switched on an off per refinement level (and for a single mglevel,
which prompted the change above).  Handle storage management with
CCTK_{In,De}creaseGroupStorage instead of
CCTK_{En,Dis}ableGroupStorage.

darcs-hash:20050201225827-891bb-eae3b6bd092ae8d6b5e49be84c6f09f0e882933e.gz

35 files changed, 734 insertions, 653 deletions

diff --git a/Carpet/Carpet/interface.ccl b/Carpet/Carpet/interface.ccl
index 414e43d0e..d6e327736 100644
--- a/Carpet/Carpet/interface.ccl
+++ b/Carpet/Carpet/interface.ccl
@@ -77,7 +77,7 @@ USES FUNCTION ConvertFromExteriorBoundary
 
 # The true prototype of the routine below:
 #   int Carpet_Regrid (const cGH * cctkGH,
-#                      gh::rexts  * bbsss,
+#                      gh::mexts  * bbsss,
 #                      gh::rbnds  * obss,
 #                      gh::rprocs * pss);
 CCTK_INT FUNCTION Carpet_Regrid (CCTK_POINTER_TO_CONST IN cctkGH,  \
diff --git a/Carpet/Carpet/src/Recompose.cc b/Carpet/Carpet/src/Recompose.cc
index 8e4618c5a..59b5ecc73 100644
--- a/Carpet/Carpet/src/Recompose.cc
+++ b/Carpet/Carpet/src/Recompose.cc
@@ -65,47 +65,50 @@ namespace Carpet {
   
   
   
-  void CheckRegions (const gh::rexts & bbsss,
+  void CheckRegions (const gh::mexts & bbsss,
                      const gh::rbnds & obss,
                      const gh::rprocs& pss)
   {
-    // At least one level
+    // At least one multigrid level
     if (bbsss.size() == 0) {
-      CCTK_WARN (0, "I cannot set up a grid hierarchy with zero refinement levels.");
+      CCTK_WARN (0, "I cannot set up a grid hierarchy with zero multigrid levels.");
     }
     assert (bbsss.size() > 0);
+    // At least one refinement level
+    if (bbsss.at(0).size() == 0) {
+      CCTK_WARN (0, "I cannot set up a grid hierarchy with zero refinement levels.");
+    }
+    assert (bbsss.at(0).size() > 0);
     // At most maxreflevels levels
-    if ((int)bbsss.size() > maxreflevels) {
+    if ((int)bbsss.at(0).size() > maxreflevels) {
       CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING,
                   "I cannot set up a grid hierarchy with more than Carpet::max_refinement_levels refinement levels.  I found Carpet::max_refinement_levels=%d, while %d levels were requested.",
-                  (int)maxreflevels, (int)bbsss.size());
-    }
-    assert ((int)bbsss.size() <= maxreflevels);
-    for (int rl=0; rl<(int)bbsss.size(); ++rl) {
-      // No empty levels
-      assert (bbsss.at(rl).size() > 0);
-      for (int c=0; c<(int)bbsss.at(rl).size(); ++c) {
-        // At least one multigrid level
-        assert (bbsss.at(rl).at(c).size() > 0);
-        for (int ml=0; ml<(int)bbsss.at(rl).at(c).size(); ++ml) {
+                  (int)maxreflevels, (int)bbsss.at(0).size());
+    }
+    assert ((int)bbsss.at(0).size() <= maxreflevels);
+    for (int ml=0; ml<(int)bbsss.size(); ++ml) {
+      for (int rl=0; rl<(int)bbsss.at(0).size(); ++rl) {
+        // No empty levels
+        assert (bbsss.at(ml).at(rl).size() > 0);
+        for (int c=0; c<(int)bbsss.at(ml).at(rl).size(); ++c) {
           // Check sizes
           // Do allow processors with zero grid points
 //           assert (all(bbsss.at(rl).at(c).at(ml).lower() <= bbsssi.at(rl).at(c).at(ml).upper()));
           // Check strides
           const int str = ipow(reffact, maxreflevels-rl-1) * ipow(mgfact, ml);
-          assert (all(bbsss.at(rl).at(c).at(ml).stride() == str));
+          assert (all(bbsss.at(ml).at(rl).at(c).stride() == str));
           // Check alignments
-          assert (all(bbsss.at(rl).at(c).at(ml).lower() % str == 0));
-          assert (all(bbsss.at(rl).at(c).at(ml).upper() % str == 0));
+          assert (all(bbsss.at(ml).at(rl).at(c).lower() % str == 0));
+          assert (all(bbsss.at(ml).at(rl).at(c).upper() % str == 0));
         }
       }
     }
     
-    assert (pss.size() == bbsss.size());
-    assert (obss.size() == bbsss.size());
-    for (int rl=0; rl<(int)bbsss.size(); ++rl) {
-      assert (obss.at(rl).size() == bbsss.at(rl).size());
-      assert (pss.at(rl).size() == bbsss.at(rl).size());
+    assert (pss.size() == bbsss.at(0).size());
+    assert (obss.size() == bbsss.at(0).size());
+    for (int rl=0; rl<(int)bbsss.at(0).size(); ++rl) {
+      assert (obss.at(rl).size() == bbsss.at(0).at(rl).size());
+      assert (pss.at(rl).size() == bbsss.at(0).at(rl).size());
     }
     
   }
@@ -128,7 +131,7 @@ namespace Carpet {
     bool did_change = false;
     BEGIN_MAP_LOOP(cgh, CCTK_GF) {
       
-      gh::rexts  bbsss = vhh.at(map)->extents();
+      gh::mexts  bbsss = vhh.at(map)->extents();
       gh::rbnds  obss  = vhh.at(map)->outer_boundaries();
       gh::rprocs pss   = vhh.at(map)->processors();
       
@@ -177,13 +180,13 @@ namespace Carpet {
   {
     CCTK_INFO ("New grid structure (grid points):");
     cout << "   Refinement level " << reflevel << ", map " << map << endl;
-    for (int rl=0; rl<hh.reflevels(); ++rl) {
-      for (int c=0; c<hh.components(rl); ++c) {
-        for (int ml=0; ml<hh.mglevels(rl,c); ++ml) {
+    for (int ml=0; ml<hh.mglevels(); ++ml) {
+      for (int rl=0; rl<hh.reflevels(); ++rl) {
+        for (int c=0; c<hh.components(rl); ++c) {
           const int convfact = ipow(mgfact, ml);
           const int levfact = ipow(reffact, rl);
-          const ivect lower = hh.extents().at(rl).at(c).at(ml).lower();
-          const ivect upper = hh.extents().at(rl).at(c).at(ml).upper();
+          const ivect lower = hh.extents().at(ml).at(rl).at(c).lower();
+          const ivect upper = hh.extents().at(ml).at(rl).at(c).upper();
           assert (all(lower * levfact % maxreflevelfact == 0));
           assert (all(upper * levfact % maxreflevelfact == 0));
           assert (all(((upper - lower) * levfact / maxreflevelfact)
@@ -205,13 +208,13 @@ namespace Carpet {
     }
     
     CCTK_INFO ("New grid structure (coordinates):");
-    for (int rl=0; rl<hh.reflevels(); ++rl) {
-      for (int c=0; c<hh.components(rl); ++c) {
-        for (int ml=0; ml<hh.mglevels(rl,c); ++ml) {
+    for (int ml=0; ml<hh.mglevels(); ++ml) {
+      for (int rl=0; rl<hh.reflevels(); ++rl) {
+        for (int c=0; c<hh.components(rl); ++c) {
           const rvect origin = origin_space.at(0);
           const rvect delta = delta_space;
-          const ivect lower = hh.extents().at(rl).at(c).at(ml).lower();
-          const ivect upper = hh.extents().at(rl).at(c).at(ml).upper();
+          const ivect lower = hh.extents().at(ml).at(rl).at(c).lower();
+          const ivect upper = hh.extents().at(ml).at(rl).at(c).upper();
           const int convfact = ipow(mgfact, ml);
           const int levfact = ipow(reffact, rl);
           cout << "   [" << ml << "][" << rl << "][" << m << "][" << c << "]"
@@ -230,7 +233,7 @@ namespace Carpet {
   
   void OutputGridStructure (const cGH * const cgh,
                             const int m,
-                            const gh::rexts & bbsss,
+                            const gh::mexts & bbsss,
                             const gh::rbnds & obss,
                             const gh::rprocs& pss)
   {
@@ -274,13 +277,13 @@ namespace Carpet {
     
     file << "iteration " << cgh->cctk_iteration << endl;
     file << "maps " << maps << endl;
-    file << m << " reflevels " << bbsss.size() << endl;
-    for (int rl=0; rl<(int)bbsss.size(); ++rl) {
-      file << m << " " << rl << " components " << bbsss.at(rl).size() << endl;
-      for (int c=0; c<(int)bbsss.at(rl).size(); ++c) {
-        file << m << " " << rl << " " << c << " mglevels " << bbsss.at(rl).at(c).size() << endl;
-        for (int ml=0; ml<(int)bbsss.at(rl).at(c).size(); ++ml) {
-          file << m << " " << rl << " " << c << " " << ml << "   " << pss.at(rl).at(c) << " " << bbsss.at(rl).at(c).at(ml) << obss.at(rl).at(c) << endl;
+    file << m << " mglevels " << bbsss.size() << endl;
+    for (int ml=0; ml<(int)bbsss.size(); ++ml) {
+      file << m << " " << ml << " reflevels " << bbsss.at(ml).size() << endl;
+      for (int rl=0; rl<(int)bbsss.at(ml).size(); ++rl) {
+        file << m << " " << ml << " " << rl << " components " << bbsss.at(ml).at(rl).size() << endl;
+        for (int c=0; c<(int)bbsss.at(ml).at(rl).size(); ++c) {
+          file << m << " " << ml << " " << rl << " " << c << "   " << pss.at(rl).at(c) << " " << bbsss.at(ml).at(rl).at(c) << obss.at(rl).at(c) << endl;
         }
       }
     }
@@ -829,19 +832,42 @@ namespace Carpet {
   }
   
   void MakeMultigridBoxes (const cGH* cgh,
-                           vector<ibbox> const & bbs,
-                           vector<bbvect> const & obs,
-                           vector<vector<ibbox> >& bbss)
+                           vector<vector<ibbox> > const & bbss,
+                           vector<vector<bbvect> > const & obss,
+                           vector<vector<vector<ibbox> > > & bbsss)
   {
-    assert (bbs.size() == obs.size());
-    ibbox base;
-    for (int c=0; c<(int)bbs.size(); ++c) {
-      base = base.expanded_containing(bbs.at(c));
+    assert (bbss.size() == obss.size());
+    for (int rl=0; rl<(int)bbss.size(); ++rl) {
+      assert (bbss.at(rl).size() == obss.at(rl).size());
     }
-    bbss.resize(bbs.size());
-    for (int c=0; c<(int)bbs.size(); ++c) {
-      MakeMultigridBoxes (cgh, base, bbs.at(c), obs.at(c), bbss.at(c));
+    
+    bbsss.resize (mglevels);
+    for (int ml=0; ml<mglevels; ++ml) {
+      bbsss.at(ml).resize (bbss.size());
+      for (int rl=0; rl<(int)bbss.size(); ++rl) {
+        bbsss.at(ml).at(rl).resize (bbss.at(rl).size());
+      }
     }
+    
+    for (int rl=0; rl<(int)bbss.size(); ++rl) {
+      
+      ibbox base;
+      for (int c=0; c<(int)bbss.at(rl).size(); ++c) {
+        base = base.expanded_containing(bbss.at(rl).at(c));
+      }
+      
+      for (int c=0; c<(int)bbss.at(rl).size(); ++c) {
+        vector<ibbox> mg_bbs;
+        MakeMultigridBoxes
+          (cgh, base, bbss.at(rl).at(c), obss.at(rl).at(c), mg_bbs);
+        assert ((int)mg_bbs.size() == mglevels);
+        for (int ml=0; ml<mglevels; ++ml) {
+          bbsss.at(ml).at(rl).at(c) = mg_bbs.at(ml);
+        }
+        
+      } // for c
+      
+    } // for rl
   }
   
 } // namespace Carpet
diff --git a/Carpet/Carpet/src/SetupGH.cc b/Carpet/Carpet/src/SetupGH.cc
index e97e8fd83..e7f06fb54 100644
--- a/Carpet/Carpet/src/SetupGH.cc
+++ b/Carpet/Carpet/src/SetupGH.cc
@@ -596,18 +596,18 @@ namespace Carpet {
     vector<int> ps;
     SplitRegions (cgh, bbs, obs, ps);
     
-    // Create all multigrid levels
-    vector<vector<ibbox> > bbss;
-    MakeMultigridBoxes (cgh, bbs, obs, bbss);
-    
     // Only one refinement level
-    vector<vector<vector<ibbox> > > bbsss(1);
+    vector<vector<ibbox> > bbss(1);
     vector<vector<bbvect> > obss(1);
     vector<vector<int> > pss(1);
-    bbsss.at(0) = bbss;
+    bbss.at(0) = bbs;
     obss.at(0) = obs;
     pss.at(0) = ps;
     
+    // Create all multigrid levels
+    vector<vector<vector<ibbox> > > bbsss;
+    MakeMultigridBoxes (cgh, bbss, obss, bbsss);
+    
     // Check the regions
     CheckRegions (bbsss, obss, pss);
       
@@ -769,11 +769,11 @@ namespace Carpet {
   }
 
   void print_grid_structure (vector<gh*> & vhh, int m) {
-    const int rl = 0;
-    for (int c=0; c<vhh.at(m)->components(rl); ++c) {
-      for (int ml=0; ml<vhh.at(m)->mglevels(rl,c); ++ml) {
-        const ivect lower = vhh.at(m)->extents().at(rl).at(c).at(ml).lower();
-        const ivect upper = vhh.at(m)->extents().at(rl).at(c).at(ml).upper();
+    for (int ml=0; ml<vhh.at(m)->mglevels(); ++ml) {
+      const int rl = 0;
+      for (int c=0; c<vhh.at(m)->components(rl); ++c) {
+        const ivect lower = vhh.at(m)->extents().at(ml).at(rl).at(c).lower();
+        const ivect upper = vhh.at(m)->extents().at(ml).at(rl).at(c).upper();
         const int convfact = ipow(mgfact, ml);
         assert (all(lower % maxreflevelfact == 0));
         assert (all(upper % maxreflevelfact == 0));
@@ -952,9 +952,17 @@ namespace Carpet {
       } else {
         SplitRegions_Automatic (cgh, bbs, obs, ps);
       }
+        
+      // Only one refinement level
+      vector<vector<ibbox> > bbss(1);
+      vector<vector<bbvect> > obss(1);
+      vector<vector<int> > pss(1);
+      bbss.at(0) = bbs;
+      obss.at(0) = obs;
+      pss.at(0) = ps;
       
       // Create all multigrid levels
-      vector<vector<ibbox> > bbss (bbs.size());
+      vector<vector<vector<ibbox> > > bbsss (mglevels);
       ivect amgfact;
       iivect aoffset;
       for (int d=0; d<dim; ++d) {
@@ -962,17 +970,17 @@ namespace Carpet {
         aoffset[d][0] = 0;
         aoffset[d][1] = convoffsets[d];
       }
-      for (size_t c=0; c<bbs.size(); ++c) {
-        bbss.at(c).resize (mglevels);
-        bbss.at(c).at(0) = bbs.at(c);
-        for (int ml=1; ml<mglevels; ++ml) {
+      bbsss.at(0) = bbss;
+      for (int ml=1; ml<mglevels; ++ml) {
+        const int rl = 0;
+        for (int c=0; c<(int)bbss.at(rl).size(); ++c) {
           // this base
           ivect const baselo = ivect(0);
           ivect const baselo1 = baselo;
           // next finer grid
-          ivect const flo = bbss.at(c).at(ml-1).lower();
-          ivect const fhi = bbss.at(c).at(ml-1).upper();
-          ivect const fstr = bbss.at(c).at(ml-1).stride();
+          ivect const flo = bbsss.at(ml-1).at(rl).at(c).lower();
+          ivect const fhi = bbsss.at(ml-1).at(rl).at(c).upper();
+          ivect const fstr = bbsss.at(ml-1).at(rl).at(c).stride();
           // this grid
           ivect const str = fstr * amgfact;
           ivect const lo = flo + xpose(obs.at(c))[0].ifthen
@@ -983,17 +991,9 @@ namespace Carpet {
                 ivect(0));
           ivect const lo1 = baselo1 + (lo - baselo1 + str - 1) / str * str;
           ivect const hi1 = lo1 + (hi - lo1) / str * str;
-          bbss.at(c).at(ml) = ibbox(lo1, hi1, str);
+          bbsss.at(ml).at(rl).at(c) = ibbox(lo1, hi1, str);
         }
       }
-        
-      // Only one refinement level
-      vector<vector<vector<ibbox> > > bbsss(1);
-      vector<vector<bbvect> > obss(1);
-      vector<vector<int> > pss(1);
-      bbsss.at(0) = bbss;
-      obss.at(0) = obs;
-      pss.at(0) = ps;
 	
       // Recompose for this map
       char * const groupname = CCTK_GroupName (group);
@@ -1027,7 +1027,7 @@ namespace Carpet {
       
       arrdata.at(group).at(m).data.resize(CCTK_NumVarsInGroupI(group));
       for (int var=0; var<(int)arrdata.at(group).at(m).data.size(); ++var) {
-        arrdata.at(group).at(m).data.at(var) = 0;
+        arrdata.at(group).at(m).data.at(var) = NULL;
       }
       
     }
diff --git a/Carpet/Carpet/src/Storage.cc b/Carpet/Carpet/src/Storage.cc
index 0e606d9bd..c0cccde03 100644
--- a/Carpet/Carpet/src/Storage.cc
+++ b/Carpet/Carpet/src/Storage.cc
@@ -35,14 +35,19 @@ namespace Carpet {
     assert (timelevels);
     for (int n=0; n<n_groups; ++n) {
       assert (groups[n] >= 0 and groups[n] < CCTK_NumGroups());
-#if 0
-      for (int nn=0; nn<n; ++nn) {
-        assert (groups[nn] != groups[n]);
-      }
-#endif
+      // TODO: timelevels[n] can also be -1; in that case, all time
+      // levels should be activated / deactivated
       assert (timelevels[n] >= 0);
     }
     
+    bool const can_do = is_meta_mode() or is_global_mode() or is_level_mode();
+    bool const all_ml = is_meta_mode();
+    int const min_ml = all_ml ? 0        : mglevel;
+    int const max_ml = all_ml ? mglevels : mglevel+1;
+    bool const all_rl = is_meta_mode() or is_global_mode();
+    int const min_rl = all_rl ? 0         : reflevel;
+    int const max_rl = all_rl ? reflevels : reflevel+1;
+    
     int total_num_timelevels = 0;
     
     for (int n=0; n<n_groups; ++n) {
@@ -68,76 +73,93 @@ namespace Carpet {
       
       // Record previous number of allocated time levels
       if (status) {
-        status[n] = groupdata.at(group).activetimelevels;
+        status[n] = groupdata.at(group).info.activetimelevels;
       }
       
       // Only do something if the number of time levels actually needs
       // to be changed -- do nothing otherwise
-      if ((inc and timelevels[n] > groupdata.at(group).activetimelevels)
-          or
-          (! inc and timelevels[n] < groupdata.at(group).activetimelevels))
-      {
+
+      const bool do_increase
+        = inc and timelevels[n] > groupdata.at(group).info.activetimelevels;
+      const bool do_decrease
+        = ! inc and timelevels[n] < groupdata.at(group).info.activetimelevels;
+      if (do_increase or do_decrease) {
         
-        // Set the new number of active time levels
-        groupdata.at(group).activetimelevels = timelevels[n];
+        if (! can_do) {
+          char * const groupname = CCTK_GroupName (group);
+          char const * const modestring
+            = (is_meta_mode() ? "meta"
+               : is_global_mode() ? "global"
+               : is_level_mode() ? "level"
+               : is_singlemap_mode() ? "singlemap"
+               : is_local_mode() ? "local"
+               : NULL);
+          CCTK_VWarn (0, __LINE__, __FILE__, CCTK_THORNSTRING,
+                      "Cannot change storage for group \"%s\" in %s mode",
+                      groupname, modestring);
+          free (groupname);
+        }
+        assert (can_do);
         
-        // There is a difference between the Cactus time levels and
-        // the Carpet time levels.  If there are n time levels, then
-        // the Cactus time levels are numbered 0 ... n-1, with the
-        // current time level being 0.  In Carpet, the time levels are
-        // numbered -(n-1) ... 0, where the current time level is also
-        // 0.
-        const int tmin = - timelevels[n] + 1;
-        const int tmax = 0;
+        // Set the new number of active time levels
+        groupdata.at(group).info.activetimelevels = timelevels[n];
         
         // Allocate the time levels as well
-        for (int m=0; m<(int)arrdata.at(group).size(); ++m) {
-          for (int var=0; var<gp.numvars; ++var) {
-            const int vectorindex
-              = (gp.contiguous
-                 ? var
-                 : gp.vectorgroup ? var % gp.vectorlength : 0);
-            const int vectorlength
-              = (gp.contiguous
-                 ? gp.numvars
-                 : gp.vectorgroup ? gp.vectorlength : 1);
-            assert (vectorindex>=0 and vectorindex<gp.numvars);
-            assert (vectorlength>0 and vectorlength<=gp.numvars);
-            ggf* vectorleader
-              = (vectorindex>0
-                 ? arrdata.at(group).at(m).data.at(var - vectorindex)
-                 : NULL);
-            const int varindex = firstvarindex + var;
-            switch (gp.vartype) {
+        for (int ml=min_ml; ml<max_ml; ++ml) {
+          for (int rl=min_rl; rl<max_rl; ++rl) {
+            for (int m=0; m<(int)arrdata.at(group).size(); ++m) {
+              for (int var=0; var<gp.numvars; ++var) {
+                const int vectorindex
+                  = gp.vectorgroup ? var % gp.vectorlength : 0;
+                const int vectorlength
+                  = gp.vectorgroup ? gp.vectorlength : 1;
+                assert (vectorindex>=0 and vectorindex<gp.numvars);
+                assert (vectorlength>0 and vectorlength<=gp.numvars);
+                ggf* const vectorleader
+                  = (vectorindex>0
+                     ? arrdata.at(group).at(m).data.at(var - vectorindex)
+                     : NULL);
+                const int varindex = firstvarindex + var;
+#warning "TODO: allocate these in SetupGH, and after recomposing"
+                if (! arrdata.at(group).at(m).data.at(var)) {
+                  switch (gp.vartype) {
 #define TYPECASE(N,T)                                                   \
-              case N:                                                   \
-                arrdata.at(group).at(m).data.at(var) = new gf<T>        \
-                (varindex, groupdata.at(group).transport_operator,      \
-                 *arrdata.at(group).at(m).tt, *arrdata.at(group).at(m).dd, \
-                 tmin, tmax, prolongation_order_time,                   \
-                 vectorlength, vectorindex, (gf<T>*)vectorleader);      \
-              break;
+                    case N:                                             \
+                      arrdata.at(group).at(m).data.at(var) = new gf<T>  \
+                      (varindex,                                        \
+                       groupdata.at(group).transport_operator,          \
+                       *arrdata.at(group).at(m).tt,                     \
+                       *arrdata.at(group).at(m).dd,                     \
+                       prolongation_order_time,                         \
+                       vectorlength, vectorindex, (gf<T>*)vectorleader); \
+                    break;
 #include "typecase"
 #undef TYPECASE
-            default:
-              UnsupportedVarType (varindex);
-            } // switch gp.vartype
-            
-            // Set the data pointers for grid arrays
-            if (gp.grouptype != CCTK_GF) {
-              assert (m == 0);
-              int const c = CCTK_MyProc(cgh);
-              for (int tl=0; tl<gp.numtimelevels; ++tl) {
-                cgh->data[varindex][tl]
-                  = (tl < groupdata.at(group).activetimelevels
-                     ? ((*arrdata.at(group).at(m).data.at(var))
-                        (-tl, 0, c, 0)->storage())
-                     : NULL);
-              }
-            }
-            
-          } // for var
-        } // for m
+                  default:
+                    UnsupportedVarType (varindex);
+                  } // switch gp.vartype
+                } // if ! allocated
+              
+                arrdata.at(group).at(m).data.at(var)->set_timelevels
+                  (ml, rl, timelevels[n]);
+              
+                // Set the data pointers for grid arrays
+                if (gp.grouptype != CCTK_GF) {
+                  assert (rl==0 and m==0);
+                  int const c = CCTK_MyProc(cgh);
+                  for (int tl=0; tl<gp.numtimelevels; ++tl) {
+                    cgh->data[varindex][tl]
+                      = (tl < groupdata.at(group).info.activetimelevels
+                         ? ((*arrdata.at(group).at(m).data.at(var))
+                            (tl, 0, c, 0)->storage())
+                         : NULL);
+                  }
+                } // if grouptype != GF
+              
+              } // for var
+            } // for m
+          } // for rl
+        } // for ml
         
       } // if really change the number of active time levels
       
@@ -146,8 +168,8 @@ namespace Carpet {
       if (gp.grouptype == CCTK_GF) {
         if (max_refinement_levels > 1) {
           if (groupdata.at(group).transport_operator != op_none) {
-            if (groupdata.at(group).activetimelevels != 0
-                and (groupdata.at(group).activetimelevels
+            if (groupdata.at(group).info.activetimelevels != 0
+                and (groupdata.at(group).info.activetimelevels
                      <= prolongation_order_time))
             {
               char * const groupname = CCTK_GroupName (group);
@@ -163,7 +185,7 @@ namespace Carpet {
 #endif
       
       // Record current number of time levels
-      total_num_timelevels += groupdata.at(group).activetimelevels;
+      total_num_timelevels += groupdata.at(group).info.activetimelevels;
       
     } // for n
     
@@ -229,11 +251,8 @@ namespace Carpet {
       group = CCTK_GroupIndex(groupname);
     }
     assert (group>=0 and group<CCTK_NumGroups());
-    const int timelevels = 0;
-    int status;
-    GroupStorageIncrease (cgh, 1, &group, &timelevels, &status);
     // Return whether storage is allocated
-    return status;
+    return groupdata.at(group).info.activetimelevels > 0;
   }
   
   
diff --git a/Carpet/Carpet/src/functions.hh b/Carpet/Carpet/src/functions.hh
index e5fdf4ed2..14ba8c50d 100644
--- a/Carpet/Carpet/src/functions.hh
+++ b/Carpet/Carpet/src/functions.hh
@@ -40,7 +40,7 @@ namespace Carpet {
    
    
   // Helpers for recomposing the grid hierarchy
-  void CheckRegions (const gh::rexts & bbsss,
+  void CheckRegions (const gh::mexts & bbsss,
                      const gh::rbnds & obss,
                      const gh::rprocs& pss);
   
@@ -48,7 +48,7 @@ namespace Carpet {
   
   void OutputGridStructure (const cGH *cgh,
                             const int m,
-                            const gh::rexts & bbsss,
+                            const gh::mexts & bbsss,
                             const gh::rbnds & obss,
                             const gh::rprocs& pss);
   
@@ -66,9 +66,9 @@ namespace Carpet {
                                vector<bbvect>& obs, vector<int>& ps);
   
   void MakeMultigridBoxes (const cGH* cgh,
-                           vector<ibbox> const & bbs,
-                           vector<bbvect> const & obs,
-                           vector<vector<ibbox> > & bbss);
+                           gh::rexts const & bbss,
+                           gh::rbnds const & obss,
+                           gh::mexts & bbsss);
   
 } // namespace Carpet
 
diff --git a/Carpet/Carpet/src/modes.cc b/Carpet/Carpet/src/modes.cc
index d489cf88b..10237e2df 100644
--- a/Carpet/Carpet/src/modes.cc
+++ b/Carpet/Carpet/src/modes.cc
@@ -90,9 +90,9 @@ namespace Carpet {
         const int m = 0;
         const int c = CCTK_MyProc(cgh);
         
-        const ibbox& base = arrdata.at(group).at(m).hh->bases().at(rl).at(ml);
+        const ibbox& base = arrdata.at(group).at(m).hh->bases().at(ml).at(rl);
         const bbvect& obnds = arrdata.at(group).at(m).hh->outer_boundaries().at(rl).at(c);
-	const ibbox& ext = arrdata.at(group).at(m).dd->boxes.at(rl).at(c).at(ml).exterior;
+	const ibbox& ext = arrdata.at(group).at(m).dd->boxes.at(ml).at(rl).at(c).exterior;
         
         ivect::ref(const_cast<int*>(groupdata.at(group).info.nghostzones))
           = arrdata.at(group).at(m).dd->lghosts;
@@ -136,7 +136,7 @@ namespace Carpet {
             for (int tl=0; tl<num_tl; ++tl) {
               ggf * const ff = arrdata.at(group).at(m).data.at(var);
               if (ff) {
-                gdata * const data = (*ff) (-tl, rl, c, ml);
+                gdata * const data = (*ff) (tl, rl, c, ml);
                 assert (data);
                 cgh->data[firstvar+var][tl] = data->storage();
               } else {
@@ -288,8 +288,8 @@ namespace Carpet {
     carpetGH.map = map = m;
     
     // Set grid shape
-    const ibbox& coarseext = vdd.at(map)->bases.at(0       ).at(mglevel).exterior;
-    const ibbox& baseext   = vdd.at(map)->bases.at(reflevel).at(mglevel).exterior;
+    const ibbox& coarseext = vdd.at(map)->bases.at(mglevel).at(0       ).exterior;
+    const ibbox& baseext   = vdd.at(map)->bases.at(mglevel).at(reflevel).exterior;
     assert (all (baseext.lower() % baseext.stride() == 0));
     assert (all ((baseext.lower() - coarseext.lower()) % baseext.stride() == 0));
     ivect::ref(cgh->cctk_levoff) = (baseext.lower() - coarseext.lower()) / baseext.stride();
@@ -353,9 +353,9 @@ namespace Carpet {
     component = c;
     
     // Set cGH fields
-    const ibbox& baseext = vdd.at(map)->bases.at(reflevel).at(mglevel).exterior;
+    const ibbox& baseext = vdd.at(map)->bases.at(mglevel).at(reflevel).exterior;
     const bbvect& obnds = vhh.at(map)->outer_boundaries().at(reflevel).at(component);
-    const ibbox& ext = vdd.at(map)->boxes.at(reflevel).at(component).at(mglevel).exterior;
+    const ibbox& ext = vdd.at(map)->boxes.at(mglevel).at(reflevel).at(component).exterior;
     
     ivect::ref(cgh->cctk_lsh) = ext.shape() / ext.stride();
     ivect::ref(cgh->cctk_lbnd)
@@ -416,7 +416,7 @@ namespace Carpet {
             for (int tl=0; tl<num_tl; ++tl) {
               ggf * const ff = arrdata.at(group).at(map).data.at(var);
               if (ff) {
-                gdata * const data = (*ff) (-tl, reflevel, component, mglevel);
+                gdata * const data = (*ff) (tl, reflevel, component, mglevel);
                 assert (data);
                 cgh->data[firstvar+var][tl] = data->storage();
               } else {
diff --git a/Carpet/Carpet/src/variables.hh b/Carpet/Carpet/src/variables.hh
index 821b77ae3..ae5ade6df 100644
--- a/Carpet/Carpet/src/variables.hh
+++ b/Carpet/Carpet/src/variables.hh
@@ -115,7 +115,6 @@ namespace Carpet {
   // Data for the groups
   struct groupdesc {
     cGroupDynamicData info;
-    int activetimelevels;
     operator_type transport_operator; // prolongation and restriction
   };
   extern vector<groupdesc> groupdata; // [group]
diff --git a/Carpet/CarpetIOASCII/src/ioascii.cc b/Carpet/CarpetIOASCII/src/ioascii.cc
index 296957809..b5a8e2121 100644
--- a/Carpet/CarpetIOASCII/src/ioascii.cc
+++ b/Carpet/CarpetIOASCII/src/ioascii.cc
@@ -550,7 +550,8 @@ namespace CarpetIOASCII {
 
                 ivect offset1;
                 if (grouptype == CCTK_GF) {
-                  const ibbox& baseext = vdd.at(Carpet::map)->bases.at(0).at(mglevel).exterior;
+                  const ibbox& baseext
+                    = vdd.at(Carpet::map)->bases.at(mglevel).at(0).exterior;
                   offset1 = baseext.lower() + offset * ext.stride();
                 } else {
                   offset1 = offset * ext.stride();
diff --git a/Carpet/CarpetIOHDF5/src/Output.cc b/Carpet/CarpetIOHDF5/src/Output.cc
index f5351b748..8f5601dce 100644
--- a/Carpet/CarpetIOHDF5/src/Output.cc
+++ b/Carpet/CarpetIOHDF5/src/Output.cc
@@ -277,7 +277,7 @@ int WriteVar (const cGH* const cctkGH, const hid_t writer,
     {
       // Using "interior" removes ghost zones and refinement boundaries.
       bboxes += arrdata.at(gindex).at(Carpet::map).dd->
-                boxes.at(refinementlevel).at(component).at(mglevel).interior;
+                boxes.at(mglevel).at(refinementlevel).at(component).interior;
     } END_COMPONENT_LOOP;
 
     // Normalise the set, i.e., try to represent the set with fewer bboxes.
@@ -372,7 +372,7 @@ int WriteVar (const cGH* const cctkGH, const hid_t writer,
         // (use either the interior or exterior here, as we did above)
         ibbox const overlap = *bbox &
           arrdata.at(gindex).at(Carpet::map).dd->
-          boxes.at(refinementlevel).at(component).at(mglevel).interior;
+          boxes.at(mglevel).at(refinementlevel).at(component).interior;
 
         // Continue if this component is not part of this combination
         if (overlap.empty())
@@ -837,7 +837,7 @@ static void AddAttributes (const cGH *const cctkGH, const char *fullname,
                              coord_handles, "COORDINATES") >= 0)
   {
     const ibbox& baseext = 
-    vdd.at(Carpet::map)->bases.at(reflevel).at(mglevel).exterior;
+    vdd.at(Carpet::map)->bases.at(mglevel).at(reflevel).exterior;
 
     const ivect pos = (bbox->lower() - baseext.lower()) / bbox->stride();
 
diff --git a/Carpet/CarpetIOHDF5/src/Recover.cc b/Carpet/CarpetIOHDF5/src/Recover.cc
index 4809c1470..854970ec2 100644
--- a/Carpet/CarpetIOHDF5/src/Recover.cc
+++ b/Carpet/CarpetIOHDF5/src/Recover.cc
@@ -761,7 +761,7 @@ static int InputVarAs (const cGH* const cctkGH, const int vindex,
       ibset all_exterior;
       for (size_t c=0; c<thedd.boxes.at(rl).size(); ++c)
       {
-        all_exterior |= thedd.boxes.at(rl).at(c).at(mglevel).exterior;
+        all_exterior |= thedd.boxes.at(mglevel).at(rl).at(c).exterior;
       }
       if (regions_read.at(m) != all_exterior)
       {
diff --git a/Carpet/CarpetInterp/src/interp.cc b/Carpet/CarpetInterp/src/interp.cc
index 0d8b14705..4340e77f9 100644
--- a/Carpet/CarpetInterp/src/interp.cc
+++ b/Carpet/CarpetInterp/src/interp.cc
@@ -272,11 +272,11 @@ namespace CarpetInterp {
           
           int const fact = maxreflevelfact / ipow(reffact, rl) * ipow(mgfact, mglevel);
           ivect const ipos = ivect(floor((pos - lower) / (delta * fact) + 0.5)) * fact;
-          assert (all(ipos % vhh.at(m)->bases().at(rl).at(ml).stride() == 0));
+          assert (all(ipos % vhh.at(m)->bases().at(ml).at(rl).stride() == 0));
           
           // TODO: use something faster than a linear search
           for (int c=0; c<vhh.at(m)->components(rl); ++c) {
-            if (vhh.at(m)->extents().at(rl).at(c).at(ml).contains(ipos)) {
+            if (vhh.at(m)->extents().at(ml).at(rl).at(c).contains(ipos)) {
               rlev.at(n) = rl;
               home.at(n) = c;
               goto found;
diff --git a/Carpet/CarpetLib/src/dh.cc b/Carpet/CarpetLib/src/dh.cc
index fff52c35c..dc3354204 100644
--- a/Carpet/CarpetLib/src/dh.cc
+++ b/Carpet/CarpetLib/src/dh.cc
@@ -79,15 +79,16 @@ void dh::recompose (const bool do_prolongate) {
   }
 }
 
-void dh::allocate_bboxes() {
-  boxes.resize(h.reflevels());
-  for (int rl=0; rl<h.reflevels(); ++rl) {
-    boxes.at(rl).resize(h.components(rl));
-    for (int c=0; c<h.components(rl); ++c) {
-      boxes.at(rl).at(c).resize(h.mglevels(rl,c));
-      for (int ml=0; ml<h.mglevels(rl,c); ++ml) {
-        const ibbox intr = h.extents().at(rl).at(c).at(ml);
-        dboxes & b = boxes.at(rl).at(c).at(ml);
+void dh::allocate_bboxes ()
+{
+  boxes.resize(h.mglevels());
+  for (int ml=0; ml<h.mglevels(); ++ml) {
+    boxes.at(ml).resize(h.reflevels());
+    for (int rl=0; rl<h.reflevels(); ++rl) {
+      boxes.at(ml).at(rl).resize(h.components(rl));
+      for (int c=0; c<h.components(rl); ++c) {
+        const ibbox intr = h.extents().at(ml).at(rl).at(c);
+        dboxes & b = boxes.at(ml).at(rl).at(c);
 
         // Interior
         // (the interior of the grid has the extent as specified by
@@ -109,21 +110,21 @@ void dh::allocate_bboxes() {
         // (interior + boundaries = exterior)
         b.boundaries = b.exterior - intr;
 
-      } // for ml
-    } // for c
-  } // for rl
+      } // for c
+    } // for rl
+  } // for ml
 }
 
-// Loops over each refinement level, each component, and each multigrid
-// level, executing the "boxesop" member function argument on the corresponding
-// element of the "boxes" member 
-void dh::foreach_reflevel_component_mglevel (dh::boxesop op) {
-
-  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) {
-        dboxes & box = boxes.at(rl).at(c).at(ml);
-        (this->*op)( box, rl, c, ml ); // evaluate member function
+// Loops over each multigrid level, each refinement level, and each
+// component, executing the "boxesop" member function argument on the
+// corresponding element of the "boxes" member
+void dh::foreach_reflevel_component_mglevel (dh::boxesop op)
+{
+  for (int ml=0; ml<h.mglevels(); ++ml) {
+    for (int rl=0; rl<h.reflevels(); ++rl) {
+      for (int c=0; c<h.components(rl); ++c) {
+        dboxes & box = boxes.at(ml).at(rl).at(c);
+        (this->*op)(box, rl, c, ml); // evaluate member function
       }
     }
   }
@@ -157,8 +158,7 @@ void dh::intersect_sync_with_interior( dh::dboxes & box, int rl, int c, int ml )
 
   // Sync boxes
   for (int cc=0; cc<h.components(rl); ++cc) {
-    dboxes & box1 = boxes.at(rl).at(cc).at(ml);
-    assert (ml<h.mglevels(rl,cc));
+    dboxes & box1 = boxes.at(ml).at(rl).at(cc);
     // intersect boundaries with interior of that component
     ibset ovlp = bnds & box1.interior;
     ovlp.normalize();
@@ -176,7 +176,7 @@ void dh::setup_multigrid_boxes( dh::dboxes & box, int rl, int c, int ml )
 
   // Multigrid boxes
   if (ml>0) {
-    dboxes & bbox = boxes.at(rl).at(c).at(ml-1);
+    dboxes & bbox = boxes.at(ml-1).at(rl).at(c);
     const ibbox intrf = bbox.interior;
     const ibbox extrf = bbox.exterior;
     // Restriction (interior)
@@ -217,7 +217,7 @@ void dh::setup_refinement_interior_boxes( dh::dboxes & box, int rl, int c, int m
   // Refinement boxes
   if (rl<h.reflevels()-1) {
     for (int cc=0; cc<h.components(rl+1); ++cc) {
-      dboxes & box1 = boxes.at(rl+1).at(cc).at(ml);
+      dboxes & box1 = boxes.at(ml).at(rl+1).at(cc);
       const ibbox intrf = box1.interior;
       // Prolongation (interior)
       // TODO: prefer boxes from the same processor
@@ -257,7 +257,7 @@ void dh::setup_refinement_exterior_boxes( dh::dboxes & box, int rl, int c, int m
   // Refinement boxes
   if (rl<h.reflevels()-1) {
     for (int cc=0; cc<h.components(rl+1); ++cc) {
-      dboxes & box1 = boxes.at(rl+1).at(cc).at(ml);
+      dboxes & box1 = boxes.at(ml).at(rl+1).at(cc);
       const ibbox intrf = box1.interior;
       const ibbox& extrf = box1.exterior;
       const ibset& bndsf = box1.boundaries;
@@ -332,7 +332,7 @@ void dh::setup_restrict_interior_boxes( dh::dboxes & box, int rl, int c, int ml
   // Refinement boxes
   if (rl<h.reflevels()-1) {
     for (int cc=0; cc<h.components(rl+1); ++cc) {
-      dboxes & box1 = boxes.at(rl+1).at(cc).at(ml);
+      dboxes & box1 = boxes.at(ml).at(rl+1).at(cc);
       const ibbox intrf = box1.interior;
       // Restriction (interior)
       {
@@ -422,7 +422,7 @@ void dh::assert_assert_assert( dh::dboxes & box, int rl, int c, int ml )
     if (rl==0) {
       const iblistvect& recv_ref_coarse = box.recv_ref_coarse;
       assert (recv_ref_coarse.empty());
-    } else {              // rl!=0
+    } else {                    // rl!=0
       const iblistvect& recv_ref_coarse = box.recv_ref_coarse;
       ibset intr = box.interior;
       ibset received;
@@ -504,23 +504,19 @@ void dh::assert_assert_assert( dh::dboxes & box, int rl, int c, int ml )
 
 void dh::calculate_bases () {
   // Calculate bases
-  bases.resize(h.reflevels());
-  for (int rl=0; rl<h.reflevels(); ++rl) {
-    if (h.components(rl)==0) {
-      bases.at(rl).resize(0);
-    } else {
-      bases.at(rl).resize(h.mglevels(rl,0));
-      for (int ml=0; ml<h.mglevels(rl,0); ++ml) {
-        dbases & base = bases.at(rl).at(ml);
-        base.exterior = ibbox();
-        base.interior = ibbox();
-        for (int c=0; c<h.components(rl); ++c) {
-          dboxes & box = boxes.at(rl).at(c).at(ml);
-          base.exterior = (base.exterior.expanded_containing(box.exterior));
-          base.interior = (base.interior.expanded_containing(box.interior));
-        }
-        base.boundaries = base.exterior - base.interior;
+  bases.resize(h.mglevels());
+  for (int ml=0; ml<h.mglevels(); ++ml) {
+    bases.at(ml).resize(h.reflevels());
+    for (int rl=0; rl<h.reflevels(); ++rl) {
+      dbases & base = bases.at(ml).at(rl);
+      base.exterior = ibbox();
+      base.interior = ibbox();
+      for (int c=0; c<h.components(rl); ++c) {
+        dboxes & box = boxes.at(ml).at(rl).at(c);
+        base.exterior = (base.exterior.expanded_containing(box.exterior));
+        base.interior = (base.interior.expanded_containing(box.interior));
       }
+      base.boundaries = base.exterior - base.interior;
     }
   }
 }
@@ -569,8 +565,8 @@ void dh::save_memory ( bool do_prolongate ) {
         assert (! vectorleader);
         (*f)->recompose_free (rl);
       } else {
-        // treat vector groups specially: delete the leader only
-        // after all other elements have been deleted
+        // treat vector groups specially: delete the leader only after
+        // all other elements have been deleted
         if ((*f)->vectorindex == 0) {
           // first element: delete nothing
           if (rl == 0) {
@@ -632,7 +628,7 @@ void dh::do_output_bboxes( dh::dboxes & box, int rl, int c, int ml )
 {
   cout << endl;
   cout << "dh bboxes:" << endl;
-  cout << "rl=" << rl << " c=" << c << " ml=" << ml << endl;
+  cout << "ml=" << ml << " rl=" << rl << " c=" << c << endl;
   cout << "exterior=" << box.exterior << endl;
   cout << "interior=" << box.interior << endl;
   cout << "send_mg_fine=" << box.send_mg_fine << endl;
@@ -652,14 +648,15 @@ void dh::do_output_bboxes( dh::dboxes & box, int rl, int c, int ml )
   cout << "recv_not=" << box.recv_not << endl;
 }
 
-void dh::output_bases () {
-  for (int rl=0; rl<h.reflevels(); ++rl) {
-    if (h.components(rl)>0) {
-      for (int ml=0; ml<h.mglevels(rl,0); ++ml) {
-        dbases & base = bases.at(rl).at(ml);
+void dh::output_bases ()
+{
+  for (int ml=0; ml<h.mglevels(); ++ml) {
+    for (int rl=0; rl<h.reflevels(); ++rl) {
+      if (h.components(rl)>0) {
+        dbases & base = bases.at(ml).at(rl);
         cout << endl;
         cout << "dh bases:" << endl;
-        cout << "rl=" << rl << " ml=" << ml << endl;
+        cout << "ml=" << ml << " rl=" << rl << endl;
         cout << "exterior=" << base.exterior << endl;
         cout << "interior=" << base.interior << endl;
         cout << "boundaries=" << base.boundaries << endl;
diff --git a/Carpet/CarpetLib/src/dh.hh b/Carpet/CarpetLib/src/dh.hh
index a0e8b8e0d..8fa11f152 100644
--- a/Carpet/CarpetLib/src/dh.hh
+++ b/Carpet/CarpetLib/src/dh.hh
@@ -55,7 +55,7 @@ public:
     iblistvect send_ref_coarse;
     iblistvect recv_ref_fine;
     iblistvect recv_ref_coarse;
-    iblistvect send_sync;        // send while syncing
+    iblistvect send_sync;       // send while syncing
     iblistvect send_ref_bnd_fine; // sent to finer grids
     
     ibset boundaries;           // boundaries
@@ -63,14 +63,6 @@ public:
     iblistvect recv_ref_bnd_coarse; // received from coarser grids
     ibset sync_not;             // not received while syncing (outer boundary of that level)
     ibset recv_not;             // not received while syncing or prolongating (globally outer boundary)
-    
-#if 0
-    // after regridding:
-    iblistvect prev_send;       // sent from previous dh
-    iblistvect recv_prev;       // received from previous dh
-    iblistvect send_prev_fine;  // sent to finer
-    iblistvect recv_prev_coarse; // received from coarser
-#endif
   };
   
 private:
@@ -81,19 +73,19 @@ private:
     ibset boundaries;           // boundaries
   };
   
-  typedef vector<dboxes> mboxes; // ... for each multigrid level
-  typedef vector<mboxes> cboxes; // ... for each component
+  typedef vector<dboxes> cboxes; // ... for each component
   typedef vector<cboxes> rboxes; // ... for each refinement level
+  typedef vector<rboxes> mboxes; // ... for each multigrid level
   
-  typedef vector<dbases> mbases; // ... for each multigrid level
-  typedef vector<mbases> rbases; // ... for each refinement level
+  typedef vector<dbases> rbases; // ... for each refinement level
+  typedef vector<rbases> mbases; // ... for each multigrid level
  
-  void allocate_bboxes();
-                                 // generic member function taking a dboxes,
-                                 // a refinement level, a component, and a 
-                                 // multigrid level
-  typedef void    (dh::*boxesop)( dboxes &, int rl, int c, int ml ); 
-  void foreach_reflevel_component_mglevel ( boxesop op );
+  void allocate_bboxes ();
+
+  // generic member function taking a dboxes, a refinement level, a
+  // component, and a multigrid level
+  typedef void (dh::*boxesop) (dboxes &, int rl, int c, int ml); 
+  void foreach_reflevel_component_mglevel (boxesop op);
 
                                  // these all of form 'boxesop'
   void setup_sync_and_refine_boxes( dboxes & b, int rl, int c, int ml );
@@ -120,8 +112,8 @@ public:                         // should be readonly
   int prolongation_order_space; // order of spatial prolongation operator
   int buffer_width;             // buffer inside refined grids
   
-  rboxes boxes;
-  rbases bases;
+  mboxes boxes;
+  mbases bases;
   
   list<ggf*> gfs;               // list of all grid functions
   
diff --git a/Carpet/CarpetLib/src/gf.cc b/Carpet/CarpetLib/src/gf.cc
index e493f9491..e3fde8dd0 100644
--- a/Carpet/CarpetLib/src/gf.cc
+++ b/Carpet/CarpetLib/src/gf.cc
@@ -14,12 +14,11 @@ using namespace std;
 template<typename T>
 gf<T>::gf (const int varindex_, const operator_type transport_operator_,
            th& t_, dh& d_,
-           const int tmin_, const int tmax_,
            const int prolongation_order_time_,
            const int vectorlength_, const int vectorindex_,
            gf* const vectorleader_)
   : ggf(varindex_, transport_operator_,
-        t_, d_, tmin_, tmax_, prolongation_order_time_,
+        t_, d_, prolongation_order_time_,
         vectorlength_, vectorindex_, vectorleader_)
 {
   // recompose ();
@@ -51,21 +50,23 @@ gf<T>::~gf () { }
 
 // Access to the data
 template<typename T>
-const data<T>* gf<T>::operator() (int tl, int rl, int c, int ml) const {
-  assert (tl>=tmin && tl<=tmax);
-  assert (rl>=0 && rl<h.reflevels());
-  assert (c>=0 && c<h.components(rl));
-  assert (ml>=0 && ml<h.mglevels(rl,c));
-  return (const data<T>*)storage.at(tl-tmin).at(rl).at(c).at(ml);
+const data<T>* gf<T>::operator() (int tl, int rl, int c, int ml) const
+{
+  assert (tl>=0 and tl<timelevels());
+  assert (rl>=0 and rl<h.reflevels());
+  assert (c>=0 and c<h.components(rl));
+  assert (ml>=0 and ml<h.mglevels());
+  return (const data<T>*)storage.at(ml).at(rl).at(c).at(tl);
 }
 
 template<typename T>
-data<T>* gf<T>::operator() (int tl, int rl, int c, int ml) {
-  assert (tl>=tmin && tl<=tmax);
-  assert (rl>=0 && rl<h.reflevels());
-  assert (c>=0 && c<h.components(rl));
-  assert (ml>=0 && ml<h.mglevels(rl,c));
-  return (data<T>*)storage.at(tl-tmin).at(rl).at(c).at(ml);
+data<T>* gf<T>::operator() (int tl, int rl, int c, int ml)
+{
+  assert (tl>=0 and tl<timelevels());
+  assert (rl>=0 and rl<h.reflevels());
+  assert (c>=0 and c<h.components(rl));
+  assert (ml>=0 and ml<h.mglevels());
+  return (data<T>*)storage.at(ml).at(rl).at(c).at(tl);
 }
 
 
@@ -76,7 +77,7 @@ ostream& gf<T>::output (ostream& os) const {
   T Tdummy;
   os << "gf<" << typestring(Tdummy) << ">:"
      << varindex << "[" << CCTK_VarName(varindex) << "],"
-     << "dt=[" << tmin << ":" << tmax<< "]";
+     << "tls=" << timelevels();
   return os;
 }
 
diff --git a/Carpet/CarpetLib/src/gf.hh b/Carpet/CarpetLib/src/gf.hh
index f3436d43b..fcfb82b3e 100644
--- a/Carpet/CarpetLib/src/gf.hh
+++ b/Carpet/CarpetLib/src/gf.hh
@@ -35,7 +35,7 @@ public:
   // Constructors
   gf (const int varindex, const operator_type transport_operator,
       th& t, dh& d,
-      const int tmin, const int tmax, const int prolongation_order_time,
+      const int prolongation_order_time,
       const int vectorlength, const int vectorindex,
       gf* const vectorleader);
   
diff --git a/Carpet/CarpetLib/src/ggf.cc b/Carpet/CarpetLib/src/ggf.cc
index c778353b8..7387dd298 100644
--- a/Carpet/CarpetLib/src/ggf.cc
+++ b/Carpet/CarpetLib/src/ggf.cc
@@ -19,24 +19,23 @@ using namespace std;
 // Constructors
 ggf::ggf (const int varindex_, const operator_type transport_operator_,
           th& t_, dh& d_,
-          const int tmin_, const int tmax_,
           const int prolongation_order_time_,
           const int vectorlength_, const int vectorindex_,
           ggf* const vectorleader_)
   : varindex(varindex_), transport_operator(transport_operator_), t(t_),
-    tmin(tmin_), tmax(tmax_),
     prolongation_order_time(prolongation_order_time_),
     h(d_.h), d(d_),
-    storage(tmax-tmin+1),
+    timelevels_(0),
+    storage(h.mglevels()),
     vectorlength(vectorlength_), vectorindex(vectorindex_),
     vectorleader(vectorleader_)
 {
   assert (&t.h == &d.h);
   
   assert (vectorlength >= 1);
-  assert (vectorindex >= 0 && vectorindex < vectorlength);
-  assert ((vectorindex==0 && !vectorleader)
-          || (vectorindex!=0 && vectorleader));
+  assert (vectorindex >= 0 and vectorindex < vectorlength);
+  assert ((vectorindex==0 and !vectorleader)
+          or (vectorindex!=0 and vectorleader));
   
   d.add(this);
 }
@@ -54,20 +53,53 @@ bool ggf::operator== (const ggf& f) const {
 
 
 // Modifiers
+void ggf::set_timelevels (const int ml, const int rl, const int new_timelevels)
+{
+  assert (ml>=0 and ml<(int)storage.size());
+  assert (rl>=0 and rl<(int)storage.at(ml).size());
+  
+  assert (new_timelevels >= 0);
+  
+  if (new_timelevels < timelevels()) {
+    
+    for (int c=0; c<(int)storage.at(ml).at(rl).size(); ++c) {
+      for (int tl=new_timelevels; tl<timelevels(); ++tl) {
+        delete storage.at(ml).at(rl).at(c).at(tl);
+      }
+      storage.at(ml).at(rl).at(c).resize (new_timelevels);
+    } // for c
+    
+  } else if (new_timelevels > timelevels()) {
+    
+    for (int c=0; c<(int)storage.at(ml).at(rl).size(); ++c) {
+      storage.at(ml).at(rl).at(c).resize (new_timelevels);
+      for (int tl=timelevels(); tl<new_timelevels; ++tl) {
+        storage.at(ml).at(rl).at(c).at(tl) = typed_data(tl,rl,c,ml);
+        storage.at(ml).at(rl).at(c).at(tl)->allocate
+          (d.boxes.at(ml).at(rl).at(c).exterior, h.proc(rl,c));
+      } // for tl
+    } // for c
+    
+  }
+  
+  timelevels_ = new_timelevels;
+}
+
+
+
 void ggf::recompose_crop ()
 {
   // Free storage that will not be needed
-  storage.resize(tmax-tmin+1);
-  for (int tl=tmin; tl<=tmax; ++tl) {
-    for (int rl=h.reflevels(); rl<(int)storage.at(tl-tmin).size(); ++rl) {
-      for (int c=0; c<(int)storage.at(tl-tmin).at(rl).size(); ++c) {
-        for (int ml=0; ml<(int)storage.at(tl-tmin).at(rl).at(c).size(); ++ml) {
-          delete storage.at(tl-tmin).at(rl).at(c).at(ml);
-        } // for ml
+  for (int ml=0; ml<h.mglevels(); ++ml) {
+    for (int rl=h.reflevels(); rl<(int)storage.at(ml).size(); ++rl) {
+      for (int c=0; c<(int)storage.at(ml).at(rl).size(); ++c) {
+        for (int tl=0; tl<(int)storage.at(ml).at(rl).at(c).size(); ++tl) {
+          delete storage.at(ml).at(rl).at(c).at(tl);
+        } // for tl
       } // for c
     } // for rl
-    storage.at(tl-tmin).resize(h.reflevels());
-  } // for tl
+    storage.at(ml).resize(h.reflevels());
+  } // for mrl
 }
 
 void ggf::recompose_allocate (const int rl)
@@ -75,58 +107,55 @@ void ggf::recompose_allocate (const int rl)
   // TODO: restructure storage only when needed
   
   // Retain storage that might be needed
-  oldstorage.resize(tmax-tmin+1);
-  for (int tl=tmin; tl<=tmax; ++tl) {
-    oldstorage.at(tl-tmin).resize(h.reflevels());
-    assert (oldstorage.at(tl-tmin).at(rl).size() == 0);
-    oldstorage.at(tl-tmin).at(rl) = storage.at(tl-tmin).at(rl);
-    storage.at(tl-tmin).at(rl).resize(0);
+  oldstorage.resize(storage.size());
+  for (int ml=0; ml<(int)storage.size(); ++ml) {
+    oldstorage.at(ml).resize(storage.at(ml).size());
+    oldstorage.at(ml).at(rl) = storage.at(ml).at(rl);
+    storage.at(ml).at(rl).resize(0);
   }
   
   // Resize structure and allocate storage
-  storage.resize(tmax-tmin+1);
-  for (int tl=tmin; tl<=tmax; ++tl) {
-    storage.at(tl-tmin).resize(h.reflevels());
-    storage.at(tl-tmin).at(rl).resize(h.components(rl));
+  storage.resize(h.mglevels());
+  for (int ml=0; ml<h.mglevels(); ++ml) {
+    storage.at(ml).resize(h.reflevels());
+    storage.at(ml).at(rl).resize(h.components(rl));
     for (int c=0; c<h.components(rl); ++c) {
-      storage.at(tl-tmin).at(rl).at(c).resize(h.mglevels(rl,c));
-      for (int ml=0; ml<h.mglevels(rl,c); ++ml) {
-        storage.at(tl-tmin).at(rl).at(c).at(ml) = typed_data(tl,rl,c,ml);
-        storage.at(tl-tmin).at(rl).at(c).at(ml)->allocate
-          (d.boxes.at(rl).at(c).at(ml).exterior, h.proc(rl,c));
-      } // for ml
+      storage.at(ml).at(rl).at(c).resize(timelevels());
+      for (int tl=0; tl<timelevels(); ++tl) {
+        storage.at(ml).at(rl).at(c).at(tl) = typed_data(tl,rl,c,ml);
+        storage.at(ml).at(rl).at(c).at(tl)->allocate
+          (d.boxes.at(ml).at(rl).at(c).exterior, h.proc(rl,c));
+      } // for tl
     } // for c
-  } // for tl
+  } // for ml
 }
 
 void ggf::recompose_fill (comm_state& state, const int rl,
                           const bool do_prolongate)
 {
   // Initialise the new storage
-  for (int c=0; c<h.components(rl); ++c) {
-    for (int ml=0; ml<h.mglevels(rl,c); ++ml) {
-      for (int tl=tmin; tl<=tmax; ++tl) {
+  for (int ml=0; ml<h.mglevels(); ++ml) {
+    for (int c=0; c<h.components(rl); ++c) {
+      for (int tl=0; tl<timelevels(); ++tl) {
         
         // Find out which regions need to be prolongated
         // (Copy the exterior because some variables are not prolongated)
         // TODO: do this once in the dh instead of for each variable here
-        ibset work (d.boxes.at(rl).at(c).at(ml).exterior);
+        ibset work (d.boxes.at(ml).at(rl).at(c).exterior);
         
         // Copy from old storage, if possible
         // TODO: copy only from interior regions?
-        if (rl<(int)oldstorage.at(tl-tmin).size()) {
-          for (int cc=0; cc<(int)oldstorage.at(tl-tmin).at(rl).size(); ++cc) {
-            if (ml<(int)oldstorage.at(tl-tmin).at(rl).at(cc).size()) {
-              // TODO: prefer same processor, etc., see dh.cc
-              ibset ovlp
-                = work & oldstorage.at(tl-tmin).at(rl).at(cc).at(ml)->extent();
-              ovlp.normalize();
-              work -= ovlp;
-              for (ibset::const_iterator r=ovlp.begin(); r!=ovlp.end(); ++r) {
-                storage.at(tl-tmin).at(rl).at(c).at(ml)->copy_from
-                  (state, oldstorage.at(tl-tmin).at(rl).at(cc).at(ml), *r);
-              }
-            } // if ml
+        if (rl<(int)oldstorage.at(ml).size()) {
+          for (int cc=0; cc<(int)oldstorage.at(ml).at(rl).size(); ++cc) {
+            // TODO: prefer same processor, etc., see dh.cc
+            ibset ovlp
+              = work & oldstorage.at(ml).at(rl).at(cc).at(tl)->extent();
+            ovlp.normalize();
+            work -= ovlp;
+            for (ibset::const_iterator r=ovlp.begin(); r!=ovlp.end(); ++r) {
+              storage.at(ml).at(rl).at(c).at(tl)->copy_from
+                (state, oldstorage.at(ml).at(rl).at(cc).at(tl), *r);
+            }
           } // for cc
         } // if rl
         
@@ -135,18 +164,17 @@ void ggf::recompose_fill (comm_state& state, const int rl,
           if (rl>0) {
             if (transport_operator != op_none) {
               const int numtl = prolongation_order_time+1;
-              assert (tmax-tmin+1 >= numtl);
+              assert (timelevels() >= numtl);
               vector<int> tls(numtl);
               vector<CCTK_REAL> times(numtl);
               for (int i=0; i<numtl; ++i) {
-                tls.at(i) = tmax - i;
+                tls.at(i) = i;
                 times.at(i) = t.time(tls.at(i),rl-1,ml);
               }
-              for (int cc=0; cc<(int)storage.at(tl-tmin).at(rl-1).size(); ++cc) {
+              for (int cc=0; cc<(int)storage.at(ml).at(rl-1).size(); ++cc) {
                 vector<const gdata*> gsrcs(numtl);
                 for (int i=0; i<numtl; ++i) {
-                  gsrcs.at(i)
-                    = storage.at(tls.at(i)-tmin).at(rl-1).at(cc).at(ml);
+                  gsrcs.at(i) = storage.at(ml).at(rl-1).at(cc).at(tls.at(i));
                   assert (gsrcs.at(i)->extent() == gsrcs.at(0)->extent());
                 }
                 const CCTK_REAL time = t.time(tl,rl,ml);
@@ -154,13 +182,17 @@ void ggf::recompose_fill (comm_state& state, const int rl,
                 // TODO: choose larger regions first
                 // TODO: prefer regions from the same processor
                 const iblist& list
-                  = d.boxes.at(rl).at(c).at(ml).recv_ref_coarse.at(cc);
-                for (iblist::const_iterator iter=list.begin(); iter!=list.end(); ++iter) {
+                  = d.boxes.at(ml).at(rl).at(c).recv_ref_coarse.at(cc);
+                for (iblist::const_iterator iter=list.begin();
+                     iter!=list.end(); ++iter)
+                {
                   ibset ovlp = work & *iter;
                   ovlp.normalize();
                   work -= ovlp;
-                  for (ibset::const_iterator r=ovlp.begin(); r!=ovlp.end(); ++r) {
-                    storage.at(tl-tmin).at(rl).at(c).at(ml)->interpolate_from
+                  for (ibset::const_iterator r=ovlp.begin();
+                       r!=ovlp.end(); ++r)
+                  {
+                    storage.at(ml).at(rl).at(c).at(tl)->interpolate_from
                       (state, gsrcs, times, *r, time,
                        d.prolongation_order_space, prolongation_order_time);
                   } // for r
@@ -176,58 +208,58 @@ void ggf::recompose_fill (comm_state& state, const int rl,
         // TODO: check this.
         
       } // for tl
-    } // for ml
-  } // for c
+    } // for c
+  } // for ml
 }
 
 void ggf::recompose_free (const int rl)
 {
   // Delete old storage
-  for (int tl=tmin; tl<=tmax; ++tl) {
-    for (int c=0; c<(int)oldstorage.at(tl-tmin).at(rl).size(); ++c) {
-      for (int ml=0; ml<(int)oldstorage.at(tl-tmin).at(rl).at(c).size(); ++ml) {
-        delete oldstorage.at(tl-tmin).at(rl).at(c).at(ml);
-      } // for ml
+  for (int ml=0; ml<(int)oldstorage.size(); ++ml) {
+    for (int c=0; c<(int)oldstorage.at(ml).at(rl).size(); ++c) {
+      for (int tl=0; tl<timelevels(); ++tl) {
+        delete oldstorage.at(ml).at(rl).at(c).at(tl);
+      } // for tl
     } // for c
-    oldstorage.at(tl-tmin).at(rl).resize(0);
-  } // for tl
+    oldstorage.at(ml).at(rl).resize(0);
+  } // for ml
 }
 
 void ggf::recompose_bnd_prolongate (comm_state& state, const int rl,
-                                       const bool do_prolongate)
+                                    const bool do_prolongate)
 {
   if (do_prolongate) {
     // Set boundaries
     if (rl>0) {
-      for (int c=0; c<h.components(rl); ++c) {
-        for (int ml=0; ml<h.mglevels(rl,c); ++ml) {
-          for (int tl=tmin; tl<=tmax; ++tl) {
+      for (int ml=0; ml<h.mglevels(); ++ml) {
+        for (int c=0; c<h.components(rl); ++c) {
+          for (int tl=0; tl<timelevels(); ++tl) {
             
             // TODO: assert that reflevel 0 boundaries are copied
             const CCTK_REAL time = t.time(tl,rl,ml);
             ref_bnd_prolongate (state,tl,rl,c,ml,time);
             
           } // for tl
-        } // for ml
-      } // for c
+        } // for c
+      } // for ml
     } // if rl
   } // if do_prolongate
 }
 
 void ggf::recompose_sync (comm_state& state, const int rl,
-                             const bool do_prolongate)
+                          const bool do_prolongate)
 {
   if (do_prolongate) {
     // Set boundaries
-    for (int c=0; c<h.components(rl); ++c) {
-      for (int ml=0; ml<h.mglevels(rl,c); ++ml) {
-        for (int tl=tmin; tl<=tmax; ++tl) {
+    for (int ml=0; ml<h.mglevels(); ++ml) {
+      for (int c=0; c<h.components(rl); ++c) {
+        for (int tl=0; tl<timelevels(); ++tl) {
           
           sync (state,tl,rl,c,ml);
           
         } // for tl
-      } // for ml
-    } // for c
+      } // for c
+    } // for ml
   } // if do_prolongate
 }
 
@@ -235,28 +267,28 @@ void ggf::recompose_sync (comm_state& state, const int rl,
 
 // Cycle the time levels by rotating the data sets
 void ggf::cycle (int rl, int c, int ml) {
-  assert (rl>=0 && rl<h.reflevels());
-  assert (c>=0 && c<h.components(rl));
-  assert (ml>=0 && ml<h.mglevels(rl,c));
-  gdata* tmpdata = storage.at(tmin-tmin).at(rl).at(c).at(ml);
-  for (int tl=tmin; tl<=tmax-1; ++tl) {
-    storage.at(tl-tmin).at(rl).at(c).at(ml) = storage.at(tl+1-tmin).at(rl).at(c).at(ml);
+  assert (rl>=0 and rl<h.reflevels());
+  assert (c>=0 and c<h.components(rl));
+  assert (ml>=0 and ml<h.mglevels());
+  gdata* tmpdata = storage.at(ml).at(rl).at(c).at(timelevels()-1);
+  for (int tl=timelevels()-1; tl>0; --tl) {
+    storage.at(ml).at(rl).at(c).at(tl) = storage.at(ml).at(rl).at(c).at(tl-1);
   }
-  storage.at(tmax-tmin).at(rl).at(c).at(ml) = tmpdata;
+  storage.at(ml).at(rl).at(c).at(0) = tmpdata;
 }
 
 // Flip the time levels by exchanging the data sets
 void ggf::flip (int rl, int c, int ml) {
-  assert (rl>=0 && rl<h.reflevels());
-  assert (c>=0 && c<h.components(rl));
-  assert (ml>=0 && ml<h.mglevels(rl,c));
-  for (int tl=0; tl<(tmax-tmin)/2; ++tl) {
-    const int tl1 = tmin + tl;
-    const int tl2 = tmax - tl;
+  assert (rl>=0 and rl<h.reflevels());
+  assert (c>=0 and c<h.components(rl));
+  assert (ml>=0 and ml<h.mglevels());
+  for (int tl=0; tl<(timelevels()-1)/2; ++tl) {
+    const int tl1 =                  tl;
+    const int tl2 = timelevels()-1 - tl;
     assert (tl1 < tl2);
-    gdata* tmpdata = storage.at(tl1-tmin).at(rl).at(c).at(ml);
-    storage.at(tl1-tmin).at(rl).at(c).at(ml) = storage.at(tl2-tmin).at(rl).at(c).at(ml);
-    storage.at(tl2-tmin).at(rl).at(c).at(ml) = tmpdata;
+    gdata* tmpdata = storage.at(ml).at(rl).at(c).at(tl1);
+    storage.at(ml).at(rl).at(c).at(tl1) = storage.at(ml).at(rl).at(c).at(tl2);
+    storage.at(ml).at(rl).at(c).at(tl2) = tmpdata;
   }
 }
 
@@ -271,21 +303,21 @@ void ggf::copycat (comm_state& state,
                    int tl2, int rl2, int ml2,
                    const ibbox dh::dboxes::* send_box)
 {
-  assert (tl1>=tmin && tl1<=tmax);
-  assert (rl1>=0 && rl1<h.reflevels());
-  assert (c1>=0 && c1<h.components(rl1));
-  assert (ml1>=0 && ml1<h.mglevels(rl1,c1));
-  assert (tl2>=tmin && tl2<=tmax);
-  assert (rl2>=0 && rl2<h.reflevels());
+  assert (tl1>=0 and tl1<timelevels());
+  assert (rl1>=0 and rl1<h.reflevels());
+  assert (c1>=0 and c1<h.components(rl1));
+  assert (ml1>=0 and ml1<h.mglevels());
+  assert (tl2>=0 and tl2<timelevels());
+  assert (rl2>=0 and rl2<h.reflevels());
   const int c2=c1;
-  assert (ml2<h.mglevels(rl2,c2));
-  const ibbox recv = d.boxes.at(rl1).at(c1).at(ml1).*recv_box;
-  const ibbox send = d.boxes.at(rl2).at(c2).at(ml2).*send_box;
+  assert (ml2<h.mglevels());
+  const ibbox recv = d.boxes.at(ml1).at(rl1).at(c1).*recv_box;
+  const ibbox send = d.boxes.at(ml2).at(rl2).at(c2).*send_box;
   assert (all(recv.shape()==send.shape()));
   // copy the content
   assert (recv==send);
-  storage.at(tl1-tmin).at(rl1).at(c1).at(ml1)->copy_from
-    (state, storage.at(tl2-tmin).at(rl2).at(c2).at(ml2), recv);
+  storage.at(ml1).at(rl1).at(c1).at(tl1)->copy_from
+    (state, storage.at(ml2).at(rl2).at(c2).at(tl2), recv);
 }
 
 // Copy regions
@@ -295,24 +327,24 @@ void ggf::copycat (comm_state& state,
                    int tl2, int rl2, int ml2,
                    const iblist dh::dboxes::* send_list)
 {
-  assert (tl1>=tmin && tl1<=tmax);
-  assert (rl1>=0 && rl1<h.reflevels());
-  assert (c1>=0 && c1<h.components(rl1));
-  assert (ml1>=0 && ml1<h.mglevels(rl1,c1));
-  assert (tl2>=tmin && tl2<=tmax);
-  assert (rl2>=0 && rl2<h.reflevels());
+  assert (tl1>=0 and tl1<timelevels());
+  assert (rl1>=0 and rl1<h.reflevels());
+  assert (c1>=0 and c1<h.components(rl1));
+  assert (ml1>=0 and ml1<h.mglevels());
+  assert (tl2>=0 and tl2<timelevels());
+  assert (rl2>=0 and rl2<h.reflevels());
   const int c2=c1;
-  assert (ml2<h.mglevels(rl2,c2));
-  const iblist recv = d.boxes.at(rl1).at(c1).at(ml1).*recv_list;
-  const iblist send = d.boxes.at(rl2).at(c2).at(ml2).*send_list;
+  assert (ml2<h.mglevels());
+  const iblist recv = d.boxes.at(ml1).at(rl1).at(c1).*recv_list;
+  const iblist send = d.boxes.at(ml2).at(rl2).at(c2).*send_list;
   assert (recv.size()==send.size());
   // walk all boxes
   for (iblist::const_iterator r=recv.begin(), s=send.begin();
        r!=recv.end(); ++r, ++s) {
     // (use the send boxes for communication)
     // copy the content
-    storage.at(tl1-tmin).at(rl1).at(c1).at(ml1)->copy_from
-      (state, storage.at(tl2-tmin).at(rl2).at(c2).at(ml2), *r);
+    storage.at(ml1).at(rl1).at(c1).at(tl1)->copy_from
+      (state, storage.at(ml2).at(rl2).at(c2).at(tl2), *r);
   }
 }
 
@@ -323,25 +355,25 @@ void ggf::copycat (comm_state& state,
                    int tl2, int rl2, int ml2,
                    const iblistvect dh::dboxes::* send_listvect)
 {
-  assert (tl1>=tmin && tl1<=tmax);
-  assert (rl1>=0 && rl1<h.reflevels());
-  assert (c1>=0 && c1<h.components(rl1));
-  assert (ml1>=0 && ml1<h.mglevels(rl1,c1));
-  assert (tl2>=tmin && tl2<=tmax);
-  assert (rl2>=0 && rl2<h.reflevels());
+  assert (tl1>=0 and tl1<timelevels());
+  assert (rl1>=0 and rl1<h.reflevels());
+  assert (c1>=0 and c1<h.components(rl1));
+  assert (ml1>=0 and ml1<h.mglevels());
+  assert (tl2>=0 and tl2<timelevels());
+  assert (rl2>=0 and rl2<h.reflevels());
   // walk all components
   for (int c2=0; c2<h.components(rl2); ++c2) {
-    assert (ml2<h.mglevels(rl2,c2));
-    const iblist recv = (d.boxes.at(rl1).at(c1).at(ml1).*recv_listvect).at(c2);
-    const iblist send = (d.boxes.at(rl2).at(c2).at(ml2).*send_listvect).at(c1);
+    assert (ml2<h.mglevels());
+    const iblist recv = (d.boxes.at(ml1).at(rl1).at(c1).*recv_listvect).at(c2);
+    const iblist send = (d.boxes.at(ml2).at(rl2).at(c2).*send_listvect).at(c1);
     assert (recv.size()==send.size());
     // walk all boxes
     for (iblist::const_iterator r=recv.begin(), s=send.begin();
       	 r!=recv.end(); ++r, ++s) {
       // (use the send boxes for communication)
       // copy the content
-      storage.at(tl1-tmin).at(rl1).at(c1).at(ml1)->copy_from
-      	(state, storage.at(tl2-tmin).at(rl2).at(c2).at(ml2), *r);
+      storage.at(ml1).at(rl1).at(c1).at(tl1)->copy_from
+      	(state, storage.at(ml2).at(rl2).at(c2).at(tl2), *r);
     }
   }
 }
@@ -354,33 +386,30 @@ void ggf::intercat (comm_state& state,
                     const ibbox dh::dboxes::* send_list,
                     CCTK_REAL time)
 {
-  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 (tl1>=0 and tl1<timelevels());
+  assert (rl1>=0 and rl1<h.reflevels());
+  assert (c1>=0 and c1<h.components(rl1));
+  assert (ml1>=0 and ml1<h.mglevels());
   for (int i=0; i<(int)tl2s.size(); ++i) {
-    assert (tl2s.at(i)>=tmin && tl2s.at(i)<=tmax);
+    assert (tl2s.at(i)>=0 and tl2s.at(i)<timelevels());
   }
   assert (rl2>=0 and rl2<h.reflevels());
   const int c2=c1;
-  assert (ml2>=0 && ml2<h.mglevels(rl2,c2));
+  assert (ml2>=0 and ml2<h.mglevels());
   
   vector<const gdata*> gsrcs(tl2s.size());
   vector<CCTK_REAL> times(tl2s.size());
   for (int i=0; i<(int)gsrcs.size(); ++i) {
-    assert (rl2<(int)storage.at(tl2s.at(i)-tmin).size());
-    assert (c2<(int)storage.at(tl2s.at(i)-tmin).at(rl2).size());
-    assert (ml2<(int)storage.at(tl2s.at(i)-tmin).at(rl2).at(c2).size());
-    gsrcs.at(i) = storage.at(tl2s.at(i)-tmin).at(rl2).at(c2).at(ml2);
+    gsrcs.at(i) = storage.at(ml2).at(rl2).at(c2).at(tl2s.at(i));
     times.at(i) = t.time(tl2s.at(i),rl2,ml2);
   }
   
-  const ibbox recv = d.boxes.at(rl1).at(c1).at(ml1).*recv_list;
-  const ibbox send = d.boxes.at(rl2).at(c2).at(ml2).*send_list;
+  const ibbox recv = d.boxes.at(ml1).at(rl1).at(c1).*recv_list;
+  const ibbox send = d.boxes.at(ml2).at(rl2).at(c2).*send_list;
   assert (all(recv.shape()==send.shape()));
   // interpolate the content
   assert (recv==send);
-  storage.at(tl1-tmin).at(rl1).at(c1).at(ml1)->interpolate_from
+  storage.at(ml1).at(rl1).at(c1).at(tl1)->interpolate_from
     (state, gsrcs, times, recv, time,
      d.prolongation_order_space, prolongation_order_time);
 }
@@ -393,36 +422,33 @@ void ggf::intercat (comm_state& state,
                     const iblist dh::dboxes::* send_list,
                     const CCTK_REAL time)
 {
-  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 (tl1>=0 and tl1<timelevels());
+  assert (rl1>=0 and rl1<h.reflevels());
+  assert (c1>=0 and c1<h.components(rl1));
+  assert (ml1>=0 and ml1<h.mglevels());
   for (int i=0; i<(int)tl2s.size(); ++i) {
-    assert (tl2s.at(i)>=tmin && tl2s.at(i)<=tmax);
+    assert (tl2s.at(i)>=0 and tl2s.at(i)<timelevels());
   }
   assert (rl2>=0 and rl2<h.reflevels());
   const int c2=c1;
-  assert (ml2>=0 && ml2<h.mglevels(rl2,c2));
+  assert (ml2>=0 and ml2<h.mglevels());
   
   vector<const gdata*> gsrcs(tl2s.size());
   vector<CCTK_REAL> times(tl2s.size());
   for (int i=0; i<(int)gsrcs.size(); ++i) {
-    assert (rl2<(int)storage.at(tl2s.at(i)-tmin).size());
-    assert (c2<(int)storage.at(tl2s.at(i)-tmin).at(rl2).size());
-    assert (ml2<(int)storage.at(tl2s.at(i)-tmin).at(rl2).at(c2).size());
-    gsrcs.at(i) = storage.at(tl2s.at(i)-tmin).at(rl2).at(c2).at(ml2);
+    gsrcs.at(i) = storage.at(ml2).at(rl2).at(c2).at(tl2s.at(i));
     times.at(i) = t.time(tl2s.at(i),rl2,ml2);
   }
   
-  const iblist recv = d.boxes.at(rl1).at(c1).at(ml1).*recv_list;
-  const iblist send = d.boxes.at(rl2).at(c2).at(ml2).*send_list;
+  const iblist recv = d.boxes.at(ml1).at(rl1).at(c1).*recv_list;
+  const iblist send = d.boxes.at(ml2).at(rl2).at(c2).*send_list;
   assert (recv.size()==send.size());
   // walk all boxes
   for (iblist::const_iterator r=recv.begin(), s=send.begin();
        r!=recv.end(); ++r, ++s) {
     // (use the send boxes for communication)
     // interpolate the content
-    storage.at(tl1-tmin).at(rl1).at(c1).at(ml1)->interpolate_from
+    storage.at(ml1).at(rl1).at(c1).at(tl1)->interpolate_from
       (state, gsrcs, times, *r, time,
        d.prolongation_order_space, prolongation_order_time);
   }
@@ -436,37 +462,34 @@ void ggf::intercat (comm_state& state,
                     const iblistvect dh::dboxes::* send_listvect,
                     const CCTK_REAL time)
 {
-  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 (tl1>=0 and tl1<timelevels());
+  assert (rl1>=0 and rl1<h.reflevels());
+  assert (c1>=0 and c1<h.components(rl1));
+  assert (ml1>=0 and ml1<h.mglevels());
   for (int i=0; i<(int)tl2s.size(); ++i) {
-    assert (tl2s.at(i)>=tmin && tl2s.at(i)<=tmax);
+    assert (tl2s.at(i)>=0 and tl2s.at(i)<timelevels());
   }
   assert (rl2>=0 and rl2<h.reflevels());
   // walk all components
   for (int c2=0; c2<h.components(rl2); ++c2) {
-    assert (ml2>=0 && ml2<h.mglevels(rl2,c2));
+    assert (ml2>=0 and ml2<h.mglevels());
     
     vector<const gdata*> gsrcs(tl2s.size());
     vector<CCTK_REAL> times(tl2s.size());
     for (int i=0; i<(int)gsrcs.size(); ++i) {
-      assert (rl2<(int)storage.at(tl2s.at(i)-tmin).size());
-      assert (c2<(int)storage.at(tl2s.at(i)-tmin).at(rl2).size());
-      assert (ml2<(int)storage.at(tl2s.at(i)-tmin).at(rl2).at(c2).size());
-      gsrcs.at(i) = storage.at(tl2s.at(i)-tmin).at(rl2).at(c2).at(ml2);
+      gsrcs.at(i) = storage.at(ml2).at(rl2).at(c2).at(tl2s.at(i));
       times.at(i) = t.time(tl2s.at(i),rl2,ml2);
     }
     
-    const iblist recv = (d.boxes.at(rl1).at(c1).at(ml1).*recv_listvect).at(c2);
-    const iblist send = (d.boxes.at(rl2).at(c2).at(ml2).*send_listvect).at(c1);
+    const iblist recv = (d.boxes.at(ml1).at(rl1).at(c1).*recv_listvect).at(c2);
+    const iblist send = (d.boxes.at(ml2).at(rl2).at(c2).*send_listvect).at(c1);
     assert (recv.size()==send.size());
     // walk all boxes
     for (iblist::const_iterator r=recv.begin(), s=send.begin();
       	 r!=recv.end(); ++r, ++s) {
       // (use the send boxes for communication)
       // interpolate the content
-      storage.at(tl1-tmin).at(rl1).at(c1).at(ml1)->interpolate_from
+      storage.at(ml1).at(rl1).at(c1).at(tl1)->interpolate_from
       	(state, gsrcs, times, *r, time,
 	 d.prolongation_order_space, prolongation_order_time);
     }
@@ -503,9 +526,9 @@ void ggf::ref_bnd_prolongate (comm_state& state,
   if (transport_operator == op_none) return;
   vector<int> tl2s;
   // Interpolation in time
-  assert (tmax-tmin+1 >= prolongation_order_time+1);
+  assert (timelevels() >= prolongation_order_time+1);
   tl2s.resize(prolongation_order_time+1);
-  for (int i=0; i<=prolongation_order_time; ++i) tl2s.at(i) = tmax - i;
+  for (int i=0; i<=prolongation_order_time; ++i) tl2s.at(i) = i;
   intercat (state,
             tl  ,rl  ,c,ml, &dh::dboxes::recv_ref_bnd_coarse,
 	    tl2s,rl-1,  ml, &dh::dboxes::send_ref_bnd_fine,
@@ -567,9 +590,9 @@ void ggf::ref_prolongate (comm_state& state,
   if (transport_operator == op_none) return;
   vector<int> tl2s;
   // Interpolation in time
-  assert (tmax-tmin+1 >= prolongation_order_time+1);
+  assert (timelevels() >= prolongation_order_time+1);
   tl2s.resize(prolongation_order_time+1);
-  for (int i=0; i<=prolongation_order_time; ++i) tl2s.at(i) = tmax - i;
+  for (int i=0; i<=prolongation_order_time; ++i) tl2s.at(i) = i;
   intercat (state,
             tl  ,rl  ,c,ml, &dh::dboxes::recv_ref_coarse,
 	    tl2s,rl-1,  ml, &dh::dboxes::send_ref_fine,
diff --git a/Carpet/CarpetLib/src/ggf.hh b/Carpet/CarpetLib/src/ggf.hh
index f1ac54baa..0fcaf960f 100644
--- a/Carpet/CarpetLib/src/ggf.hh
+++ b/Carpet/CarpetLib/src/ggf.hh
@@ -34,10 +34,10 @@ class ggf {
   typedef vector<iblist> iblistvect;
   
   typedef gdata*        tdata;  // data ...
-  typedef vector<tdata> mdata;  // ... for each multigrid level
-  typedef vector<mdata> cdata;  // ... for each component
+  typedef vector<tdata> fdata;  // ... for each time level
+  typedef vector<fdata> cdata;  // ... for each component
   typedef vector<cdata> rdata;  // ... for each refinement level
-  typedef vector<rdata> fdata;  // ... for each time level
+  typedef vector<rdata> mdata;  // ... for each multigrid level
   
 public:				// should be readonly
   
@@ -45,15 +45,17 @@ public:				// should be readonly
   const int varindex;           // Cactus variable index
   const operator_type transport_operator;
   
-  const th &t;               // time hierarchy
-  const int tmin, tmax;		// timelevels
+  const th &t;                  // time hierarchy
   const int prolongation_order_time; // order of temporal prolongation operator
   
-  const gh &h;               // grid hierarchy
+  const gh &h;                  // grid hierarchy
   dh &d;			// data hierarchy
 
+private:
+  int timelevels_;              // time levels
+
 protected:
-  fdata storage;		// storage
+  mdata storage;		// storage
   
 public:
   const int vectorlength;       // vector length
@@ -61,14 +63,13 @@ public:
   const ggf* vectorleader;      // first vector element
   
 private:
-  fdata oldstorage;
+  mdata oldstorage;            // temporary storage
   
 public:
 
   // Constructors
   ggf (const int varindex, const operator_type transport_operator,
        th& t, dh& d,
-       const int tmin, const int tmax,
        const int prolongation_order_time,
        const int vectorlength, const int vectorindex,
        ggf* const vectorleader);
@@ -78,11 +79,18 @@ public:
 
   // Comparison
   bool operator== (const ggf& f) const;
-
-
-
+  
+  // Querying
+  int timelevels () const
+  {
+    return timelevels_;
+  }
+  
+  
+  
   // Modifiers
-  // void recompose ();
+  void set_timelevels (int ml, int rl, int new_timelevels);
+
   void recompose_crop ();
   void recompose_allocate (int rl);
   void recompose_fill (comm_state& state, int rl, bool do_prolongate);
diff --git a/Carpet/CarpetLib/src/gh.cc b/Carpet/CarpetLib/src/gh.cc
index c73bc9dd3..7b5f97c81 100644
--- a/Carpet/CarpetLib/src/gh.cc
+++ b/Carpet/CarpetLib/src/gh.cc
@@ -29,7 +29,7 @@ gh::gh (const int reffact_, const centering refcent_,
 gh::~gh () { }
 
 // Modifiers
-void gh::recompose (const rexts& exts,
+void gh::recompose (const mexts& exts,
                     const rbnds& outer_bounds,
                     const rprocs& procs,
                     const bool do_prolongate)
@@ -70,45 +70,47 @@ void gh::recompose (const rexts& exts,
 
 void gh::check_processor_number_consistency () {
   for (int rl=0; rl<reflevels(); ++rl) {
-    assert (processors().size() == extents().size());
-    assert (outer_boundaries().size() == extents().size());
+    assert (processors().size() == extents().at(0).size());
+    assert (outer_boundaries().size() == extents().at(0).size());
     for (int c=0; c<components(rl); ++c) {
-      assert (processors().at(rl).size() == extents().at(rl).size());
-      assert (outer_boundaries().at(rl).size() == extents().at(rl).size());
+      assert (processors().at(rl).size() == extents().at(0).at(rl).size());
+      assert (outer_boundaries().at(rl).size() == extents().at(0).at(rl).size());
     }
   }
 }
   
-void gh::check_multigrid_consistency () {
-  for (int rl=0; rl<reflevels(); ++rl) {
-    for (int c=0; c<components(rl); ++c) {
-      assert (mglevels(rl,c)>0);
-      for (int ml=1; ml<mglevels(rl,c); ++ml) {
-	assert (all(extents().at(rl).at(c).at(ml).stride()
-		    == ivect(mgfact) * extents().at(rl).at(c).at(ml-1).stride()));
+void gh::check_multigrid_consistency ()
+{
+  assert (mglevels()>0);
+  for (int ml=1; ml<mglevels(); ++ml) {
+    for (int rl=0; rl<reflevels(); ++rl) {
+      for (int c=0; c<components(rl); ++c) {
+	assert (all(extents().at(ml).at(rl).at(c).stride()
+		    == ivect(mgfact) * extents().at(ml-1).at(rl).at(c).stride()));
         // TODO: put the check back in, taking outer boundaries into
         // account
 #if 0
-	assert (extents().at(rl).at(c).at(ml)
-		.contracted_for(extents().at(rl).at(c).at(ml-1))
-		.is_contained_in(extents().at(rl).at(c).at(ml-1)));
+	assert (extents().at(ml).at(rl).at(c)
+		.contracted_for(extents().at(ml-1).at(rl).at(c))
+		.is_contained_in(extents().at(ml-1).at(rl).at(c)));
 #endif
       }
     }
   }
 }
   
-void gh::check_component_consistency () {
-  for (int rl=0; rl<reflevels(); ++rl) {
-    assert (components(rl)>0);
-    for (int c=0; c<components(rl); ++c) {
-      for (int ml=0; ml<mglevels(rl,c); ++ml) {
-        const ibbox &b  = extents().at(rl).at(c).at(ml);
-        const ibbox &b0 = extents().at(rl).at(0).at(ml);
+void gh::check_component_consistency ()
+{
+  for (int ml=0; ml<mglevels(); ++ml) {
+    for (int rl=0; rl<reflevels(); ++rl) {
+      assert (components(rl)>0);
+      for (int c=0; c<components(rl); ++c) {
+        const ibbox &b  = extents().at(ml).at(rl).at(c);
+        const ibbox &b0 = extents().at(ml).at(rl).at(0);
 	assert (all(b.stride() == b0.stride()));
 	assert (b.is_aligned_with(b0));
         for (int cc=c+1; cc<components(rl); ++cc) {
-          assert ((b & extents().at(rl).at(cc).at(ml)).empty());
+          assert ((b & extents().at(ml).at(rl).at(cc)).empty());
         }
       }
     }
@@ -127,40 +129,40 @@ void gh::check_base_grid_extent () {
   }
 }
   
-void gh::check_refinement_levels () {
-  for (int rl=1; rl<reflevels(); ++rl) {
-    assert (all(extents().at(rl-1).at(0).at(0).stride()
-		== ivect(reffact) * extents().at(rl).at(0).at(0).stride()));
-    // Check contained-ness:
-    // first take all coarse grids ...
-    ibset all;
-    for (int c=0; c<components(rl-1); ++c) {
-      all |= extents().at(rl-1).at(c).at(0);
-    }
-    // ... remember their size ...
-    const int sz = all.size();
-    // ... then add the coarsified fine grids ...
-    for (int c=0; c<components(rl); ++c) {
-      all |= extents().at(rl).at(c).at(0).contracted_for(extents().at(rl-1).at(0).at(0));
+void gh::check_refinement_levels ()
+{
+  for (int ml=0; ml<mglevels(); ++ml) {
+    for (int rl=1; rl<reflevels(); ++rl) {
+      assert (all(extents().at(ml).at(rl-1).at(0).stride()
+                  == ivect(reffact) * extents().at(ml).at(rl).at(0).stride()));
+      // Check contained-ness:
+      // first take all coarse grids ...
+      ibset all;
+      for (int c=0; c<components(rl-1); ++c) {
+        all |= extents().at(ml).at(rl-1).at(c);
+      }
+      // ... remember their size ...
+      const int sz = all.size();
+      // ... then add the coarsified fine grids ...
+      for (int c=0; c<components(rl); ++c) {
+        all |= extents().at(ml).at(rl).at(c).contracted_for(extents().at(ml).at(rl-1).at(0));
+      }
+      // ... and then check the sizes:
+      assert (all.size() == sz);
     }
-    // ... and then check the sizes:
-    assert (all.size() == sz);
   }
 }
 
-void gh::calculate_base_extents_of_all_levels () {
-  _bases.resize(reflevels());
-  for (int rl=0; rl<reflevels(); ++rl) {
-    if (components(rl)==0) {
-      _bases.at(rl).resize(0);
-    } else {
-      _bases.at(rl).resize(mglevels(rl,0));
-      for (int ml=0; ml<mglevels(rl,0); ++ml) {
-        _bases.at(rl).at(ml) = ibbox();
-        ibbox &bb = _bases.at(rl).at(ml);
-	for (int c=0; c<components(rl); ++c) {
-	  bb = bb.expanded_containing(extents().at(rl).at(c).at(ml));
-	}
+void gh::calculate_base_extents_of_all_levels ()
+{
+  _bases.resize(mglevels());
+  for (int ml=0; ml<mglevels(); ++ml) {
+    _bases.at(ml).resize(reflevels());
+    for (int rl=0; rl<reflevels(); ++rl) {
+      _bases.at(ml).at(rl) = ibbox();
+      ibbox &bb = _bases.at(ml).at(rl);
+      for (int c=0; c<components(rl); ++c) {
+        bb = bb.expanded_containing(extents().at(ml).at(rl).at(c));
       }
     }
   }
@@ -198,14 +200,15 @@ void gh::remove (dh* d) {
 }
 
 
-void gh::do_output_bboxes (ostream& os) const {
-  for (int rl=0; rl<reflevels(); ++rl) {
-    for (int c=0; c<components(rl); ++c) {
-      for (int ml=0; ml<mglevels(rl,c); ++ml) {
+void gh::do_output_bboxes (ostream& os) const
+{
+  for (int ml=0; ml<mglevels(); ++ml) {
+    for (int rl=0; rl<reflevels(); ++rl) {
+      for (int c=0; c<components(rl); ++c) {
         os << endl;
         os << "gh bboxes:" << endl;
-        os << "rl=" << rl << " c=" << c << " ml=" << ml << endl;
-        os << "extent=" << extents().at(rl).at(c).at(ml) << endl;
+        os << "ml=" << ml << " rl=" << rl << " c=" << c << endl;
+        os << "extent=" << extents().at(ml).at(rl).at(c) << endl;
         os << "outer_boundary=" << outer_boundaries().at(rl).at(c) << endl;
         os << "processor=" << processors().at(rl).at(c) << endl;
       }
@@ -213,14 +216,15 @@ void gh::do_output_bboxes (ostream& os) const {
   }
 }
 
-void gh::do_output_bases (ostream& os) const {
-  for (int rl=0; rl<reflevels(); ++rl) {
-    if (components(rl)>0) {
-      for (int ml=0; ml<mglevels(rl,0); ++ml) {
+void gh::do_output_bases (ostream& os) const
+{
+  for (int ml=0; ml<mglevels(); ++ml) {
+    for (int rl=0; rl<reflevels(); ++rl) {
+      if (components(rl)>0) {
         os << endl;
         os << "gh bases:" << endl;
-        os << "rl=" << rl << " ml=" << ml << endl;
-        os << "base=" << bases().at(rl).at(ml) << endl;
+        os << "ml=" << ml << " rl=" << rl << endl;
+        os << "base=" << bases().at(ml).at(rl) << endl;
       }
     }
   }
diff --git a/Carpet/CarpetLib/src/gh.hh b/Carpet/CarpetLib/src/gh.hh
index 6a7bb2e83..3b7f0a2c8 100644
--- a/Carpet/CarpetLib/src/gh.hh
+++ b/Carpet/CarpetLib/src/gh.hh
@@ -33,9 +33,9 @@ class gh {
 public:
   
   // Types
-  typedef vector<ibbox> mexts;	// ... for each multigrid level
-  typedef vector<mexts> cexts;	// ... for each component
+  typedef vector<ibbox> cexts;	// ... for each component
   typedef vector<cexts> rexts;	// ... for each refinement level
+  typedef vector<rexts> mexts;	// ... for each multigrid level
   
   typedef vector<bbvect> cbnds;	// ... for each component
   typedef vector<cbnds> rbnds;	// ... for each refinement level
@@ -56,9 +56,9 @@ public:				// should be readonly
   
   
 private:
-  vector<vector<ibbox> > _bases; // [rl][ml]
+  vector<vector<ibbox> > _bases; // [ml][rl]
   // TODO: invent structure for this
-  rexts _extents;		// extents of all grids
+  mexts _extents;		// extents of all grids
   rbnds _outer_boundaries;	// boundary descriptions of all grids
   rprocs _processors;		// processor numbers of all grids
 
@@ -76,11 +76,12 @@ public:
   virtual ~gh ();
   
   // Modifiers
-  void recompose (const rexts& exts, const rbnds& outer_bounds,
+  void recompose (const mexts& exts, const rbnds& outer_bounds,
 		  const rprocs& procs,
                   const bool do_prolongate);
 
-  const rexts & extents() const {
+  const mexts & extents() const
+  {
     return _extents;
   }
 
@@ -97,16 +98,20 @@ public:
   }
   
   // Accessors
-  int reflevels () const {
+  int mglevels () const
+  {
     return (int)_extents.size();
   }
   
-  int components (const int rl) const {
-    return (int)_extents.at(rl).size();
+  int reflevels () const
+  {
+    if (mglevels() == 0) return 0;
+    return (int)_extents.at(0).size();
   }
   
-  int mglevels (const int rl, const int c) const {
-    return (int)_extents.at(rl).at(c).size();
+  int components (const int rl) const
+  {
+    return (int)_extents.at(0).at(rl).size();
   }
   
   bbvect outer_boundary (const int rl, const int c) const {
diff --git a/Carpet/CarpetLib/src/th.cc b/Carpet/CarpetLib/src/th.cc
index 18559ae75..77415fffc 100644
--- a/Carpet/CarpetLib/src/th.cc
+++ b/Carpet/CarpetLib/src/th.cc
@@ -25,29 +25,29 @@ th::~th () {
 }
 
 // Modifiers
-void th::recompose () {
-  times.resize(h.reflevels());
-  deltas.resize(h.reflevels());
-  for (int rl=0; rl<h.reflevels(); ++rl) {
-    const int old_mglevels = times.at(rl).size();
-    CCTK_REAL mgtime;
-    // Select default time
-    if (old_mglevels==0 && rl==0) {
-      mgtime = 0;
-    } else if (old_mglevels==0) {
-      mgtime = times.at(rl-1).at(0);
-    } else {
-      mgtime = times.at(rl).at(old_mglevels-1);
-    }
-    times.at(rl).resize(h.mglevels(rl,0), mgtime);
-    deltas.at(rl).resize(h.mglevels(rl,0));
-    for (int ml=0; ml<h.mglevels(rl,0); ++ml) {
-      if (rl==0 && ml==0) {
-	deltas.at(rl).at(ml) = delta;
+void th::recompose ()
+{
+  const int old_mglevels = times.size();
+  times.resize(h.mglevels());
+  deltas.resize(h.mglevels());
+  for (int ml=0; ml<h.mglevels(); ++ml) {
+    const int old_reflevels = times.at(ml).size();
+    times.at(ml).resize(h.reflevels());
+    deltas.at(ml).resize(h.reflevels());
+    for (int rl=0; rl<h.reflevels(); ++rl) {
+      if (old_mglevels==0) {
+        times.at(ml).at(rl) = 0;
+      } else if (rl < old_reflevels) {
+        // do nothing
+      } else {
+        times.at(ml).at(rl) = times.at(ml).at(rl-1);
+      }
+      if (ml==0 and rl==0) {
+	deltas.at(ml).at(rl) = delta;
       } else if (ml==0) {
-	deltas.at(rl).at(ml) = deltas.at(rl-1).at(ml) / h.reffact;
+	deltas.at(ml).at(rl) = deltas.at(ml).at(rl-1) / h.reffact;
       } else {
-	deltas.at(rl).at(ml) = deltas.at(rl).at(ml-1) * h.mgfact;
+	deltas.at(ml).at(rl) = deltas.at(ml-1).at(rl) * h.mgfact;
       }
     }
   }
@@ -59,11 +59,13 @@ void th::recompose () {
 void th::output (ostream& os) const {
   os << "th:"
      << "times={";
-  for (int rl=0; rl<h.reflevels(); ++rl) {
-    for (int ml=0; ml<h.mglevels(rl,0); ++ml) {
-      if (!(rl==0 && ml==0)) os << ",";
-      os << rl << ":" << ml << ":"
-	 << times.at(rl).at(ml) << "(" << deltas.at(rl).at(ml) << ")";
+  const char * sep = "";
+  for (int ml=0; ml<h.mglevels(); ++ml) {
+    for (int rl=0; rl<h.reflevels(); ++rl) {
+      os << sep
+         << ml << ":" << rl << ":"
+	 << times.at(ml).at(rl) << "(" << deltas.at(ml).at(rl) << ")";
+      sep = ",";
     }
   }
   os << "}";
diff --git a/Carpet/CarpetLib/src/th.hh b/Carpet/CarpetLib/src/th.hh
index 6ef1308aa..42bd6504e 100644
--- a/Carpet/CarpetLib/src/th.hh
+++ b/Carpet/CarpetLib/src/th.hh
@@ -33,8 +33,8 @@ public:				// should be readonly
 private:
   
   CCTK_REAL delta;		// time step
-  vector<vector<CCTK_REAL> > times; // current times
-  vector<vector<CCTK_REAL> > deltas; // time steps
+  vector<vector<CCTK_REAL> > times; // current times [ml][rl]
+  vector<vector<CCTK_REAL> > deltas; // time steps [ml][rl]
   
 public:
   
@@ -48,38 +48,44 @@ public:
   void recompose ();
   
   // Time management
-  CCTK_REAL get_time (const int rl, const int ml) const {
-    assert (rl>=0 && rl<h.reflevels());
-    assert (ml>=0 && ml<h.mglevels(rl,0));
-    return times.at(rl).at(ml);
+  CCTK_REAL get_time (const int rl, const int ml) const
+  {
+    assert (rl>=0 and rl<h.reflevels());
+    assert (ml>=0 and ml<h.mglevels());
+    return times.at(ml).at(rl);
   }
   
-  void set_time (const int rl, const int ml, const CCTK_REAL t) {
-    assert (rl>=0 && rl<h.reflevels());
-    assert (ml>=0 && ml<h.mglevels(rl,0));
-    times.at(rl).at(ml) = t;
+  void set_time (const int rl, const int ml, const CCTK_REAL t)
+  {
+    assert (rl>=0 and rl<h.reflevels());
+    assert (ml>=0 and ml<h.mglevels());
+    times.at(ml).at(rl) = t;
   }
   
-  void advance_time (const int rl, const int ml) {
+  void advance_time (const int rl, const int ml)
+  {
     set_time(rl,ml, get_time(rl,ml) + get_delta(rl,ml));
   }
   
-  CCTK_REAL get_delta (const int rl, const int ml) const {
-    assert (rl>=0 && rl<h.reflevels());
-    assert (ml>=0 && ml<h.mglevels(rl,0));
-    return deltas.at(rl).at(ml);
+  CCTK_REAL get_delta (const int rl, const int ml) const
+  {
+    assert (rl>=0 and rl<h.reflevels());
+    assert (ml>=0 and ml<h.mglevels());
+    return deltas.at(ml).at(rl);
   }
   
-  void set_delta (const int rl, const int ml, const CCTK_REAL dt) {
-    assert (rl>=0 && rl<h.reflevels());
-    assert (ml>=0 && ml<h.mglevels(rl,0));
-    deltas.at(rl).at(ml) = dt;
+  void set_delta (const int rl, const int ml, const CCTK_REAL dt)
+  {
+    assert (rl>=0 and rl<h.reflevels());
+    assert (ml>=0 and ml<h.mglevels());
+    deltas.at(ml).at(rl) = dt;
   }
   
-  CCTK_REAL time (const int tl, const int rl, const int ml) const {
-    assert (rl>=0 && rl<h.reflevels());
-    assert (ml>=0 && ml<h.mglevels(rl,0));
-    return get_time(rl, ml) + tl * get_delta(rl, ml);
+  CCTK_REAL time (const int tl, const int rl, const int ml) const
+  {
+    assert (rl>=0 and rl<h.reflevels());
+    assert (ml>=0 and ml<h.mglevels());
+    return get_time(rl, ml) - tl * get_delta(rl, ml);
   }
   
   // Output
diff --git a/Carpet/CarpetReduce/src/mask_carpet.cc b/Carpet/CarpetReduce/src/mask_carpet.cc
index c068a2fde..342c3901f 100644
--- a/Carpet/CarpetReduce/src/mask_carpet.cc
+++ b/Carpet/CarpetReduce/src/mask_carpet.cc
@@ -45,22 +45,22 @@ namespace CarpetMask {
       gh const & hh = *vhh.at(Carpet::map);
       dh const & dd = *vdd.at(Carpet::map);
       
-      ibbox const & base = hh.bases().at(reflevel).at(mglevel);
+      ibbox const & base = hh.bases().at(mglevel).at(reflevel);
       
       
       
       // Calculate the union of all refined regions
       ibset refined;
       for (int c=0; c<hh.components(reflevel); ++c) {
-        refined |= hh.extents().at(reflevel).at(c).at(mglevel);
+        refined |= hh.extents().at(mglevel).at(reflevel).at(c);
       }
       refined.normalize();
       
       // Calculate the union of all coarse regions
       ibset parent;
       for (int c=0; c<hh.components(reflevel-1); ++c) {
-//         parent |= hh.extents().at(reflevel-1).at(c).at(mglevel).expanded_for(base);
-        parent |= hh.extents().at(reflevel-1).at(c).at(mglevel).expand(ivect(reffact-1),ivect(reffact-1)).contracted_for(base);
+//         parent |= hh.extents().at(mglevel).at(reflevel-1).at(c).expanded_for(base);
+        parent |= hh.extents().at(mglevel).at(reflevel-1).at(c).expand(ivect(reffact-1),ivect(reffact-1)).contracted_for(base);
       }
       parent.normalize();
       
@@ -102,7 +102,7 @@ namespace CarpetMask {
           DECLARE_CCTK_ARGUMENTS;
           
           ibbox const & ext
-            = dd.boxes.at(reflevel).at(component).at(mglevel).exterior;
+            = dd.boxes.at(mglevel).at(reflevel).at(component).exterior;
           
           for (int d=0; d<dim; ++d) {
             for (ibset::const_iterator bi = boundaries[d].begin();
@@ -164,7 +164,7 @@ namespace CarpetMask {
           DECLARE_CCTK_ARGUMENTS;
           
           ibbox const & ext
-            = dd.boxes.at(reflevel).at(component).at(mglevel).exterior;
+            = dd.boxes.at(mglevel).at(reflevel).at(component).exterior;
           
           for (ibset::const_iterator bi = refined.begin();
                bi != refined.end();
diff --git a/Carpet/CarpetRegrid/interface.ccl b/Carpet/CarpetRegrid/interface.ccl
index 7c38c0efb..6d8643c63 100644
--- a/Carpet/CarpetRegrid/interface.ccl
+++ b/Carpet/CarpetRegrid/interface.ccl
@@ -52,7 +52,7 @@ USES FUNCTION ConvertFromPhysicalBoundary
 
 # The true prototype of the routine below:
 #   int Carpet_Regrid (const cGH * cctkGH,
-#                      gh::rexts  * bbsss,
+#                      gh::mexts  * bbsss,
 #                      gh::rbnds  * obss,
 #                      gh::rprocs * pss);
 CCTK_INT FUNCTION Carpet_Regrid (CCTK_POINTER_TO_CONST IN cctkGH,  \
diff --git a/Carpet/CarpetRegrid/src/automatic.cc b/Carpet/CarpetRegrid/src/automatic.cc
index 2e29195e5..abe685d47 100644
--- a/Carpet/CarpetRegrid/src/automatic.cc
+++ b/Carpet/CarpetRegrid/src/automatic.cc
@@ -26,7 +26,7 @@ namespace CarpetRegrid {
   
   int Automatic (cGH const * const cctkGH,
                  gh const & hh,
-                 gh::rexts  & bbsss,
+                 gh::mexts  & bbsss,
                  gh::rbnds  & obss,
                  gh::rprocs & pss)
   {
@@ -35,6 +35,7 @@ namespace CarpetRegrid {
     assert (refinement_levels >= 1);
     
     assert (bbsss.size() >= 1);
+    vector<vector<ibbox> > bbss = bbsss.at(0);
     
     
     
@@ -68,32 +69,29 @@ namespace CarpetRegrid {
     gh::cprocs ps;
     SplitRegions (cctkGH, bbs, obs, ps);
     
-    // make multigrid aware
-    vector<vector<ibbox> > bbss;
-    MakeMultigridBoxes (cctkGH, bbs, obs, bbss);
-    
-    
-    
     if (bbss.size() == 0) {
       // remove all finer levels
-      bbsss.resize(reflevel+1);
+      bbss.resize(reflevel+1);
       obss.resize(reflevel+1);
       pss.resize(reflevel+1);
     } else {
-      assert (reflevel < (int)bbsss.size());
-      if (reflevel+1 == (int)bbsss.size()) {
+      assert (reflevel < (int)bbss.size());
+      if (reflevel+1 == (int)bbss.size()) {
 	// add a finer level
-	bbsss.push_back (bbss);
+	bbss.push_back (bbs);
 	obss.push_back (obs);
 	pss.push_back (ps);
       } else {
 	// change a finer level
-	bbsss.at(reflevel+1) = bbss;
+	bbss.at(reflevel+1) = bbs;
 	obss.at(reflevel+1) = obs;
 	pss.at(reflevel+1) = ps;
       }
     }
     
+    // make multigrid aware
+    MakeMultigridBoxes (cctkGH, bbss, obss, bbsss);
+    
     return 1;
   }
   
diff --git a/Carpet/CarpetRegrid/src/baselevel.cc b/Carpet/CarpetRegrid/src/baselevel.cc
index ce4aee1c2..a2579d8ec 100644
--- a/Carpet/CarpetRegrid/src/baselevel.cc
+++ b/Carpet/CarpetRegrid/src/baselevel.cc
@@ -19,7 +19,7 @@ namespace CarpetRegrid {
   
   int BaseLevel (cGH const * const cctkGH,
                  gh const & hh,
-                 gh::rexts  & bbsss,
+                 gh::mexts  & bbsss,
                  gh::rbnds  & obss,
                  gh::rprocs & pss)
   {
diff --git a/Carpet/CarpetRegrid/src/centre.cc b/Carpet/CarpetRegrid/src/centre.cc
index 2c5dd6ce0..2bcbda6eb 100644
--- a/Carpet/CarpetRegrid/src/centre.cc
+++ b/Carpet/CarpetRegrid/src/centre.cc
@@ -19,7 +19,7 @@ namespace CarpetRegrid {
   
   int Centre (cGH const * const cctkGH,
               gh const & hh,
-              gh::rexts  & bbsss,
+              gh::mexts  & bbsss,
               gh::rbnds  & obss,
               gh::rprocs & pss)
   {
@@ -31,8 +31,9 @@ namespace CarpetRegrid {
     if (reflevel == refinement_levels) return 0;
     
     assert (bbsss.size() >= 1);
+    vector<vector<ibbox> > bbss = bbsss.at(0);
     
-    bbsss.resize (refinement_levels);
+    bbss.resize (refinement_levels);
     obss.resize (refinement_levels);
     pss.resize (refinement_levels);
     
@@ -45,7 +46,7 @@ namespace CarpetRegrid {
     
     assert (! smart_outer_boundaries);
     
-    for (size_t rl=1; rl<bbsss.size(); ++rl) {
+    for (size_t rl=1; rl<bbss.size(); ++rl) {
       
       // save old values
       ivect const oldrlb = rlb;
@@ -74,14 +75,13 @@ namespace CarpetRegrid {
       gh::cprocs ps;
       SplitRegions (cctkGH, bbs, obs, ps);
       
-      // make multigrid aware
-      vector<vector<ibbox> > bbss;
-      MakeMultigridBoxes (cctkGH, bbs, obs, bbss);
-      
-      bbsss.at(rl) = bbss;
+      bbss.at(rl) = bbs;
       obss.at(rl) = obs;
       pss.at(rl) = ps;
       
+      // make multigrid aware
+      MakeMultigridBoxes (cctkGH, bbss, obss, bbsss);
+      
     } // for rl
     
     return 1;
diff --git a/Carpet/CarpetRegrid/src/manualcoordinatelist.cc b/Carpet/CarpetRegrid/src/manualcoordinatelist.cc
index 04dffd5e5..2467bae24 100644
--- a/Carpet/CarpetRegrid/src/manualcoordinatelist.cc
+++ b/Carpet/CarpetRegrid/src/manualcoordinatelist.cc
@@ -23,7 +23,7 @@ namespace CarpetRegrid {
   
   int ManualCoordinateList (cGH const * const cctkGH,
                             gh const & hh,
-                            gh::rexts  & bbsss,
+                            gh::mexts  & bbsss,
                             gh::rbnds  & obss,
                             gh::rprocs & pss)
   {
@@ -36,6 +36,7 @@ namespace CarpetRegrid {
     if (reflevel == refinement_levels) return 0;
     
     assert (bbsss.size() >= 1);
+    vector<vector<ibbox> > bbss = bbsss.at(0);
     
     jjvect nboundaryzones, is_internal, is_staggered, shiftout;
     ierr = GetBoundarySpecification
@@ -52,7 +53,7 @@ namespace CarpetRegrid {
        &exterior_min[0], &exterior_max[0], &base_spacing[0]);
     assert (!ierr);
     
-    bbsss.resize (refinement_levels);
+    bbss.resize (refinement_levels);
     obss.resize (refinement_levels);
     pss.resize (refinement_levels);
     
@@ -183,14 +184,13 @@ namespace CarpetRegrid {
       gh::cprocs ps;
       SplitRegions (cctkGH, bbs, obs, ps);
       
-      // make multigrid aware
-      vector<vector<ibbox> > bbss;
-      MakeMultigridBoxes (cctkGH, bbs, obs, bbss);
-      
-      bbsss.at(rl) = bbss;
+      bbss.at(rl) = bbs;
       obss.at(rl) = obs;
       pss.at(rl) = ps;
       
+      // make multigrid aware
+      MakeMultigridBoxes (cctkGH, bbss, obss, bbsss);
+      
     } // for rl
     
     return 1;
diff --git a/Carpet/CarpetRegrid/src/manualcoordinates.cc b/Carpet/CarpetRegrid/src/manualcoordinates.cc
index cc757fdce..4495f7ac9 100644
--- a/Carpet/CarpetRegrid/src/manualcoordinates.cc
+++ b/Carpet/CarpetRegrid/src/manualcoordinates.cc
@@ -21,7 +21,7 @@ namespace CarpetRegrid {
   
   int ManualCoordinates (cGH const * const cctkGH,
                          gh const & hh,
-                         gh::rexts  & bbsss,
+                         gh::mexts  & bbsss,
                          gh::rbnds  & obss,
                          gh::rprocs & pss)
   {
@@ -36,8 +36,9 @@ namespace CarpetRegrid {
     if (reflevel == refinement_levels) return 0;
     
     assert (bbsss.size() >= 1);
+    vector<vector<ibbox> > bbss = bbsss.at(0);
     
-    bbsss.resize (refinement_levels);
+    bbss.resize (refinement_levels);
     obss.resize (refinement_levels);
     pss.resize (refinement_levels);
     
@@ -51,7 +52,7 @@ namespace CarpetRegrid {
     
     assert (! smart_outer_boundaries);
     
-    for (size_t rl=1; rl<bbsss.size(); ++rl) {
+    for (size_t rl=1; rl<bbss.size(); ++rl) {
       
       bbvect const ob (false);
       
@@ -66,14 +67,13 @@ namespace CarpetRegrid {
       gh::cprocs ps;
       SplitRegions (cctkGH, bbs, obs, ps);
       
-      // make multigrid aware
-      vector<vector<ibbox> > bbss;
-      MakeMultigridBoxes (cctkGH, bbs, obs, bbss);
-      
-      bbsss.at(rl) = bbss;
+      bbss.at(rl) = bbs;
       obss.at(rl) = obs;
       pss.at(rl) = ps;
       
+      // make multigrid aware
+      MakeMultigridBoxes (cctkGH, bbss, obss, bbsss);
+      
     } // for rl
     
     return 1;
diff --git a/Carpet/CarpetRegrid/src/manualgridpointlist.cc b/Carpet/CarpetRegrid/src/manualgridpointlist.cc
index c4f534fd1..93dfe6563 100644
--- a/Carpet/CarpetRegrid/src/manualgridpointlist.cc
+++ b/Carpet/CarpetRegrid/src/manualgridpointlist.cc
@@ -23,7 +23,7 @@ namespace CarpetRegrid {
   
   int ManualGridpointList (cGH const * const cctkGH,
                            gh const & hh,
-                           gh::rexts  & bbsss,
+                           gh::mexts  & bbsss,
                            gh::rbnds  & obss,
                            gh::rprocs & pss)
   {
@@ -35,8 +35,9 @@ namespace CarpetRegrid {
     if (reflevel == refinement_levels) return 0;
     
     assert (bbsss.size() >= 1);
+    vector<vector<ibbox> > bbss = bbsss.at(0);
     
-    bbsss.resize (refinement_levels);
+    bbss.resize (refinement_levels);
     obss.resize (refinement_levels);
     pss.resize (refinement_levels);
     
@@ -104,14 +105,13 @@ namespace CarpetRegrid {
       gh::cprocs ps;
       SplitRegions (cctkGH, bbs, obs, ps);
       
-      // make multigrid aware
-      vector<vector<ibbox> > bbss;
-      MakeMultigridBoxes (cctkGH, bbs, obs, bbss);
-      
-      bbsss.at(rl) = bbss;
+      bbss.at(rl) = bbs;
       obss.at(rl) = obs;
       pss.at(rl) = ps;
       
+      // make multigrid aware
+      MakeMultigridBoxes (cctkGH, bbss, obss, bbsss);
+      
     } // for rl
     
     return 1;
diff --git a/Carpet/CarpetRegrid/src/manualgridpoints.cc b/Carpet/CarpetRegrid/src/manualgridpoints.cc
index b4cb85a17..92432686c 100644
--- a/Carpet/CarpetRegrid/src/manualgridpoints.cc
+++ b/Carpet/CarpetRegrid/src/manualgridpoints.cc
@@ -22,7 +22,7 @@ namespace CarpetRegrid {
   
   int ManualGridpoints (cGH const * const cctkGH,
                         gh const & hh,
-                        gh::rexts  & bbsss,
+                        gh::mexts  & bbsss,
                         gh::rbnds  & obss,
                         gh::rprocs & pss)
   {
@@ -37,8 +37,9 @@ namespace CarpetRegrid {
     if (reflevel == refinement_levels) return 0;
     
     assert (bbsss.size() >= 1);
+    vector<vector<ibbox> > bbss = bbsss.at(0);
     
-    bbsss.resize (refinement_levels);
+    bbss.resize (refinement_levels);
     obss.resize (refinement_levels);
     pss.resize (refinement_levels);
     
@@ -52,7 +53,7 @@ namespace CarpetRegrid {
     
     assert (! smart_outer_boundaries);
     
-    for (size_t rl=1; rl<bbsss.size(); ++rl) {
+    for (size_t rl=1; rl<bbss.size(); ++rl) {
       
       bbvect const ob (false);
       
@@ -67,14 +68,13 @@ namespace CarpetRegrid {
       gh::cprocs ps;
       SplitRegions (cctkGH, bbs, obs, ps);
       
-      // make multigrid aware
-      vector<vector<ibbox> > bbss;
-      MakeMultigridBoxes (cctkGH, bbs, obs, bbss);
-      
-      bbsss.at(rl) = bbss;
+      bbss.at(rl) = bbs;
       obss.at(rl) = obs;
       pss.at(rl) = ps;
       
+      // make multigrid aware
+      MakeMultigridBoxes (cctkGH, bbss, obss, bbsss);
+      
     } // for rl
     
     return 1;
diff --git a/Carpet/CarpetRegrid/src/moving.cc b/Carpet/CarpetRegrid/src/moving.cc
index 0985399d5..7936865c3 100644
--- a/Carpet/CarpetRegrid/src/moving.cc
+++ b/Carpet/CarpetRegrid/src/moving.cc
@@ -20,7 +20,7 @@ namespace CarpetRegrid {
   
   int Moving (cGH const * const cctkGH,
               gh const & hh,
-              gh::rexts  & bbsss,
+              gh::mexts  & bbsss,
               gh::rbnds  & obss,
               gh::rprocs & pss)
   {
@@ -30,8 +30,9 @@ namespace CarpetRegrid {
     assert (refinement_levels >= 1);
     
     assert (bbsss.size() >= 1);
+    vector<vector<ibbox> > bbss = bbsss.at(0);
     
-    bbsss.resize (refinement_levels);
+    bbss.resize (refinement_levels);
     obss.resize (refinement_levels);
     pss.resize (refinement_levels);
     
@@ -40,7 +41,7 @@ namespace CarpetRegrid {
     
     assert (! smart_outer_boundaries);
     
-    for (size_t rl=1; rl<bbsss.size(); ++rl) {
+    for (size_t rl=1; rl<bbss.size(); ++rl) {
       
       // calculate new extent
       CCTK_REAL const argument = 2*M_PI * moving_circle_frequency * cctk_time;
@@ -63,14 +64,13 @@ namespace CarpetRegrid {
       gh::cprocs ps;
       SplitRegions (cctkGH, bbs, obs, ps);
       
-      // make multigrid aware
-      vector<vector<ibbox> > bbss;
-      MakeMultigridBoxes (cctkGH, bbs, obs, bbss);
-      
-      bbsss.at(rl) = bbss;
+      bbss.at(rl) = bbs;
       obss.at(rl) = obs;
       pss.at(rl) = ps;
       
+      // make multigrid aware
+      MakeMultigridBoxes (cctkGH, bbss, obss, bbsss);
+      
     } // for rl
     
     return 1;
diff --git a/Carpet/CarpetRegrid/src/regrid.cc b/Carpet/CarpetRegrid/src/regrid.cc
index 94b281513..92b234d24 100644
--- a/Carpet/CarpetRegrid/src/regrid.cc
+++ b/Carpet/CarpetRegrid/src/regrid.cc
@@ -28,7 +28,7 @@ namespace CarpetRegrid {
     
     const cGH * const cctkGH = (const cGH *) cctkGH_;
     
-    gh::rexts  & bbsss = * (gh::rexts  *) bbsss_;
+    gh::mexts  & bbsss = * (gh::mexts  *) bbsss_;
     gh::rbnds  & obss  = * (gh::rbnds  *) obss_;
     gh::rprocs & pss   = * (gh::rprocs *) pss_;
     
diff --git a/Carpet/CarpetRegrid/src/regrid.hh b/Carpet/CarpetRegrid/src/regrid.hh
index fea47e073..30038d8b2 100644
--- a/Carpet/CarpetRegrid/src/regrid.hh
+++ b/Carpet/CarpetRegrid/src/regrid.hh
@@ -30,7 +30,7 @@ namespace CarpetRegrid {
     
     /* Aliased functions */
 //     CCTK_INT CarpetRegrid_Regrid (const cGH * const cctkGH,
-//                                   gh<dim>::rexts  * bbsss,
+//                                   gh<dim>::mexts  * bbsss,
 //                                   gh<dim>::rbnds  * obss,
 //                                   gh<dim>::rprocs * pss);
     CCTK_INT CarpetRegrid_Regrid (CCTK_POINTER_TO_CONST const cctkGH_,
@@ -44,19 +44,19 @@ namespace CarpetRegrid {
   
   int BaseLevel (cGH const * const cctkGH,
                  gh const & hh,
-                 gh::rexts  & bbsss,
+                 gh::mexts  & bbsss,
                  gh::rbnds  & obss,
                  gh::rprocs & pss);
 
   int Centre (cGH const * const cctkGH,
               gh const & hh,
-              gh::rexts  & bbsss,
+              gh::mexts  & bbsss,
               gh::rbnds  & obss,
               gh::rprocs & pss);
 
   int ManualGridpoints (cGH const * const cctkGH,
                         gh const & hh,
-                        gh::rexts  & bbsss,
+                        gh::mexts  & bbsss,
                         gh::rbnds  & obss,
                         gh::rprocs & pss);
   
@@ -72,7 +72,7 @@ namespace CarpetRegrid {
 
   int ManualCoordinates (cGH const * const cctkGH,
                          gh const & hh,
-                         gh::rexts  & bbsss,
+                         gh::mexts  & bbsss,
                          gh::rbnds  & obss,
                          gh::rprocs & pss);
   
@@ -97,25 +97,25 @@ namespace CarpetRegrid {
 
   int ManualGridpointList (cGH const * const cctkGH,
                            gh const & hh,
-                           gh::rexts  & bbsss,
+                           gh::mexts  & bbsss,
                            gh::rbnds  & obss,
                            gh::rprocs & pss);
   
   int ManualCoordinateList (cGH const * const cctkGH,
                             gh const & hh,
-                            gh::rexts  & bbsss,
+                            gh::mexts  & bbsss,
                             gh::rbnds  & obss,
                             gh::rprocs & pss);
 
   int Moving (cGH const * const cctkGH,
               gh const & hh,
-              gh::rexts  & bbsss,
+              gh::mexts  & bbsss,
               gh::rbnds  & obss,
               gh::rprocs & pss);
   
   int Automatic (cGH const * const cctkGH,
                  gh const & hh,
-                 gh::rexts  & bbsss,
+                 gh::mexts  & bbsss,
                  gh::rbnds  & obss,
                  gh::rprocs & pss);
   
diff --git a/Carpet/CarpetSlab/src/GetHyperslab.cc b/Carpet/CarpetSlab/src/GetHyperslab.cc
index c22e0b1b4..1c9e687cc 100644
--- a/Carpet/CarpetSlab/src/GetHyperslab.cc
+++ b/Carpet/CarpetSlab/src/GetHyperslab.cc
@@ -187,8 +187,8 @@ namespace CarpetSlab {
         
         // Calculate overlapping extents
         const bboxset<int,dim> myextents
-          = ((mydd->boxes.at(rl).at(component).at(mglevel).sync_not
-              | mydd->boxes.at(rl).at(component).at(mglevel).interior)
+          = ((mydd->boxes.at(mglevel).at(rl).at(component).sync_not
+              | mydd->boxes.at(mglevel).at(rl).at(component).interior)
              & hextent);
         
         // Loop over overlapping extents
diff --git a/Carpet/CarpetSlab/src/slab.cc b/Carpet/CarpetSlab/src/slab.cc
index 9ca04ae6a..ac8ea0ef9 100644
--- a/Carpet/CarpetSlab/src/slab.cc
+++ b/Carpet/CarpetSlab/src/slab.cc
@@ -219,8 +219,8 @@ namespace CarpetSlab {
         
         // Calculate overlapping extents
         const bboxset<int,dim> myextents
-          = ((mydd->boxes.at(rl).at(component).at(mglevel).sync_not
-              | mydd->boxes.at(rl).at(component).at(mglevel).interior)
+          = ((mydd->boxes.at(mglevel).at(rl).at(component).sync_not
+              | mydd->boxes.at(mglevel).at(rl).at(component).interior)
              & hextent);
         
         // Loop over overlapping extents