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

16 files changed, 2220 insertions, 0 deletions

diff --git a/CarpetDev/CarpetAdaptiveRegrid/README b/CarpetDev/CarpetAdaptiveRegrid/README
new file mode 100644
index 000000000..43bc0dae1
--- /dev/null
+++ b/CarpetDev/CarpetAdaptiveRegrid/README
@@ -0,0 +1,13 @@
+CVS info   : $Header$
+
+Cactus Code Thorn CarpetAdaptiveRegrid
+Thorn Author(s)     : Ian Hawke <hawke@aei.mpg.de>
+Thorn Maintainer(s) : Ian Hawke <hawke@aei.mpg.de>
+--------------------------------------------------------------------------
+
+Purpose of the thorn:
+
+This is a Berger-Rigoutsos style clusterer for Carpet.
+
+Given a GF containing an error estimate it will return to Carpet a
+list of boxes for regridding purposes.
\ No newline at end of file
diff --git a/CarpetDev/CarpetAdaptiveRegrid/doc/documentation.tex b/CarpetDev/CarpetAdaptiveRegrid/doc/documentation.tex
new file mode 100644
index 000000000..4c947626c
--- /dev/null
+++ b/CarpetDev/CarpetAdaptiveRegrid/doc/documentation.tex
@@ -0,0 +1,144 @@
+% *======================================================================*
+%  Cactus Thorn template for ThornGuide documentation
+%  Author: Ian Kelley
+%  Date: Sun Jun 02, 2002
+%  $Header: /cactusdevcvs/Cactus/doc/ThornGuide/template.tex,v 1.12 2004/01/07 20:12:39 rideout Exp $
+%
+%  Thorn documentation in the latex file doc/documentation.tex
+%  will be included in ThornGuides built with the Cactus make system.
+%  The scripts employed by the make system automatically include
+%  pages about variables, parameters and scheduling parsed from the
+%  relevant thorn CCL files.
+%
+%  This template contains guidelines which help to assure that your
+%  documentation will be correctly added to ThornGuides. More
+%  information is available in the Cactus UsersGuide.
+%
+%  Guidelines:
+%   - Do not change anything before the line
+%       % START CACTUS THORNGUIDE",
+%     except for filling in the title, author, date, etc. fields.
+%        - Each of these fields should only be on ONE line.
+%        - Author names should be separated with a \\ or a comma.
+%   - You can define your own macros, but they must appear after
+%     the START CACTUS THORNGUIDE line, and must not redefine standard
+%     latex commands.
+%   - To avoid name clashes with other thorns, 'labels', 'citations',
+%     'references', and 'image' names should conform to the following
+%     convention:
+%       ARRANGEMENT_THORN_LABEL
+%     For example, an image wave.eps in the arrangement CactusWave and
+%     thorn WaveToyC should be renamed to CactusWave_WaveToyC_wave.eps
+%   - Graphics should only be included using the graphicx package.
+%     More specifically, with the "\includegraphics" command.  Do
+%     not specify any graphic file extensions in your .tex file. This
+%     will allow us to create a PDF version of the ThornGuide
+%     via pdflatex.
+%   - References should be included with the latex "\bibitem" command.
+%   - Use \begin{abstract}...\end{abstract} instead of \abstract{...}
+%   - Do not use \appendix, instead include any appendices you need as
+%     standard sections.
+%   - For the benefit of our Perl scripts, and for future extensions,
+%     please use simple latex.
+%
+% *======================================================================*
+%
+% Example of including a graphic image:
+%    \begin{figure}[ht]
+% 	\begin{center}
+%    	   \includegraphics[width=6cm]{MyArrangement_MyThorn_MyFigure}
+% 	\end{center}
+% 	\caption{Illustration of this and that}
+% 	\label{MyArrangement_MyThorn_MyLabel}
+%    \end{figure}
+%
+% Example of using a label:
+%   \label{MyArrangement_MyThorn_MyLabel}
+%
+% Example of a citation:
+%    \cite{MyArrangement_MyThorn_Author99}
+%
+% Example of including a reference
+%   \bibitem{MyArrangement_MyThorn_Author99}
+%   {J. Author, {\em The Title of the Book, Journal, or periodical}, 1 (1999),
+%   1--16. {\tt http://www.nowhere.com/}}
+%
+% *======================================================================*
+
+% If you are using CVS use this line to give version information
+% $Header: /cactusdevcvs/Cactus/doc/ThornGuide/template.tex,v 1.12 2004/01/07 20:12:39 rideout Exp $
+
+\documentclass{article}
+
+% Use the Cactus ThornGuide style file
+% (Automatically used from Cactus distribution, if you have a
+%  thorn without the Cactus Flesh download this from the Cactus
+%  homepage at www.cactuscode.org)
+\usepackage{../../../../doc/latex/cactus}
+
+\begin{document}
+
+% The author of the documentation
+\author{Ian Hawke \textless hawke@aei.mpg.de\textgreater}
+
+% The title of the document (not necessarily the name of the Thorn)
+\title{}
+
+% the date your document was last changed, if your document is in CVS,
+% please use:
+%    \date{$ $Date: 2004/01/07 20:12:39 $ $}
+\date{November 05 2004}
+
+\maketitle
+
+% Do not delete next line
+% START CACTUS THORNGUIDE
+
+% Add all definitions used in this documentation here
+%   \def\mydef etc
+
+% Add an abstract for this thorn's documentation
+\begin{abstract}
+
+\end{abstract}
+
+% The following sections are suggestive only.
+% Remove them or add your own.
+
+\section{Introduction}
+
+\section{Physical System}
+
+\section{Numerical Implementation}
+
+\section{Using This Thorn}
+
+\subsection{Obtaining This Thorn}
+
+\subsection{Basic Usage}
+
+\subsection{Special Behaviour}
+
+\subsection{Interaction With Other Thorns}
+
+\subsection{Examples}
+
+\subsection{Support and Feedback}
+
+\section{History}
+
+\subsection{Thorn Source Code}
+
+\subsection{Thorn Documentation}
+
+\subsection{Acknowledgements}
+
+
+\begin{thebibliography}{9}
+
+\end{thebibliography}
+
+% Do not delete next line
+% END CACTUS THORNGUIDE
+
+\end{document}
diff --git a/CarpetDev/CarpetAdaptiveRegrid/interface.ccl b/CarpetDev/CarpetAdaptiveRegrid/interface.ccl
new file mode 100644
index 000000000..fe14102e7
--- /dev/null
+++ b/CarpetDev/CarpetAdaptiveRegrid/interface.ccl
@@ -0,0 +1,64 @@
+# Interface definition for thorn CarpetAdaptiveRegrid
+# $Header$
+
+implements: CarpetAdaptiveRegrid
+inherits: Carpet CarpetLib
+
+uses include header: carpet.hh
+
+uses include header: defs.hh
+
+uses include header: bbox.hh
+uses include header: bboxset.hh
+uses include header: vect.hh
+
+uses include header: gf.hh
+uses include header: gh.hh
+
+# The location of the boundary points
+CCTK_INT FUNCTION GetBoundarySpecification \
+  (CCTK_INT IN size, \
+   CCTK_INT OUT ARRAY nboundaryzones, \
+   CCTK_INT OUT ARRAY is_internal, \
+   CCTK_INT OUT ARRAY is_staggered, \
+   CCTK_INT OUT ARRAY shiftout)
+USES FUNCTION GetBoundarySpecification
+
+# The overall size of the domain
+CCTK_INT FUNCTION GetDomainSpecification \
+  (CCTK_INT IN size, \
+   CCTK_REAL OUT ARRAY physical_min, \
+   CCTK_REAL OUT ARRAY physical_max, \
+   CCTK_REAL OUT ARRAY interior_min, \
+   CCTK_REAL OUT ARRAY interior_max, \
+   CCTK_REAL OUT ARRAY exterior_min, \
+   CCTK_REAL OUT ARRAY exterior_max, \
+   CCTK_REAL OUT ARRAY spacing)
+USES FUNCTION GetDomainSpecification
+
+# Convert between boundaries types
+CCTK_INT FUNCTION ConvertFromPhysicalBoundary \
+  (CCTK_INT IN size, \
+   CCTK_REAL IN  ARRAY physical_min, \
+   CCTK_REAL IN  ARRAY physical_max, \
+   CCTK_REAL OUT ARRAY interior_min, \
+   CCTK_REAL OUT ARRAY interior_max, \
+   CCTK_REAL OUT ARRAY exterior_min, \
+   CCTK_REAL OUT ARRAY exterior_max, \
+   CCTK_REAL IN  ARRAY spacing)
+USES FUNCTION ConvertFromPhysicalBoundary
+
+
+
+# The true prototype of the routine below:
+#   int Carpet_Regrid (const cGH * cctkGH,
+#                      gh<dim>::rexts  * bbsss,
+#                      gh<dim>::rbnds  * obss,
+#                      gh<dim>::rprocs * pss);
+CCTK_INT FUNCTION Carpet_Regrid (CCTK_POINTER_TO_CONST IN cctkGH,  \
+                                 CCTK_POINTER IN bbsss,            \
+                                 CCTK_POINTER IN obss,             \
+                                 CCTK_POINTER IN pss,		   \
+				 CCTK_INT IN force)
+PROVIDES FUNCTION Carpet_Regrid WITH CarpetAdaptiveRegrid_Regrid LANGUAGE C
+
diff --git a/CarpetDev/CarpetAdaptiveRegrid/param.ccl b/CarpetDev/CarpetAdaptiveRegrid/param.ccl
new file mode 100644
index 000000000..d853de204
--- /dev/null
+++ b/CarpetDev/CarpetAdaptiveRegrid/param.ccl
@@ -0,0 +1,30 @@
+# Parameter definitions for thorn CarpetAdaptiveRegrid
+# $Header$
+
+# Refinement criteria for automatic refining
+
+CCTK_INT minwidth "Minimum width of refined region" STEERABLE=always
+{
+  1:* :: "must be positive"
+} 8
+
+CCTK_REAL minfraction "Minimum fraction of points in need of refinement in a refined region" STEERABLE=always
+{
+  0:1 :: "must be positive and less than one"
+} 0.75
+
+CCTK_REAL maxerror "Maximum allowed error for non-refined grid points" STEERABLE=always
+{
+  *:* :: "everything goes"
+} 1.0
+
+CCTK_STRING errorvar "Name of grid function that contains the error" STEERABLE=always
+{
+  ".*" :: "must be the name of a grid function"
+} ""
+
+shares: driver
+
+USES INT ghost_size ""
+{
+}
diff --git a/CarpetDev/CarpetAdaptiveRegrid/schedule.ccl b/CarpetDev/CarpetAdaptiveRegrid/schedule.ccl
new file mode 100644
index 000000000..3ce3e0440
--- /dev/null
+++ b/CarpetDev/CarpetAdaptiveRegrid/schedule.ccl
@@ -0,0 +1,8 @@
+# Schedule definitions for thorn CarpetAdaptiveRegrid
+# $Header$
+
+schedule CarpetAdaptiveRegridParamcheck at PARAMCHECK
+{
+  LANG: C
+  OPTIONS: global
+} "Check Parameters"
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