CarpetAdaptiveRegrid

The first part of CarpetAdaptiveRegrid. At the moment the only piece
of "working" code is the F90 that does the box manipulation. This can
be tested using the standalone makefile in src/test.

darcs-hash:20041113152820-5d2e4-7bbcb3804e5643c7753c9ab663df4dc7d357d7e9.gz

11 files changed, 1961 insertions, 0 deletions

diff --git a/CarpetDev/CarpetAdaptiveRegrid/src/CAR.cc b/CarpetDev/CarpetAdaptiveRegrid/src/CAR.cc
new file mode 100644
index 000000000..e76d20716
--- /dev/null
+++ b/CarpetDev/CarpetAdaptiveRegrid/src/CAR.cc
@@ -0,0 +1,138 @@
+#include <assert.h>
+
+#include <sstream>
+#include <string>
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+
+#include "gh.hh"
+#include "vect.hh"
+
+#include "carpet.hh"
+#include "CAR.hh"
+
+extern "C" {
+  static const char* rcsid = "$Header:$";
+  CCTK_FILEVERSION(Carpet_CarpetAdaptiveregrid_regrid_cc);
+}
+
+
+
+namespace CarpetAdaptiveRegrid {
+  
+  using namespace std;
+  using namespace Carpet;
+
+  extern "C" {
+    void CCTK_FCALL CCTK_FNAME(copy_mask)
+      (const CCTK_INT& snx, const CCTK_INT& sny, const CCTK_INT& snz,
+       const CCTK_INT* smask, const CCTK_INT sbbox[3][3],
+       const CCTK_INT& dnx, const CCTK_INT& dny, const CCTK_INT& dnz,
+       CCTK_INT* dmask, const CCTK_INT dbbox[3][3]);
+    void CCTK_FCALL CCTK_FNAME(check_box)
+      (const CCTK_INT& nx, const CCTK_INT& ny, const CCTK_INT& nz,
+       const CCTK_INT* mask,
+       CCTK_INT* sum_x, CCTK_INT* sum_y, CCTK_INT* sum_z,
+       CCTK_INT* sig_x, CCTK_INT* sig_y, CCTK_INT* sig_z,
+       const CCTK_INT bbox[3][3],
+       CCTK_INT newbbox1[3][3], CCTK_INT newbbox2[3][3],
+       const CCTK_INT& min_width, const CCTK_REAL& min_density,
+       CCTK_INT& didit);
+  }
+  
+  CCTK_INT CarpetAdaptiveRegrid_Regrid (CCTK_POINTER_TO_CONST const cctkGH_,
+                                CCTK_POINTER const bbsss_,
+                                CCTK_POINTER const obss_,
+                                CCTK_POINTER const pss_,
+				CCTK_INT force)
+  {
+    DECLARE_CCTK_PARAMETERS;
+    
+    const cGH * const cctkGH = (const cGH *) cctkGH_;
+    
+    gh<dim>::rexts  & bbsss = * (gh<dim>::rexts  *) bbsss_;
+    gh<dim>::rbnds  & obss  = * (gh<dim>::rbnds  *) obss_;
+    gh<dim>::rprocs & pss   = * (gh<dim>::rprocs *) pss_;
+    
+    gh<dim> const & hh = *vhh.at(Carpet::map);
+    
+    assert (is_singlemap_mode());
+    
+    // In force mode (force == true) we do not check the
+    // CarpetAdaptiveregrid parameters
+
+    if (!force) {
+
+      assert (regrid_every == -1 || regrid_every == 0
+	      || regrid_every % maxmglevelfact == 0);
+    
+      // Return if no regridding is desired
+      if (regrid_every == -1) return 0;
+      
+      // Return if we want to regrid during initial data only, and this
+      // is not the time for initial data
+      if (regrid_every == 0 && cctkGH->cctk_iteration != 0) return 0;
+
+      // Return if we want to regrid regularly, but not at this time
+      if (regrid_every > 0 && cctkGH->cctk_iteration != 0
+	  && (cctkGH->cctk_iteration-1) % regrid_every != 0)
+      {
+	return 0;
+      }
+
+    }    
+
+    int do_recompose;
+
+    // Make a cboxlist (or set) from the boxes on the level.
+
+    // cboxlist is list of cbox's.
+
+    // cbox is a bbox 
+    // plus a dim-D array of booleans mask(i,j,k)
+    // plus 3 sum arrays (sigma_i = sum_{j,k} mask(i,:,:) etc)
+    // plus 3 signature arrays (1D laplacian of sigma arrays)
+
+    // First thing to do; make a cbox containing all active boxes
+    // on this level.
+    // cbox list contains this cbox.
+    // Create the mask for this cbox by looking at the error GF.
+
+    // Then recursively for every cbox in the list until no changes:
+
+    //   Prune the cbox (remove all zeros from the ends)
+
+    //   Split the cbox (find zero in sigma_:, split along that line)
+
+    //   if density too low then find zero crossings of signature and
+    //   split there
+
+    // if none of these steps are taken the cbox is finished.
+
+    // Most of the actual work is done by the utility functions in
+    // CAR_utils.F77.
+    // What we have to do is
+
+    // loop over every box in the list
+    //   for each box call check_box
+    //   if it returns zero the box is done
+    //   if it returns one then the box needs shrinking
+    //   if it returns two then the box needs splitting
+
+    // so, the method will be:
+
+    // Create a cboxlist with one member
+
+    // loop over the list:
+    //   call check_box
+    //   if the return is two then create a new box from newbbox2 and add it
+    //     to the end of the list
+    //   if the return value is one or two shrink the current box to the
+    //     values given in newbbox1
+    //   if the return value is one or two go redo this loop again.
+    
+    return do_recompose;
+  }
+  
+} // namespace CarpetAdaptiveRegrid
diff --git a/CarpetDev/CarpetAdaptiveRegrid/src/CAR.h b/CarpetDev/CarpetAdaptiveRegrid/src/CAR.h
new file mode 100644
index 000000000..cf48c4073
--- /dev/null
+++ b/CarpetDev/CarpetAdaptiveRegrid/src/CAR.h
@@ -0,0 +1,21 @@
+/* $Header:$ */
+
+#ifndef CARPETADAPTIVEREGRID_H
+#define CARPETADAPTIVEREGRID_H
+
+#include "cctk_Arguments.h"
+
+#ifdef __cplusplus
+namespace CarpetAdaptiveRegrid {
+  extern "C" {
+#endif
+    
+    /* Scheduled functions */
+    int CarpetAdaptiveregridParamcheck (CCTK_ARGUMENTS);
+
+#ifdef __cplusplus
+  } /* extern "C" */
+} /* namespace CarpetAdaptiveregrid */
+#endif
+
+#endif /* !defined(CARPETADAPTIVEREGRID_H) */
diff --git a/CarpetDev/CarpetAdaptiveRegrid/src/CAR.hh b/CarpetDev/CarpetAdaptiveRegrid/src/CAR.hh
new file mode 100644
index 000000000..d5371068c
--- /dev/null
+++ b/CarpetDev/CarpetAdaptiveRegrid/src/CAR.hh
@@ -0,0 +1,46 @@
+// $Header:$
+
+#ifndef CARPETADAPTIVEREGRID_HH
+#define CARPETADAPTIVEREGRID_HH
+
+#include <list>
+
+#include "cctk.h"
+#include "cctk_Arguments.h"
+
+#include "bbox.hh"
+#include "gf.hh"
+#include "gh.hh"
+#include "vect.hh"
+
+#include "carpet.hh"
+
+
+
+namespace CarpetAdaptiveRegrid {
+  
+  using namespace std;
+  using namespace Carpet;
+  
+  
+  
+  extern "C" {
+    
+    /* Scheduled functions */
+    int CarpetAdaptiveRegridParamcheck (CCTK_ARGUMENTS);
+    
+    /* Aliased functions */
+//     CCTK_INT CarpetAdaptiveRegrid_Regrid (const cGH * const cctkGH,
+//                                   gh<dim>::rexts  * bbsss,
+//                                   gh<dim>::rbnds  * obss,
+//                                   gh<dim>::rprocs * pss);
+    CCTK_INT CarpetAdaptiveRegrid_Regrid (CCTK_POINTER_TO_CONST const cctkGH_,
+                                  CCTK_POINTER const bbsss_,
+                                  CCTK_POINTER const obss_,
+                                  CCTK_POINTER const pss_,
+				  CCTK_INT force);
+  }
+
+} // namespace CarpetAdaptiveregrid
+
+#endif // !defined(CARPETADAPTIVEREGRID_HH)
diff --git a/CarpetDev/CarpetAdaptiveRegrid/src/CAR_Paramcheck.cc b/CarpetDev/CarpetAdaptiveRegrid/src/CAR_Paramcheck.cc
new file mode 100644
index 000000000..5a6970565
--- /dev/null
+++ b/CarpetDev/CarpetAdaptiveRegrid/src/CAR_Paramcheck.cc
@@ -0,0 +1,36 @@
+#include <assert.h>
+
+#include "cctk.h"
+#include "cctk_Arguments.h"
+#include "cctk_Parameters.h"
+
+#include "carpet.hh"
+#include "CAR.hh"
+
+static const char* rcsid = "$Header:$";
+
+CCTK_FILEVERSION(CAR_Paramcheck_cc)
+
+namespace CarpetAdaptiveRegrid {
+  
+  using namespace std;
+  using namespace Carpet;
+  
+  int CarpetAdaptiveRegridParamcheck (CCTK_ARGUMENTS)
+  {
+    DECLARE_CCTK_ARGUMENTS;
+    DECLARE_CCTK_PARAMETERS;
+    
+    int type;
+    const CCTK_INT * const domain_from_coordbase
+      = (const CCTK_INT *) CCTK_ParameterGet ("domain_from_coordbase", "Carpet", &type);
+    assert (domain_from_coordbase);
+    assert (type == PARAMETER_BOOLEAN);
+    if (! *domain_from_coordbase) {
+      CCTK_PARAMWARN ("CarpetAdaptiveRegrid requires that Carpet::domain_from_coordbase=yes");
+    }
+    
+    return 0;
+  }
+  
+}
diff --git a/CarpetDev/CarpetAdaptiveRegrid/src/CAR_utils.F90 b/CarpetDev/CarpetAdaptiveRegrid/src/CAR_utils.F90
new file mode 100644
index 000000000..64018c999
--- /dev/null
+++ b/CarpetDev/CarpetAdaptiveRegrid/src/CAR_utils.F90
@@ -0,0 +1,776 @@
+!!$  $Header:$
+      
+#include "cctk.h"
+
+!!$ Given the mask set the 1D sums
+
+subroutine count_points(nx, ny, nz, mask, sum_x, sum_y, sum_z)
+
+  implicit none
+
+  CCTK_INT :: nx, ny, nz
+  CCTK_INT :: mask(nx, ny, nz)
+  CCTK_INT :: sum_x(nx), sum_y(ny), sum_z(nz)
+  
+  CCTK_INT :: i, j, k
+  
+  do i = 1, nx
+    
+    sum_x(i) = 0
+    
+  end do
+  
+  do j = 1, ny
+    
+    sum_y(j) = 0
+    
+  end do
+  
+  do k = 1, nz
+    
+    sum_z(k) = 0
+    
+  end do
+  
+  do k = 1, nz
+    do j = 1, ny
+      do i = 1, nx
+        
+        sum_x(i) = sum_x(i) + mask(i, j, k)
+        sum_y(j) = sum_y(j) + mask(i, j, k)
+        sum_z(k) = sum_z(k) + mask(i, j, k)
+        
+      end do
+    end do
+  end do
+  
+end subroutine count_points
+
+!!$ Given a sum compute the signature
+
+subroutine signature(nx, lsum, sig)
+
+  implicit none
+  
+  CCTK_INT :: nx
+  CCTK_INT :: lsum(nx), sig(nx)
+  
+  CCTK_INT :: i
+  
+  sig(1)  = 0
+  sig(nx) = 0
+  
+  do i = 2, nx - 1
+    
+    sig(i) = lsum(i - 1) - 2 * lsum(i) + lsum(i + 1)
+    
+  end do
+  
+end subroutine signature
+
+!!$ Given a sum prune the zeros.
+!!$ That is, find the minimum and maximum non-zero indices
+
+subroutine prune(nx, lsum, ilo, ihi)
+
+  implicit none
+  
+  CCTK_INT :: nx
+  CCTK_INT :: lsum(nx)
+  CCTK_INT :: ilo, ihi
+  
+  CCTK_INT :: i
+  
+  ilo = 0
+  ihi = 0
+  
+  i = 0
+  
+11 if (ilo .eq. 0) then
+    i = i + 1
+    if (i > nx) then
+      call CCTK_WARN(0, "Error in prune; sum is all zero!")
+    end if
+    if (lsum(i) > 0) then
+      ilo = i
+    end if
+    goto 11
+  end if
+  
+  i = nx + 1
+12 if (ihi .eq. 0) then
+    i = i - 1
+!!$    write(*,*) "prune:",nx,i,lsum(i)
+    if (i < 1) then
+      call CCTK_WARN(0, "Error in prune; sum is all zero!")
+    end if
+    if (lsum(i) > 0) then
+      ihi = i
+    end if
+    goto 12
+  end if
+  
+end subroutine prune
+
+!!$ Prune box.
+!!$ didit is 0 is nothing was done, otherwise 1
+!!$ newbbox is always set.
+
+subroutine prune_box(nx, ny, nz, sum_x, sum_y, sum_z, &
+                     bbox, newbbox, didit)
+
+  implicit none
+  
+  CCTK_INT :: nx, ny, nz
+  CCTK_INT :: sum_x(nx), sum_y(ny), sum_z(nz)
+  CCTK_INT :: bbox(3,3), newbbox(3,3)
+  CCTK_INT :: didit
+  
+  CCTK_INT :: lo, hi
+  
+  didit = 0
+
+  newbbox(:,3) = bbox(:,3)
+
+!!$  x direction
+
+  call prune(nx, sum_x, lo, hi)
+  
+  newbbox(1,1) = bbox(1,1) + lo - 1
+  newbbox(1,2) = bbox(1,2) + hi - nx
+
+!!$  write(*,*) "x dirn",lo,hi,nx,bbox(1,:)
+  
+  if ( (lo > 1) .or. (hi < nx) ) then
+    didit = 1
+  end if
+  
+!!$  y direction
+
+  call prune(ny, sum_y, lo, hi)
+  
+  newbbox(2,1) = bbox(2,1) + lo - 1
+  newbbox(2,2) = bbox(2,2) + hi - ny
+
+!!$  write(*,*) "y dirn",lo,hi,ny,bbox(2,:)
+  
+  if ( (lo > 1) .or. (hi < ny) ) then
+    didit = 1
+  end if
+  
+!!$  z direction
+
+  call prune(nz, sum_z, lo, hi)
+  
+  newbbox(3,1) = bbox(3,1) + lo - 1
+  newbbox(3,2) = bbox(3,2) + hi - nz
+
+!!$  write(*,*) "z dirn",lo,hi,nz,bbox(3,:)
+  
+  if ( (lo > 1) .or. (hi < nz) ) then
+    didit = 1
+  end if
+  
+end subroutine prune_box
+
+!!$ Find holes.
+!!$ That is, find the first location within the grid 
+!!$ (taking into account the minimum width parameters)
+!!$ that contains a zero
+!!$
+!!$ ihole is set to zero if there is no hole
+
+subroutine find_hole(nx, lsum, min_width, ihole)
+  
+  implicit none
+  
+  CCTK_INT :: nx
+  CCTK_INT :: lsum(nx)
+  CCTK_INT :: min_width, ihole
+  
+  CCTK_INT :: i
+  
+  ihole = 0
+  
+  if (nx < 2 * min_width + 1) return
+  
+!!$  write(*,*) "find hole",nx,min_width
+
+  i = min_width 
+13 if ( (ihole .eq. 0) .and. (i < nx - min_width + 1) ) then
+    i = i + 1
+!!$    write(*,*) "fh",i,lsum(i)
+    if (lsum(i) .eq. 0) then
+      ihole = i
+    end if
+    goto 13
+  end if
+
+!!$  write(*,*) "Done find hole",ihole
+
+end subroutine find_hole
+
+!!$ Split a box at a hole if one exists
+!!$ didit is zero if nothing is done and 1 otherwise
+!!$ newbbox1 and 2 are always set, but newbbox2 is 
+!!$ meaningless if nothing is done.
+
+subroutine split_box_at_hole(nx, ny, nz, sum_x, sum_y, sum_z, &
+                             bbox, newbbox1, newbbox2, min_width, &
+                             didit)
+
+  implicit none
+
+  CCTK_INT :: nx, ny, nz
+  CCTK_INT :: sum_x(nx), sum_y(ny), sum_z(nz)
+  CCTK_INT :: bbox(3,3), newbbox1(3,3), newbbox2(3,3)
+  CCTK_INT :: min_width
+  CCTK_INT :: didit
+  
+  CCTK_INT :: ihole
+  
+  didit = 0
+  
+  newbbox1(1,1) = bbox(1,1)
+  newbbox1(2,1) = bbox(2,1)
+  newbbox1(3,1) = bbox(3,1)
+  newbbox1(1,2) = bbox(1,2)
+  newbbox1(2,2) = bbox(2,2)
+  newbbox1(3,2) = bbox(3,2)
+  
+  newbbox2(1,1) = bbox(1,1)
+  newbbox2(2,1) = bbox(2,1)
+  newbbox2(3,1) = bbox(3,1)
+  newbbox2(1,2) = bbox(1,2)
+  newbbox2(2,2) = bbox(2,2)
+  newbbox2(3,2) = bbox(3,2)
+  
+  newbbox1(:,3) = bbox(:,3)
+  newbbox2(:,3) = bbox(:,3)
+
+!!$  write(*,*) nx,ny,nz
+
+!!$  Try splitting in z first
+
+  call find_hole(nz, sum_z, min_width, ihole)
+
+!!$  write(*,*) "Find hole in z:", ihole
+
+  if (ihole > 0) then
+
+    didit = 1
+    
+    newbbox1(3,2) = bbox(3,1) + ihole - 2
+    newbbox2(3,1) = newbbox1(3,2) + 1
+    
+    return
+    
+  end if
+  
+!!$  Then try splitting in y next
+
+  call find_hole(ny, sum_y, min_width, ihole)
+
+!!$  write(*,*) "Find hole in y:", ihole
+
+  if (ihole > 0) then
+    
+    didit = 1
+    
+    newbbox1(2,2) = bbox(2,1) + ihole - 2
+    newbbox2(2,1) = newbbox1(2,2) + 1
+    
+    return
+    
+  end if
+
+!!$  Then try splitting in x last
+
+  call find_hole(nx, sum_x, min_width, ihole)
+
+!!$  write(*,*) "Find hole in x:", ihole
+
+  if (ihole > 0) then
+    
+    didit = 1
+    
+    newbbox1(1,2) = bbox(1,1) + ihole - 2
+    newbbox2(1,1) = newbbox1(1,2) + 1
+    
+    return
+    
+  end if
+  
+!!$  We should only reach here if we did not find any holes
+!!$  So set newbbox2 to dummy values
+
+  newbbox2(1,1) = 0
+  newbbox2(2,1) = 0
+  newbbox2(3,1) = 0
+  newbbox2(1,2) = 0
+  newbbox2(2,2) = 0
+  newbbox2(3,2) = 0
+  
+end subroutine split_box_at_hole
+
+!!$ Find the density.
+!!$ That is, find the ratio of marked points to total points.
+
+subroutine find_density(nx, ny, nz, mask, density)
+
+  implicit none
+
+  CCTK_INT :: nx, ny, nz
+  CCTK_INT :: mask(nx, ny, nz)
+  CCTK_REAL :: density
+  
+  CCTK_INT :: marked, i, j, k
+  
+  marked = 0
+  
+  do k = 1, nz
+    do j = 1, ny
+      do i = 1, nx
+        
+        marked = marked + mask(i, j, k)
+        
+      end do
+    end do
+  end do
+  
+  density = dble(marked) / dble(nx * ny * nz)
+  
+end subroutine find_density
+
+!!$ Split along a direction.
+!!$ That is, find the zero crossing in the signature
+!!$ (taking into account the minimum width parameters)
+!!$ with the largest first difference in the signature.
+!!$
+!!$ If the grid should not be split, isplit is zero
+!!$ max_jump is passed in so that we can try all directions.
+
+subroutine split(nx, signature, min_width, isplit, max_jump)
+
+  implicit none
+  
+  CCTK_INT :: nx
+  CCTK_INT :: signature(nx)
+  CCTK_INT :: min_width, isplit, max_jump
+  
+  CCTK_INT :: i
+  CCTK_INT :: jump
+  
+  isplit = 0
+  
+  if (nx < 2 * min_width + 1) return
+  
+  do i = min_width + 1, nx - min_width
+    
+!!$    write(*,*) "split:",i,signature(i),signature(i-1)
+
+    if (signature(i) * signature(i - 1) < 0) then
+      
+      jump = abs( signature(i) - signature(i - 1) )
+      
+!!$      write(*,*) "jump at",i,"is",jump
+      
+      if (jump > max_jump) then
+        
+        isplit = i
+        max_jump = jump
+        
+      end if
+      
+    end if
+    
+  end do
+  
+end subroutine split
+
+!!$ Split the box at turning points in the signature
+!!$ Arguments and general behaviour are identical to split_box_at_hole
+
+subroutine split_box_at_sig(nx, ny, nz, sig_x, sig_y, sig_z, &
+                            bbox, newbbox1, newbbox2, min_width, &
+                            didit)
+
+  implicit none
+
+  CCTK_INT :: nx, ny, nz
+  CCTK_INT :: sig_x(nx), sig_y(ny), sig_z(nz)
+  CCTK_INT :: bbox(3,3), newbbox1(3,3), newbbox2(3,3)
+  CCTK_INT :: min_width
+  CCTK_INT :: didit
+  
+  CCTK_INT :: isplit, max_jump
+  
+  didit = 0
+  max_jump = 0
+  
+  newbbox1(1,1) = bbox(1,1)
+  newbbox1(2,1) = bbox(2,1)
+  newbbox1(3,1) = bbox(3,1)
+  newbbox1(1,2) = bbox(1,2)
+  newbbox1(2,2) = bbox(2,2)
+  newbbox1(3,2) = bbox(3,2)
+  
+  newbbox2(1,1) = bbox(1,1)
+  newbbox2(2,1) = bbox(2,1)
+  newbbox2(3,1) = bbox(3,1)
+  newbbox2(1,2) = bbox(1,2)
+  newbbox2(2,2) = bbox(2,2)
+  newbbox2(3,2) = bbox(3,2)
+  
+  newbbox1(:,3) = bbox(:,3)
+  newbbox2(:,3) = bbox(:,3)
+
+!!$  Try splitting in z first
+
+  call split(nz, sig_z, min_width, isplit, max_jump)
+  
+!!$  write(*,*) "Split in z at sig:",isplit
+  
+  if (isplit > 0) then
+    
+    didit = 1
+    
+    newbbox1(3,2) = bbox(3,1) + isplit - 2
+    newbbox2(3,1) = newbbox1(3,2) + 1
+    
+  end if
+  
+!!$  Then try splitting in y next
+
+  call split(ny, sig_y, min_width, isplit, max_jump)
+  
+!!$  write(*,*) "Split in y at sig:",isplit
+  
+  if (isplit > 0) then
+    
+    didit = 1
+    
+    newbbox1(3,2) = bbox(3,2)
+    newbbox2(3,1) = bbox(3,1)
+    newbbox1(2,2) = bbox(2,1) + isplit - 2
+    newbbox2(2,1) = newbbox1(2,2) + 1
+    
+  end if
+  
+!!$  Then try splitting in x last
+
+  call split(nx, sig_x, min_width, isplit, max_jump)
+  
+!!$  write(*,*) "Split in x at sig:",isplit
+  
+  if (isplit > 0) then
+    
+    didit = 1
+    
+    newbbox1(3,2) = bbox(3,2)
+    newbbox2(3,1) = bbox(3,1)
+    newbbox1(2,2) = bbox(2,2)
+    newbbox2(2,1) = bbox(2,1)
+    newbbox1(1,2) = bbox(1,1) + isplit - 2
+    newbbox2(1,1) = newbbox1(1,2) + 1
+    
+  end if
+  
+  if (didit > 0) return
+
+!!$  We should only reach here if we did not find any holes
+!!$  So set newbbox2 to dummy values
+
+  newbbox2(1,1) = 0
+  newbbox2(2,1) = 0
+  newbbox2(3,1) = 0
+  newbbox2(1,2) = 0
+  newbbox2(2,2) = 0
+  newbbox2(3,2) = 0
+  
+end subroutine split_box_at_sig
+
+!!$ Prune a new box
+!!$ This is to be done after the grid is split
+
+subroutine prune_new_box(nx, ny, nz, mask, bbox, newbbox)
+
+  implicit none
+
+  CCTK_INT :: nx, ny, nz
+  CCTK_INT :: mask(nx, ny, nz)
+  CCTK_INT :: bbox(3,3), newbbox(3,3)
+
+  CCTK_INT, dimension(:,:,:), allocatable :: tmp_mask
+  CCTK_INT, dimension(:), allocatable :: tmp_sum_x, tmp_sum_y, tmp_sum_z
+  CCTK_INT, dimension(:), allocatable :: tmp_sig_x, tmp_sig_y, tmp_sig_z
+  CCTK_INT :: tmp_nx, tmp_ny, tmp_nz
+  CCTK_INT :: tmp_bbox(3,3)
+  CCTK_INT :: tmp_didit
+  CCTK_INT :: ierr
+
+  tmp_nx = newbbox(1,2) - newbbox(1,1) + 1
+  tmp_ny = newbbox(2,2) - newbbox(2,1) + 1
+  tmp_nz = newbbox(3,2) - newbbox(3,1) + 1
+
+  newbbox(:,3) = bbox(:,3)
+  
+  allocate(tmp_mask(tmp_nx,tmp_ny,tmp_nz),&
+           tmp_sum_x(tmp_nx),tmp_sum_y(tmp_ny),tmp_sum_z(tmp_nz),&
+           tmp_sig_x(tmp_nx),tmp_sig_y(tmp_ny),tmp_sig_z(tmp_nz),&
+           STAT=ierr)
+  if (ierr .ne. 0) call CCTK_WARN(0, "Allocation error!")
+  
+  call copy_mask(nx, ny, nz, mask, bbox, &
+                 tmp_nx, tmp_ny, tmp_nz, tmp_mask, newbbox)
+
+  call count_points(tmp_nx, tmp_ny, tmp_nz, tmp_mask, &
+                    tmp_sum_x, tmp_sum_y, tmp_sum_z)
+
+  call prune_box(tmp_nx, tmp_ny, tmp_nz, &
+                 tmp_sum_x, tmp_sum_y, tmp_sum_z, &
+                 newbbox, tmp_bbox, tmp_didit)
+
+!!$  write(*,*) "Pruning split box. Did it:", tmp_didit
+
+  if (tmp_didit > 0) then
+    newbbox = tmp_bbox
+  end if
+  
+  deallocate(tmp_mask, tmp_sum_x, tmp_sum_y, tmp_sum_z, &
+             tmp_sig_x, tmp_sig_y, tmp_sig_z, &
+             STAT=ierr)
+  if (ierr .ne. 0) call CCTK_WARN(0, "Deallocation error!")
+
+end subroutine prune_new_box
+
+!!$ Check a box
+!!$ This function runs through the entire clustering algorithm for a 
+!!$ single box.
+!!$
+!!$ The return value didit has 3 states:
+!!$
+!!$     0: nothing happened so this box is now fine
+!!$     1: the box requires pruning
+!!$     2: the box requires splitting
+!!$
+!!$ As the new box would require memory, pass it back up to the C++ code.
+
+subroutine check_box(nx, ny, nz, mask, sum_x, sum_y, sum_z, &
+                                       sig_x, sig_y, sig_z, &
+                     bbox, newbbox1, newbbox2, &
+                     min_width, min_density, &
+                     didit)
+
+  implicit none
+
+  CCTK_INT :: nx, ny, nz
+  CCTK_INT :: mask(nx, ny, nz)
+  CCTK_INT :: sum_x(nx), sum_y(ny), sum_z(nz)
+  CCTK_INT :: sig_x(nx), sig_y(ny), sig_z(nz)
+  CCTK_INT :: bbox(3,3), newbbox1(3,3), newbbox2(3,3)
+  CCTK_INT :: min_width
+  CCTK_REAL :: min_density
+  CCTK_INT :: didit
+  
+  CCTK_INT :: i
+
+  CCTK_REAL :: density
+
+!!$  write(*,*) nx, ny, nz
+
+!!$  First set up the sums
+
+  call count_points(nx, ny, nz, mask, sum_x, sum_y, sum_z)
+
+!!$  Then prune the box
+
+  call prune_box(nx, ny, nz, sum_x, sum_y, sum_z, bbox, newbbox1, &
+                 didit)
+
+!!$  If it needs pruning go back to C++.
+
+!!$  write(*,*) "Pruned. Did it:", didit
+
+  if (didit > 0) return
+
+!!$  Otherwise the box doesn't need pruning. Try finding a hole.
+
+  call split_box_at_hole(nx, ny, nz, sum_x, sum_y, sum_z, &
+                         bbox, newbbox1, newbbox2, min_width, didit)
+
+!!$  If it needs splitting then go back to C++
+
+!!$  write(*,*) "Split box at hole. Did it:", didit
+  
+  if (didit > 0) then
+
+!!$    Prune the split boxes
+
+!!$    Box 1
+
+    call prune_new_box(nx, ny, nz, mask, bbox, newbbox1)
+
+!!$    Box 2
+
+    call prune_new_box(nx, ny, nz, mask, bbox, newbbox2)
+
+    didit = 2
+    return
+  end if
+
+!!$  Otherwise there are no holes. Check the density.
+
+!!$  write(*,*) "Finding density"
+      
+  call find_density(nx, ny, nz, mask, density)
+  
+!!$  write(*,*) "Found density:", density
+  
+!!$  If the density is sufficient then this box is done.
+
+  if (density > min_density) then
+    didit = 0
+    return
+  end if
+
+!!$  Otherwise we try and split the box.
+
+!!$  Set up the signature arrays.
+
+!!$  write(*,*) "Setting up signatures"
+
+  call signature(nx, sum_x, sig_x)
+  call signature(ny, sum_y, sig_y)
+  call signature(nz, sum_z, sig_z)
+
+!!$  write(*,*) "Sums:"
+!!$
+!!$  do i = ny, 1, -1
+!!$    write(*,'(i3)') sum_y(i)
+!!$  end do
+!!$  
+!!$  do i = 1, nx
+!!$    write(*,'("   ",i3," ")',advance="no") sum_x(i)
+!!$  end do
+!!$  write(*,*)
+!!$
+!!$  write(*,*) "Sigs:"
+!!$
+!!$  do i = ny, 1, -1
+!!$    write(*,'(i3)') sig_y(i)
+!!$  end do
+!!$  
+!!$  do i = 1, nx
+!!$    write(*,'("   ",i3," ")',advance="no") sig_x(i)
+!!$  end do
+!!$  write(*,*)
+
+!!$  write(*,*) "Sum x:"
+!!$  write(*,*) sum_x
+!!$  write(*,*) "Sig x:"
+!!$  write(*,*) sig_x
+!!$  write(*,*) "Sum y:"
+!!$  write(*,*) sum_y
+!!$  write(*,*) "Sig y:"
+!!$  write(*,*) sig_y
+!!$  write(*,*) "Sum z:"
+!!$  write(*,*) sum_z
+!!$  write(*,*) "Sig z:"
+!!$  write(*,*) sig_z
+
+!!$  Try splitting by the signature
+
+!!$  write(*,*) "Splitting at signatures"
+
+  call split_box_at_sig(nx, ny, nz, sig_x, sig_y, sig_z, &
+                        bbox, newbbox1, newbbox2, min_width, &
+                        didit)
+
+!!$  This is the last trick up our sleeve. So we just check if
+!!$  it succeeded so that we can correct the value of didit.
+
+!!$  write(*,*) "Split box at signature. Did it:", didit
+
+  if (didit > 0) then
+
+!!$    Prune the split boxes
+
+!!$    Box 1
+
+    call prune_new_box(nx, ny, nz, mask, bbox, newbbox1)
+
+!!$    Box 2
+
+    call prune_new_box(nx, ny, nz, mask, bbox, newbbox2)
+
+    didit = 2
+  end if
+
+!!$  Then we are done
+
+end subroutine check_box
+
+!!$ Copy the mask.
+!!$ Given two masks, one strictly contained within the other, copy
+!!$ the intersection to the smaller mask.
+!!$ 
+!!$ All "s" quantities are the "source" mask - the larger one
+!!$ All "d" quantities are the "destination" mask
+!!$ The location of the masks in "grid space" is denoted by 
+!!$ the bboxes.
+!!$ ?bbox(:,1) is the lower boundary
+!!$ ?bbox(:,2) is the upper boundary
+!!$ ?bbox(:,3) is the stride and is irrelevant
+
+subroutine copy_mask(snx, sny, snz, smask, sbbox, &
+                     dnx, dny, dnz, dmask, dbbox)
+
+  implicit none
+
+  CCTK_INT :: snx, sny, snz
+  CCTK_INT :: smask(snx, sny, snz)
+  CCTK_INT :: sbbox(3, 2)
+  CCTK_INT :: dnx, dny, dnz
+  CCTK_INT :: dmask(dnx, dny, dnz)
+  CCTK_INT :: dbbox(3, 2)
+  
+  CCTK_INT :: i, j, k, si, sj, sk, di, dj, dk
+  
+  if ( (dbbox(1,1) < sbbox(1,1)) .or. &
+       (dbbox(2,1) < sbbox(2,1)) .or. &
+       (dbbox(3,1) < sbbox(3,1)) .or. &
+       (dbbox(1,2) > sbbox(1,2)) .or. &
+       (dbbox(2,2) > sbbox(2,2)) .or. &
+       (dbbox(3,2) > sbbox(3,2)) ) then
+    call CCTK_WARN(0, &
+        "The destination mask is not contained in the source mask!")
+  end if
+
+  do k = dbbox(3,1), dbbox(3,2)
+    do j = dbbox(2,1), dbbox(2,2)
+      do i = dbbox(1,1), dbbox(1,2)
+        
+        si = 1 + i - sbbox(1,1)
+        sj = 1 + j - sbbox(2,1)
+        sk = 1 + k - sbbox(3,1)
+        
+        di = 1 + i - dbbox(1,1)
+        dj = 1 + j - dbbox(2,1)
+        dk = 1 + k - dbbox(3,1)
+        
+        dmask(di, dj, dk) = smask(si, sj, sk)
+
+!!$        write(*,*) "copying mask: loc",i,j,k
+!!$        write(*,*) "copying mask:d", di, dj, dk
+!!$        write(*,'("(",i2,":",i2,"):(",i2,":",i2,"):(",i2,":",i2,")")') &
+!!$             dbbox
+!!$        write(*,*) "copying mask:s", si, sj, sk
+!!$        write(*,'("(",i2,":",i2,"):(",i2,":",i2,"):(",i2,":",i2,")")') &
+!!$             sbbox
+        
+      end do
+    end do
+  end do
+
+end subroutine copy_mask
diff --git a/CarpetDev/CarpetAdaptiveRegrid/src/make.code.defn b/CarpetDev/CarpetAdaptiveRegrid/src/make.code.defn
new file mode 100644
index 000000000..049bdfb27
--- /dev/null
+++ b/CarpetDev/CarpetAdaptiveRegrid/src/make.code.defn
@@ -0,0 +1,10 @@
+# Main make.code.defn file for thorn CarpetAdaptiveRegrid
+# $Header$
+
+# Source files in this directory
+SRCS = 	CAR_Paramcheck.cc \\
+	CAR.cc \\
+	CAR_utils.F90
+
+# Subdirectories containing source files
+SUBDIRS = 
diff --git a/CarpetDev/CarpetAdaptiveRegrid/src/test/CAR_test.F90 b/CarpetDev/CarpetAdaptiveRegrid/src/test/CAR_test.F90
new file mode 100644
index 000000000..00640fa30
--- /dev/null
+++ b/CarpetDev/CarpetAdaptiveRegrid/src/test/CAR_test.F90
@@ -0,0 +1,134 @@
+#include "cctk.h"
+
+program CAR_test
+
+  CCTK_INT :: nx, ny, nz
+  CCTK_INT, dimension(:,:,:), allocatable :: mask
+  CCTK_INT, dimension(:), allocatable :: sum_x, sum_y, sum_z, &
+       sig_x, sig_y, sig_z
+  CCTK_INT, dimension(3,2) :: bbox, newbbox1, newbbox2
+  CCTK_INT :: min_width
+  CCTK_REAL :: min_density
+  CCTK_INT :: didit
+
+  CCTK_INT :: ierr, i, j
+
+  nx = 10
+  ny = 12
+  nz = 8
+!!$  nx = 9
+!!$  ny = 8
+!!$  nz = 1
+
+  min_width = 3
+!!$  min_density = 0.8
+  min_density = 0.9
+
+  didit = -1
+
+  bbox(:,1) = 1
+  bbox(1,2) = nx
+  bbox(2,2) = ny
+  bbox(3,2) = nz
+
+  newbbox1 = 0
+  newbbox2 = 0
+
+  allocate(mask(nx,ny,nz),&
+       sum_x(nx),sig_x(nx),&
+       sum_y(ny),sig_y(ny),&
+       sum_z(nz),sig_z(nz),STAT=ierr)
+  if (ierr .ne. 0) call CCTK_WARN(0, "Allocation error")
+
+  sum_x = 0
+  sum_y = 0
+  sum_z = 0
+  sig_x = 0
+  sig_y = 0
+  sig_z = 0
+
+  mask = 0
+
+!!$  Pruning check:
+!!$  mask(3:7,4:8,1:6) = 1
+
+!!$  Split at hole check (x)
+!!$  mask(:4,:,:) = 1
+!!$  mask(7:,:,:) = 1
+
+!!$  Split at hole check (y)
+!!$  mask(:,:5,:) = 1
+!!$  mask(:,8:,:) = 1
+
+!!$  Split at hole check (z)
+!!$  mask(:,:,:3) = 1
+!!$  mask(:,:,6:) = 1
+
+!!$  Density check
+!!$  mask = 1
+!!$  mask(4:5,7:9,3:5) = 0
+
+!!$  Signature check
+!!$  This should have a density of 2/3
+  mask(:4,:8,:) = 1
+  mask(5:,5:,:) = 1
+
+!!$  Signature check
+!!$  This should have a density of 3/4
+!!$  mask(:5,:,:) = 1
+!!$  mask(6:,5:,:) = 1
+
+!!$  Signature check
+!!$  This should have a density of just under 0.8
+!!$  mask(:5,:,:) = 1
+!!$  mask(6:,6:,:) = 1
+
+!!$  Test from
+!!$ http://www-lmc.imag.fr/IDOPT/AGRIF/documentation/doc_AGRIF3DUser/node9.html
+!!$  Requires
+!!$  nx = 9
+!!$  ny = 8
+!!$  min_density > 0.681
+!!$  This one is a bit pathological; I don't think the website description
+!!$  is correct
+!!$  mask(1:2,:,:) = 1
+!!$  mask(3,1:2,:) = 1; mask(3,4:5,:) = 1; mask(3,7:8,:) = 1
+!  mask(3,1:5,:) = 1; mask(3,7:8,:) = 1
+!  mask(3,:,:)   = 1
+!!$  mask(4,:,:)   = 1
+!!$  mask(5,3:6,:) = 1
+!  mask(6,3,:)   = 1; mask(6,5:6,:) = 1
+!!$  mask(6,3:6,:) = 1
+!!$  mask(7:9,3:6,:) = 1
+
+  do j = ny, 1, -1
+    do i = 1, nx
+      write(*,'(i1," ")',advance="no") mask(i,j,1)
+    end do
+    write(*,*) " "
+  end do
+
+  call check_box(nx, ny, nz, mask, sum_x, sum_y, sum_z, &
+       sig_x, sig_y, sig_z, bbox, newbbox1, newbbox2, &
+       min_width, min_density, didit)
+  
+  write(*,*) "Done. Did it:", didit
+
+  select case (didit)
+    case(0)
+      write(*,*) "Unchanged"
+    case(1)
+      write(*,*) "Single new bbox"
+      write(*,'(a10,"(",i2,":",i2,"):(",i2,":",i2,"):(",i2,":",i2,")")') &
+           "new bbox: ", newbbox1(1,:),newbbox1(2,:),newbbox1(3,:)
+    case(2)
+      write(*,*) "Two new bboxes"
+      write(*,'(a11,"(",i2,":",i2,"):(",i2,":",i2,"):(",i2,":",i2,")")') &
+           "new bbox1: ", newbbox1(1,:),newbbox1(2,:),newbbox1(3,:)
+      write(*,'(a11,"(",i2,":",i2,"):(",i2,":",i2,"):(",i2,":",i2,")")') &
+           "new bbox2: ", newbbox2(1,:),newbbox2(2,:),newbbox2(3,:)
+    case default
+      write(*,*) "Error in return; didit is",didit
+  end select
+
+end program CAR_test
diff --git a/CarpetDev/CarpetAdaptiveRegrid/src/test/CAR_utils.F90 b/CarpetDev/CarpetAdaptiveRegrid/src/test/CAR_utils.F90
new file mode 100644
index 000000000..64018c999
--- /dev/null
+++ b/CarpetDev/CarpetAdaptiveRegrid/src/test/CAR_utils.F90
@@ -0,0 +1,776 @@
+!!$  $Header:$
+      
+#include "cctk.h"
+
+!!$ Given the mask set the 1D sums
+
+subroutine count_points(nx, ny, nz, mask, sum_x, sum_y, sum_z)
+
+  implicit none
+
+  CCTK_INT :: nx, ny, nz
+  CCTK_INT :: mask(nx, ny, nz)
+  CCTK_INT :: sum_x(nx), sum_y(ny), sum_z(nz)
+  
+  CCTK_INT :: i, j, k
+  
+  do i = 1, nx
+    
+    sum_x(i) = 0
+    
+  end do
+  
+  do j = 1, ny
+    
+    sum_y(j) = 0
+    
+  end do
+  
+  do k = 1, nz
+    
+    sum_z(k) = 0
+    
+  end do
+  
+  do k = 1, nz
+    do j = 1, ny
+      do i = 1, nx
+        
+        sum_x(i) = sum_x(i) + mask(i, j, k)
+        sum_y(j) = sum_y(j) + mask(i, j, k)
+        sum_z(k) = sum_z(k) + mask(i, j, k)
+        
+      end do
+    end do
+  end do
+  
+end subroutine count_points
+
+!!$ Given a sum compute the signature
+
+subroutine signature(nx, lsum, sig)
+
+  implicit none
+  
+  CCTK_INT :: nx
+  CCTK_INT :: lsum(nx), sig(nx)
+  
+  CCTK_INT :: i
+  
+  sig(1)  = 0
+  sig(nx) = 0
+  
+  do i = 2, nx - 1
+    
+    sig(i) = lsum(i - 1) - 2 * lsum(i) + lsum(i + 1)
+    
+  end do
+  
+end subroutine signature
+
+!!$ Given a sum prune the zeros.
+!!$ That is, find the minimum and maximum non-zero indices
+
+subroutine prune(nx, lsum, ilo, ihi)
+
+  implicit none
+  
+  CCTK_INT :: nx
+  CCTK_INT :: lsum(nx)
+  CCTK_INT :: ilo, ihi
+  
+  CCTK_INT :: i
+  
+  ilo = 0
+  ihi = 0
+  
+  i = 0
+  
+11 if (ilo .eq. 0) then
+    i = i + 1
+    if (i > nx) then
+      call CCTK_WARN(0, "Error in prune; sum is all zero!")
+    end if
+    if (lsum(i) > 0) then
+      ilo = i
+    end if
+    goto 11
+  end if
+  
+  i = nx + 1
+12 if (ihi .eq. 0) then
+    i = i - 1
+!!$    write(*,*) "prune:",nx,i,lsum(i)
+    if (i < 1) then
+      call CCTK_WARN(0, "Error in prune; sum is all zero!")
+    end if
+    if (lsum(i) > 0) then
+      ihi = i
+    end if
+    goto 12
+  end if
+  
+end subroutine prune
+
+!!$ Prune box.
+!!$ didit is 0 is nothing was done, otherwise 1
+!!$ newbbox is always set.
+
+subroutine prune_box(nx, ny, nz, sum_x, sum_y, sum_z, &
+                     bbox, newbbox, didit)
+
+  implicit none
+  
+  CCTK_INT :: nx, ny, nz
+  CCTK_INT :: sum_x(nx), sum_y(ny), sum_z(nz)
+  CCTK_INT :: bbox(3,3), newbbox(3,3)
+  CCTK_INT :: didit
+  
+  CCTK_INT :: lo, hi
+  
+  didit = 0
+
+  newbbox(:,3) = bbox(:,3)
+
+!!$  x direction
+
+  call prune(nx, sum_x, lo, hi)
+  
+  newbbox(1,1) = bbox(1,1) + lo - 1
+  newbbox(1,2) = bbox(1,2) + hi - nx
+
+!!$  write(*,*) "x dirn",lo,hi,nx,bbox(1,:)
+  
+  if ( (lo > 1) .or. (hi < nx) ) then
+    didit = 1
+  end if
+  
+!!$  y direction
+
+  call prune(ny, sum_y, lo, hi)
+  
+  newbbox(2,1) = bbox(2,1) + lo - 1
+  newbbox(2,2) = bbox(2,2) + hi - ny
+
+!!$  write(*,*) "y dirn",lo,hi,ny,bbox(2,:)
+  
+  if ( (lo > 1) .or. (hi < ny) ) then
+    didit = 1
+  end if
+  
+!!$  z direction
+
+  call prune(nz, sum_z, lo, hi)
+  
+  newbbox(3,1) = bbox(3,1) + lo - 1
+  newbbox(3,2) = bbox(3,2) + hi - nz
+
+!!$  write(*,*) "z dirn",lo,hi,nz,bbox(3,:)
+  
+  if ( (lo > 1) .or. (hi < nz) ) then
+    didit = 1
+  end if
+  
+end subroutine prune_box
+
+!!$ Find holes.
+!!$ That is, find the first location within the grid 
+!!$ (taking into account the minimum width parameters)
+!!$ that contains a zero
+!!$
+!!$ ihole is set to zero if there is no hole
+
+subroutine find_hole(nx, lsum, min_width, ihole)
+  
+  implicit none
+  
+  CCTK_INT :: nx
+  CCTK_INT :: lsum(nx)
+  CCTK_INT :: min_width, ihole
+  
+  CCTK_INT :: i
+  
+  ihole = 0
+  
+  if (nx < 2 * min_width + 1) return
+  
+!!$  write(*,*) "find hole",nx,min_width
+
+  i = min_width 
+13 if ( (ihole .eq. 0) .and. (i < nx - min_width + 1) ) then
+    i = i + 1
+!!$    write(*,*) "fh",i,lsum(i)
+    if (lsum(i) .eq. 0) then
+      ihole = i
+    end if
+    goto 13
+  end if
+
+!!$  write(*,*) "Done find hole",ihole
+
+end subroutine find_hole
+
+!!$ Split a box at a hole if one exists
+!!$ didit is zero if nothing is done and 1 otherwise
+!!$ newbbox1 and 2 are always set, but newbbox2 is 
+!!$ meaningless if nothing is done.
+
+subroutine split_box_at_hole(nx, ny, nz, sum_x, sum_y, sum_z, &
+                             bbox, newbbox1, newbbox2, min_width, &
+                             didit)
+
+  implicit none
+
+  CCTK_INT :: nx, ny, nz
+  CCTK_INT :: sum_x(nx), sum_y(ny), sum_z(nz)
+  CCTK_INT :: bbox(3,3), newbbox1(3,3), newbbox2(3,3)
+  CCTK_INT :: min_width
+  CCTK_INT :: didit
+  
+  CCTK_INT :: ihole
+  
+  didit = 0
+  
+  newbbox1(1,1) = bbox(1,1)
+  newbbox1(2,1) = bbox(2,1)
+  newbbox1(3,1) = bbox(3,1)
+  newbbox1(1,2) = bbox(1,2)
+  newbbox1(2,2) = bbox(2,2)
+  newbbox1(3,2) = bbox(3,2)
+  
+  newbbox2(1,1) = bbox(1,1)
+  newbbox2(2,1) = bbox(2,1)
+  newbbox2(3,1) = bbox(3,1)
+  newbbox2(1,2) = bbox(1,2)
+  newbbox2(2,2) = bbox(2,2)
+  newbbox2(3,2) = bbox(3,2)
+  
+  newbbox1(:,3) = bbox(:,3)
+  newbbox2(:,3) = bbox(:,3)
+
+!!$  write(*,*) nx,ny,nz
+
+!!$  Try splitting in z first
+
+  call find_hole(nz, sum_z, min_width, ihole)
+
+!!$  write(*,*) "Find hole in z:", ihole
+
+  if (ihole > 0) then
+
+    didit = 1
+    
+    newbbox1(3,2) = bbox(3,1) + ihole - 2
+    newbbox2(3,1) = newbbox1(3,2) + 1
+    
+    return
+    
+  end if
+  
+!!$  Then try splitting in y next
+
+  call find_hole(ny, sum_y, min_width, ihole)
+
+!!$  write(*,*) "Find hole in y:", ihole
+
+  if (ihole > 0) then
+    
+    didit = 1
+    
+    newbbox1(2,2) = bbox(2,1) + ihole - 2
+    newbbox2(2,1) = newbbox1(2,2) + 1
+    
+    return
+    
+  end if
+
+!!$  Then try splitting in x last
+
+  call find_hole(nx, sum_x, min_width, ihole)
+
+!!$  write(*,*) "Find hole in x:", ihole
+
+  if (ihole > 0) then
+    
+    didit = 1
+    
+    newbbox1(1,2) = bbox(1,1) + ihole - 2
+    newbbox2(1,1) = newbbox1(1,2) + 1
+    
+    return
+    
+  end if
+  
+!!$  We should only reach here if we did not find any holes
+!!$  So set newbbox2 to dummy values
+
+  newbbox2(1,1) = 0
+  newbbox2(2,1) = 0
+  newbbox2(3,1) = 0
+  newbbox2(1,2) = 0
+  newbbox2(2,2) = 0
+  newbbox2(3,2) = 0
+  
+end subroutine split_box_at_hole
+
+!!$ Find the density.
+!!$ That is, find the ratio of marked points to total points.
+
+subroutine find_density(nx, ny, nz, mask, density)
+
+  implicit none
+
+  CCTK_INT :: nx, ny, nz
+  CCTK_INT :: mask(nx, ny, nz)
+  CCTK_REAL :: density
+  
+  CCTK_INT :: marked, i, j, k
+  
+  marked = 0
+  
+  do k = 1, nz
+    do j = 1, ny
+      do i = 1, nx
+        
+        marked = marked + mask(i, j, k)
+        
+      end do
+    end do
+  end do
+  
+  density = dble(marked) / dble(nx * ny * nz)
+  
+end subroutine find_density
+
+!!$ Split along a direction.
+!!$ That is, find the zero crossing in the signature
+!!$ (taking into account the minimum width parameters)
+!!$ with the largest first difference in the signature.
+!!$
+!!$ If the grid should not be split, isplit is zero
+!!$ max_jump is passed in so that we can try all directions.
+
+subroutine split(nx, signature, min_width, isplit, max_jump)
+
+  implicit none
+  
+  CCTK_INT :: nx
+  CCTK_INT :: signature(nx)
+  CCTK_INT :: min_width, isplit, max_jump
+  
+  CCTK_INT :: i
+  CCTK_INT :: jump
+  
+  isplit = 0
+  
+  if (nx < 2 * min_width + 1) return
+  
+  do i = min_width + 1, nx - min_width
+    
+!!$    write(*,*) "split:",i,signature(i),signature(i-1)
+
+    if (signature(i) * signature(i - 1) < 0) then
+      
+      jump = abs( signature(i) - signature(i - 1) )
+      
+!!$      write(*,*) "jump at",i,"is",jump
+      
+      if (jump > max_jump) then
+        
+        isplit = i
+        max_jump = jump
+        
+      end if
+      
+    end if
+    
+  end do
+  
+end subroutine split
+
+!!$ Split the box at turning points in the signature
+!!$ Arguments and general behaviour are identical to split_box_at_hole
+
+subroutine split_box_at_sig(nx, ny, nz, sig_x, sig_y, sig_z, &
+                            bbox, newbbox1, newbbox2, min_width, &
+                            didit)
+
+  implicit none
+
+  CCTK_INT :: nx, ny, nz
+  CCTK_INT :: sig_x(nx), sig_y(ny), sig_z(nz)
+  CCTK_INT :: bbox(3,3), newbbox1(3,3), newbbox2(3,3)
+  CCTK_INT :: min_width
+  CCTK_INT :: didit
+  
+  CCTK_INT :: isplit, max_jump
+  
+  didit = 0
+  max_jump = 0
+  
+  newbbox1(1,1) = bbox(1,1)
+  newbbox1(2,1) = bbox(2,1)
+  newbbox1(3,1) = bbox(3,1)
+  newbbox1(1,2) = bbox(1,2)
+  newbbox1(2,2) = bbox(2,2)
+  newbbox1(3,2) = bbox(3,2)
+  
+  newbbox2(1,1) = bbox(1,1)
+  newbbox2(2,1) = bbox(2,1)
+  newbbox2(3,1) = bbox(3,1)
+  newbbox2(1,2) = bbox(1,2)
+  newbbox2(2,2) = bbox(2,2)
+  newbbox2(3,2) = bbox(3,2)
+  
+  newbbox1(:,3) = bbox(:,3)
+  newbbox2(:,3) = bbox(:,3)
+
+!!$  Try splitting in z first
+
+  call split(nz, sig_z, min_width, isplit, max_jump)
+  
+!!$  write(*,*) "Split in z at sig:",isplit
+  
+  if (isplit > 0) then
+    
+    didit = 1
+    
+    newbbox1(3,2) = bbox(3,1) + isplit - 2
+    newbbox2(3,1) = newbbox1(3,2) + 1
+    
+  end if
+  
+!!$  Then try splitting in y next
+
+  call split(ny, sig_y, min_width, isplit, max_jump)
+  
+!!$  write(*,*) "Split in y at sig:",isplit
+  
+  if (isplit > 0) then
+    
+    didit = 1
+    
+    newbbox1(3,2) = bbox(3,2)
+    newbbox2(3,1) = bbox(3,1)
+    newbbox1(2,2) = bbox(2,1) + isplit - 2
+    newbbox2(2,1) = newbbox1(2,2) + 1
+    
+  end if
+  
+!!$  Then try splitting in x last
+
+  call split(nx, sig_x, min_width, isplit, max_jump)
+  
+!!$  write(*,*) "Split in x at sig:",isplit
+  
+  if (isplit > 0) then
+    
+    didit = 1
+    
+    newbbox1(3,2) = bbox(3,2)
+    newbbox2(3,1) = bbox(3,1)
+    newbbox1(2,2) = bbox(2,2)
+    newbbox2(2,1) = bbox(2,1)
+    newbbox1(1,2) = bbox(1,1) + isplit - 2
+    newbbox2(1,1) = newbbox1(1,2) + 1
+    
+  end if
+  
+  if (didit > 0) return
+
+!!$  We should only reach here if we did not find any holes
+!!$  So set newbbox2 to dummy values
+
+  newbbox2(1,1) = 0
+  newbbox2(2,1) = 0
+  newbbox2(3,1) = 0
+  newbbox2(1,2) = 0
+  newbbox2(2,2) = 0
+  newbbox2(3,2) = 0
+  
+end subroutine split_box_at_sig
+
+!!$ Prune a new box
+!!$ This is to be done after the grid is split
+
+subroutine prune_new_box(nx, ny, nz, mask, bbox, newbbox)
+
+  implicit none
+
+  CCTK_INT :: nx, ny, nz
+  CCTK_INT :: mask(nx, ny, nz)
+  CCTK_INT :: bbox(3,3), newbbox(3,3)
+
+  CCTK_INT, dimension(:,:,:), allocatable :: tmp_mask
+  CCTK_INT, dimension(:), allocatable :: tmp_sum_x, tmp_sum_y, tmp_sum_z
+  CCTK_INT, dimension(:), allocatable :: tmp_sig_x, tmp_sig_y, tmp_sig_z
+  CCTK_INT :: tmp_nx, tmp_ny, tmp_nz
+  CCTK_INT :: tmp_bbox(3,3)
+  CCTK_INT :: tmp_didit
+  CCTK_INT :: ierr
+
+  tmp_nx = newbbox(1,2) - newbbox(1,1) + 1
+  tmp_ny = newbbox(2,2) - newbbox(2,1) + 1
+  tmp_nz = newbbox(3,2) - newbbox(3,1) + 1
+
+  newbbox(:,3) = bbox(:,3)
+  
+  allocate(tmp_mask(tmp_nx,tmp_ny,tmp_nz),&
+           tmp_sum_x(tmp_nx),tmp_sum_y(tmp_ny),tmp_sum_z(tmp_nz),&
+           tmp_sig_x(tmp_nx),tmp_sig_y(tmp_ny),tmp_sig_z(tmp_nz),&
+           STAT=ierr)
+  if (ierr .ne. 0) call CCTK_WARN(0, "Allocation error!")
+  
+  call copy_mask(nx, ny, nz, mask, bbox, &
+                 tmp_nx, tmp_ny, tmp_nz, tmp_mask, newbbox)
+
+  call count_points(tmp_nx, tmp_ny, tmp_nz, tmp_mask, &
+                    tmp_sum_x, tmp_sum_y, tmp_sum_z)
+
+  call prune_box(tmp_nx, tmp_ny, tmp_nz, &
+                 tmp_sum_x, tmp_sum_y, tmp_sum_z, &
+                 newbbox, tmp_bbox, tmp_didit)
+
+!!$  write(*,*) "Pruning split box. Did it:", tmp_didit
+
+  if (tmp_didit > 0) then
+    newbbox = tmp_bbox
+  end if
+  
+  deallocate(tmp_mask, tmp_sum_x, tmp_sum_y, tmp_sum_z, &
+             tmp_sig_x, tmp_sig_y, tmp_sig_z, &
+             STAT=ierr)
+  if (ierr .ne. 0) call CCTK_WARN(0, "Deallocation error!")
+
+end subroutine prune_new_box
+
+!!$ Check a box
+!!$ This function runs through the entire clustering algorithm for a 
+!!$ single box.
+!!$
+!!$ The return value didit has 3 states:
+!!$
+!!$     0: nothing happened so this box is now fine
+!!$     1: the box requires pruning
+!!$     2: the box requires splitting
+!!$
+!!$ As the new box would require memory, pass it back up to the C++ code.
+
+subroutine check_box(nx, ny, nz, mask, sum_x, sum_y, sum_z, &
+                                       sig_x, sig_y, sig_z, &
+                     bbox, newbbox1, newbbox2, &
+                     min_width, min_density, &
+                     didit)
+
+  implicit none
+
+  CCTK_INT :: nx, ny, nz
+  CCTK_INT :: mask(nx, ny, nz)
+  CCTK_INT :: sum_x(nx), sum_y(ny), sum_z(nz)
+  CCTK_INT :: sig_x(nx), sig_y(ny), sig_z(nz)
+  CCTK_INT :: bbox(3,3), newbbox1(3,3), newbbox2(3,3)
+  CCTK_INT :: min_width
+  CCTK_REAL :: min_density
+  CCTK_INT :: didit
+  
+  CCTK_INT :: i
+
+  CCTK_REAL :: density
+
+!!$  write(*,*) nx, ny, nz
+
+!!$  First set up the sums
+
+  call count_points(nx, ny, nz, mask, sum_x, sum_y, sum_z)
+
+!!$  Then prune the box
+
+  call prune_box(nx, ny, nz, sum_x, sum_y, sum_z, bbox, newbbox1, &
+                 didit)
+
+!!$  If it needs pruning go back to C++.
+
+!!$  write(*,*) "Pruned. Did it:", didit
+
+  if (didit > 0) return
+
+!!$  Otherwise the box doesn't need pruning. Try finding a hole.
+
+  call split_box_at_hole(nx, ny, nz, sum_x, sum_y, sum_z, &
+                         bbox, newbbox1, newbbox2, min_width, didit)
+
+!!$  If it needs splitting then go back to C++
+
+!!$  write(*,*) "Split box at hole. Did it:", didit
+  
+  if (didit > 0) then
+
+!!$    Prune the split boxes
+
+!!$    Box 1
+
+    call prune_new_box(nx, ny, nz, mask, bbox, newbbox1)
+
+!!$    Box 2
+
+    call prune_new_box(nx, ny, nz, mask, bbox, newbbox2)
+
+    didit = 2
+    return
+  end if
+
+!!$  Otherwise there are no holes. Check the density.
+
+!!$  write(*,*) "Finding density"
+      
+  call find_density(nx, ny, nz, mask, density)
+  
+!!$  write(*,*) "Found density:", density
+  
+!!$  If the density is sufficient then this box is done.
+
+  if (density > min_density) then
+    didit = 0
+    return
+  end if
+
+!!$  Otherwise we try and split the box.
+
+!!$  Set up the signature arrays.
+
+!!$  write(*,*) "Setting up signatures"
+
+  call signature(nx, sum_x, sig_x)
+  call signature(ny, sum_y, sig_y)
+  call signature(nz, sum_z, sig_z)
+
+!!$  write(*,*) "Sums:"
+!!$
+!!$  do i = ny, 1, -1
+!!$    write(*,'(i3)') sum_y(i)
+!!$  end do
+!!$  
+!!$  do i = 1, nx
+!!$    write(*,'("   ",i3," ")',advance="no") sum_x(i)
+!!$  end do
+!!$  write(*,*)
+!!$
+!!$  write(*,*) "Sigs:"
+!!$
+!!$  do i = ny, 1, -1
+!!$    write(*,'(i3)') sig_y(i)
+!!$  end do
+!!$  
+!!$  do i = 1, nx
+!!$    write(*,'("   ",i3," ")',advance="no") sig_x(i)
+!!$  end do
+!!$  write(*,*)
+
+!!$  write(*,*) "Sum x:"
+!!$  write(*,*) sum_x
+!!$  write(*,*) "Sig x:"
+!!$  write(*,*) sig_x
+!!$  write(*,*) "Sum y:"
+!!$  write(*,*) sum_y
+!!$  write(*,*) "Sig y:"
+!!$  write(*,*) sig_y
+!!$  write(*,*) "Sum z:"
+!!$  write(*,*) sum_z
+!!$  write(*,*) "Sig z:"
+!!$  write(*,*) sig_z
+
+!!$  Try splitting by the signature
+
+!!$  write(*,*) "Splitting at signatures"
+
+  call split_box_at_sig(nx, ny, nz, sig_x, sig_y, sig_z, &
+                        bbox, newbbox1, newbbox2, min_width, &
+                        didit)
+
+!!$  This is the last trick up our sleeve. So we just check if
+!!$  it succeeded so that we can correct the value of didit.
+
+!!$  write(*,*) "Split box at signature. Did it:", didit
+
+  if (didit > 0) then
+
+!!$    Prune the split boxes
+
+!!$    Box 1
+
+    call prune_new_box(nx, ny, nz, mask, bbox, newbbox1)
+
+!!$    Box 2
+
+    call prune_new_box(nx, ny, nz, mask, bbox, newbbox2)
+
+    didit = 2
+  end if
+
+!!$  Then we are done
+
+end subroutine check_box
+
+!!$ Copy the mask.
+!!$ Given two masks, one strictly contained within the other, copy
+!!$ the intersection to the smaller mask.
+!!$ 
+!!$ All "s" quantities are the "source" mask - the larger one
+!!$ All "d" quantities are the "destination" mask
+!!$ The location of the masks in "grid space" is denoted by 
+!!$ the bboxes.
+!!$ ?bbox(:,1) is the lower boundary
+!!$ ?bbox(:,2) is the upper boundary
+!!$ ?bbox(:,3) is the stride and is irrelevant
+
+subroutine copy_mask(snx, sny, snz, smask, sbbox, &
+                     dnx, dny, dnz, dmask, dbbox)
+
+  implicit none
+
+  CCTK_INT :: snx, sny, snz
+  CCTK_INT :: smask(snx, sny, snz)
+  CCTK_INT :: sbbox(3, 2)
+  CCTK_INT :: dnx, dny, dnz
+  CCTK_INT :: dmask(dnx, dny, dnz)
+  CCTK_INT :: dbbox(3, 2)
+  
+  CCTK_INT :: i, j, k, si, sj, sk, di, dj, dk
+  
+  if ( (dbbox(1,1) < sbbox(1,1)) .or. &
+       (dbbox(2,1) < sbbox(2,1)) .or. &
+       (dbbox(3,1) < sbbox(3,1)) .or. &
+       (dbbox(1,2) > sbbox(1,2)) .or. &
+       (dbbox(2,2) > sbbox(2,2)) .or. &
+       (dbbox(3,2) > sbbox(3,2)) ) then
+    call CCTK_WARN(0, &
+        "The destination mask is not contained in the source mask!")
+  end if
+
+  do k = dbbox(3,1), dbbox(3,2)
+    do j = dbbox(2,1), dbbox(2,2)
+      do i = dbbox(1,1), dbbox(1,2)
+        
+        si = 1 + i - sbbox(1,1)
+        sj = 1 + j - sbbox(2,1)
+        sk = 1 + k - sbbox(3,1)
+        
+        di = 1 + i - dbbox(1,1)
+        dj = 1 + j - dbbox(2,1)
+        dk = 1 + k - dbbox(3,1)
+        
+        dmask(di, dj, dk) = smask(si, sj, sk)
+
+!!$        write(*,*) "copying mask: loc",i,j,k
+!!$        write(*,*) "copying mask:d", di, dj, dk
+!!$        write(*,'("(",i2,":",i2,"):(",i2,":",i2,"):(",i2,":",i2,")")') &
+!!$             dbbox
+!!$        write(*,*) "copying mask:s", si, sj, sk
+!!$        write(*,'("(",i2,":",i2,"):(",i2,":",i2,"):(",i2,":",i2,")")') &
+!!$             sbbox
+        
+      end do
+    end do
+  end do
+
+end subroutine copy_mask
diff --git a/CarpetDev/CarpetAdaptiveRegrid/src/test/Makefile b/CarpetDev/CarpetAdaptiveRegrid/src/test/Makefile
new file mode 100644
index 000000000..0fa2fd877
--- /dev/null
+++ b/CarpetDev/CarpetAdaptiveRegrid/src/test/Makefile
@@ -0,0 +1,11 @@
+CARtest: CAR_test.F90 ../CAR_utils.F90 cctk.f cctk.h
+	ifc -c -g cctk.f
+	cp ../CAR_utils.F90 .
+	icc -E CAR_utils.F90 > CAR_utils.f90
+	ifc -c -g CAR_utils.f90
+	icc -E CAR_test.F90 > CAR_test.f90
+	ifc -c -g CAR_test.f90
+	ifc -g -o CARtest CAR_test.o CAR_utils.o cctk.o
+
+clean:
+	rm -f *.o *.f90 *~ CARtest
diff --git a/CarpetDev/CarpetAdaptiveRegrid/src/test/cctk.f b/CarpetDev/CarpetAdaptiveRegrid/src/test/cctk.f
new file mode 100644
index 000000000..ef5c17d4f
--- /dev/null
+++ b/CarpetDev/CarpetAdaptiveRegrid/src/test/cctk.f
@@ -0,0 +1,11 @@
+      subroutine CCTK_WARN(a,b)
+
+      implicit none
+
+      integer a
+      character*100 b
+      
+      write(*,*) b
+      stop
+
+      end subroutine CCTK_WARN
diff --git a/CarpetDev/CarpetAdaptiveRegrid/src/test/cctk.h b/CarpetDev/CarpetAdaptiveRegrid/src/test/cctk.h
new file mode 100644
index 000000000..25d4243c4
--- /dev/null
+++ b/CarpetDev/CarpetAdaptiveRegrid/src/test/cctk.h
@@ -0,0 +1,2 @@
+#define CCTK_INT integer
+#define CCTK_REAL real