[[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

Removed files:
  exactblendbound.F		renamed to  blended_boundary.F
  exactboundary.F		renamed to  boundary.F
  exactcartblendbound.F		renamed to  xyz_blended_boundary.F
  exactdata.F			renamed to  Bona_Masso_data.F
  exactgauge.F			renamed to  gauge.F
  exactinitialize.F		renamed to  initialize.F
  exactmetric.F			renamed to  metric.F
  linextraponebound.F		renamed to  linear_extrap_one_bndry.F


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

8 files changed, 0 insertions, 1304 deletions

diff --git a/src/exactblendbound.F b/src/exactblendbound.F
deleted file mode 100644
index 17db5c1..0000000
--- a/src/exactblendbound.F
+++ /dev/null
@@ -1,192 +0,0 @@
-#include "cctk.h"
-#include "cctk_Parameters.h"
-#include "cctk_Arguments.h"
-
-      subroutine exactblendboundary(CCTK_ARGUMENTS)
-
-      implicit none
-
-      DECLARE_CCTK_ARGUMENTS
-      DECLARE_CCTK_PARAMETERS
-
-      logical doKij, doGij, doLapse, doShift
-
-      integer i,j,k
-      integer nx,ny,nz
-      integer CCTK_Equals
-
-      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 exactdata(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/exactboundary.F b/src/exactboundary.F
deleted file mode 100644
index ab0edd8..0000000
--- a/src/exactboundary.F
+++ /dev/null
@@ -1,123 +0,0 @@
-
-#include "cctk.h"
-#include "cctk_Parameters.h"
-#include "cctk_Arguments.h"
-
-      subroutine exactboundary(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 exactgauge.
-
-#define EXACTDATAPOINT \
-               call exactdata(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/exactcartblendbound.F b/src/exactcartblendbound.F
deleted file mode 100644
index 7fe41f2..0000000
--- a/src/exactcartblendbound.F
+++ /dev/null
@@ -1,206 +0,0 @@
-#include "cctk.h"
-#include "cctk_Parameters.h"
-#include "cctk_Arguments.h"
-
-      subroutine exactcartblendboundary(CCTK_ARGUMENTS)
-
-      implicit none
-
-      DECLARE_CCTK_ARGUMENTS
-      DECLARE_CCTK_PARAMETERS
-
-      logical doKij, doGij, doDs, doVs, doLapse, doShift
-
-      integer i,j,k
-      integer nx,ny,nz
-      integer ninterps
-      integer CCTK_Equals
-
-      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 exactdata(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 
diff --git a/src/exactdata.F b/src/exactdata.F
deleted file mode 100644
index ad8a289..0000000
--- a/src/exactdata.F
+++ /dev/null
@@ -1,202 +0,0 @@
-
-C     This routine calculates BM initial data from a subroutine
-C     exactmetric that calculates the spacetime metric and its 
-C     inverse.
-
-#include "cctk.h"
-#include "cctk_Parameters.h"
-#include "cctk_Arguments.h"
-
-      subroutine  exactdata(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_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 exactmetric(
-     $     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 exactmetric(
-     $     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 exactmetric(
-     $     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 exactmetric(
-     $     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 exactmetric(
-     $     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 exactmetric(
-     $     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 exactmetric(
-     $     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 exactmetric(
-     $     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 exactmetric(
-     $     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/exactgauge.F b/src/exactgauge.F
deleted file mode 100644
index e45d686..0000000
--- a/src/exactgauge.F
+++ /dev/null
@@ -1,145 +0,0 @@
-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. 
-
-#include "cctk.h"
-#include "cctk_Parameters.h"
-#include "cctk_Arguments.h"
-
-      subroutine exactgauge(CCTK_ARGUMENTS)
-
-      implicit none
-
-      DECLARE_CCTK_ARGUMENTS
-      DECLARE_CCTK_PARAMETERS
-
-      integer i,j,k
-      integer nx,ny,nz
-      integer CCTK_Equals
-
-      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 exactdata(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 exactdata(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 exactdata(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,'exactgauge has been called without doing anything')
-      end if
-
-      return
-      end
diff --git a/src/exactinitialize.F b/src/exactinitialize.F
deleted file mode 100644
index 4cc4be1..0000000
--- a/src/exactinitialize.F
+++ /dev/null
@@ -1,73 +0,0 @@
-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.
-      
-#include "cctk.h"
-#include "cctk_Parameters.h"
-#include "cctk_Arguments.h"
-
-      subroutine exactinitialize(CCTK_ARGUMENTS)
-
-      implicit none
-
-      DECLARE_CCTK_ARGUMENTS
-      DECLARE_CCTK_PARAMETERS
-
-      integer i,j,k
-      integer nx,ny,nz
-      integer CCTK_Equals
-
-      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 exactdata(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/exactmetric.F b/src/exactmetric.F
deleted file mode 100644
index 16bbee2..0000000
--- a/src/exactmetric.F
+++ /dev/null
@@ -1,223 +0,0 @@
-#include "cctk.h"
-#include "cctk_Parameters.h"
-#include "cctk_Arguments.h"
-
-C     This is just a wrapper for different routines in place of exactmetric.
-
-      subroutine exactmetric(
-     $     x, y, z, t,
-     $     gdtt, gdtx, gdty, gdtz, 
-     $     gdxx, gdyy, gdzz, gdxy, gdyz, gdxz,
-     $     gutt, gutx, guty, gutz, 
-     $     guxx, guyy, guzz, guxy, guyz, guxz)
-
-      implicit none
-
-      DECLARE_CCTK_PARAMETERS
-
-      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
-
-      logical firstcall
-
-      integer CCTK_Equals
-      
-C     Call the corresponding routine.
-
-      if (CCTK_Equals(exactmodel,"boostrot").eq.1) then
-
-         call boostrotmetric(
-     $        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 (CCTK_Equals(exactmodel,"finkelstein").eq.1) then
-
-         call finkelstein(
-     $        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 (CCTK_Equals(exactmodel,"kerrschild").eq.1) then
-
-         call 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 (CCTK_Equals(exactmodel,"flatschwarz").eq.1) then
-
-         call flatschwarz(
-     $        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 (CCTK_Equals(exactmodel,"starschwarz").eq.1) then
-
-         call starschwarz(
-     $        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 (CCTK_Equals(exactmodel,"novikov").eq.1) then
-
-         call 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 (CCTK_Equals(exactmodel,"bowl").eq.1) then
-
-         call 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 (CCTK_Equals(exactmodel,"minkowski").eq.1) then
-
-         call 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 (CCTK_Equals(exactmodel,"fakebinary").eq.1) then
-
-         call 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)
-
-      elseif (CCTK_Equals(exactmodel,"multiBH").eq.1) then
-
-         call multi_BHs(
-     $        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 (CCTK_Equals(exactmodel,"flatfunny").eq.1) then
-
-         call flatfunny(
-     $        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 (CCTK_Equals(exactmodel,"flatshift").eq.1) then
-
-         call flatshift(
-     $        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 (CCTK_Equals(exactmodel,"bianchiI").eq.1) then
-
-         call bianchiI(
-     $        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 (CCTK_Equals(exactmodel,"rob-wal").eq.1) then
-
-         call rob_wal(
-     $        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 (CCTK_Equals(exactmodel,"godel").eq.1) then
-
-         call godel(
-     $        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 (CCTK_Equals(exactmodel,"kerr").eq.1) then
-
-         call kerr(
-     $        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 (CCTK_Equals(exactmodel,"desitter").eq.1) then
-
-         call desitter(
-     $        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 (CCTK_Equals(exactmodel,"kasner").eq.1) then
-
-         call kasner(
-     $        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 (CCTK_Equals(exactmodel,"milne").eq.1) then
-
-         call 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)
-
-
-      elseif (CCTK_Equals(exactmodel,"gauge_wave").eq.1) then
-
-         call 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)
-
-
-      else
-
-
-         call CCTK_WARN(0,"Unknown value of parameter exactmodel")
-
-      end if
-
-      return
-      end
diff --git a/src/linextraponebound.F b/src/linextraponebound.F
deleted file mode 100644
index 271aa0c..0000000
--- a/src/linextraponebound.F
+++ /dev/null
@@ -1,140 +0,0 @@
-#include "cctk.h"
-#include "cctk_Arguments.h"
-
-      subroutine linextraponebound(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