[[This is a redo of my "cvs import" of 2002/06/11, this time using proper

cvs operations (commit/delete/add) to preserve the full CVS history of this
thorn.]]

This is a major cleanup/revision of AEIThorns/Exact.

Major user-visible changes:
* major expansion of doc/documentation.tex
* major expansion of documentation in param.ccl file
* rename all parameters, systematize spacetime/coordinate/parameter names
  (there is a perl script in par/convert-pars.pl to convert old parameter
  files to the new names)
* [from Mitica Vulcanov] many additions and fixes to
  cosmological solutions and Schwarzschild-Lemaitre
* fix stress-energy tensor computations so they work -- before they were
  all disabled in CVS (INCLUDES lines were commented out in interface.ccl)
  due to requiring excessive friendship with evolution thorns
  and/or public parameters; new code copies parameters to restricted
  grid scalars, which Cactus automagically "pushes" to friends
* added some more tests to testsuite, though these don't yet work fully

Additional internal changes:
* rename many Fortran subroutines (and a few C ones too)
  so their names start with the thorn name
  to reduce the chances of name collisions with other thorns
* move all metrics to subdirectory so the main source directory isn't
  so cluttered
* move two files containing subroutines which were never called
  (they didn't belong in this thorn, but somehow got into cvs by accident)
  into new archive/ directory
* some (small) improvements in efficiency -- the exact_model parameter
  is now decoded from a keyword (string) to an integer once at INITIAL,
  and that integer tested by the stress-energy tensor code,
  rather than requiring a separate series of string tests at each grid
  point (!) like the old stress-energy tensor code did

Modified Files:
  ParamCheck.c			added a check to make sure we don't try
				to set the shift if it doesn't have storage
  Startup.c
  make.code.defn
  slice_data.F
  slice_evolve.F
  slice_initialize.F

Added Files:
  Bona_Masso_data.F		moved from old  exactdata.F
  blended_boundary.F		moved from old  exactblendbound.F
  boundary.F			moved from old  exactboundary.F
  decode_pars.F			new file to decode exact_model into integer,
				copy parameters to grid scalars for Calc_Tmunu
				code
  gauge.F			moved from old  exactgauge.F
  initialize.F			moved from old  exactinitialize.F
  linear_extrap_one_bndry.F	moved from old  linextraponebound.F
  metric.F			moved from old  exactmetric.F
  xyz_blended_boundary.F	moved from old  exactcartblendbound.F


git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/Exact/trunk@101 e296648e-0e4f-0410-bd07-d597d9acff87

15 files changed, 1695 insertions, 89 deletions

diff --git a/src/Bona_Masso_data.F b/src/Bona_Masso_data.F
new file mode 100644
index 0000000..9805e0a
--- /dev/null
+++ b/src/Bona_Masso_data.F
@@ -0,0 +1,214 @@
+c This routine calculates Bona-Masso initial data, making use of the
+c subroutine Exact__metric() to calculate the spacetime metric and its
+c inverse.  Note that this use of the Bona-Masso variables is independent
+c of how (or even if) we are evolving the Einstein equations -- here
+c the Bona-Masso variables are "just" used as intermediate variables.
+
+c $Header$
+
+#include "cctk.h"
+
+      subroutine  Exact__Bona_Masso_data(
+     $     decoded_exact_model,
+     $     x, y, z, t,
+     $     gxx, gyy, gzz, gxy, gyz, gxz,
+     $     hxx, hyy, hzz, hxy, hyz, hxz,
+     $     dxgxx, dxgyy, dxgzz, dxgxy, dxgyz, dxgxz,
+     $     dygxx, dygyy, dygzz, dygxy, dygyz, dygxz,
+     $     dzgxx, dzgyy, dzgzz, dzgxy, dzgyz, dzgxz,
+     $     alp, ax, ay, az, betax, betay, betaz,
+     $     bxx, bxy, bxz, byx, byy, byz, bzx, bzy, bzz)
+
+      implicit none
+      CCTK_INT decoded_exact_model
+      CCTK_REAL x, y, z, t,
+     $       gxx, gyy, gzz, gxy, gyz, gxz,
+     $       hxx, hyy, hzz, hxy, hyz, hxz,
+     $       dxgxx, dxgyy, dxgzz, dxgxy, dxgyz, dxgxz,
+     $       dygxx, dygyy, dygzz, dygxy, dygyz, dygxz,
+     $       dzgxx, dzgyy, dzgzz, dzgxy, dzgyz, dzgxz,
+     $       alp, ax, ay, az, betax, betay, betaz,
+     $       bxx, bxy, bxz, byx, byy, byz, bzx, bzy, bzz
+
+C     gxx is g_xx etc.
+C     hxx is K_xx etc.
+C     dxgyy is (/2) dg_yy / dx (sic!)
+C     alp is N, betax is N^x etc.
+C     bxy is (/2) dN^y / dx (sic and sic!)
+C     ax is dN / dx / N (sic!)
+
+      CCTK_REAL eps, 
+     $       gdtt, gdtx, gdty, gdtz, 
+     $       gutt, gutx, guty, gutz, 
+     $       guxx, guyy, guzz, guxy, guyz, guxz,
+     $       gdtt_p, gdtx_p, gdty_p, gdtz_p, 
+     $       gdxx_p, gdyy_p, gdzz_p, gdxy_p, gdyz_p, gdxz_p,
+     $       gutt_p, gutx_p, guty_p, gutz_p, 
+     $       guxx_p, guyy_p, guzz_p, guxy_p, guyz_p, guxz_p,
+     $       gdtt_m, gdtx_m, gdty_m, gdtz_m, 
+     $       gdxx_m, gdyy_m, gdzz_m, gdxy_m, gdyz_m, gdxz_m,
+     $       gutt_m, gutx_m, guty_m, gutz_m, 
+     $       guxx_m, guyy_m, guzz_m, guxy_m, guyz_m, guxz_m
+
+      parameter (eps=1.d-6)
+
+C     Get the spacetime metric and its inverse at the base point.
+
+      call Exact__metric(
+     $     decoded_exact_model,
+     $     x, y, z, t,
+     $     gdtt, gdtx, gdty, gdtz, 
+     $     gxx, gyy, gzz, gxy, gyz, gxz,
+     $     gutt, gutx, guty, gutz, 
+     $     guxx, guyy, guzz, guxy, guyz, guxz)
+
+C     Calculate lapse and shift from the upper metric.
+
+      alp = 1.d0 / sqrt(-gutt)
+
+      betax = - gutx / gutt
+      betay = - guty / gutt
+      betaz = - gutz / gutt
+
+C     In order to calculate the derivatives of the lapse and shift from
+C     the contravariant metric, use g^tt = -1/N^2 and g^i0 = N^i / N^2
+C     Note that ax is dN/dx / N and that bxy is 1/2 dN^y / dx.
+
+C     Calculate x-derivatives. 
+
+      call Exact__metric(
+     $     decoded_exact_model,
+     $     x+eps, y, z, t,
+     $     gdtt_p, gdtx_p, gdty_p, gdtz_p,
+     $     gdxx_p, gdyy_p, gdzz_p, gdxy_p, gdyz_p, gdxz_p,
+     $     gutt_p, gutx_p, guty_p, gutz_p,
+     $     guxx_p, guyy_p, guzz_p, guxy_p, guyz_p, guxz_p)
+      call Exact__metric(
+     $     decoded_exact_model,
+     $     x-eps, y, z, t,
+     $     gdtt_m, gdtx_m, gdty_m, gdtz_m,
+     $     gdxx_m, gdyy_m, gdzz_m, gdxy_m, gdyz_m, gdxz_m,
+     $     gutt_m, gutx_m, guty_m, gutz_m,
+     $     guxx_m, guyy_m, guzz_m, guxy_m, guyz_m, guxz_m)
+
+      dxgxx = (gdxx_p  - gdxx_m)  / 4.d0 / eps
+      dxgyy = (gdyy_p  - gdyy_m)  / 4.d0 / eps
+      dxgzz = (gdzz_p  - gdzz_m)  / 4.d0 / eps
+      dxgxy = (gdxy_p  - gdxy_m)  / 4.d0 / eps
+      dxgyz = (gdyz_p  - gdyz_m)  / 4.d0 / eps
+      dxgxz = (gdxz_p  - gdxz_m)  / 4.d0 / eps
+
+      ax = - 0.5d0 * (gutt_p - gutt_m) / 2.d0 / eps / gutt
+
+      bxx = ((- gutx_p / gutt_p) - (- gutx_m / gutt_m)) / 4.d0 / eps
+      bxy = ((- guty_p / gutt_p) - (- guty_m / gutt_m)) / 4.d0 / eps
+      bxz = ((- gutz_p / gutt_p) - (- gutz_m / gutt_m)) / 4.d0 / eps
+
+C     Calculate y-derivatives.
+
+      call Exact__metric(
+     $     decoded_exact_model,
+     $     x, y+eps, z, t,
+     $     gdtt_p, gdtx_p, gdty_p, gdtz_p,
+     $     gdxx_p, gdyy_p, gdzz_p, gdxy_p, gdyz_p, gdxz_p,
+     $     gutt_p, gutx_p, guty_p, gutz_p,
+     $     guxx_p, guyy_p, guzz_p, guxy_p, guyz_p, guxz_p)
+      call Exact__metric(
+     $     decoded_exact_model,
+     $     x, y-eps, z, t,
+     $     gdtt_m, gdtx_m, gdty_m, gdtz_m,
+     $     gdxx_m, gdyy_m, gdzz_m, gdxy_m, gdyz_m, gdxz_m,
+     $     gutt_m, gutx_m, guty_m, gutz_m,
+     $     guxx_m, guyy_m, guzz_m, guxy_m, guyz_m, guxz_m)
+
+      dygxx = (gdxx_p  - gdxx_m)  / 4.d0 / eps
+      dygyy = (gdyy_p  - gdyy_m)  / 4.d0 / eps
+      dygzz = (gdzz_p  - gdzz_m)  / 4.d0 / eps
+      dygxy = (gdxy_p  - gdxy_m)  / 4.d0 / eps
+      dygyz = (gdyz_p  - gdyz_m)  / 4.d0 / eps
+      dygxz = (gdxz_p  - gdxz_m)  / 4.d0 / eps
+
+      ay = - 0.5d0 * (gutt_p - gutt_m) / 2.d0 / eps / gutt
+
+      byx = ((- gutx_p / gutt_p) - (- gutx_m / gutt_m)) / 4.d0 / eps
+      byy = ((- guty_p / gutt_p) - (- guty_m / gutt_m)) / 4.d0 / eps
+      byz = ((- gutz_p / gutt_p) - (- gutz_m / gutt_m)) / 4.d0 / eps
+
+C     Calculate z-derivatives.
+
+      call Exact__metric(
+     $     decoded_exact_model,
+     $     x, y, z+eps, t,
+     $     gdtt_p, gdtx_p, gdty_p, gdtz_p,
+     $     gdxx_p, gdyy_p, gdzz_p, gdxy_p, gdyz_p, gdxz_p,
+     $     gutt_p, gutx_p, guty_p, gutz_p,
+     $     guxx_p, guyy_p, guzz_p, guxy_p, guyz_p, guxz_p)
+      call Exact__metric(
+     $     decoded_exact_model,
+     $     x, y, z-eps, t,
+     $     gdtt_m, gdtx_m, gdty_m, gdtz_m,
+     $     gdxx_m, gdyy_m, gdzz_m, gdxy_m, gdyz_m, gdxz_m,
+     $     gutt_m, gutx_m, guty_m, gutz_m,
+     $     guxx_m, guyy_m, guzz_m, guxy_m, guyz_m, guxz_m)
+
+      dzgxx = (gdxx_p  - gdxx_m)  / 4.d0 / eps
+      dzgyy = (gdyy_p  - gdyy_m)  / 4.d0 / eps
+      dzgzz = (gdzz_p  - gdzz_m)  / 4.d0 / eps
+      dzgxy = (gdxy_p  - gdxy_m)  / 4.d0 / eps
+      dzgyz = (gdyz_p  - gdyz_m)  / 4.d0 / eps
+      dzgxz = (gdxz_p  - gdxz_m)  / 4.d0 / eps
+
+      az = - 0.5d0 * (gutt_p - gutt_m) / 2.d0 / eps / gutt
+
+      bzx = ((- gutx_p / gutt_p) - (- gutx_m / gutt_m)) / 4.d0 / eps
+      bzy = ((- guty_p / gutt_p) - (- guty_m / gutt_m)) / 4.d0 / eps
+      bzz = ((- gutz_p / gutt_p) - (- gutz_m / gutt_m)) / 4.d0 / eps
+
+C     Calculate t-derivatives, and extrinsic curvature.
+
+      call Exact__metric(
+     $     decoded_exact_model,
+     $     x, y, z, t+eps,
+     $     gdtt_p, gdtx_p, gdty_p, gdtz_p,
+     $     gdxx_p, gdyy_p, gdzz_p, gdxy_p, gdyz_p, gdxz_p,
+     $     gutt_p, gutx_p, guty_p, gutz_p,
+     $     guxx_p, guyy_p, guzz_p, guxy_p, guyz_p, guxz_p)
+      call Exact__metric(
+     $     decoded_exact_model,
+     $     x, y, z, t-eps,
+     $     gdtt_m, gdtx_m, gdty_m, gdtz_m,
+     $     gdxx_m, gdyy_m, gdzz_m, gdxy_m, gdyz_m, gdxz_m,
+     $     gutt_m, gutx_m, guty_m, gutz_m,
+     $     guxx_m, guyy_m, guzz_m, guxy_m, guyz_m, guxz_m)
+
+      hxx = - (gdxx_p - gdxx_m) / 4.d0 / eps / alp
+     $     + (dxgxx * betax + dygxx * betay + dzgxx * betaz
+     $     + 2.d0 * (bxx * gxx + bxy * gxy + bxz * gxz)) / alp
+
+      hyy = - (gdyy_p - gdyy_m) / 4.d0 / eps / alp
+     $     + (dxgyy * betax + dygyy * betay + dzgyy * betaz
+     $     + 2.d0 * (byx * gxy + byy * gyy + byz * gyz)) / alp
+
+      hzz = - (gdzz_p - gdzz_m) / 4.d0 / eps / alp
+     $     + (dxgzz * betax + dygzz * betay + dzgzz * betaz
+     $     + 2.d0 * (bzx * gxz + bzy * gyz + bzz * gzz)) / alp
+
+      hxy = - (gdxy_p - gdxy_m) / 4.d0 / eps / alp
+     $     + (dxgxy * betax + dygxy * betay + dzgxy * betaz
+     $     + bxx * gxy + bxy * gyy + bxz * gyz
+     $     + byx * gxx + byy * gxy + byz * gxz) / alp
+
+      hyz = - (gdyz_p - gdyz_m) / 4.d0 / eps / alp
+     $     + (dxgyz * betax + dygyz * betay + dzgyz * betaz
+     $     + byx * gxz + byy * gyz + byz * gzz
+     $     + bzx * gxy + bzy * gyy + bzz * gyz) / alp
+
+      hxz = - (gdxz_p - gdxz_m) / 4.d0 / eps / alp
+     $     + (dxgxz * betax + dygxz * betay + dzgxz * betaz      
+     $     + bxx * gxz + bxy * gyz + bxz * gzz
+     $     + bzx * gxx + bzy * gxy + bzz * gxz) / alp
+
+      return
+      end
+      
+
diff --git a/src/ParamCheck.c b/src/ParamCheck.c
index 24e5c62..27c10a4 100644
--- a/src/ParamCheck.c
+++ b/src/ParamCheck.c
@@ -53,13 +53,16 @@ void Exact_ParamCheck(CCTK_ARGUMENTS);
    @calls     
    @calledby   
    @history 
- 
+   @hdate      Tue Jun 11 18:25:49 CEST 2002
+   @hauthor    Jonathan Thornburg
+   @desc       Add a test that we're not trying to set the shift
+               with storage *not* present for it
+   @enddesc
    @endhistory 
  
  @@*/
 void Exact_ParamCheck(CCTK_ARGUMENTS)
 {
-  DECLARE_CCTK_ARGUMENTS;
   DECLARE_CCTK_PARAMETERS;
 
   if(! CCTK_EQUALS(metric_type, "physical"))
@@ -67,11 +70,11 @@ void Exact_ParamCheck(CCTK_ARGUMENTS)
     CCTK_PARAMWARN("Unknown ADMBase::metric_type - known types are \"physical\"");
   }
 
-  if(CCTK_EQUALS(shift_evolution_method, "exact") && ! CCTK_EQUALS(initial_shift, "exact"))
+  if ( CCTK_EQUALS(shift_evolution_method, "exact")
+       && CCTK_EQUALS(initial_shift, "none") )
   {
-    CCTK_PARAMWARN("Exact shift evolution requires exact shift initial data");
+    CCTK_PARAMWARN("can't set the shift if there's no storage for it!");
   }
-
 }
 
 /********************************************************************
diff --git a/src/Startup.c b/src/Startup.c
index 7f7d0d2..664e44e 100644
--- a/src/Startup.c
+++ b/src/Startup.c
@@ -1,9 +1,11 @@
+/* startup routine for Exact thorn */
+/* $Header$ */
 
 #include "cctk.h"
 #include "cctk_Parameters.h"
 #include "CactusEinstein/CoordGauge/src/Slicing.h"
 
-void Exact_RegisterSlicing(void) 
+void Exact__RegisterSlicing(void) 
 {
   int handle;
   handle=Einstein_RegisterSlicing("exact");
diff --git a/src/blended_boundary.F b/src/blended_boundary.F
new file mode 100644
index 0000000..f74beb3
--- /dev/null
+++ b/src/blended_boundary.F
@@ -0,0 +1,196 @@
+C $Header$
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+#include "cctk_Arguments.h"
+
+      subroutine Exact__blended_boundary(CCTK_ARGUMENTS)
+
+      implicit none
+
+      DECLARE_CCTK_ARGUMENTS
+      DECLARE_CCTK_PARAMETERS
+      DECLARE_CCTK_FUNCTIONS
+
+      logical doKij, doGij, doLapse, doShift
+
+      integer i,j,k
+      integer nx,ny,nz
+
+      CCTK_REAL router, rinner, frac, onemfrac
+      CCTK_REAL gxxe, gyye, gzze, gxye, gyze, gxze
+      CCTK_REAL kxxe, kyye, kzze, kxye, kyze, kxze
+      CCTK_REAL dxgxxe, dxgyye, dxgzze, dxgxye, dxgyze, dxgxze
+      CCTK_REAL dygxxe, dygyye, dygzze, dygxye, dygyze, dygxze
+      CCTK_REAL dzgxxe, dzgyye, dzgzze, dzgxye, dzgyze, dzgxze
+
+      CCTK_REAL alpe, axe, aye, aze
+      CCTK_REAL betaxe,betaye,betaze
+      CCTK_REAL bxxe,bxye,bxze,byxe,byye,byze,bzxe,bzye,bzze
+
+      CCTK_REAL det, uxx, uxy, uxz, uyy, uyz, uzz
+      CCTK_REAL vxe,vye,vze,sav
+
+      CCTK_REAL :: dx,dy,dz,time
+      CCTK_REAL :: ierr,xmx,xmn,ymx,ymn,zmx,zmn,rmx
+
+C     Grid parameters.
+
+      nx = cctk_lsh(1)
+      ny = cctk_lsh(2)
+      nz = cctk_lsh(3)
+
+      dx = cctk_delta_space(1)
+      dy = cctk_delta_space(2)
+      dz = cctk_delta_space(3)
+
+      time = cctk_time
+
+C     Other parameters.
+
+      doKij = (exblend_Ks.eq.1)
+      doGij = (exblend_gs.eq.1)
+
+      doLapse  = ((exblend_gauge.eq.1).and.
+     $            (CCTK_Equals(lapse_evolution_method,"exact").ne.0))
+      doShift  = ((exblend_gauge.eq.1).and.
+     $            (CCTK_Equals(shift_evolution_method,"exact").ne.0))
+
+      call CCTK_CoordRange(ierr,cctkGH,xmn,xmx,-1,"x","cart3d")
+      call CCTK_CoordRange(ierr,cctkGH,ymn,ymx,-1,"y","cart3d")
+      call CCTK_CoordRange(ierr,cctkGH,zmn,zmx,-1,"z","cart3d")
+
+      rmx = min(xmx,ymx,zmx)
+
+      if (exblend_rout.lt.0) then 
+         router = rmx - 2.0d0*dx
+      else
+         router = exblend_rout
+      endif
+
+      if (exblend_width.lt.0) then 
+         rinner = router + exblend_width*dx
+      else
+         rinner = router - exblend_width
+      endif
+
+      do k=1,nz
+         do j=1,ny
+            do i=1,nx
+
+c     We only do anything if r >= rinner so only evaluate exact data 
+c     there.
+
+               if (r(i,j,k) .ge. rinner) then
+
+                  call Exact__Bona_Masso_data(
+     $                 decoded_exact_model,
+     $                 x(i,j,k), y(i,j,k), z(i,j,k), time,
+     $                 gxxe, gyye, gzze, gxye, gyze, gxze,
+     $                 kxxe, kyye, kzze, kxye, kyze, kxze,
+     $                 dxgxxe, dxgyye, dxgzze, dxgxye, dxgyze, dxgxze,
+     $                 dygxxe, dygyye, dygzze, dygxye, dygyze, dygxze,
+     $                 dzgxxe, dzgyye, dzgzze, dzgxye, dzgyze, dzgxze,
+     $                 alpe, axe, aye, aze, betaxe, betaye, betaze,
+     $                 bxxe, bxye, bxze, byxe, 
+     $                 byye, byze, bzxe, bzye, bzze)
+
+c     This sucks, but we want the exact vs so we can blend them also.
+
+                  det = -(gxze**2*gyye) 
+     &                 + 2.d0*gxye*gxze*gyze 
+     &                 - gxxe*gyze**2 
+     &                 - gxye**2*gzze 
+     &                 + gxxe*gyye*gzze
+    
+                  uxx=(-gyze**2 + gyye*gzze)/det
+                  uxy=(gxze*gyze - gxye*gzze)/det
+                  uyy=(-gxze**2 + gxxe*gzze)/det
+                  uxz=(-gxze*gyye + gxye*gyze)/det
+                  uyz=(gxye*gxze - gxxe*gyze)/det
+                  uzz=(-gxye**2 + gxxe*gyye)/det
+
+c     Outside of router we want to place exact data on our grid
+
+                  if (r(i,j,k) .gt. router) then
+
+c     This is one of those things I will invariably screw up if I type
+c     it in so let the computer do it
+
+#define exassign(q) q(i,j,k) = q e
+#define exassign_grp(p) \
+        exassign(p xx) &&\
+        exassign(p xy) &&\
+        exassign(p xz) &&\
+        exassign(p yy) &&\
+        exassign(p yz) &&\
+        exassign(p zz)
+
+c     Note this plays on the nasty trick that fortran doesnt give a
+c     damn about spaces so gxx e is the same as gxxe for the parser...
+c     Grody but effective!
+
+                    if (doGij) then 
+                       exassign_grp(g)
+                    endif
+                    if (doKij) then 
+                       exassign_grp(k)
+                    endif
+
+                    if (doLapse) then 
+                       exassign(alp)
+                    endif
+                    
+                    if (doShift.and.(shift_state.ne.0)) then 
+                       exassign(betax)
+                       exassign(betay)
+                       exassign(betaz)
+                    endif
+
+c     OK so we dont want to blend so use a goto to jump.
+                 else
+
+c     Evaluate the linear weighting fraction. Obvious...
+
+                    frac = (r(i,j,k) - rinner) / (router - rinner)
+                    onemfrac = 1.0D0 - frac
+
+c     Once again some c-preprocessor tricks based on the whole fortran
+c     space thing...
+
+#define INTPOINT(f,v) f(i,j,k) = frac * v + onemfrac * f(i,j,k)
+#define intone(f) INTPOINT(f, f e)
+#define int_grp(p) \
+        intone(p xx) &&\
+        intone(p xy) &&\
+        intone(p xz) &&\
+        intone(p yy) &&\
+        intone(p yz) &&\
+        intone(p zz)
+
+                    if (doGij) then 
+                       int_grp(g)
+                    endif
+                    if (doKij) then
+                       int_grp(k)
+                    endif
+
+                    if (doLapse) then 
+                       intone(alp)
+                    endif
+
+                    if (doShift.and.(shift_state.ne.0)) then 
+                       intone(betax)
+                       intone(betay)
+                       intone(betaz)
+                    endif
+
+                 endif          ! r > router else
+              endif             ! r > rinner
+
+            enddo
+         enddo
+      enddo
+
+      return
+      end 
diff --git a/src/boundary.F b/src/boundary.F
new file mode 100644
index 0000000..1862e43
--- /dev/null
+++ b/src/boundary.F
@@ -0,0 +1,125 @@
+C $Header$
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+#include "cctk_Arguments.h"
+
+      subroutine Exact__boundary(CCTK_ARGUMENTS)
+
+      implicit none
+
+      DECLARE_CCTK_ARGUMENTS
+
+      integer i,j,k
+      integer nx,ny,nz
+
+      CCTK_REAL tplusone
+      CCTK_REAL 
+     $     dxgxxjunk, dxgyyjunk, dxgzzjunk, 
+     $     dxgxyjunk, dxgyzjunk, dxgxzjunk,
+     $     dygxxjunk, dygyyjunk, dygzzjunk, 
+     $     dygxyjunk, dygyzjunk, dygxzjunk,
+     $     dzgxxjunk, dzgyyjunk, dzgzzjunk, 
+     $     dzgxyjunk, dzgyzjunk, dzgxzjunk,
+     $     axjunk, ayjunk, azjunk, 
+     $     bxxjunk, bxyjunk, bxzjunk, 
+     $     byxjunk, byyjunk, byzjunk, 
+     $     bzxjunk, bzyjunk, bzzjunk
+
+C     Grid parameters.
+
+      nx = cctk_lsh(1)
+      ny = cctk_lsh(2)
+      nz = cctk_lsh(3)
+
+C     Set all initial data including dijk and vi on all points which 
+C     are on the boundary of the domain if it really is the boundary
+C     of the complete grid. Treat all six sides of the grid cube this way.
+
+c     Set t = time + dt.  This is necessary here because by the time
+c     we reach this point the geometry has been evolved one time step
+c     but the variable `time' still hasn't been updated. 
+
+      tplusone = cctk_time + cctk_delta_time
+
+C     Note we also always set the lapse and shift at the boundaries at
+C     time t+1. This is to provide boundary conditions for testing
+C     elliptic gauge conditions. If they are not used, they will be
+C     overwritten by Exact__gauge.
+
+#define EXACTDATAPOINT                                          \
+               call Exact__Bona_Masso_data(                     \
+                   decoded_exact_model,                         \
+                   x(i,j,k), y(i,j,k), z(i,j,k), tplusone,      \
+                   gxx(i,j,k), gyy(i,j,k), gzz(i,j,k),          \
+                   gxy(i,j,k), gyz(i,j,k), gxz(i,j,k),          \
+                   kxx(i,j,k), kyy(i,j,k), kzz(i,j,k),          \
+                   kxy(i,j,k), kyz(i,j,k), kxz(i,j,k),          \
+                   dxgxxjunk, dxgyyjunk, dxgzzjunk,             \
+                   dxgxyjunk, dxgyzjunk, dxgxzjunk,             \
+                   dygxxjunk, dygyyjunk, dygzzjunk,             \
+                   dygxyjunk, dygyzjunk, dygxzjunk,             \
+                   dzgxxjunk, dzgyyjunk, dzgzzjunk,             \
+                   dzgxyjunk, dzgyzjunk, dzgxzjunk,             \
+                   alp(i,j,k), axjunk, ayjunk, azjunk,          \
+                   betax(i,j,k), betay(i,j,k), betaz(i,j,k),    \
+                   bxxjunk, bxyjunk, bxzjunk,                   \
+                   byxjunk, byyjunk, byzjunk,                   \
+                   bzxjunk, bzyjunk, bzzjunk)
+
+      if (cctk_bbox(1) .eq. 1) then 
+         i=1
+         do j=1,ny
+            do k=1,nz
+               EXACTDATAPOINT
+            end do
+         end do
+      end if
+
+      if (cctk_bbox(2) .eq. 1) then 
+         i=nx
+         do j=1,ny
+            do k=1,nz
+               EXACTDATAPOINT
+            end do
+         end do
+      end if
+
+      if (cctk_bbox(3) .eq. 1) then 
+         j=1
+         do i=1,nx
+            do k=1,nz
+               EXACTDATAPOINT
+            end do
+         end do
+      end if
+
+      if (cctk_bbox(4) .eq. 1) then 
+         j=ny
+         do i=1,nx
+            do k=1,nz
+               EXACTDATAPOINT
+           end do
+         end do
+      end if
+
+      if (cctk_bbox(5) .eq. 1) then 
+         k=1
+         do j=1,ny
+            do i=1,nx
+               EXACTDATAPOINT
+            end do
+         end do
+      end if
+
+      if (cctk_bbox(6) .eq. 1) then 
+         k=nz
+         do j=1,ny
+            do i=1,nx
+               EXACTDATAPOINT
+            end do
+         end do
+      end if
+
+      return
+      end
diff --git a/src/decode_pars.F b/src/decode_pars.F
new file mode 100644
index 0000000..e81fd8a
--- /dev/null
+++ b/src/decode_pars.F
@@ -0,0 +1,201 @@
+c/*@@
+c @file		decode_pars.F
+c @date		Fri Jun  7 19:47:46 CEST 2002
+c @author	Jonathan Thornburg <jthorn@aei.mpg.de>
+c @desc
+c	Decode/copy parameters for this thorn into grid scalars
+c	so we can share this with friends
+c	so the Calc_Tmunu code in ../include/Scalar_CalcTmunu.inc
+c	can use them in computing the stress-energy tensor
+c
+c	Actually we only have to copy those parameters which are
+c	used by the Calc_Tmunu code.  For simplicity we decode
+c	exact_model into the integer decoded_exact_model, then
+c	copy only the per-model parameters for those models which
+c	have stress-energy tensor code in ../include/Scalar_CalcTmunu.inc .
+c
+c @enddesc
+c @version	$Header$
+c@@*/
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+#include "cctk_Arguments.h"
+
+#include "param_defs.inc"
+
+c/*@@
+c @routine	Exact__decode_pars
+c @date		Fri Jun  7 19:47:46 CEST 2002
+c @author	Jonathan Thornburg <jthorn@aei.mpg.de>
+c @desc
+c	Decode/copy parameters for this thorn into grid scalars, so
+c	we can share this with friends for the use of the Calc_Tmunu
+c	code in ../include/Scalar_CalcTmunu.inc in computing the
+c	stress-energy tensor.
+c @enddesc
+c @version	$Header$
+c@@*/
+	subroutine Exact__decode_pars(CCTK_ARGUMENTS)
+	implicit none
+	DECLARE_CCTK_ARGUMENTS
+	DECLARE_CCTK_PARAMETERS
+	DECLARE_CCTK_FUNCTIONS
+
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+c
+c decode  exact_model  into the integer  decoded_exact_model
+c
+
+c Minkowski spacetime
+	if      (CCTK_Equals(exact_model, "Minkowski") .ne. 0) then
+	     decoded_exact_model = EXACT__Minkowski
+	elseif (CCTK_Equals(exact_model, "Minkowski/shift") .ne. 0) then
+	     decoded_exact_model = EXACT__Minkowski_shift
+	elseif (CCTK_Equals(exact_model, "Minkowski/funny") .ne. 0) then
+	     decoded_exact_model = EXACT__Minkowski_funny
+	elseif (CCTK_Equals(exact_model, "Minkowski/gauge wave") .ne. 0) then
+	     decoded_exact_model = EXACT__Minkowski_gauge_wave
+
+c black hole spacetimes
+	elseif (CCTK_Equals(exact_model, "Schwarzschild/EF") .ne. 0) then
+	     decoded_exact_model = EXACT__Schwarzschild_EF
+	elseif (CCTK_Equals(exact_model, "Schwarzschild/PG") .ne. 0) then
+	     decoded_exact_model = EXACT__Schwarzschild_PG
+	elseif (CCTK_Equals(exact_model, "Schwarzschild/Novikov") .ne. 0) then
+	     decoded_exact_model = EXACT__Schwarzschild_Novikov
+	elseif (CCTK_Equals(exact_model, "Kerr/Boyer-Lindquist") .ne. 0) then
+	     decoded_exact_model = EXACT__Kerr_BoyerLindquist
+	elseif (CCTK_Equals(exact_model, "Kerr/Kerr-Schild") .ne. 0) then
+	     decoded_exact_model = EXACT__Kerr_KerrSchild
+	elseif (CCTK_Equals(exact_model, "Schwarzschild-Lemaitre") .ne. 0) then
+	     decoded_exact_model = EXACT__Schwarzschild_Lemaitre
+	elseif (CCTK_Equals(exact_model, "multi-BH") .ne. 0) then
+	     decoded_exact_model = EXACT__multi_BH
+	elseif (CCTK_Equals(exact_model, "Alvi") .ne. 0) then
+	     decoded_exact_model = EXACT__Alvi
+	elseif (CCTK_Equals(exact_model, "Thorne-fakebinary") .ne. 0) then
+	     decoded_exact_model = EXACT__Thorne_fakebinary
+
+c cosmological spacetimes
+	elseif (CCTK_Equals(exact_model, "Lemaitre") .ne. 0) then
+	     decoded_exact_model = EXACT__Lemaitre
+	elseif (CCTK_Equals(exact_model, "Robertson-Walker") .ne. 0) then
+	     decoded_exact_model = EXACT__Robertson_Walker
+	elseif (CCTK_Equals(exact_model, "de Sitter") .ne. 0) then
+	     decoded_exact_model = EXACT__de_Sitter
+	elseif (CCTK_Equals(exact_model, "de Sitter+Lambda") .ne. 0) then
+	     decoded_exact_model = EXACT__de_Sitter_Lambda
+	elseif (CCTK_Equals(exact_model, "anti-de Sitter+Lambda") .ne. 0) then
+	     decoded_exact_model = EXACT__anti_de_Sitter_Lambda
+	elseif (CCTK_Equals(exact_model, "Bianchi I") .ne. 0) then
+	     decoded_exact_model = EXACT__Bianchi_I
+	elseif (CCTK_Equals(exact_model, "Goedel") .ne. 0) then
+	     decoded_exact_model = EXACT__Goedel
+	elseif (CCTK_Equals(exact_model, "Bertotti") .ne. 0) then
+	     decoded_exact_model = EXACT__Bertotti
+	elseif (CCTK_Equals(exact_model, "Kasner-like") .ne. 0) then
+	     decoded_exact_model = EXACT__Kasner_like
+	elseif (CCTK_Equals(exact_model, "Kasner-axisymmetric") .ne. 0) then
+	     decoded_exact_model = EXACT__Kasner_axisymmetric
+	elseif (CCTK_Equals(exact_model, "Kasner-generalized") .ne. 0) then
+	     decoded_exact_model = EXACT__Kasner_generalized
+	elseif (CCTK_Equals(exact_model, "Milne") .ne. 0) then
+	     decoded_exact_model = EXACT__Milne
+
+c miscellaneous spacetimes
+	elseif (CCTK_Equals(exact_model, "boost-rotation symmetric") .ne. 0) then
+	     decoded_exact_model = EXACT__boost_rotation_symmetric
+	elseif (CCTK_Equals(exact_model, "bowl") .ne. 0) then
+	     decoded_exact_model = EXACT__bowl
+	elseif (CCTK_Equals(exact_model, "constant density star") .ne. 0) then
+	     decoded_exact_model = EXACT__constant_density_star
+	else
+	     call CCTK_WARN(0, "invalid exact_model parameter!")
+	endif
+
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+c
+c parameters for Schwarzschild-Lemaitre metric
+c (Schwarzschiled black hole with cosmological constant)
+c
+	Schwarzschild_Lemaitre___Lambda = Schwarzschild_Lemaitre__Lambda
+	Schwarzschild_Lemaitre___mass   = Schwarzschild_Lemaitre__mass
+
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+c
+c parameters for Lemaitre-type spacetime
+c
+	Lemaitre___kappa    = Lemaitre__kappa
+	Lemaitre___Lambda   = Lemaitre__Lambda
+	Lemaitre___epsilon0 = Lemaitre__epsilon0
+	Lemaitre___R0       = Lemaitre__R0
+
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+c
+c parameters for Robertson-Walker spacetime
+c
+	Robertson_Walker___R0       = Robertson_Walker__R0
+	Robertson_Walker___rho      = Robertson_Walker__rho
+	Robertson_Walker___k        = Robertson_Walker__k
+	Robertson_Walker___pressure = Robertson_Walker__pressure
+
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+c
+c parameters for de Sitter spacetime
+c
+	de_Sitter___scale = de_Sitter__scale
+
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+c
+c parameters for de Sitter spacetime with cosmological constant
+c
+	de_Sitter_Lambda___scale = de_Sitter_Lambda__scale
+
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+c
+c parameters for anti-de Sitter spacetime with cosmological constant
+c
+	anti_de_Sitter_Lambda___scale = anti_de_Sitter_Lambda__scale
+
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+c
+c parameters for Bertotti spacetime
+c
+	Bertotti___Lambda = Bertotti__Lambda
+
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+c
+c parameters for Kasner-like spacetime
+c
+	Kasner_like___q = Kasner_like__q
+
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+c
+c parameters for generalized Kasner spacetime
+c
+	Kasner_generalized___p1 = Kasner_generalized__p1
+	Kasner_generalized___p2 = Kasner_generalized__p2
+
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+c
+c parameters for constant density (Schwarzschild) star
+c
+	constant_density_star___mass   = constant_density_star__mass
+	constant_density_star___radius = constant_density_star__radius
+
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+	return
+	end
diff --git a/src/gauge.F b/src/gauge.F
new file mode 100644
index 0000000..f5c2e88
--- /dev/null
+++ b/src/gauge.F
@@ -0,0 +1,152 @@
+C     This routine sets the lapse and/or shift by calling a routine
+C     that does it pointwise. Note that it could be easily modified
+C     to set the Bona-Masso variables B_xx etc. 
+C $Header$
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+#include "cctk_Arguments.h"
+
+      subroutine Exact__gauge(CCTK_ARGUMENTS)
+
+      implicit none
+
+      DECLARE_CCTK_ARGUMENTS
+      DECLARE_CCTK_PARAMETERS
+      DECLARE_CCTK_FUNCTIONS
+
+      integer i,j,k
+      integer nx,ny,nz
+
+      CCTK_REAL tplushalf
+      CCTK_REAL gxxjunk, gyyjunk, gzzjunk, 
+     $     gxyjunk, gyzjunk, gxzjunk,
+     $     hxxjunk, hyyjunk, hzzjunk, 
+     $     hxyjunk, hyzjunk, hxzjunk,
+     $     dxgxxjunk, dxgyyjunk, dxgzzjunk, 
+     $     dxgxyjunk, dxgyzjunk, dxgxzjunk,
+     $     dygxxjunk, dygyyjunk, dygzzjunk, 
+     $     dygxyjunk, dygyzjunk, dygxzjunk,
+     $     dzgxxjunk, dzgyyjunk, dzgzzjunk, 
+     $     dzgxyjunk, dzgyzjunk, dzgxzjunk,
+     $     alpjunk, axjunk, ayjunk, azjunk, 
+     $     betaxjunk, betayjunk, betazjunk,
+     $     bxxjunk, bxyjunk, bxzjunk, 
+     $     byxjunk, byyjunk, byzjunk, 
+     $     bzxjunk, bzyjunk, bzzjunk
+
+C     Grid parameters.
+
+      nx = cctk_lsh(1)
+      ny = cctk_lsh(2)
+      nz = cctk_lsh(3)
+
+C     Set t = time + dt/2.  This is for consistency with a second
+C     order numerical scheme.
+
+      tplushalf = cctk_time + 0.5D0*cctk_delta_time
+
+C     Set both lapse and shift.
+
+      if (( (CCTK_Equals(lapse_evolution_method,"exact").ne.0) .and.
+     $      (CCTK_Equals(shift_evolution_method,"exact").ne.0) )
+     $    .or.
+     $    ( (CCTK_Equals(initial_lapse,"exact").ne.0) .and.
+     $      (CCTK_Equals(initial_shift,"exact").ne.0) )) then
+
+         do k=1,nz
+            do j=1,ny
+               do i=1,nx
+
+                  call Exact__Bona_Masso_data(
+     $                 decoded_exact_model,
+     $                 x(i,j,k), y(i,j,k), z(i,j,k), tplushalf,
+     $                 gxxjunk, gyyjunk, gzzjunk, 
+     $                 gxyjunk, gyzjunk, gxzjunk,
+     $                 hxxjunk, hyyjunk, hzzjunk, 
+     $                 hxyjunk, hyzjunk, hxzjunk,
+     $                 dxgxxjunk, dxgyyjunk, dxgzzjunk, 
+     $                 dxgxyjunk, dxgyzjunk, dxgxzjunk,
+     $                 dygxxjunk, dygyyjunk, dygzzjunk, 
+     $                 dygxyjunk, dygyzjunk, dygxzjunk,
+     $                 dzgxxjunk, dzgyyjunk, dzgzzjunk, 
+     $                 dzgxyjunk, dzgyzjunk, dzgxzjunk,
+     $                 alp(i,j,k), axjunk, ayjunk, azjunk, 
+     $                 betax(i,j,k), betay(i,j,k), betaz(i,j,k),
+     $                 bxxjunk, bxyjunk, bxzjunk, 
+     $                 byxjunk, byyjunk, byzjunk, 
+     $                 bzxjunk, bzyjunk, bzzjunk)
+               end do
+            end do
+         end do
+
+C     Set lapse only.
+
+      else if ((CCTK_Equals(lapse_evolution_method,"exact").ne.0)
+     $        .or. (CCTK_Equals(initial_lapse,"exact").ne.0)) then
+
+         do k=1,nz
+            do j=1,ny
+               do i=1,nx
+
+                  call Exact__Bona_Masso_data(
+     $                 decoded_exact_model,
+     $                 x(i,j,k), y(i,j,k), z(i,j,k), tplushalf,
+     $                 gxxjunk, gyyjunk, gzzjunk, 
+     $                 gxyjunk, gyzjunk, gxzjunk,
+     $                 hxxjunk, hyyjunk, hzzjunk, 
+     $                 hxyjunk, hyzjunk, hxzjunk,
+     $                 dxgxxjunk, dxgyyjunk, dxgzzjunk, 
+     $                 dxgxyjunk, dxgyzjunk, dxgxzjunk,
+     $                 dygxxjunk, dygyyjunk, dygzzjunk, 
+     $                 dygxyjunk, dygyzjunk, dygxzjunk,
+     $                 dzgxxjunk, dzgyyjunk, dzgzzjunk, 
+     $                 dzgxyjunk, dzgyzjunk, dzgxzjunk,
+     $                 alp(i,j,k), axjunk, ayjunk, azjunk, 
+     $                 betaxjunk, betayjunk, betazjunk,
+     $                 bxxjunk, bxyjunk, bxzjunk, 
+     $                 byxjunk, byyjunk, byzjunk, 
+     $                 bzxjunk, bzyjunk, bzzjunk)
+
+               end do
+            end do
+         end do
+
+C     Set shift only.
+
+      else if ((CCTK_Equals(shift_evolution_method,"exact").ne.0) 
+     $        .or. (CCTK_Equals(initial_shift,"exact").ne.0)) then
+
+         do k=1,nz
+            do j=1,ny
+               do i=1,nx
+
+                  call Exact__Bona_Masso_data(
+     $                 decoded_exact_model,
+     $                 x(i,j,k), y(i,j,k), z(i,j,k), tplushalf,
+     $                 gxxjunk, gyyjunk, gzzjunk, 
+     $                 gxyjunk, gyzjunk, gxzjunk,
+     $                 hxxjunk, hyyjunk, hzzjunk, 
+     $                 hxyjunk, hyzjunk, hxzjunk,
+     $                 dxgxxjunk, dxgyyjunk, dxgzzjunk, 
+     $                 dxgxyjunk, dxgyzjunk, dxgxzjunk,
+     $                 dygxxjunk, dygyyjunk, dygzzjunk, 
+     $                 dygxyjunk, dygyzjunk, dygxzjunk,
+     $                 dzgxxjunk, dzgyyjunk, dzgzzjunk, 
+     $                 dzgxyjunk, dzgyzjunk, dzgxzjunk,
+     $                 alpjunk, axjunk, ayjunk, azjunk, 
+     $                 betax(i,j,k), betay(i,j,k), betaz(i,j,k),
+     $                 bxxjunk, bxyjunk, bxzjunk, 
+     $                 byxjunk, byyjunk, byzjunk, 
+     $                 bzxjunk, bzyjunk, bzzjunk)
+
+               end do
+            end do
+         end do
+
+      else
+         call CCTK_WARN(1,'Exact__gauge has been called without doing anything')
+      end if
+
+      return
+      end
diff --git a/src/initialize.F b/src/initialize.F
new file mode 100644
index 0000000..69a4c75
--- /dev/null
+++ b/src/initialize.F
@@ -0,0 +1,76 @@
+C     Wrapper for boostrotdata. Calls it and vectorini.
+C     Sets Cauchy data, lapse and shift, and what else is needed
+C     in the Bona-Masso formalism, at an initial time.
+C $Header$
+      
+#include "cctk.h"
+#include "cctk_Parameters.h"
+#include "cctk_Arguments.h"
+
+      subroutine Exact__initialize(CCTK_ARGUMENTS)
+
+      implicit none
+
+      DECLARE_CCTK_ARGUMENTS
+      DECLARE_CCTK_PARAMETERS
+      DECLARE_CCTK_FUNCTIONS
+
+      integer i,j,k
+      integer nx,ny,nz
+
+      CCTK_REAL alpjunk, axjunk, ayjunk, azjunk, 
+     $     betaxjunk, betayjunk, betazjunk,
+     $     bxxjunk, bxyjunk, bxzjunk, 
+     $     byxjunk, byyjunk, byzjunk, 
+     $     bzxjunk, bzyjunk, bzzjunk
+      CCTK_REAL 
+     $     dxgxxjunk, dxgyyjunk, dxgzzjunk, 
+     $     dxgxyjunk, dxgyzjunk, dxgxzjunk,
+     $     dygxxjunk, dygyyjunk, dygzzjunk, 
+     $     dygxyjunk, dygyzjunk, dygxzjunk,
+     $     dzgxxjunk, dzgyyjunk, dzgzzjunk, 
+     $     dzgxyjunk, dzgyzjunk, dzgxzjunk
+
+C     Note I assume time has been initialized to physical time.
+C     Set data pointwise.
+
+      nx = cctk_lsh(1)
+      ny = cctk_lsh(2)
+      nz = cctk_lsh(3)
+
+      do k=1,nz
+         do j=1,ny
+            do i=1,nx
+               
+               call Exact__Bona_Masso_data(
+     $              decoded_exact_model,
+     $              x(i,j,k), y(i,j,k), z(i,j,k), cctk_time,
+     $              gxx(i,j,k), gyy(i,j,k), gzz(i,j,k), 
+     $              gxy(i,j,k), gyz(i,j,k), gxz(i,j,k),
+     $              kxx(i,j,k), kyy(i,j,k), kzz(i,j,k), 
+     $              kxy(i,j,k), kyz(i,j,k), kxz(i,j,k),
+     $              dxgxxjunk, dxgyyjunk, dxgzzjunk, 
+     $              dxgxyjunk, dxgyzjunk, dxgxzjunk,
+     $              dygxxjunk, dygyyjunk, dygzzjunk, 
+     $              dygxyjunk, dygyzjunk, dygxzjunk,
+     $              dzgxxjunk, dzgyyjunk, dzgzzjunk, 
+     $              dzgxyjunk, dzgyzjunk, dzgxzjunk,
+     $              alpjunk, axjunk, ayjunk, azjunk, 
+     $              betaxjunk, betayjunk, betazjunk,
+     $              bxxjunk, bxyjunk, bxzjunk,
+     $              byxjunk, byyjunk, byzjunk,
+     $              bzxjunk, bzyjunk, bzzjunk)
+
+            end do
+         end do
+      end do
+      
+C     Tell the code there is no need to treat the conformal factor
+C     as a separate field. That is, we have set the physical metric here.
+c     Commented out in einstein revamp, now Exact doesnot inherit anything
+c     about the conformal factor
+c      conformal_state = 0
+c      psi = 1.0D0
+
+      return
+      end
diff --git a/src/linear_extrap_one_bndry.F b/src/linear_extrap_one_bndry.F
new file mode 100644
index 0000000..e10e770
--- /dev/null
+++ b/src/linear_extrap_one_bndry.F
@@ -0,0 +1,144 @@
+c this subroutine linearly extrapolates one Cactus variable
+c on one boundary of the Cactus grid box
+C $Header$
+
+#include "cctk.h"
+#include "cctk_Arguments.h"
+
+      subroutine Exact__linear_extrap_one_bndry(CCTK_ARGUMENTS,var)
+
+      implicit none
+
+      DECLARE_CCTK_ARGUMENTS
+
+      integer nx,ny,nz
+
+      CCTK_REAL var(cctk_lsh(1),cctk_lsh(2),cctk_lsh(3))
+
+C     Grid parameters.
+
+      nx = cctk_lsh(1)
+      ny = cctk_lsh(2)
+      nz = cctk_lsh(3)
+
+C     Linear extrapolation from the interiors to the boundaries.
+C     Does not support octant or quadrant.
+
+C     6 faces.
+
+      if (cctk_bbox(1).eq.1 .and. nx.ge.4 ) then 
+         var(1,:,:) = 2.d0 * var(2,:,:) - var(3,:,:)
+      endif
+      if (cctk_bbox(3).eq.1 .and. ny.ge.4 ) then 
+         var(:,1,:) = 2.d0 * var(:,2,:) - var(:,3,:)
+      endif
+      if (cctk_bbox(5).eq.1 .and. nz.ge.4 ) then 
+         var(:,:,1) = 2.d0 * var(:,:,2) - var(:,:,3)
+      endif
+      if (cctk_bbox(2).eq.1 .and. nx.ge.4 ) then 
+         var(nx,:,:) = 2.d0 * var(nx-1,:,:) - var(nx-2,:,:)
+      endif
+      if (cctk_bbox(4).eq.1 .and. ny.ge.4 ) then 
+         var(:,ny,:) = 2.d0 * var(:,ny-1,:) - var(:,ny-2,:)
+      endif
+      if (cctk_bbox(6).eq.1 .and. nz.ge.4 ) then 
+         var(:,:,nz) = 2.d0 * var(:,:,nz-1) - var(:,:,nz-2)
+      endif
+
+C     12 edges.
+C     4 round face x=min.
+      if ( cctk_bbox(1).eq.1 .and. nx.ge.4 
+     $     .and. cctk_bbox(3).eq.1 .and. ny.ge.4 ) then 
+         var(1,1,:) = 2.d0 * var(2,2,:) - var(3,3,:)
+      end if
+
+      if ( cctk_bbox(1).eq.1 .and. nx.ge.4 
+     $     .and. cctk_bbox(4).eq.1 .and. ny.ge.4 ) then 
+         var(1,ny,:) = 2.d0 * var(2,ny-1,:) - var(2,ny-2,:)
+      end if
+
+      if ( cctk_bbox(1).eq.1 .and. nx.ge.4 
+     $     .and. cctk_bbox(5).eq.1 .and. nz.ge.4 ) then 
+         var(1,:,1) = 2.d0 * var(2,:,2) - var(3,:,3)
+      end if
+
+      if ( cctk_bbox(1).eq.1 .and. nx.ge.4 
+     $     .and. cctk_bbox(6).eq.1 .and. nz.ge.4 ) then 
+         var(1,:,nz) = 2.d0 * var(2,:,nz-1) - var(3,:,nz-2)
+      end if
+
+C     4 around face x=max.
+      if ( cctk_bbox(2).eq.1 .and. nx.ge.4 
+     $     .and. cctk_bbox(3).eq.1 .and. ny.ge.4 ) then 
+         var(nx,1,:) = 2.d0 * var(nx-1,2,:) - var(nx-2,3,:)
+      end if
+
+      if ( cctk_bbox(2).eq.1 .and. nx.ge.4 
+     $     .and. cctk_bbox(4).eq.1 .and. ny.ge.4 ) then 
+         var(nx,ny,:) = 2.d0 * var(nx-1,ny-1,:) - var(nx-2,ny-2,:)
+      end if
+
+      if ( cctk_bbox(2).eq.1 .and. nx.ge.4 
+     $     .and. cctk_bbox(5).eq.1 .and. nz.ge.4 ) then 
+         var(nx,:,1) = 2.d0 * var(nx-1,:,2) - var(nx-2,:,3)
+      end if
+
+      if ( cctk_bbox(2).eq.1 .and. nx.ge.4 
+     $     .and. cctk_bbox(6).eq.1 .and. nz.ge.4 ) then 
+         var(nx,:,nz) = 2.d0 * var(nx-1,:,nz-1) - var(nx-2,:,nz-2)
+      end if
+
+C     Remaining 2 in y=min.
+      if ( cctk_bbox(3).eq.1 .and. ny.ge.4
+     $     .and. cctk_bbox(5).eq.1 .and. nz.ge.4 ) then
+         var(:,1,1) = 2.d0 * var(:,2,2) - var(:,3,3)
+      end if
+
+      if ( cctk_bbox(3).eq.1 .and. ny.ge.4
+     $     .and. cctk_bbox(6).eq.1 .and. nz.ge.4 ) then 
+         var(:,1,nz) = 2.d0 * var(:,2,nz-1) - var(:,2,nz-2)
+      end if
+
+C     Remaining 2 in y=ymax.      
+      if ( cctk_bbox(4).eq.1 .and. ny.ge.4
+     $     .and. cctk_bbox(5).eq.1 .and. nz.ge.4 ) then
+         var(:,ny,1) = 2.d0 * var(:,ny-1,2) - var(:,ny-2,3)  
+      end if
+
+      if ( cctk_bbox(4).eq.1 .and. ny.ge.4
+     $     .and. cctk_bbox(6).eq.1 .and. nz.ge.4 ) then 
+         var(:,ny,nz) = 2.d0 * var(:,ny-1,nz-1) - var(:,ny-2,nz-2)
+      end if
+
+C     8 corners.
+
+      if (nx.ge.4 .and. ny.ge.4 .and. nz.ge.4) then
+         if (cctk_bbox(1).eq.1 .and. cctk_bbox(3).eq.1 .and. cctk_bbox(5).eq.1) then
+            var(1,1,1) = 2.d0*var(2,2,2) - var(3,3,3)
+         end if
+         if (cctk_bbox(1).eq.1 .and. cctk_bbox(3).eq.1 .and. cctk_bbox(6).eq.1) then
+            var(1,1,nz) = 2.d0*var(2,2,nz-1) - var(3,3,nz-2)
+         end if
+         if (cctk_bbox(1).eq.1 .and. cctk_bbox(4).eq.1 .and. cctk_bbox(5).eq.1) then
+            var(1,ny,1) = 2.d0*var(2,ny-1,2) - var(3,ny-2,3)
+         end if
+         if (cctk_bbox(1).eq.1 .and. cctk_bbox(4).eq.1 .and. cctk_bbox(6).eq.1) then
+            var(1,ny,nz) = 2.d0*var(2,ny-1,nz-1) - var(3,ny-2,nz-2)
+         end if
+        if (cctk_bbox(2).eq.1 .and. cctk_bbox(3).eq.1 .and. cctk_bbox(5).eq.1) then
+            var(nx,1,1) = 2.d0*var(nx-1,2,2) - var(nx-2,3,3)
+         end if
+         if (cctk_bbox(2).eq.1 .and. cctk_bbox(3).eq.1 .and. cctk_bbox(6).eq.1) then
+            var(nx,1,nz) = 2.d0*var(nx-1,2,nz-1) - var(nx-2,3,nz-2)
+         end if
+         if (cctk_bbox(2).eq.1 .and. cctk_bbox(4).eq.1 .and. cctk_bbox(5).eq.1) then
+            var(nx,ny,1) = 2.d0*var(nx-1,ny-1,2) - var(nx-2,ny-2,3)
+         end if
+         if (cctk_bbox(2).eq.1 .and. cctk_bbox(4).eq.1 .and. cctk_bbox(6).eq.1) then
+            var(nx,ny,nz) = 2.d0*var(nx-1,ny-1,nz-1) - var(nx-2,ny-2,nz-2)
+         end if
+      end if
+
+
+      return
+      end
diff --git a/src/make.code.defn b/src/make.code.defn
index f4946f0..136b470 100644
--- a/src/make.code.defn
+++ b/src/make.code.defn
@@ -2,41 +2,26 @@
 # $Header$
 
 # Source files in this directory
-SRCS = ParamCheck.c \
-       exactinitialize.F \
-       exactgauge.F \
-       exactboundary.F \
-       exactdata.F \
-       exactmetric.F \
-       boostrotmetric.F \
-       finkelstein.F \
-       kerrschild.F \
-       flatschwarz.F \
-       starschwarz.F \
-       bowl.F \
-       novikov.F \
-       minkowski.F \
-       fakebinary.F \
-       multiBH.F \
-       flatfunny.F \
-       flatshift.F \
-       bianchiI.F \
-       rob-wal.F \
-       godel.F \
-       kerrcart.F \
-       desitter.F \
-       kasner_l.F \
-       milne.F \
-       exactblendbound.F \
-       exactcartblendbound.F \
-       slice_initialize.F \
-       slice_evolve.F \
-       slice_data.F \
-       linextraponebound.F \
-       Startup.c \
-       GaugeWave.F
-#      ComparisonSolutions.F Comparison.c
-
-# Subdirectories containing source files
-SUBDIRS = include
+SRCS = ParamCheck.c			\
+       Startup.c			\
+       decode_pars.F			\
+       initialize.F			\
+					\
+       slice_initialize.F		\
+       slice_evolve.F			\
+       slice_data.F			\
+					\
+       gauge.F				\
+       Bona_Masso_data.F		\
+       metric.F				\
+					\
+       boundary.F			\
+       blended_boundary.F		\
+       xyz_blended_boundary.F		\
+       linear_extrap_one_bndry.F
 
+# Subdirectories containing source files to be compiled
+# n.b. ./include/ contains source fragments to be included in other
+#      thorns for calculating T_mu_nu, but these source fragments
+#      do *not* go in the SUBDIRS list here
+SUBDIRS = metrics
diff --git a/src/metric.F b/src/metric.F
new file mode 100644
index 0000000..912f399
--- /dev/null
+++ b/src/metric.F
@@ -0,0 +1,288 @@
+c This subroutine calculates the 4-metric and its inverse at an event,
+c by decoding  decoded_exact_model  and calling the appropriate subroutine
+C $Header$
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+#include "cctk_Arguments.h"
+
+#include "param_defs.inc"
+
+      subroutine Exact__metric(
+     $     decoded_exact_model,
+     $     x, y, z, t,
+     $     gdtt, gdtx, gdty, gdtz,
+     $     gdxx, gdyy, gdzz, gdxy, gdyz, gdxz,
+     $     gutt, gutx, guty, gutz,
+     $     guxx, guyy, guzz, guxy, guyz, guxz,rama)
+
+      implicit none
+      DECLARE_CCTK_FUNCTIONS
+
+c arguments
+      CCTK_INT decoded_exact_model
+      CCTK_REAL x, y, z, t,
+     $     gdtt, gdtx, gdty, gdtz,
+     $     gdxx, gdyy, gdzz, gdxy, gdyz, gdxz,
+     $     gutt, gutx, guty, gutz,
+     $     guxx, guyy, guzz, guxy, guyz, guxz, rama
+
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+c
+c Minkowski spacetime
+c
+
+      if     (decoded_exact_model .eq. EXACT__Minkowski) then
+         call Exact__Minkowski(
+     $        x, y, z, t,
+     $        gdtt, gdtx, gdty, gdtz,
+     $        gdxx, gdyy, gdzz, gdxy, gdyz, gdxz,
+     $        gutt, gutx, guty, gutz,
+     $        guxx, guyy, guzz, guxy, guyz, guxz)
+
+      elseif (decoded_exact_model .eq. EXACT__Minkowski_shift) then
+         call Exact__Minkowski_shift(
+     $        x, y, z, t,
+     $        gdtt, gdtx, gdty, gdtz,
+     $        gdxx, gdyy, gdzz, gdxy, gdyz, gdxz,
+     $        gutt, gutx, guty, gutz,
+     $        guxx, guyy, guzz, guxy, guyz, guxz)
+
+      elseif (decoded_exact_model .eq. EXACT__Minkowski_funny) then
+         call Exact__Minkowski_funny(
+     $        x, y, z, t,
+     $        gdtt, gdtx, gdty, gdtz,
+     $        gdxx, gdyy, gdzz, gdxy, gdyz, gdxz,
+     $        gutt, gutx, guty, gutz,
+     $        guxx, guyy, guzz, guxy, guyz, guxz)
+
+      elseif (decoded_exact_model .eq. EXACT__Minkowski_gauge_wave) then
+         call Exact__Minkowski_gauge_wave(
+     $        x, y, z, t,
+     $        gdtt, gdtx, gdty, gdtz,
+     $        gdxx, gdyy, gdzz, gdxy, gdyz, gdxz,
+     $        gutt, gutx, guty, gutz,
+     $        guxx, guyy, guzz, guxy, guyz, guxz)
+
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+c
+c black hole spacetimes
+c
+
+      elseif (decoded_exact_model .eq. EXACT__Schwarzschild_EF) then
+         call Exact__Schwarzschild_EF(
+     $        x, y, z, t,
+     $        gdtt, gdtx, gdty, gdtz,
+     $        gdxx, gdyy, gdzz, gdxy, gdyz, gdxz,
+     $        gutt, gutx, guty, gutz,
+     $        guxx, guyy, guzz, guxy, guyz, guxz)
+
+      elseif (decoded_exact_model .eq. EXACT__Schwarzschild_PG) then
+         call Exact__Schwarzschild_PG(
+     $        x, y, z, t,
+     $        gdtt, gdtx, gdty, gdtz,
+     $        gdxx, gdyy, gdzz, gdxy, gdyz, gdxz,
+     $        gutt, gutx, guty, gutz,
+     $        guxx, guyy, guzz, guxy, guyz, guxz)
+
+      elseif (decoded_exact_model .eq. EXACT__Schwarzschild_Novikov) then
+         call Exact__Schwarzschild_Novikov(x,y,z,t,
+     $        gdtt, gdtx, gdty, gdtz,
+     $        gdxx, gdyy, gdzz, gdxy, gdyz, gdxz,
+     $        gutt, gutx, guty, gutz,
+     $        guxx, guyy, guzz, guxy, guyz, guxz)
+
+      elseif (decoded_exact_model .eq. EXACT__Kerr_BoyerLindquist) then
+         call Exact__Kerr_BoyerLindquist(
+     $        x, y, z, t,
+     $        gdtt, gdtx, gdty, gdtz,
+     $        gdxx, gdyy, gdzz, gdxy, gdyz, gdxz,
+     $        gutt, gutx, guty, gutz,
+     $        guxx, guyy, guzz, guxy, guyz, guxz)
+
+      elseif (decoded_exact_model .eq. EXACT__Kerr_KerrSchild) then
+         call Exact__Kerr_KerrSchild(
+     $        x, y, z, t,
+     $        gdtt, gdtx, gdty, gdtz,
+     $        gdxx, gdyy, gdzz, gdxy, gdyz, gdxz,
+     $        gutt, gutx, guty, gutz,
+     $        guxx, guyy, guzz, guxy, guyz, guxz)
+
+      elseif (decoded_exact_model .eq. EXACT__Schwarzschild_Lemaitre) then
+         call Exact__Schwarzschild_Lemaitre(
+     $        x, y, z, t,
+     $        gdtt, gdtx, gdty, gdtz,
+     $        gdxx, gdyy, gdzz, gdxy, gdyz, gdxz,
+     $        gutt, gutx, guty, gutz,
+     $        guxx, guyy, guzz, guxy, guyz, guxz)
+
+      elseif (decoded_exact_model .eq. EXACT__multi_BH) then
+         call Exact__multi_BH(
+     $        x, y, z, t,
+     $        gdtt, gdtx, gdty, gdtz,
+     $        gdxx, gdyy, gdzz, gdxy, gdyz, gdxz,
+     $        gutt, gutx, guty, gutz,
+     $        guxx, guyy, guzz, guxy, guyz, guxz)
+
+c
+c not fully implemented yet -- see Nina Jansen for details
+c
+c     elseif (decoded_exact_model .eq. EXACT__Alvi) then
+c        call Exact__Alvi(
+c    $        x, y, z, t,
+c    $        gdtt, gdtx, gdty, gdtz,
+c    $        gdxx, gdyy, gdzz, gdxy, gdyz, gdxz,
+c    $        gutt, gutx, guty, gutz,
+c    $        guxx, guyy, guzz, guxy, guyz, guxz)
+c
+
+      elseif (decoded_exact_model .eq. EXACT__Thorne_fakebinary) then
+         call Exact__Thorne_fakebinary(
+     $        x, y, z, t,
+     $        gdtt, gdtx, gdty, gdtz,
+     $        gdxx, gdyy, gdzz, gdxy, gdyz, gdxz,
+     $        gutt, gutx, guty, gutz,
+     $        guxx, guyy, guzz, guxy, guyz, guxz)
+
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+c
+c cosmological spacetimes
+c
+
+      elseif (decoded_exact_model .eq. EXACT__Lemaitre) then
+         call Exact__Lemaitre(
+     $        x, y, z, t,
+     $        gdtt, gdtx, gdty, gdtz,
+     $        gdxx, gdyy, gdzz, gdxy, gdyz, gdxz,
+     $        gutt, gutx, guty, gutz,
+     $        guxx, guyy, guzz, guxy, guyz, guxz)
+
+      elseif (decoded_exact_model .eq. EXACT__Robertson_Walker) then
+         call Exact__Robertson_Walker(
+     $        x, y, z, t,
+     $        gdtt, gdtx, gdty, gdtz,
+     $        gdxx, gdyy, gdzz, gdxy, gdyz, gdxz,
+     $        gutt, gutx, guty, gutz,
+     $        guxx, guyy, guzz, guxy, guyz, guxz,rama)
+
+      elseif (decoded_exact_model .eq. EXACT__de_Sitter) then
+         call Exact__de_Sitter(
+     $        x, y, z, t,
+     $        gdtt, gdtx, gdty, gdtz,
+     $        gdxx, gdyy, gdzz, gdxy, gdyz, gdxz,
+     $        gutt, gutx, guty, gutz,
+     $        guxx, guyy, guzz, guxy, guyz, guxz)
+
+      elseif (decoded_exact_model .eq. EXACT__de_Sitter_Lambda) then
+         call Exact__de_Sitter_Lambda(
+     $        x, y, z, t,
+     $        gdtt, gdtx, gdty, gdtz,
+     $        gdxx, gdyy, gdzz, gdxy, gdyz, gdxz,
+     $        gutt, gutx, guty, gutz,
+     $        guxx, guyy, guzz, guxy, guyz, guxz)
+
+      elseif (decoded_exact_model .eq. EXACT__anti_de_Sitter_Lambda) then
+         call Exact__anti_de_Sitter_Lambda(
+     $        x, y, z, t,
+     $        gdtt, gdtx, gdty, gdtz,
+     $        gdxx, gdyy, gdzz, gdxy, gdyz, gdxz,
+     $        gutt, gutx, guty, gutz,
+     $        guxx, guyy, guzz, guxy, guyz, guxz)
+
+      elseif (decoded_exact_model .eq. EXACT__Bianchi_I) then
+         call Exact__Bianchi_I(
+     $        x, y, z, t,
+     $        gdtt, gdtx, gdty, gdtz,
+     $        gdxx, gdyy, gdzz, gdxy, gdyz, gdxz,
+     $        gutt, gutx, guty, gutz,
+     $        guxx, guyy, guzz, guxy, guyz, guxz)
+
+      elseif (decoded_exact_model .eq. EXACT__Goedel) then
+         call Exact__Goedel(
+     $        x, y, z, t,
+     $        gdtt, gdtx, gdty, gdtz,
+     $        gdxx, gdyy, gdzz, gdxy, gdyz, gdxz,
+     $        gutt, gutx, guty, gutz,
+     $        guxx, guyy, guzz, guxy, guyz, guxz)
+
+      elseif (decoded_exact_model .eq. EXACT__Bertotti) then
+         call Exact__Bertotti(
+     $        x, y, z, t,
+     $        gdtt, gdtx, gdty, gdtz,
+     $        gdxx, gdyy, gdzz, gdxy, gdyz, gdxz,
+     $        gutt, gutx, guty, gutz,
+     $        guxx, guyy, guzz, guxy, guyz, guxz)
+
+      elseif (decoded_exact_model .eq. EXACT__Kasner_like) then
+         call Exact__Kasner_like(
+     $        x, y, z, t,
+     $        gdtt, gdtx, gdty, gdtz,
+     $        gdxx, gdyy, gdzz, gdxy, gdyz, gdxz,
+     $        gutt, gutx, guty, gutz,
+     $        guxx, guyy, guzz, guxy, guyz, guxz)
+
+      elseif (decoded_exact_model .eq. EXACT__Kasner_axisymmetric) then
+         call Exact__Kasner_axisymmetric(
+     $        x, y, z, t,
+     $        gdtt, gdtx, gdty, gdtz,
+     $        gdxx, gdyy, gdzz, gdxy, gdyz, gdxz,
+     $        gutt, gutx, guty, gutz,
+     $        guxx, guyy, guzz, guxy, guyz, guxz)
+
+      elseif (decoded_exact_model .eq. EXACT__Kasner_generalized) then
+         call Exact__Kasner_generalized(
+     $        x, y, z, t,
+     $        gdtt, gdtx, gdty, gdtz,
+     $        gdxx, gdyy, gdzz, gdxy, gdyz, gdxz,
+     $        gutt, gutx, guty, gutz,
+     $        guxx, guyy, guzz, guxy, guyz, guxz)
+
+      elseif (decoded_exact_model .eq. EXACT__Milne) then
+         call Exact__Milne(
+     $        x, y, z, t,
+     $        gdtt, gdtx, gdty, gdtz,
+     $        gdxx, gdyy, gdzz, gdxy, gdyz, gdxz,
+     $        gutt, gutx, guty, gutz,
+     $        guxx, guyy, guzz, guxy, guyz, guxz)
+
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+c
+c miscellaneous spacetimes
+c
+
+      elseif (decoded_exact_model .eq. EXACT__boost_rotation_symmetric) then
+         call Exact__boost_rotation_symmetric(
+     $        x, y, z, t,
+     $        gdtt, gdtx, gdty, gdtz,
+     $        gdxx, gdyy, gdzz, gdxy, gdyz, gdxz,
+     $        gutt, gutx, guty, gutz,
+     $        guxx, guyy, guzz, guxy, guyz, guxz)
+
+      elseif (decoded_exact_model .eq. EXACT__bowl) then
+         call Exact__bowl(
+     $        x, y, z, t,
+     $        gdtt, gdtx, gdty, gdtz,
+     $        gdxx, gdyy, gdzz, gdxy, gdyz, gdxz,
+     $        gutt, gutx, guty, gutz,
+     $        guxx, guyy, guzz, guxy, guyz, guxz)
+
+      elseif (decoded_exact_model .eq. EXACT__constant_density_star) then
+         call Exact__constant_density_star(
+     $        x, y, z, t,
+     $        gdtt, gdtx, gdty, gdtz,
+     $        gdxx, gdyy, gdzz, gdxy, gdyz, gdxz,
+     $        gutt, gutx, guty, gutz,
+     $        guxx, guyy, guzz, guxy, guyz, guxz)
+
+cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
+
+      else
+         call CCTK_WARN(0,"Unknown value of decoded_exact_model")
+      endif
+
+      return
+      end
diff --git a/src/slice_data.F b/src/slice_data.F
index 3ca811d..08b49cf 100644
--- a/src/slice_data.F
+++ b/src/slice_data.F
@@ -1,10 +1,11 @@
 C     Extract Cauchy data from the slice x^A(x^i) stored in slicex,
 C     slicey, slicez, slicet, and calculate dx^A/dt.
+C $Header$
 
 #include "cctk.h"
 #include "cctk_Arguments.h"
 
-      subroutine slice_data(CCTK_ARGUMENTS)
+      subroutine Exact__slice_data(CCTK_ARGUMENTS)
 
       implicit none
 
@@ -55,7 +56,8 @@ C     Calculate first derivatives of slice coordinates.
 
 C     Now we need the exact solution metric in the preferred coordinates x^A.
 
-               call exactmetric(
+               call Exact__metric(
+     $              decoded_exact_model,
      $              slicex(i,j,k), slicey(i,j,k), slicez(i,j,k), 
      $              slicet(i,j,k),
      $              gd(4,4), gd(1,4), gd(2,4), gd(3,4),
@@ -162,7 +164,8 @@ C     with an editor (hint for proofreading).
 C     Calculate g_AB,C. Need to sum explicitly over C. Do this with 
 C     the editor.
 
-               call exactmetric(
+               call Exact__metric(
+     $              decoded_exact_model,
      $              slicex(i,j,k)+ex_eps, slicey(i,j,k), slicez(i,j,k), 
      $              slicet(i,j,k),
      $              gd_p(4,4), gd_p(1,4), gd_p(2,4), gd_p(3,4),
@@ -171,7 +174,8 @@ C     the editor.
      $              gu(4,4), gu(1,4), gu(2,4), gu(3,4),
      $              gu(1,1), gu(2,2), gu(3,3), 
      $              gu(1,2), gu(2,3), gu(1,3))
-               call exactmetric(
+               call Exact__metric(
+     $              decoded_exact_model,
      $              slicex(i,j,k)-ex_eps, slicey(i,j,k), slicez(i,j,k), 
      $              slicet(i,j,k),
      $              gd_m(4,4), gd_m(1,4), gd_m(2,4), gd_m(3,4),
@@ -186,7 +190,8 @@ C     the editor.
                   end do
                end do
 
-               call exactmetric(
+               call Exact__metric(
+     $              decoded_exact_model,
      $              slicex(i,j,k), slicey(i,j,k)+ex_eps, slicez(i,j,k), 
      $              slicet(i,j,k),
      $              gd_p(4,4), gd_p(1,4), gd_p(2,4), gd_p(3,4),
@@ -195,7 +200,8 @@ C     the editor.
      $              gu(4,4), gu(1,4), gu(2,4), gu(3,4),
      $              gu(1,1), gu(2,2), gu(3,3), 
      $              gu(1,2), gu(2,3), gu(1,3))
-               call exactmetric(
+               call Exact__metric(
+     $              decoded_exact_model,
      $              slicex(i,j,k), slicey(i,j,k)-ex_eps, slicez(i,j,k), 
      $              slicet(i,j,k),
      $              gd_m(4,4), gd_m(1,4), gd_m(2,4), gd_m(3,4),
@@ -210,7 +216,8 @@ C     the editor.
                   end do
                end do
 
-               call exactmetric(
+               call Exact__metric(
+     $              decoded_exact_model,
      $              slicex(i,j,k), slicey(i,j,k), slicez(i,j,k)+ex_eps, 
      $              slicet(i,j,k),
      $              gd_p(4,4), gd_p(1,4), gd_p(2,4), gd_p(3,4),
@@ -219,7 +226,8 @@ C     the editor.
      $              gu(4,4), gu(1,4), gu(2,4), gu(3,4),
      $              gu(1,1), gu(2,2), gu(3,3), 
      $              gu(1,2), gu(2,3), gu(1,3))
-               call exactmetric(
+               call Exact__metric(
+     $              decoded_exact_model,
      $              slicex(i,j,k), slicey(i,j,k), slicez(i,j,k)-ex_eps, 
      $              slicet(i,j,k),
      $              gd_m(4,4), gd_m(1,4), gd_m(2,4), gd_m(3,4),
@@ -234,7 +242,8 @@ C     the editor.
                   end do
                end do
 
-               call exactmetric(
+               call Exact__metric(
+     $              decoded_exact_model,
      $              slicex(i,j,k), slicey(i,j,k), slicez(i,j,k), 
      $              slicet(i,j,k)+ex_eps,
      $              gd_p(4,4), gd_p(1,4), gd_p(2,4), gd_p(3,4),
@@ -243,7 +252,8 @@ C     the editor.
      $              gu(4,4), gu(1,4), gu(2,4), gu(3,4),
      $              gu(1,1), gu(2,2), gu(3,3), 
      $              gu(1,2), gu(2,3), gu(1,3))
-               call exactmetric(
+               call Exact__metric(
+     $              decoded_exact_model,
      $              slicex(i,j,k), slicey(i,j,k), slicez(i,j,k), 
      $              slicet(i,j,k)-ex_eps,
      $              gd_m(4,4), gd_m(1,4), gd_m(2,4), gd_m(3,4),
@@ -302,30 +312,30 @@ C     Calculate K_ij.
 
 C     Synchronize and bound slicetmp2, which contains dx^A/dt.
 
-      call linextraponebound(CCTK_ARGUMENTS,slicetmp2x)
-      call linextraponebound(CCTK_ARGUMENTS,slicetmp2y)
-      call linextraponebound(CCTK_ARGUMENTS,slicetmp2z)
-      call linextraponebound(CCTK_ARGUMENTS,slicetmp2t)
+      call Exact__linear_extrap_one_bndry(CCTK_ARGUMENTS,slicetmp2x)
+      call Exact__linear_extrap_one_bndry(CCTK_ARGUMENTS,slicetmp2y)
+      call Exact__linear_extrap_one_bndry(CCTK_ARGUMENTS,slicetmp2z)
+      call Exact__linear_extrap_one_bndry(CCTK_ARGUMENTS,slicetmp2t)
  
       call CCTK_SyncGroup(ierr,cctkGH,"Exact::Exact_slicetemp2")
 
 C     Bound and synchronize the 3-metric and extrinsic curvature.
 
-      call linextraponebound(CCTK_ARGUMENTS,gxx)
-      call linextraponebound(CCTK_ARGUMENTS,gxy)
-      call linextraponebound(CCTK_ARGUMENTS,gxz)
-      call linextraponebound(CCTK_ARGUMENTS,gyy)
-      call linextraponebound(CCTK_ARGUMENTS,gyz)
-      call linextraponebound(CCTK_ARGUMENTS,gzz)
+      call Exact__linear_extrap_one_bndry(CCTK_ARGUMENTS,gxx)
+      call Exact__linear_extrap_one_bndry(CCTK_ARGUMENTS,gxy)
+      call Exact__linear_extrap_one_bndry(CCTK_ARGUMENTS,gxz)
+      call Exact__linear_extrap_one_bndry(CCTK_ARGUMENTS,gyy)
+      call Exact__linear_extrap_one_bndry(CCTK_ARGUMENTS,gyz)
+      call Exact__linear_extrap_one_bndry(CCTK_ARGUMENTS,gzz)
 
       call CCTK_SyncGroup(ierr,cctkGH,"admbase::metric")
 
-      call linextraponebound(CCTK_ARGUMENTS,kxx)
-      call linextraponebound(CCTK_ARGUMENTS,kxy)
-      call linextraponebound(CCTK_ARGUMENTS,kxz)
-      call linextraponebound(CCTK_ARGUMENTS,kyy)
-      call linextraponebound(CCTK_ARGUMENTS,kyz)
-      call linextraponebound(CCTK_ARGUMENTS,kzz)
+      call Exact__linear_extrap_one_bndry(CCTK_ARGUMENTS,kxx)
+      call Exact__linear_extrap_one_bndry(CCTK_ARGUMENTS,kxy)
+      call Exact__linear_extrap_one_bndry(CCTK_ARGUMENTS,kxz)
+      call Exact__linear_extrap_one_bndry(CCTK_ARGUMENTS,kyy)
+      call Exact__linear_extrap_one_bndry(CCTK_ARGUMENTS,kyz)
+      call Exact__linear_extrap_one_bndry(CCTK_ARGUMENTS,kzz)
 
       call CCTK_SyncGroup(ierr,cctkGH,"admbase::curv")
 
diff --git a/src/slice_evolve.F b/src/slice_evolve.F
index 8380593..9c95134 100644
--- a/src/slice_evolve.F
+++ b/src/slice_evolve.F
@@ -1,9 +1,10 @@
 C     Evolve the slice in the exact spacetime.
+C $Header$
 
 #include "cctk.h"
 #include "cctk_Arguments.h"
 
-      subroutine slice_evolve(CCTK_ARGUMENTS)
+      subroutine Exact__slice_evolve(CCTK_ARGUMENTS)
 
       implicit none
 
@@ -39,10 +40,10 @@ C     by slice_data.
 
 C     Synchronize and bound slice.
 
-      call linextraponebound(CCTK_ARGUMENTS,slicetmp1x)
-      call linextraponebound(CCTK_ARGUMENTS,slicetmp1y)
-      call linextraponebound(CCTK_ARGUMENTS,slicetmp1z)
-      call linextraponebound(CCTK_ARGUMENTS,slicetmp1t)
+      call Exact__linear_extrap_one_bndry(CCTK_ARGUMENTS,slicetmp1x)
+      call Exact__linear_extrap_one_bndry(CCTK_ARGUMENTS,slicetmp1y)
+      call Exact__linear_extrap_one_bndry(CCTK_ARGUMENTS,slicetmp1z)
+      call Exact__linear_extrap_one_bndry(CCTK_ARGUMENTS,slicetmp1t)
 
       call CCTK_SyncGroup(ierr,cctkGH,"Exact::Exact_slicetemp1")
 
@@ -73,7 +74,8 @@ C     Calculate first derivatives of slice coordinates.
 
 C     Now we need the exact solution metric in the preferred coordinates
 C     x^A. 
-               call exactmetric(
+               call Exact__metric(
+     $              decoded_exact_model,
      $              slicetmp1x(i,j,k), slicetmp1y(i,j,k), slicetmp1z(i,j,k), 
      $              slicetmp1t(i,j,k),
      $              gd(4,4), gd(1,4), gd(2,4), gd(3,4),
@@ -92,10 +94,10 @@ C     Calculate n^A and dx^A/dt
 
 C     Synchronize and bound slicetmp2, which contains dx^A/dt.
 
-      call linextraponebound(CCTK_ARGUMENTS,slicetmp2x)
-      call linextraponebound(CCTK_ARGUMENTS,slicetmp2y)
-      call linextraponebound(CCTK_ARGUMENTS,slicetmp2z)
-      call linextraponebound(CCTK_ARGUMENTS,slicetmp2t)
+      call Exact__linear_extrap_one_bndry(CCTK_ARGUMENTS,slicetmp2x)
+      call Exact__linear_extrap_one_bndry(CCTK_ARGUMENTS,slicetmp2y)
+      call Exact__linear_extrap_one_bndry(CCTK_ARGUMENTS,slicetmp2z)
+      call Exact__linear_extrap_one_bndry(CCTK_ARGUMENTS,slicetmp2t)
 
       call CCTK_SyncGroup(ierr,cctkGH,"Exact::Exact_slicetemp2")
 
@@ -108,20 +110,17 @@ C Leapfrog step.
 
 C     Synchronize and bound slice.
 
-      call linextraponebound(CCTK_ARGUMENTS,slicex)
-      call linextraponebound(CCTK_ARGUMENTS,slicey)
-      call linextraponebound(CCTK_ARGUMENTS,slicez)
-      call linextraponebound(CCTK_ARGUMENTS,slicet)
+      call Exact__linear_extrap_one_bndry(CCTK_ARGUMENTS,slicex)
+      call Exact__linear_extrap_one_bndry(CCTK_ARGUMENTS,slicey)
+      call Exact__linear_extrap_one_bndry(CCTK_ARGUMENTS,slicez)
+      call Exact__linear_extrap_one_bndry(CCTK_ARGUMENTS,slicet)
 
       call CCTK_SyncGroup(ierr,cctkGH,"Exact::Exact_slice")
 
 C     Extract Cauchy data at the new position, and store dxA/dt 
 C     for use in the next Lax step.
 
-      call slice_data(CCTK_ARGUMENTS)
+      call Exact__slice_data(CCTK_ARGUMENTS)
 
       return
       end
-
-
-
diff --git a/src/slice_initialize.F b/src/slice_initialize.F
index 4493bbc..ccd9f5d 100644
--- a/src/slice_initialize.F
+++ b/src/slice_initialize.F
@@ -1,9 +1,10 @@
+C $Header$
 
 #include "cctk.h"
 #include "cctk_Parameters.h"
 #include "cctk_Arguments.h"
 
-      subroutine slice_initialize(CCTK_ARGUMENTS)
+      subroutine Exact__slice_initialize(CCTK_ARGUMENTS)
 
       implicit none
 
@@ -22,7 +23,7 @@ C     exact solution, and the exact solution time is a Gaussian.
 
 C     Calculate Cauchy data and dx^A/dt.
 
-      call slice_data(CCTK_ARGUMENTS)
+      call Exact__slice_data(CCTK_ARGUMENTS)
 
 C     This tells the code that what we have set is the physical metric,
 C     not a conformally rescaled one.
diff --git a/src/xyz_blended_boundary.F b/src/xyz_blended_boundary.F
new file mode 100644
index 0000000..98a3735
--- /dev/null
+++ b/src/xyz_blended_boundary.F
@@ -0,0 +1,210 @@
+C $Header$
+
+#include "cctk.h"
+#include "cctk_Parameters.h"
+#include "cctk_Arguments.h"
+
+      subroutine Exact__xyz_blended_boundary(CCTK_ARGUMENTS)
+
+      implicit none
+
+      DECLARE_CCTK_ARGUMENTS
+      DECLARE_CCTK_PARAMETERS
+      DECLARE_CCTK_FUNCTIONS
+
+      logical doKij, doGij, doDs, doVs, doLapse, doShift
+
+      integer i,j,k
+      integer nx,ny,nz
+      integer ninterps
+
+      CCTK_REAL xblend, yblend, zblend
+      CCTK_REAL xmin, xmax, ymin, ymax, zmin, zmax
+      CCTK_REAL xfrac, yfrac, zfrac, onemxfrac, onemyfrac, onemzfrac
+      CCTK_REAL oonints
+      CCTK_REAL sfrac, onemsfrac
+      CCTK_REAL gxxe, gyye, gzze, gxye, gyze, gxze
+      CCTK_REAL kxxe, kyye, kzze, kxye, kyze, kxze
+      CCTK_REAL dxgxxe, dxgyye, dxgzze, dxgxye, dxgyze, dxgxze
+      CCTK_REAL dygxxe, dygyye, dygzze, dygxye, dygyze, dygxze
+      CCTK_REAL dzgxxe, dzgyye, dzgzze, dzgxye, dzgyze, dzgxze
+      CCTK_REAL alpe, axe, aye, aze
+      CCTK_REAL betaxe,betaye,betaze
+      CCTK_REAL bxxe,bxye,bxze,byxe,byye,byze,bzxe,bzye,bzze
+      CCTK_REAL det, uxx, uxy, uxz, uyy, uyz, uzz
+      CCTK_REAL vxe,vye,vze,sav
+
+      CCTK_REAL :: dx,dy,dz,time,ierr
+
+C     Grid parameters.
+
+      nx = cctk_lsh(1)
+      ny = cctk_lsh(2)
+      nz = cctk_lsh(3)
+
+      dx = cctk_delta_space(1)
+      dy = cctk_delta_space(2)
+      dz = cctk_delta_space(3)
+
+      time = cctk_time
+
+C     Other parameters.
+
+      doKij = (exblend_Ks.eq.1)
+      doGij = (exblend_gs.eq.1)
+
+      doLapse  = ((exblend_gauge.eq.1).and.
+     $            (CCTK_Equals("lapse_evolution_method","exact").ne.0))
+      doShift  = ((exblend_gauge.eq.1).and.
+     $            (CCTK_Equals(shift_evolution_method,"exact").ne.0))
+
+      if (exblend_width.lt.0) then 
+         xblend = - exblend_width*dx
+         yblend = - exblend_width*dy
+         zblend = - exblend_width*dz
+      else
+         xblend = exblend_width
+         yblend = exblend_width
+         zblend = exblend_width
+      endif
+
+      call CCTK_CoordRange(ierr,cctkGH,xmin,xmax,-1,"x","cart3d")
+      call CCTK_CoordRange(ierr,cctkGH,ymin,ymax,-1,"y","cart3d")
+      call CCTK_CoordRange(ierr,cctkGH,zmin,zmax,-1,"z","cart3d")
+
+      do k=1,nz
+         do j=1,ny
+            do i=1,nx
+
+c     We only do anything if in the blending region
+
+               if (x(i,j,k) .ge. xmax - xblend .or.
+     $             x(i,j,k) .le. xmin + xblend .or.
+     $             y(i,j,k) .ge. ymax - yblend .or.
+     $             y(i,j,k) .le. ymin + yblend .or.
+     $             z(i,j,k) .ge. zmax - zblend .or.
+     $             z(i,j,k) .le. zmin + zblend) then
+
+                  call Exact__Bona_Masso_data(
+     $                 decoded_exact_model,
+     $                 x(i,j,k), y(i,j,k), z(i,j,k), time,
+     $                 gxxe, gyye, gzze, gxye, gyze, gxze,
+     $                 kxxe, kyye, kzze, kxye, kyze, kxze,
+     $                 dxgxxe, dxgyye, dxgzze, dxgxye, dxgyze, dxgxze,
+     $                 dygxxe, dygyye, dygzze, dygxye, dygyze, dygxze,
+     $                 dzgxxe, dzgyye, dzgzze, dzgxye, dzgyze, dzgxze,
+     $                 alpe, axe, aye, aze, betaxe, betaye, betaze,
+     $                 bxxe, bxye, bxze, byxe, 
+     $                 byye, byze, bzxe, bzye, bzze)
+
+c     This sucks, but we want the exact vs so we can blend them also.
+
+                  det = -(gxze**2*gyye) 
+     &                 + 2.d0*gxye*gxze*gyze 
+     &                 - gxxe*gyze**2 
+     &                 - gxye**2*gzze 
+     &                 + gxxe*gyye*gzze
+  
+                  uxx=(-gyze**2 + gyye*gzze)/det
+                  uxy=(gxze*gyze - gxye*gzze)/det
+                  uyy=(-gxze**2 + gxxe*gzze)/det
+                  uxz=(-gxze*gyye + gxye*gyze)/det
+                  uyz=(gxye*gxze - gxxe*gyze)/det
+                  uzz=(-gxye**2 + gxxe*gyye)/det
+
+c     OK so 6 blending cases. If frac = 1 we get all exact
+
+                  ninterps = 0
+
+                  xfrac     = 0.0D0
+                  onemxfrac = 0.0D0
+                  yfrac     = 0.0D0
+                  onemyfrac = 0.0D0
+                  zfrac     = 0.0D0
+                  onemzfrac = 0.0D0
+
+                  if (x(i,j,k) .le. xmin + xblend) then 
+                     xfrac = 1.0D0 - (x(i,j,k)-xmin) / xblend
+                     onemxfrac = 1.0D0 - xfrac
+                     ninterps = ninterps + 1
+                  endif
+
+                  if (x(i,j,k) .ge. xmax - xblend) then 
+                     xfrac = 1.0D0 - (xmax - x(i,j,k)) / xblend
+                     onemxfrac = 1.0D0 - xfrac
+                     ninterps = ninterps + 1
+                  endif
+
+                  if (y(i,j,k) .le. ymin + yblend) then 
+                     yfrac = 1.0D0 - (y(i,j,k)-ymin) / yblend
+                     onemyfrac = 1.0D0 - yfrac
+                     ninterps = ninterps + 1
+                  endif
+
+                  if (y(i,j,k) .ge. ymax - yblend) then 
+                     yfrac = 1.0D0 - (ymax - y(i,j,k)) / yblend
+                     onemyfrac = 1.0D0 - yfrac
+                     ninterps = ninterps + 1
+                  endif
+
+                  if (z(i,j,k) .le. zmin + zblend) then 
+                     zfrac = 1.0D0 - (z(i,j,k)-zmin) / zblend
+                     onemzfrac = 1.0D0 - zfrac
+                     ninterps = ninterps + 1
+                  endif
+
+                  if (z(i,j,k) .ge. zmax - zblend) then 
+                     zfrac = 1.0D0 - (zmax - z(i,j,k)) / zblend
+                     onemzfrac = 1.0D0 - zfrac
+                     ninterps = ninterps + 1
+                  endif
+
+                  oonints = 1.0D0 / ninterps
+
+                  if (ninterps .eq. 0 .or. ninterps .gt. 3) then 
+                     print *,"NINTERPS error", ninterps
+                     STOP
+                  endif
+                     
+                  sfrac = (xfrac + yfrac + zfrac) * oonints
+                  onemsfrac = 1.0D0 - sfrac
+
+c     Once again some c-preprocessor tricks based on the whole fortran
+c     space thing...
+
+#define INTPOINT(f,v) f(i,j,k) = sfrac * v  + onemsfrac * f(i,j,k)
+#define intone(f) INTPOINT(f, f e)
+#define int_grp(p) \
+        intone(p xx) &&\
+        intone(p xy) &&\
+        intone(p xz) &&\
+        intone(p yy) &&\
+        intone(p yz) &&\
+        intone(p zz)
+
+                  if (doGij) then 
+                     int_grp(g)
+                  endif
+
+                  if (doKij) then
+                     int_grp(k)
+                  endif
+
+                  if (doLapse) then 
+                     intone(alp)
+                  endif
+
+                  if (doShift.and.(shift_state.ne.0)) then 
+                     intone(betax)
+                     intone(betay)
+                     intone(betaz)
+                  endif
+
+               endif            ! r > rinner
+
+            enddo
+         enddo
+      enddo
+
+      return
+      end