[[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:
  GaugeWave.F		moved to  metrics/Minkowski_gauge_wave.F
  bianchiI.F		moved to  metrics/Bianchi_I.F
  boostrotmetric.F	moved to  metrics/boost_rotation_symmetric.F
  bowl.F		moved to  metrics/bowl.F
  desitter.F		moved to  metrics/de_Sitter.F
  fakebinary.F		moved to  metrics/Thorne_fakebinary.F
  finkelstein.F		moved to  metrics/Schwarzschild_EF.F
  flatfunny.F		moved to  metrics/Minkowski_funny.F
  flatschwarz.F		moved to  metrics/Schwarzschild_PG.F
  flatshift.F		moved to  metrics/Minkowski_shift.F
  godel.F		moved to  metrics/Goedel.F
  kasner_l.F		moved to  metrics/Kasner_like.F
  kerrcart.F		moved to  metrics/Kerr_BoyerLindquist.F
  kerrschild.F		moved to  metrics/Kerr_KerrSchild.F
  milne.F		moved to  metrics/Milne.F
  minkowski.F		moved to  metrics/Minkowski.F
  multiBH.F		moved to  metrics/multi_BH.F
  novikov.F		moved to  metrics/Schwarzschild_Novikov.F
  rob-wal.F		moved to  metrics/Robertson_Walker.F
  starschwarz.F		moved to  metrics/constant_density_star.F


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

20 files changed, 0 insertions, 2157 deletions

diff --git a/src/GaugeWave.F b/src/GaugeWave.F
deleted file mode 100644
index 6ba0fef..0000000
--- a/src/GaugeWave.F
+++ /dev/null
@@ -1,124 +0,0 @@
-#include "cctk.h"
-#include "cctk_Arguments.h"
-#include "cctk_Parameters.h"
-
-#define Pi (4 * atan(1.d0))
-
-
-      subroutine  gauge_wave(
-     $     x, y, z, t,
-     $     gdtt, gdtx, gdty, gdtz, 
-     $     gdxx, gdyy, gdzz, gdxy, gdyz, gdzx,
-     $     gutt, gutx, guty, gutz, 
-     $     guxx, guyy, guzz, guxy, guyz, guzx)
-
-      implicit none
-      integer CCTK_Equals
-
-      DECLARE_CCTK_PARAMETERS
-
-      CCTK_REAL x, y, z, t
-      CCTK_REAL gdtt, gdtx, gdty, gdtz, 
-     $       gdxx, gdyy, gdzz, gdxy, gdyz, gdzx,
-     $       gutt, gutx, guty, gutz, 
-     $       guxx, guyy, guzz, guxy, guyz, guzx
-
-      logical firstcall
-
-      CCTK_REAL a, o, H, d, fs
-      CCTK_REAL zero,half,one
-
-      data firstcall /.true./
-      save firstcall, a, o, half, d
-
-C     Numbers.
-
-      zero = 0.0d0
-      half = 0.5d0
-      one  = 1.0d0
-
-C     Get parameters of the exact solution.
-
-      if (firstcall) then
-
-         a = GW_a
-         o = GW_omega
-         d = GW_del
-         fs = GW_phase_shift
-
-         firstcall = .false.
-
-      end if
-
-C     How should the wave look like.
-
-      if (CCTK_Equals(GW_H,"sin").eq.1) then
-
-        d = GW_del * half / Pi
-
-        if (GW_diagonal.ne.0) then
-
-          H = one - a * sin((x-y)/d - o*t/d)
- 
-        else
-
-          H = one - a * sin((x-o*t)/d - fs)
-
-        end if
-
-      elseif (CCTK_Equals(GW_H,"expsin").eq.1) then
-
-        H = exp(a*sin(x/d)*cos(t/d))
-
-      elseif (CCTK_Equals(GW_H,"gau").eq.1) then
-
-        H = one - a*dexp(-(x-t)**2/d**2)
-
-      end if
-
-C     write metric.
-
-      if (GW_diagonal .eq. 1) then
-
-        gdxx = half * H + half
-        gdxy = -half* H + half
-        gdyy = half * H + half
-
-        guxx = half / H + half
-        guxy =-half / H + half
-        guyy = half / H + half
-
-      else
-
-        gdxx = H
-        gdxy = zero
-        gdyy = one
-
-        guxx = one/H
-        guxy = zero
-        guyy = one
-
-      end if
-      
-      gdtt = - H
-      gdtx = zero
-      gdty = zero
-      gdtz = zero
-
-      gdzx = zero
-      gdyz = zero
-      gdzz = one
-
-C     and upper metric.
-
-      gutt = - one/H
-      gutx = zero
-      guty = zero
-      gutz = zero
-
-      guyz = zero
-      guzx = zero
-      guzz = one
-
-      return
-      end
diff --git a/src/bianchiI.F b/src/bianchiI.F
deleted file mode 100644
index 6d51575..0000000
--- a/src/bianchiI.F
+++ /dev/null
@@ -1,65 +0,0 @@
-C Bianchi-I spacetime !!!! Fake one...
-C It is not a real Bianchi I spacetime, it just emulates
-C one, taking the two BianchiI functions as harmonis ones...
-C Author : D. Vulcanov (Timisoara, Romania)
-
-#include "cctk.h"
-#include "cctk_Parameters.h"
-
-      subroutine  bianchiI(
-     $     x, y, z, t,
-     $     gdtt, gdtx, gdty, gdtz, 
-     $     gdxx, gdyy, gdzz, gdxy, gdyz, gdzx,
-     $     gutt, gutx, guty, gutz, 
-     $     guxx, guyy, guzz, guxy, guyz, guzx)
-
-      implicit none
-
-      DECLARE_CCTK_PARAMETERS
-
-      CCTK_REAL x, y, z, t
-      CCTK_REAL gdtt, gdtx, gdty, gdtz, 
-     $       gdxx, gdyy, gdzz, gdxy, gdyz, gdzx,
-     $       gutt, gutx, guty, gutz, 
-     $       guxx, guyy, guzz, guxy, guyz, guzx
-
-      logical firstcall
-    
-      CCTK_REAL bx, by, arad
-
-      data firstcall /.true./
-      save firstcall
-      if (firstcall) then
-          arad = bia
-          PRINT *, firstcall, arad, bia
-          firstcall = .false.
-      end if
-
-
-      bx = arad*sin(x+t)
-      by = arad*cos(x+t)
-     
-      gdtt = - 1.d0
-      gdtx = 0.d0
-      gdty = 0.d0
-      gdtz = 0.d0
-      gdxx = bx**2 
-      gdyy = by**2 
-      gdzz = by**2
-      gdxy = 0.d0
-      gdyz = 0.d0
-      gdzx = 0.d0
-
-      gutt = - 1.d0
-      gutx = 0.d0
-      guty = 0.d0
-      gutz = 0.d0
-      guxx = bx**(-2)
-      guyy = by**(-2)
-      guzz = by**(-2)
-      guxy = 0.d0
-      guyz = 0.d0 
-      guzx = 0.d0
-
-      return
-      end
diff --git a/src/boostrotmetric.F b/src/boostrotmetric.F
deleted file mode 100644
index bd3bab3..0000000
--- a/src/boostrotmetric.F
+++ /dev/null
@@ -1,165 +0,0 @@
-
-#include "cctk.h"
-#include "cctk_Parameters.h"
-
-      subroutine  boostrotmetric(
-     $     x, y, z, t,
-     $     gdtt, gdtx, gdty, gdtz, 
-     $     gdxx, gdyy, gdzz, gdxy, gdyz, gdzx,
-     $     gutt, gutx, guty, gutz, 
-     $     guxx, guyy, guzz, guxy, guyz, guzx)
-
-      implicit none
-
-      DECLARE_CCTK_PARAMETERS
-
-      CCTK_REAL x, y, z, t,
-     $    gdtt, gdtx, gdty, gdtz, 
-     $    gdxx, gdyy, gdzz, gdxy, gdyz, gdzx,
-     $    gutt, gutx, guty, gutz, 
-     $    guxx, guyy, guzz, guxy, guyz, guzx
-      CCTK_REAL a, b, mu0, mu1, lam1, mu2, lam2,
-     $    lam3, mu4, lam4, mu5, lam5, num, div, f,
-     $    elam, emu0, delta, gfunc, tmp
-      CCTK_REAL h, d, numlim
-
-      external gfunc
-
-      logical firstcall
-
-      data firstcall /.true./
-      save firstcall, h, d, numlim
-
-C     Get parameters of the exact solution.
-
-      if (firstcall) then
-
-         h = boostrotscale
-         d = boostrotstrength
-
-         numlim = boostrotsafedistance
-
-         firstcall = .false.
-
-      end if
-
-C     Intermediate quantities.
-
-      a = x**2 + y**2
-      b = z**2 - t**2
-
-      num = (0.5d0 * (a + b) - h)**2 + 2.d0 * h * a
-
-C     Make sure we are not sitting on one of the two source wordlines,
-C     given by x = y = 0, z = +/- sqrt(h^2 + t^2)
-
-      if (num / h**4 .le. numlim) stop 'too close to source wordline'
-
-      div = 1.d0 / sqrt(num**3)
-      f = d**2  * ((0.25d0 * (a + b)**2 - h**2)**2 
-     $     - 0.5d0 * h**2 * a * b) / num**4
-
-      mu0 = - d * div * (0.5d0 * a**2 + h * a)
-      mu1 = - d * div * (0.5d0 * b + a - h)
-      lam1 = d * div * (0.5d0 * b - h) - a * f
-      mu2 = gfunc(b, mu1)
-      lam2 =  gfunc(b, lam1)
-
-      lam3 = d * div * (0.5d0 * b**2 - h * b)
-      lam4 = - d * div * (0.5d0 * a + h) - b * f
-      mu4 = - d * div * (0.5d0 * a + b + h) 
-      mu5 = gfunc(a, - mu4)
-      lam5 = gfunc(a, lam4)
-
-      elam = exp(lam3 + a * lam4)
-      emu0 = exp(mu0)
-      delta = exp(lam3) * (mu5 - lam5)
-
-C     All nonvanishing metric coefficients (downstairs).
-
-      gdxx = elam + y**2 * Delta
-      gdyy = elam + x**2 * Delta
-      gdxy = - x * y * Delta
-      gdzz = emu0 * (1.d0 + lam2 * z**2 - mu2 * t**2)
-      gdtz = - emu0 * z * t * (lam2 - mu2)
-      gdtt = - emu0 * (1.d0 + mu2 * z**2 - lam2 * t**2)
-
-C     Others.
-
-      gdzx = 0.d0
-      gdyz = 0.d0
-      gdtx = 0.d0
-      gdty = 0.d0
-
-C     It is clear that the 3-metric is always spacelike in the xy plane. So
-C     we only need to check that gdzz is positive.
-
-      if (gdzz .le. 0.d0) then
-         write(*,*) 'WARNING 3-metric not spacelike in boostrot at'
-         write(*,*) 't =', t, 'z =', z
-         write(*,*) 'x =', x, 'y =', y
-         pause
-      end if
-
-C     Calculate inverse metric. That is not too difficult as it is
-c     in block-diagonal form.
-
-      tmp = gdtt * gdzz - gdtz**2
-
-      if (tmp .eq. 0.d0) then
-         write(*,*) 'boostrot metric inverse failed in tz plane'
-         STOP
-      end if
-
-      gutt = gdzz / tmp
-      guzz = gdtt / tmp
-      gutz = - gdtz / tmp
-      
-      tmp = gdxx * gdyy - gdxy**2
-
-      if (tmp .eq. 0.d0) then
-         write(*,*) 'boostrot metric inverse failed in xy plane'
-         STOP
-      end if
-
-      guxx = gdyy / tmp
-      guyy = gdxx / tmp
-      guxy = - gdxy / tmp
-
-      guzx = 0.d0
-      guyz = 0.d0
-      gutx = 0.d0
-      guty = 0.d0
-
-      return
-      end
-
-ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
-
-C     Calculates g = [exp (x f) - 1] / x as a power series for small x, 
-C     so that the expression is regular at x = 0.
-
-      function gfunc(x, f)
-    
-      implicit none
-
-      integer n
-
-      CCTK_REAL x, f, gfunc
-      CCTK_REAL sum, tmp
-
-      if (abs(x*f) .ge. 1.d-6) then
-         gfunc = (exp(x*f) - 1.d0) / x
-      else
-         sum = 0.d0
-         tmp = f
-         do n=1,10
-            tmp = tmp / dble(n) 
-            sum = sum + tmp
-            tmp = tmp * x * f
-         end do
-         gfunc = sum
-      end if
-
-      return
-      end
diff --git a/src/bowl.F b/src/bowl.F
deleted file mode 100644
index 124269b..0000000
--- a/src/bowl.F
+++ /dev/null
@@ -1,253 +0,0 @@
-#include "cctk.h"
-#include "cctk_Parameters.h"
-
-      subroutine bowl(
-     $     x, y, z, t,
-     $     gdtt, gdtx, gdty, gdtz, 
-     $     gdxx, gdyy, gdzz, gdxy, gdyz, gdzx,
-     $     gutt, gutx, guty, gutz, 
-     $     guxx, guyy, guzz, guxy, guyz, guzx)
-
-c     The metric given here is not a solution of
-c     Einsteins equations!  It is nevertheless
-c     useful for tests since it has a particularly
-c     nice geometry.  It is a static, spherically
-c     symmetric metric with no shift.
-c
-c     In spherical coordinates, the metric has the
-c     form:
-c
-c       2      2        2        2
-c     ds  =  dr  +  R(r)  d Omega
-c
-c     Clearly, r measures radial proper distance, and R(r)
-c     is the areal (Schwarzschild) radius.
-c
-c     I choose a form of R(r) such that:
-c
-c     R  -->   r     r<<1, r>>1
-c
-c     So close in, and far away we have a flat metric.
-c     In the middle region, I take R to be smaller than
-c     r, but still larger than zero.  This deficit in
-c     areal radius produces the geometry of a "bag of gold".
-c
-c     The size of the deviation from a flat geometry
-c     is controled by the parameter "bowl_a".
-c     If this parameter is 0, we are in flat space.
-c     The width of the curved region is controled by
-c     the paramter "bowl_s", and the place where the
-c     curvature becomes significant (the center of the
-c     deformation) is controled by "bowl_c".
-c
-c     The specific form of the function R(r) is:
-c
-c     R(r) = ( r - a f(r) )
-c
-c     and the form of thr function f(r) depends on the
-c     parameter bowl_type:
-c                                            2   2
-c     bowl_type = "gauss"   f(r) = exp(-(r-c) / s )
-c
-c     bowl_type = "fermi"   f(r) = 1 / ( 1 + exp(-s(r-c)) )
-c
-c     There are three extra paramters (bowl_dx,bowl_dy,bowl_dz)
-c     that scale the (x,y,z) axis respectively.  Their default
-c     value is 1.  These parameters are useful to hide the
-c     spherical symmetry of the metric.
-
-      implicit none
-
-      DECLARE_CCTK_PARAMETERS
-
-c     Input.
-      CCTK_REAL x, y, z, t
-
-c     Output.
-      CCTK_REAL gdtt, gdtx, gdty, gdtz, 
-     $          gdxx, gdyy, gdzz, gdxy, gdyz, gdzx,
-     $          gutt, gutx, guty, gutz, 
-     $          guxx, guyy, guzz, guxy, guyz, guzx
-
-      logical firstcall,evolve
-
-      integer type
-      integer CCTK_Equals
-
-      CCTK_REAL a,c,s
-      CCTK_REAL dx,dy,dz
-      CCTK_REAL r,r2,rr2
-      CCTK_REAL xr,yr,zr,xr2,yr2,zr2
-      CCTK_REAL fac,det
-      CCTK_REAL tfac,st,t0
-      CCTK_REAL zero,one,two
-
-      data firstcall /.true./
-      save firstcall,evolve,type,a,c,s,dx,dy,dz,t0,st
-
-c     Find {zero,one,two}.
-
-      zero = 0.0d0
-      one  = 1.0d0
-      two  = 2.0d0
-
-c     Get parameters of the metric.
-
-      if (firstcall) then
-
-         a = bowl_a
-         c = bowl_c
-         s = bowl_s
-
-         dx = bowl_dx
-         dy = bowl_dy
-         dz = bowl_dz
-
-         print *, a,c,s,dx,dy,dz
-
-         if (dx.le.zero) then
-            call CCTK_WARN(0,"The parameter bowl_dx must be greater than 0")
-         end if
-
-         if (dy.le.zero) then
-            call CCTK_WARN(0,"The parameter bowl_dy must be greater than 0")
-         end if
-
-         if (dz.le.zero) then
-            call CCTK_WARN(0,"The parameter bowl_dz must be greater than 0")
-         end if
-
-         if (CCTK_Equals(bowl_type,"gauss").ne.0) then
-            type = 1
-         else if (CCTK_Equals(bowl_type,"fermi").ne.0) then
-            type = 2
-         else
-            call CCTK_WARN(0,"Unknown type of bowl metric")
-         end if
-
-         if (bowl_evolve.eq.1) then
-            evolve = .true.
-            t0 = bowl_t0
-            st = bowl_st
-         else
-            evolve = .false.
-         end if
-
-         firstcall = .false.
-
-      end if
-
-c     Multiply the bowl strength "a" with a time evolution factor.
-c     The time evolution factor is taken to be a Fermi step centered
-c     in `t0' and with a width `st'.  The size of this step is always
-c     1 so that far in the past we will always have flat space, and
-c     far in the future we will have the static bowl.
-
-      if (evolve) then
-         tfac = one/(one + dexp(-st*(t-t0)))
-      else
-         tfac = one
-      end if
-
-      a = a*tfac
-
-c     Find {r2,r}.
-
-      r2 = (x/dx)**2 + (y/dy)**2 + (z/dz)**2
-      r  = dsqrt(r2)
-
-c     Find the form function rr2
-c
-c                      2    2                 2
-c     rr2  =  (r - a f)  / r  =  (1 - a f / r)
-
-      if (type.eq.1) then
-
-c        Gaussian bowl:
-c                                2  2       2
-c        rr2 = [ 1 - a exp(-(r-c) /s ) / r ]
-c
-c        Notice that this really does not go to 1 at the
-c        origin.  To fix this, I multiply the gaussian
-c        with the factor:
-c
-c        fac  =  1  -  sech(4r)
-c
-c        This goes smoothly to 0 at the origin, and climbs
-c        fast to a limiting value of 1 (at r=1 it is already
-c        equal to 0.96).
-
-         fac = one - two/(dexp(4.0d0*r) + dexp(-4.0d0*r))
-         rr2 = (one - a*fac*dexp(-((r-c)/s)**2)/r)**2
-
-      else if (type.eq.2) then
-
-c        Fermi bowl:
-c                                                  2
-c        rr2 = [ 1 - 1 / ( 1 + exp(-s(r-c)) ) / r ]
-c
-c        Again, this doesnt really go to 1 at the origin, so
-c        I use the same trick as above.
-
-         fac = one - two/(dexp(4.0d0*r) + dexp(-4.0d0*r))
-         rr2 = (one - a*fac/(one + dexp(-s*(r-c)))/r)**2
-
-      end if
-
-c     Give metric components.
-
-      gdtt = - one
-      gdtx = zero
-      gdty = zero
-      gdtz = zero
-
-      if (r.ne.0) then
-
-         xr = (x/dx)/r
-         yr = (y/dy)/r
-         zr = (z/dz)/r
-
-         xr2 = xr**2
-         yr2 = yr**2
-         zr2 = zr**2
-
-         gdxx = (xr2 + rr2*(yr2 + zr2))/dx**2
-         gdyy = (yr2 + rr2*(xr2 + zr2))/dy**2
-         gdzz = (zr2 + rr2*(xr2 + yr2))/dz**2
-
-         gdxy = xr*yr*(one - rr2)/(dx*dy)
-         gdyz = yr*zr*(one - rr2)/(dy*dz)
-         gdzx = xr*zr*(one - rr2)/(dx*dz)
-
-      else
-
-         gdxx = one
-         gdyy = one
-         gdzz = one
-
-         gdxy = zero
-         gdyz = zero
-         gdzx = zero
-
-      end if
-
-c     Find inverse metric.
-
-      gutt = - one
-      gutx = zero
-      guty = zero
-      gutz = zero
-
-      det = gdxx*gdyy*gdzz + two*gdxy*gdzx*gdyz
-     .    - gdxx*gdyz**2 - gdyy*gdzx**2 - gdzz*gdxy**2
-
-      guxx = (gdyy*gdzz - gdyz**2)/det
-      guyy = (gdxx*gdzz - gdzx**2)/det
-      guzz = (gdxx*gdyy - gdxy**2)/det
-
-      guxy = (gdzx*gdyz - gdzz*gdxy)/det
-      guyz = (gdxy*gdzx - gdxx*gdyz)/det
-      guzx = (gdxy*gdyz - gdyy*gdzx)/det
-
-      return
-      end
diff --git a/src/desitter.F b/src/desitter.F
deleted file mode 100644
index 4af0bfb..0000000
--- a/src/desitter.F
+++ /dev/null
@@ -1,73 +0,0 @@
-C Einstein-DeSitter metric spacetime !!!!
-C It emulates the Robertson-Walker universe
-C near t=0, with zero pressure, and k=0
-C See :J.N. Islam, An Introduction to 
-C Mathematical Cosmology, Cambridge, 1992 and 
-C S. Hawking, G.F.R. Ellis, The Large Scale
-C Structure of space-time, Cambridge, 1973
-C Author : D. Vulcanov (Timisoara, Romania)
-
-#include "cctk.h"
-#include "cctk_Parameters.h"
-
-      subroutine  desitter(
-     $     x, y, z, t,
-     $     gdtt, gdtx, gdty, gdtz, 
-     $     gdxx, gdyy, gdzz, gdxy, gdyz, gdzx,
-     $     gutt, gutx, guty, gutz, 
-     $     guxx, guyy, guzz, guxy, guyz, guzx)
-
-      implicit none
-
-      DECLARE_CCTK_PARAMETERS
-
-      CCTK_REAL x, y, z, t
-      CCTK_REAL gdtt, gdtx, gdty, gdtz, 
-     $       gdxx, gdyy, gdzz, gdxy, gdyz, gdzx,
-     $       gutt, gutx, guty, gutz, 
-     $       guxx, guyy, guzz, guxy, guyz, guzx
-
-      logical firstcall
-    
-      CCTK_REAL arad,r2,x2,y2,z2,am
-
-      data firstcall /.true./
-      save firstcall
-      if (firstcall) then
-          arad = desitt_a
-          PRINT *, firstcall, arad, desitt_a
-          firstcall = .false.
-      end if
- 
-      x2=x*x
-      y2=y*y
-      z2=z*z
-      r2 = x2+y2+z2
-      am=arad*t**(4/3)     
-
-      gdtt = -1.d0 
-      gdtx = 0.d0 
-      gdty = 0.d0
-      gdtz = 0.d0
-      gdxx = 1.d0+(-1.d0 +am)*x2/r2 
-      gdyy = 1.d0+(-1.d0 +am)*y2/r2
-      gdzz = 1.d0+(-1.d0 +am)*z2/r2
-      gdxy = (-1.d0   +am)*x*y/r2
-      gdyz = (-1.d0   +am)*y*z/r2
-      gdzx = (-1.d0   +am)*x*z/r2
-
-      gutt = -1.d0
-      gutx = 0.d0
-      guty = 0.d0
-      gutz = 0.d0
-      guxx = 1.d0+(-1.d0 +1.d0/am)*x2/r2 
-      guyy = 1.d0+(-1.d0 +1.d0/am)*y2/r2
-      guzz = 1.d0+(-1.d0 +1.d0/am)*z2/r2
-      guxy = (-1.d0  +1.d0/am)*x*y/r2
-      guyz = (-1.d0  +1.d0/am)*y*z/r2
-      guzx = (-1.d0  +1.d0/am)*x*z/r2
- 
-
-
-      return
-      end
diff --git a/src/fakebinary.F b/src/fakebinary.F
deleted file mode 100644
index 0e496c2..0000000
--- a/src/fakebinary.F
+++ /dev/null
@@ -1,173 +0,0 @@
-C     fakebinary.F
-C     Bernd Bruegmann, 6/98
-C
-C     Compute Thorne s four metric that, although not a solution to the 
-C     Einstein equations, has several characteristic features of a binary
-C     star system. See gr-qc/98......
-
-#include "cctk.h"
-#include "cctk_Parameters.h"
-
-      subroutine 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)
-
-      implicit none
-
-      DECLARE_CCTK_PARAMETERS
-
-C     input
-      CCTK_REAL x, y, z, t
-
-C     output
-      CCTK_REAL gdtt, gdtx, gdty, gdtz, 
-     $       gdxx, gdyy, gdzz, gdxy, gdyz, gdxz,
-     $       gutt, gutx, guty, gutz, 
-     $       guxx, guyy, guzz, guxy, guyz, guxz
-
-C     static
-
-      logical firstcall
-
-      integer CCTK_Equals
-
-      CCTK_REAL eps, m, a0, omega0, bround, atype, aretarded
-
-      data firstcall /.true./
-      save firstcall, eps, m, a0, omega0, bround, atype, aretarded
-
-C     temps
-      CCTK_REAL a, omega, tau, f
-      CCTK_REAL c, c0, c1, c2, c3 
-      CCTK_REAL rho, r, sinp, cosp, phi, sint, cost, tx, ty, tz, px, py, pz
-      CCTK_REAL a2, b2, bx, by, bz, detgd
-
-C     get parameters of the exact solution.
-
-      if (firstcall) then
-
-         firstcall = .false.
-
-         eps = fakebinary_eps
-         m   = fakebinary_m
-         a0  = fakebinary_a0
-
-         omega0 = fakebinary_omega0
-         bround = fakebinary_bround
-
-         bround = max(bround, eps)
-
-         if (CCTK_Equals(fakebinary_atype,"constant").ne.0) then
-            atype = 0.d0
-         elseif (CCTK_Equals(fakebinary_atype,"quadrupole").ne.0) then
-            atype = 1.d0
-         else
-            call CCTK_Warn(0,"Unknown value of parameter fakebinary_atype")
-         endif
-
-         if (fakebinary_retarded.ne.0) then
-            aretarded = 1.d0
-         else
-            aretarded = 0.d0
-         endif
-
-      end if
-
-C     spherical coordinates
-
-      rho = max(sqrt(x**2 + y**2), eps)
-      r = sqrt(rho**2 + z**2)
-      sinp = y / rho
-      cosp = x / rho
-      phi = acos(cosp)
-      sint = rho / r
-      cost = z / r
-      tx = cost*cosp
-      ty = cost*sinp
-      tz = sint
-      px = - sinp
-      py = cosp
-      pz = 0
-
-C     distance function a(T-R)
-
-      tau = 5.d0/128.d0 * a0**4 / m**3
-      a = a0 * (1.d0 - atype * 4.d0*(t - aretarded*r)/tau)**(0.25d0)
-
-C     orbital frequency omega(T-R)
-
-      omega = 0.5d0*(m/a**3)**2
-
-C     1/r type potential f
-
-      c = y**2 + z**2 + bround**2;
-      f = ((x-a)**2 + c)**(-0.5d0) + ((x+a)**2 + c)**(-0.5d0)
-
-C     the three metric, tt part
-
-      c3 = 2.d0*(phi + omega*r)
-      c0 =  - 4.d0 * m * a**2 * omega**3 * (omega*r)**4 
-     .     / (1 + (omega*r)**2)**(2.5d0)
-      c1 = (1 + cost**2) * cos(c3) * c0
-      c2 = - 2.d0 * cost * sin(c3) * c0
-      gdxx = c1 * (tx*tx - px*px) +  c2 * (tx*px + px*tx)
-      gdxy = c1 * (tx*ty - px*py) +  c2 * (tx*py + px*ty)
-      gdxz = c1 * (tx*tz - px*pz) +  c2 * (tx*pz + px*tz)
-      gdyy = c1 * (ty*ty - py*py) +  c2 * (ty*py + py*ty)
-      gdyz = c1 * (ty*tz - py*pz) +  c2 * (ty*pz + py*tz)
-      gdzz = c1 * (tz*tz - pz*pz) +  c2 * (tz*pz + pz*tz)
-
-C     the three metric, add conformally flat part
-
-      c = (1.d0 + m * f)**2
-      gdxx = gdxx + c
-      gdyy = gdyy + c
-      gdzz = gdzz + c
-
-C     the shift vector and covector
-
-      c = (1.d0 -  2*m*a**2/(r**2+a**2) * f) * omega * rho
-      bx = c * px
-      by = c * py
-      bz = c * pz
-      gdtx = gdxx*bx + gdxy*by + gdxz*bz
-      gdty = gdxy*bx + gdyy*by + gdyz*bz
-      gdtz = gdxz*bx + gdyz*by + gdzz*bz
-      b2 = gdtx*bx + gdty*by + gdtz*bz
-
-C     lapse squard and time-time component of the four metric
-
-      a2 = (1.d0 - m * f)**2
-      gdtt = b2 - a2
-
-C     done with metric, find its inverse
-C     inverse three metric
-
-      detgd = -(gdxz**2*gdyy) + 2*gdxy*gdxz*gdyz - gdxx*gdyz**2
-     .        - gdxy**2*gdzz -  gdxx*gdyy*gdzz
-      guxx = (-gdyz**2 + gdyy*gdzz) / detgd
-      guxy = (gdxz*gdyz - gdxy*gdzz) / detgd
-      guxz = (-(gdxz*gdyy) + gdxy*gdyz) / detgd
-      guyy = (-gdxz**2 + gdxx*gdzz) / detgd 
-      guyz = (gdxy*gdxz - gdxx*gdyz) / detgd 
-      guzz = (-gdxy**2 + gdxx*gdyy) / detgd 
-
-C     inverse four metric
-
-      gutt = - 1.d0/a2
-      gutx = bx/a2
-      guty = by/a2
-      gutz = bz/a2
-      guxx = guxx - bx*bx/a2
-      guxy = guxy - bx*by/a2
-      guxz = guxz - bx*bz/a2
-      guyy = guyy - by*by/a2
-      guyz = guyz - by*bz/a2
-      guzz = guzz - bz*bz/a2
-
-C     done!
-      return
-      end
diff --git a/src/finkelstein.F b/src/finkelstein.F
deleted file mode 100644
index 33f6737..0000000
--- a/src/finkelstein.F
+++ /dev/null
@@ -1,66 +0,0 @@
-#include "cctk.h"
-#include "cctk_Parameters.h"
-
-      subroutine  finkelstein(
-     $     x, y, z, t,
-     $     gdtt, gdtx, gdty, gdtz, 
-     $     gdxx, gdyy, gdzz, gdxy, gdyz, gdzx,
-     $     gutt, gutx, guty, gutz, 
-     $     guxx, guyy, guzz, guxy, guyz, guzx)
-
-      implicit none
-
-      DECLARE_CCTK_PARAMETERS
-
-      CCTK_REAL x, y, z, t
-      CCTK_REAL gdtt, gdtx, gdty, gdtz, 
-     $       gdxx, gdyy, gdzz, gdxy, gdyz, gdzx,
-     $       gutt, gutx, guty, gutz, 
-     $       guxx, guyy, guzz, guxy, guyz, guzx
-
-      logical firstcall
-
-      CCTK_REAL eps, m
-
-      data firstcall /.true./
-      save firstcall, eps, m
-
-      CCTK_REAL r
-
-C     Get parameters of the exact solution.
-
-      if (firstcall) then
-
-         eps = kerrschild_eps
-         m   = kerrschild_m
-
-         firstcall = .false.
-
-      end if
-
-      r = max(sqrt(x**2 + y**2 + z**2), eps)
-
-      gdtt = - (1.d0 - 2.d0 * m / r)
-      gdtx = 2.d0 * m * x / r**2
-      gdty = 2.d0 * m * y / r**2
-      gdtz = 2.d0 * m * z / r**2
-      gdxx = 1.d0 + 2.d0 * m * x**2 / r**3
-      gdyy = 1.d0 + 2.d0 * m * y**2 / r**3
-      gdzz = 1.d0 + 2.d0 * m * z**2 / r**3
-      gdxy = 2.d0 * m * x * y / r**3
-      gdyz = 2.d0 * m * y * z / r**3
-      gdzx = 2.d0 * m * z * x / r**3
-
-      gutt = - (1.d0 + 2.d0 * m / r)
-      gutx = 2.d0 * m * x / r**2
-      guty = 2.d0 * m * y / r**2
-      gutz = 2.d0 * m * z / r**2
-      guxx = 1.d0 - 2.d0 * m * x**2 / r**3
-      guyy = 1.d0 - 2.d0 * m * y**2 / r**3
-      guzz = 1.d0 - 2.d0 * m * z**2 / r**3
-      guxy = - 2.d0 * m * x * y / r**3
-      guyz = - 2.d0 * m * y * z / r**3
-      guzx = - 2.d0 * m * z * x / r**3
-
-      return
-      end
diff --git a/src/flatfunny.F b/src/flatfunny.F
deleted file mode 100644
index 48e4864..0000000
--- a/src/flatfunny.F
+++ /dev/null
@@ -1,105 +0,0 @@
-#include "cctk.h"
-#include "cctk_Parameters.h"
-
-      subroutine flatfunny(
-     $     x, y, z, t,
-     $     gdtt, gdtx, gdty, gdtz, 
-     $     gdxx, gdyy, gdzz, gdxy, gdyz, gdzx,
-     $     gutt, gutx, guty, gutz, 
-     $     guxx, guyy, guzz, guxy, guyz, guzx)
-
-c     The metric given here corresponds to that of flat spacetime
-c     but with non-trivial spatial coordinates.  Basically, I take
-c     the flat metric in spherical coordinates, define a new
-c     radial coordinate such that:
-c
-c     r  =  r    ( 1  -  a f(r   ) )
-c            new              new
-c
-c     where f(r) is a gaussian centered at r=0 with amplitude 1.
-c     Finally, I transform back to cartesian coordinates.
-c     For  0 <= a < 1, the transformation above is monotonic.
-
-      implicit none
-
-      DECLARE_CCTK_PARAMETERS
-
-c     Input.
-
-      CCTK_REAL x,y,z,t
-
-c     Output.
-
-      CCTK_REAL gdtt, gdtx, gdty, gdtz, 
-     $          gdxx, gdyy, gdzz, gdxy, gdyz, gdzx,
-     $          gutt, gutx, guty, gutz, 
-     $          guxx, guyy, guzz, guxy, guyz, guzx
-
-c     Internal.
-
-      CCTK_REAL a,s
-      CCTK_REAL r,r2,x2,y2,z2
-      CCTK_REAL f,fp,g11,g22
-      CCTK_REAL det,zero,one,two
-
-c     Define numbers.
-
-      zero = 0.0d0
-      one  = 1.0d0
-      two  = 2.0d0
-
-c     Read parameters
-
-      a = flatfunny_a
-      s = flatfunny_s
-
-c     Find transformation function.
-
-      x2 = x**2
-      y2 = y**2
-      z2 = z**2
-
-      r2 = x2 + y2 + z2
-      r  = dsqrt(r2)
-
-      f  = dexp(-r2/s**2)
-      fp = - two*r/s**2*f
-
-c     Give metric components.
-
-      g11 = (one - a*(f + r*fp))**2
-      g22 = (one - a*f)**2
-
-      gdtt = - one
-      gdtx = zero
-      gdty = zero
-      gdtz = zero
-
-      gdxx = (x2*g11 + (y2 + z2)*g22)/r2
-      gdyy = (y2*g11 + (x2 + z2)*g22)/r2
-      gdzz = (z2*g11 + (x2 + y2)*g22)/r2
-
-      gdxy = x*y*(g11 - g22)/r2
-      gdzx = x*z*(g11 - g22)/r2
-      gdyz = y*z*(g11 - g22)/r2
-
-c     Find inverse metric.
-
-      gutt = - one
-      gutx = zero
-      guty = zero
-      gutz = zero
-
-      det = gdxx*gdyy*gdzz + two*gdxy*gdzx*gdyz
-     .    - gdxx*gdyz**2 - gdyy*gdzx**2 - gdzz*gdxy**2
-
-      guxx = (gdyy*gdzz - gdyz**2)/det
-      guyy = (gdxx*gdzz - gdzx**2)/det
-      guzz = (gdxx*gdyy - gdxy**2)/det
-
-      guxy = (gdzx*gdyz - gdzz*gdxy)/det
-      guyz = (gdxy*gdzx - gdxx*gdyz)/det
-      guzx = (gdxy*gdyz - gdyy*gdzx)/det
-
-      return
-      end
diff --git a/src/flatschwarz.F b/src/flatschwarz.F
deleted file mode 100644
index a228741..0000000
--- a/src/flatschwarz.F
+++ /dev/null
@@ -1,70 +0,0 @@
-#include "cctk.h"
-#include "cctk_Parameters.h"
-
-      subroutine  flatschwarz(
-     $     x, y, z, t,
-     $     gdtt, gdtx, gdty, gdtz, 
-     $     gdxx, gdyy, gdzz, gdxy, gdyz, gdzx,
-     $     gutt, gutx, guty, gutz, 
-     $     guxx, guyy, guzz, guxy, guyz, guzx)
-
-      implicit none
-
-      DECLARE_CCTK_PARAMETERS
-
-      CCTK_REAL x, y, z, t
-      CCTK_REAL gdtt, gdtx, gdty, gdtz, 
-     $       gdxx, gdyy, gdzz, gdxy, gdyz, gdzx,
-     $       gutt, gutx, guty, gutz, 
-     $       guxx, guyy, guzz, guxy, guyz, guzx
-
-      logical firstcall
-
-      CCTK_REAL eps, m
-
-      data firstcall /.true./
-      save firstcall, eps, m
-
-      CCTK_REAL r, bx, by, bz, b2
-
-C     Get parameters of the exact solution.
-
-      if (firstcall) then
-
-         eps = kerrschild_eps
-         m   = kerrschild_m
-
-         firstcall = .false.
-
-      end if
-
-      r = max(sqrt(x**2 + y**2 + z**2), eps)
-      bx = sqrt(2.d0 * m / r) * x / r
-      by = sqrt(2.d0 * m / r) * y / r
-      bz = sqrt(2.d0 * m / r) * z / r
-      b2 = 2.d0 * m / r
-
-      gdtt = - 1.d0 + b2
-      gdtx = bx
-      gdty = by
-      gdtz = bz
-      gdxx = 1.d0 
-      gdyy = 1.d0 
-      gdzz = 1.d0 
-      gdxy = 0.d0
-      gdyz = 0.d0
-      gdzx = 0.d0
-
-      gutt = - 1.d0
-      gutx = bx
-      guty = by
-      gutz = bz
-      guxx = 1.d0 - bx**2
-      guyy = 1.d0 - by**2
-      guzz = 1.d0 - bz**2
-      guxy = - bx * by
-      guyz = - by * bz
-      guzx = - bz * bx
-
-      return
-      end
diff --git a/src/flatshift.F b/src/flatshift.F
deleted file mode 100644
index b6feccd..0000000
--- a/src/flatshift.F
+++ /dev/null
@@ -1,107 +0,0 @@
-#include "cctk.h"
-#include "cctk_Parameters.h"
-
-      subroutine flatshift(
-     $     x, y, z, t,
-     $     gdtt, gdtx, gdty, gdtz, 
-     $     gdxx, gdyy, gdzz, gdxy, gdyz, gdzx,
-     $     gutt, gutx, guty, gutz, 
-     $     guxx, guyy, guzz, guxy, guyz, guzx)
-
-c     The metric given here corresponds to that of flat spacetime
-c     but with non-trivial slicing and a shift vector such that
-c     the resulting metric is still time independent.  I take
-c     the flat metric in spherical coordinates and define a new
-c     time coordinate as:
-c
-c     t    =   t   -  a f(r)
-c      new
-c
-c     where f(r) is a gaussian centered at r=0 with amplitude 1.
-c     Finally, I transform back to cartesian coordinates.
-c     For  -1 < fp < 1, the transformation above results in spatial
-c     slices.
-
-      implicit none
-
-      DECLARE_CCTK_PARAMETERS
-
-c     Input.
-
-      CCTK_REAL x,y,z,t
-
-c     Output.
-
-      CCTK_REAL gdtt, gdtx, gdty, gdtz, 
-     $          gdxx, gdyy, gdzz, gdxy, gdyz, gdzx,
-     $          gutt, gutx, guty, gutz, 
-     $          guxx, guyy, guzz, guxy, guyz, guzx
-
-
-c     Internal.
-
-      CCTK_REAL a,s
-      CCTK_REAL r,r2,x2,y2,z2
-      CCTK_REAL f,fp,fpr,fpr2
-      CCTK_REAL zero,one,two
-
-c     Define numbers.
-
-      zero = 0.0d0
-      one  = 1.0d0
-      two  = 2.0d0
-
-c     Read parameters
-
-      a = flatshift_a
-      s = flatshift_s
-
-c     Find transformation function.
-
-      x2 = x**2
-      y2 = y**2
-      z2 = z**2
-
-      r2 = x2 + y2 + z2
-      r  = dsqrt(r2)
-
-      f    = a*dexp(-r2/s**2)
-      fp   = - two*f*r/s**2
-      fpr  = fp/r
-      fpr2 = fpr**2
-
-c     Give metric components.
-
-      gdtt = - one
-
-      gdtx = - x*fpr
-      gdty = - y*fpr
-      gdtz = - z*fpr
-
-      gdxx = one - x2*fpr2
-      gdyy = one - y2*fpr2
-      gdzz = one - z2*fpr2
-
-      gdxy = - x*y*fpr2
-      gdzx = - x*z*fpr2
-      gdyz = - y*z*fpr2
-
-c     Inverse metric. And yes, it is this simple, simpler
-c     than the metric itself.  I tripled checked!
-
-      gutt = - one + fp**2
-
-      gutx = - x*fpr
-      guty = - y*fpr
-      gutz = - z*fpr
-
-      guxx = one
-      guyy = one
-      guzz = one
-
-      guxy = zero
-      guzx = zero
-      guyz = zero
-
-      return
-      end
diff --git a/src/godel.F b/src/godel.F
deleted file mode 100644
index 811d816..0000000
--- a/src/godel.F
+++ /dev/null
@@ -1,63 +0,0 @@
-C Godel spacetime !!!!
-C See: S. Hawking, G.F.R. Ellis, The Large
-C Scale Structure of space-time, Cambridge, 1973
-C Author : D. Vulcanov (Timsioara, Romania)
-
-#include "cctk.h"
-#include "cctk_Parameters.h"
-
-      subroutine  godel(
-     $     x, y, z, t,
-     $     gdtt, gdtx, gdty, gdtz, 
-     $     gdxx, gdyy, gdzz, gdxy, gdyz, gdzx,
-     $     gutt, gutx, guty, gutz, 
-     $     guxx, guyy, guzz, guxy, guyz, guzx)
-
-      implicit none
-
-      DECLARE_CCTK_PARAMETERS
-
-      CCTK_REAL x, y, z, t
-      CCTK_REAL gdtt, gdtx, gdty, gdtz, 
-     $       gdxx, gdyy, gdzz, gdxy, gdyz, gdzx,
-     $       gutt, gutx, guty, gutz, 
-     $       guxx, guyy, guzz, guxy, guyz, guzx
-
-      logical firstcall
-    
-      CCTK_REAL arad
-
-      data firstcall /.true./
-      save firstcall
-      if (firstcall) then
-          arad = godel_a
-          PRINT *, firstcall, arad, godel_a
-          firstcall = .false.
-      end if
-
-
-     
-      gdtt = -arad*arad
-      gdtx = 0.d0
-      gdty = 0.d0
-      gdtz = 0.d0
-      gdxx = (1/2)*arad*arad*exp(x)*exp(x)
-      gdyy = -arad*arad
-      gdzz =  arad*arad*exp(x)
-      gdxy = 0.d0
-      gdyz = 0.d0
-      gdzx = arad*arad*exp(x)
-
-      gutt = -1/(arad*arad)
-      gutx = 0.d0
-      guty = 0.d0
-      gutz = 0.d0
-      guxx = -2/(arad*arad*exp(x)*exp(x))
-      guyy = -1/(arad*arad)
-      guzz = -1/(arad*arad*exp(x))
-      guxy = 0.d0
-      guyz = 0.d0 
-      guzx = 2/(arad*arad*exp(x))
-
-      return
-      end
diff --git a/src/kasner_l.F b/src/kasner_l.F
deleted file mode 100644
index 136cc0d..0000000
--- a/src/kasner_l.F
+++ /dev/null
@@ -1,68 +0,0 @@
-C Kasner_like metric spacetime !!!!
-C See: L. Pimentel, Int. Journ. of Theor. Physics,
-C, vol. 32, No. 6, 1993, p. 979 (and the references
-C cited here), and S. Gotlober, et. al. Early Evolution
-C of the Universe and Formation Structure, Akad. Verlag, 1990
-C Author : D. Vulcanov
-
-#include "cctk.h"
-#include "cctk_Parameters.h"
-
-      subroutine  kasner(
-     $     x, y, z, t,
-     $     gdtt, gdtx, gdty, gdtz, 
-     $     gdxx, gdyy, gdzz, gdxy, gdyz, gdzx,
-     $     gutt, gutx, guty, gutz, 
-     $     guxx, guyy, guzz, guxy, guyz, guzx)
-
-      implicit none
-
-      DECLARE_CCTK_PARAMETERS
-
-      CCTK_REAL x, y, z, t
-      CCTK_REAL gdtt, gdtx, gdty, gdtz, 
-     $       gdxx, gdyy, gdzz, gdxy, gdyz, gdzx,
-     $       gutt, gutx, guty, gutz, 
-     $       guxx, guyy, guzz, guxy, guyz, guzx
-
-      logical firstcall
-    
-      CCTK_REAL qq,a1,a3
-
-      data firstcall /.true./
-      save firstcall
-      if (firstcall) then
-          qq = kasner_q
-          PRINT *, firstcall, qq, kasner_q
-          firstcall = .false.
-      end if
- 
-      a1= t**qq
-      a3= t**(1-2*qq)    
-
-      gdtt = -1.d0 
-      gdtx = 0.d0 
-      gdty = 0.d0
-      gdtz = 0.d0
-      gdxx = a1**2 
-      gdyy = a1**2
-      gdzz = a3**2
-      gdxy = 0.d0
-      gdyz = 0.d0
-      gdzx = 0.d0
-
-      gutt = -1.d0
-      gutx = 0.d0
-      guty = 0.d0
-      gutz = 0.d0
-      guxx = a1**(-2)
-      guyy = a1**(-2)
-      guzz = a3**(-2)
-      guxy = 0.d0
-      guyz = 0.d0
-      guzx = 0.d0
- 
-
-
-      return
-      end
diff --git a/src/kerrcart.F b/src/kerrcart.F
deleted file mode 100644
index c2db4da..0000000
--- a/src/kerrcart.F
+++ /dev/null
@@ -1,65 +0,0 @@
-C Kerr metric in cartesian coordinates spacetime !!!!
-C See: S. Hawking, G.F.R. Ellis, The Large
-C Scale Structure of space-time, Cambridge, 1973 and 
-C "The phone book" - Misner, Thorne, Wheeler 
-C Author : D. Vulcanov (Timisoara, Romania)
-
-#include "cctk.h"
-#include "cctk_Parameters.h"
-
-      subroutine  kerr(
-     $     x, y, z, t,
-     $     gdtt, gdtx, gdty, gdtz, 
-     $     gdxx, gdyy, gdzz, gdxy, gdyz, gdzx,
-     $     gutt, gutx, guty, gutz, 
-     $     guxx, guyy, guzz, guxy, guyz, guzx)
-
-      implicit none
-
-      DECLARE_CCTK_PARAMETERS
-
-      CCTK_REAL x, y, z, t
-      CCTK_REAL gdtt, gdtx, gdty, gdtz, 
-     $       gdxx, gdyy, gdzz, gdxy, gdyz, gdzx,
-     $       gutt, gutx, guty, gutz, 
-     $       guxx, guyy, guzz, guxy, guyz, guzx
-
-      logical firstcall
-    
-      CCTK_REAL arad, marad
-
-      data firstcall /.true./
-      save firstcall
-      if (firstcall) then
-          arad = kerrc_a
-         marad = kerrc_m
-          PRINT *, firstcall, arad, kerrc_a, kerrc_m
-          firstcall = .false.
-      end if
-
-
-     
-      gdtt = -(y**2*arad**2+x**2-2*marad*x)/(x**2+y**2*arad**2)
-      gdtx = 2*(arad*marad*x*(y**2-1))/(x**2+y**2*arad**2)
-      gdty = 0.d0
-      gdtz = 0.d0
-      gdxx = -(x**4+x**2*arad**2+2*arad**2*marad*x+arad**2*y**2*x**2 - 2*arad**2*y**2*marad*x+arad**4*y**2)*(y**2-1)/(x**2+arad**2*y**2)   
-      gdyy = (x**2+y**2*arad**2)/(x**2-2*marad*x+arad**2) 
-      gdzz =  -(x**2+y**2*arad**2)/(y**2-1)
-      gdxy = 0.d0
-      gdyz = 0.d0
-      gdzx = 0.d0
-
-      gutt = -(-x**4-x**2*arad**2-2*arad**2*marad*x-arad**2*y**2*x**2 +2*arad**2*y**2*marad*x-arad**4*y**2)/((x**2+arad**2*y**2)*(-x**2+2*marad*x-arad**2))
-      gutx = 2*(arad*marad*x)/((x**2+arad**2*y**2)*(-x**2+2*marad*x-arad**2))
-      guty = 0.d0
-      gutz = 0.d0
-      guxx = -(-arad**2*y**2-x**2+2*marad*x)/((x**2+arad**2*y**2)*(y**2-1)*(-x**2+2*marad*x-arad**2))
-      guyy = -(-x**2+2*marad*x-arad**2)/(x**2+arad**2*y**2)
-      guzz = -(y**2-1)/(x**2+arad**2*y**2)
-      guxy = 0.d0
-      guyz = 0.d0
-      guzx = 0.d0
-
-      return
-      end
diff --git a/src/kerrschild.F b/src/kerrschild.F
deleted file mode 100644
index d1d5e7c..0000000
--- a/src/kerrschild.F
+++ /dev/null
@@ -1,127 +0,0 @@
-#include "cctk.h"
-#include "cctk_Parameters.h"
-
-C     Kerr-Schild form of boosted rotating black hole.
-C     Program g_ab = eta_ab + H l_a l_b, g^ab = eta^ab - H l^a l^b.
-C     Here eta_ab is Minkowski in Cartesian coordinates, H is a scalar,
-C     and l is a null vector.
-
-      subroutine kerrschild(
-     $     x, y, z, t,
-     $     gdtt, gdtx, gdty, gdtz, 
-     $     gdxx, gdyy, gdzz, gdxy, gdyz, gdzx,
-     $     gutt, gutx, guty, gutz, 
-     $     guxx, guyy, guzz, guxy, guyz, guzx)
-
-      implicit none
-
-      DECLARE_CCTK_PARAMETERS
-
-      CCTK_REAL x, y, z, t
-      CCTK_REAL gdtt, gdtx, gdty, gdtz, 
-     $       gdxx, gdyy, gdzz, gdxy, gdyz, gdzx,
-     $       gutt, gutx, guty, gutz, 
-     $       guxx, guyy, guzz, guxy, guyz, guzx
-
-      CCTK_REAL gamma, t0, z0, x0, y0, rho02, r02, r0, costheta0, 
-     $     lt0, lx0, ly0, lz0, hh, lt, lx, ly, lz
-
-      logical firstcall
-      CCTK_REAL boostv, eps, m, a
-      data firstcall /.true./
-      save firstcall, boostv, eps, m, a
-
-
-C     Get parameters of the exact solution.
-
-      if (firstcall) then
-
-         boostv = kerrschild_boostv
-         eps    = kerrschild_eps
-
-         m = kerrschild_m
-         a = kerrschild_a
-
-         firstcall = .false.
-
-      end if
-
-C     Boost factor.
-
-      gamma = 1.d0 / sqrt(1.d0 - boostv**2)
-
-C     Lorentz transform t,x,y,z -> t0,x0,y0,z0. 
-C     t0 is never used, but is here for illustration, and we introduce
-C     x0 and y0 also only for clarity.
-C     Note that z0 = 0 means z = vt for the BH.
-
-      t0 = gamma * (t - boostv * z)
-      z0 = gamma * (z - boostv * t)
-      x0 = x
-      y0 = y
-
-C     Coordinate distance to center of black hole. Note it moves!
-
-      rho02 = x0**2 + y0**2 + z0**2
-
-C     Spherical auxiliary coordinate r and angle theta in BH rest frame.
-
-      r02 = 0.5d0 * (rho02 - a**2) 
-     $     + sqrt(0.25d0 * (rho02 - a**2)**2 + a**2 * z0**2)
-      if (r02 .lt. eps) then
-         write(6,*) 'r02 too small in kerrschild'
-         write(6,*) x0, y0, z0
-         write(6,*) rho02, a**2, r02
-      end if
-      r02 = max(r02,eps)
-      r0 = sqrt(r02)
-      costheta0 = z0 / r0
-      
-C     Coefficient H. Note this transforms as a scalar, so does not carry
-C     the suffix 0.
-      hh = m * r0 / (r0**2 + a**2 * costheta0**2)
-      
-C     Components of l_a in rest frame. Note indices down.
-      lt0 = 1.d0
-      lx0 = (r0 * x0 + a * y0) / (r0**2 + a**2)
-      ly0 = (r0 * y0 - a * x0) / (r0**2 + a**2)
-      lz0 = z0 / r0
-
-C     Now boost it to coordinates x, y, z, t.
-C     This is the reverse Lorentz transformation, but applied
-C     to a one-form, so the sign of boostv is the same as the forward
-C     Lorentz transformation applied to the coordinates.
-
-      lt = gamma * (lt0 - boostv * lz0) 
-      lz = gamma * (lz0 - boostv * lt0) 
-      lx = lx0
-      ly = ly0
-      
-C     Down metric. g_ab = flat_ab + H l_a l_b
-
-      gdtt = - 1.d0 + 2.d0 * hh * lt * lt
-      gdtx = 2.d0 * hh * lt * lx
-      gdty = 2.d0 * hh * lt * ly
-      gdtz = 2.d0 * hh * lt * lz
-      gdxx = 1.d0 + 2.d0 * hh * lx * lx
-      gdyy = 1.d0 + 2.d0 * hh * ly * ly
-      gdzz = 1.d0 + 2.d0 * hh * lz * lz
-      gdxy = 2.d0 * hh * lx * ly
-      gdyz = 2.d0 * hh * ly * lz
-      gdzx = 2.d0 * hh * lz * lx
-
-C     Up metric. g^ab = flat^ab - H l^a l^b.
-C     Notice that g^ab = g_ab and l^i = l_i and l^0 = - l_0 in flat spacetime.
-      gutt = - 1.d0 - 2.d0 * hh * lt * lt 
-      gutx = 2.d0 * hh * lt * lx
-      guty = 2.d0 * hh * lt * ly
-      gutz = 2.d0 * hh * lt * lz
-      guxx = 1.d0 - 2.d0 * hh * lx * lx
-      guyy = 1.d0 - 2.d0 * hh * ly * ly
-      guzz = 1.d0 - 2.d0 * hh * lz * lz
-      guxy = - 2.d0 * hh * lx * ly
-      guyz = - 2.d0 * hh * ly * lz
-      guzx = - 2.d0 * hh * lz * lx
-
-      return
-      end
diff --git a/src/milne.F b/src/milne.F
deleted file mode 100644
index 0d78fa3..0000000
--- a/src/milne.F
+++ /dev/null
@@ -1,56 +0,0 @@
-C De Milne spacetime metric 
-C Suggested by Matteo Rossi and E. Onofri (Univ. di Parma, Italy)
-C They inted to use thsi metric for simulating an Pre-Big-Bang
-C Cosmology, as proposed by Veneziano some year ago 
-C Author : D. Vulcanov (Timsoara, Romania)
-#include "cctk.h"
-
-      subroutine milne(
-     $     x, y, z, t,
-     $     gdtt, gdtx, gdty, gdtz, 
-     $     gdxx, gdyy, gdzz, gdxy, gdyz, gdzx,
-     $     gutt, gutx, guty, gutz, 
-     $     guxx, guyy, guzz, guxy, guyz, guzx)
-
-      implicit none
-
-c     Input.
-
-      CCTK_REAL x, y, z, t
-
-c     Output.
-
-      CCTK_REAL gdtt, gdtx, gdty, gdtz, 
-     $       gdxx, gdyy, gdzz, gdxy, gdyz, gdzx,
-     $       gutt, gutx, guty, gutz, 
-     $       guxx, guyy, guzz, guxy, guyz, guzx,
-     $       coef
-
-
-      coef= t**2/(1+ x**2 +y**2+ z**2)
-
-      gdtt = -1.d0 
-      gdtx = 0.d0  
-      gdty = 0.d0 
-      gdtz = 0.d0 
-      gdxx =  coef*(1+y**2+z**2) 
-      gdyy =  coef*(1+x**2+z**2) 
-      gdzz =  coef*(1+x**2+y**2)
-      gdxy = -coef*x*y 
-      gdyz = -coef*y*z 
-      gdzx = -coef*x*z
-
-      gutt = -1.d0 
-      gutx = 0.d0 
-      guty = 0.d0 
-      gutz = 0.d0 
-      guxx = (1+x**2)/t**2
-      guyy = (1+y**2)/t**2
-      guzz = (1+z**2)/t**2
-      guxy = x*y/t**2
-      guyz = y*z/t**2 
-      guzx = x*z/t**2
-
-
-      return
-      end
diff --git a/src/minkowski.F b/src/minkowski.F
deleted file mode 100644
index ed76b17..0000000
--- a/src/minkowski.F
+++ /dev/null
@@ -1,46 +0,0 @@
-#include "cctk.h"
-
-      subroutine minkowski(
-     $     x, y, z, t,
-     $     gdtt, gdtx, gdty, gdtz, 
-     $     gdxx, gdyy, gdzz, gdxy, gdyz, gdzx,
-     $     gutt, gutx, guty, gutz, 
-     $     guxx, guyy, guzz, guxy, guyz, guzx)
-
-      implicit none
-
-c     Input.
-
-      CCTK_REAL x, y, z, t
-
-c     Output.
-
-      CCTK_REAL gdtt, gdtx, gdty, gdtz, 
-     $       gdxx, gdyy, gdzz, gdxy, gdyz, gdzx,
-     $       gutt, gutx, guty, gutz, 
-     $       guxx, guyy, guzz, guxy, guyz, guzx
-
-      gdtt = -1.d0 
-      gdtx = 0.d0 
-      gdty = 0.d0 
-      gdtz = 0.d0 
-      gdxx = 1.d0 
-      gdyy = 1.d0 
-      gdzz = 1.d0 
-      gdxy = 0.d0 
-      gdyz = 0.d0 
-      gdzx = 0.d0
-
-      gutt = -1.d0 
-      gutx = 0.d0 
-      guty = 0.d0 
-      gutz = 0.d0 
-      guxx = 1.d0 
-      guyy = 1.d0 
-      guzz = 1.d0 
-      guxy = 0.d0 
-      guyz = 0.d0 
-      guzx = 0.d0
-
-      return
-      end
diff --git a/src/multiBH.F b/src/multiBH.F
deleted file mode 100644
index 8125200..0000000
--- a/src/multiBH.F
+++ /dev/null
@@ -1,133 +0,0 @@
-c=======================================================================
-c23456789012345678901234567890123456789012345678901234567890123456789012
-c-----------------------------------------------------------------------
-c by Hisaaki Shinkai   shinkai@wurel.wustl.edu  19980603
-c-----------------------------------------------------------------------
-c  This is for maximally charged multi BH solutions such as 
-c  Majumdar-Papapetrou (1947) or Kastor-Traschen (1993) solution.
-c  See also doc/KTsol.tex for brief review of this solution. 
-c-----------------------------------------------------------------------
-c  For usage:  in your par file
-c     exactmodel = "multiBH"
-c     kt_hubble       = 0.0   # 0.0 means MP solution
-c     kt_nBH          = 2     # number of BHs (upto 4, currently)
-c     m_bh1,co_bh1x,co_bh1y,co_bh1z =  1.0, -2.0,0.0,0.0  # masses and 
-c     m_bh2,co_bh2x,co_bh2y,co_bh2z =  1.0,  2.0,0.0,0.0  #   locations
-c-----------------------------------------------------------------------
-
-#include "cctk.h"
-#include "cctk_Parameters.h"
-
-      subroutine  multi_BHs(
-     $     x, y, z, t,
-     $     gdtt, gdtx, gdty, gdtz, 
-     $     gdxx, gdyy, gdzz, gdxy, gdyz, gdzx,
-     $     gutt, gutx, guty, gutz, 
-     $     guxx, guyy, guzz, guxy, guyz, guzx)
-
-      implicit none
-
-      DECLARE_CCTK_PARAMETERS
-
-c     Input.
-      CCTK_REAL x, y, z, t
-
-c     Output.
-      CCTK_REAL gdtt, gdtx, gdty, gdtz, 
-     $       gdxx, gdyy, gdzz, gdxy, gdyz, gdzx,
-     $       gutt, gutx, guty, gutz, 
-     $       guxx, guyy, guzz, guxy, guyz, guzx
-
-      logical firstcall
-
-      CCTK_REAL kt_r, kt_aa, kt_omega
-      CCTK_REAL one,zero
-      CCTK_REAL kt_xbh(10),kt_ybh(10),kt_zbh(10),kt_mbh(10)
-
-      integer i
-
-      data firstcall /.true./
-      save firstcall,kt_xbh,kt_ybh,kt_zbh,kt_mbh
-
-c     Get parameters of the exact solution.
-
-      if (firstcall) then
-
-         write(*,*) ' welcome to Kastor-Traschen (Majumdar-Papapetrou)'
-
-         if(kt_nBH.ge.1) then
-            kt_xbh(1) = co_bh1x
-            kt_ybh(1) = co_bh1y
-            kt_zbh(1) = co_bh1z
-            kt_mbh(1) = m_bh1
-         endif
-
-         if(kt_nBH.ge.2) then
-            kt_xbh(2) = co_bh2x
-            kt_ybh(2) = co_bh2y
-            kt_zbh(2) = co_bh2z
-            kt_mbh(2) = m_bh2
-         endif
-
-         if(kt_nBH.ge.3) then
-            kt_xbh(3) = co_bh3x
-            kt_ybh(3) = co_bh3y
-            kt_zbh(3) = co_bh3z
-            kt_mbh(3) = m_bh3
-         endif
-
-         if(kt_nBH.ge.4) then
-            kt_xbh(4) = co_bh4x
-            kt_ybh(4) = co_bh4y
-            kt_zbh(4) = co_bh4z
-            kt_mbh(4) = m_bh4
-         endif
-
-         write(*,*) 'time=',t
-         write(*,*) '   mass    BH(x,y,z) '
-
-         do i=1,kt_nBH
-            write(*,'(4e12.3)') kt_mbh(i),kt_xbh(i),kt_ybh(i),kt_zbh(i)
-         enddo
-
-         firstcall = .false.
-
-      end if
-
-      one =1.0
-      zero=0.0 
-
-      kt_aa=exp(kt_hubble*t)
-      kt_omega=1.0
-
-      do i=1,kt_nBH
-         kt_r=sqrt((x-kt_xbh(i))**2+(y-kt_ybh(i))**2+(z-kt_zbh(i))**2)
-         kt_omega= kt_omega+kt_mbh(i)/kt_aa/kt_r
-      enddo
-      
-c     write(*,*) kt_omega,kt_aa
-
-      gdtt = -1.0/kt_omega**2
-      gdtx = zero
-      gdty = zero
-      gdtz = zero
-      gdxx = (kt_aa*kt_omega)**2
-      gdyy = (kt_aa*kt_omega)**2
-      gdzz = (kt_aa*kt_omega)**2
-      gdxy = zero
-      gdyz = zero
-      gdzx = zero
-
-      gutt = one/gdtt 
-      gutx = zero
-      guty = zero
-      gutz = zero
-      guxx = one/gdxx
-      guyy = one/gdyy
-      guzz = one/gdzz
-      guxy = zero
-      guyz = zero
-      guzx = zero
-
-      return
-      end
diff --git a/src/novikov.F b/src/novikov.F
deleted file mode 100644
index c0034f8..0000000
--- a/src/novikov.F
+++ /dev/null
@@ -1,201 +0,0 @@
-#include "cctk.h"
-#include "cctk_Parameters.h"
-
-      subroutine novikov(
-     $     x, y, z, t,
-     $     gdtt, gdtx, gdty, gdtz, 
-     $     gdxx, gdyy, gdzz, gdxy, gdyz, gdzx,
-     $     gutt, gutx, guty, gutz, 
-     $     guxx, guyy, guzz, guxy, guyz, guzx)
-
-c     The metric given here corresponds to the novikov solution
-c     in isotropic coordinates, as presented first in Bruegman96 
-c     then in correct form in Cactus paper 1. This code is the code
-c     which was used for the comparisons in cactus paper 1, and is written
-c     by PW with input from BB.
-
-      implicit none
-
-      DECLARE_CCTK_PARAMETERS
-
-c     Input.
-
-      CCTK_REAL x, y, z, t
-
-c     Output.
-
-      CCTK_REAL gdtt, gdtx, gdty, gdtz, 
-     $       gdxx, gdyy, gdzz, gdxy, gdyz, gdzx,
-     $       gutt, gutx, guty, gutz, 
-     $       guxx, guyy, guzz, guxy, guyz, guzx
-
-c     Internal.
-
-      logical firstcall
-
-      CCTK_REAL mass,radius
-      CCTK_REAL r,c,psi4
-      CCTK_REAL zero,one,two
-
-      CCTK_REAL nov_dr_drmax, nov_rmax, nov_r
-      CCTK_REAL grr, gqq, detg
-
-      CCTK_REAL psi4_o_r2
-
-      data firstcall /.true./
-      save firstcall
-
-c     Find {zero,one,two}.
-
-      zero = 0.0D0
-      one  = 1.0D0
-      two  = 2.0D0
-
-c     Get parameters of the metric. (just the mass)
-      mass = 1.0D0
-
-c     Find r.
-      r = dsqrt(x**2 + y**2 + z**2)
-
-c     Find conformal factor.
-      c = mass/(two*r)
-      psi4 = (one + c)**4
-
-c     Evaluate novikov stuff. Note abs(t) since the data is time
-c     symmetric (the metric is, at least...)
-      grr = nov_dr_drmax(abs(t),abs(r))
-      gqq = nov_r(abs(t),abs(r))
-
-      grr = grr **2
-      gqq = gqq**2 / nov_rmax(abs(r))**2
-
-
-c     Find metric components.
-      psi4_o_r2 = psi4 / r**2
-
-      gdtt = - 1.0D0
-
-      gdtx = zero
-      gdty = zero
-      gdtz = zero
-
-c     This is just straightforward spherical -> cartesian I hope... ;-)
-c     Note at t=0 (grr = gqq = 1) this gives the expected result
-c     (namely diagonal psi^4, since psi4_o_r2 = psi^4 / r^2)
-      gdxx = (grr * x**2 + gqq * (y**2 + z**2)) * psi4_o_r2
-      gdyy = (grr * y**2 + gqq * (x**2 + z**2)) * psi4_o_r2
-      gdzz = (grr * z**2 + gqq * (x**2 + y**2)) * psi4_o_r2
- 
-      gdxy = (grr - gqq) * x * y * psi4_o_r2
-      gdzx = (grr - gqq) * x * z * psi4_o_r2
-      gdyz = (grr - gqq) * y * z * psi4_o_r2
-
-c     Find inverse metric.
-      gutt = one/gdtt
-      gutx = zero
-      guty = zero
-      gutz = zero
-      detg = gdtt*gdxx*gdyy*gdzz-gdtt*gdxx*gdyz**2/4.D0-gdtt*gdxy**2*gdz
-     $z/4.D0+gdtt*gdxy*gdzx*gdyz/4.D0-gdtt*gdzx**2*gdyy/4.D0
-      guxx = -gdtt*(-4.D0*gdyy*gdzz+gdyz**2)/(4.D0 * detg)
-      guxy = -gdtt*(2.D0*gdxy*gdzz-gdzx*gdyz)/(4.D0 * detg)
-      guzx = gdtt*(gdxy*gdyz-2.D0*gdzx*gdyy)/(4.D0 * detg)
-      guyy = gdtt*(4.D0*gdxx*gdzz-gdzx**2)/(4.D0 * detg)
-      guyz = -gdtt*(2.D0*gdxx*gdyz-gdxy*gdzx)/(4.D0 * detg)
-      guzz = gdtt*(4.D0*gdxx*gdyy-gdxy**2)/(4.D0 * detg)
-
-      guxx = one/psi4
-      guyy = one/psi4
-      guzz = one/psi4
-
-      guxy = zero
-      guyz = zero
-      guzx = zero
-
-      return
-      end
-
-c     These are functions which evaluate the novikov stuff.
-
-c     dr/drmax
-      CCTK_REAL function nov_dr_drmax(tauin,rbarin)
-
-      implicit none
-
-      CCTK_REAL rbarin, tauin
-      CCTK_REAL rt, nov_r, rmaxt, nov_rmax
-
-      rt = nov_r(tauin, rbarin)
-      rmaxt = nov_rmax(rbarin)
-      nov_dr_drmax = 1.5D0 - rt / (2.0D0 * rmaxt) + 
-     $     1.5D0 * sqrt(rmaxt / rt - 1.0D0) *
-     $     acos(sqrt(rt/rmaxt))
-
-      return
-      end
-
-
-
-
-c     Bisection to invert the function below. This is pretty crappy
-c     but it works.
-      CCTK_REAL function nov_r(tauin, rbarin)
-      implicit none
-c     input
-      CCTK_REAL tauin, rbarin
-
-c     funtions
-      CCTK_REAL nov_rmax, nov_tau
-
-c     temps
-      CCTK_REAL rg, drg, delt, ttmp, rmt
-      CCTK_REAL eps
-      integer nit
-      nit = 0
-      delt = 1000.0D0
-      rmt = nov_rmax(rbarin)
-      rg  = rmt
-      drg = rg / 2.0D0
-      eps = 1.d-6 * rmt
-      do while (delt .gt. eps .and. nit .lt. 100)
-         ttmp = nov_tau(rg, rmt)
-         delt = abs(tauin - ttmp)
-         if (delt .gt. eps) then 
-            if (ttmp .gt. tauin .or. rg .lt. drg) then 
-               rg = rg + drg
-c     Enforce upper bound
-               if (rg .gt. rmt) rg = rmt
-               drg = drg / 2.0D0
-            else
-               rg = rg - drg
-            endif
-         endif
-c         write (*,*) rg, ttmp, tauin
-         nit = nit + 1
-      enddo
-      if (nit .ge. 100) then
-         write (*,*) "Novikov: inversion did not converge"
-      endif
-      nov_r = rg
-      return
-      end
-
-c     Evaluate tau as a function of r and rmax
-      CCTK_REAL function nov_tau(r, rmax)
-      implicit none
-      CCTK_REAL r, rmax
-
-      nov_tau=  rmax * dsqrt(0.5D0 * r * (1.0D0 - r / rmax)) +
-     $     2.0D0 * (rmax / 2)**(3.0/2.0) *
-     $     acos (dsqrt(r/rmax))
-
-      return
-      end
-
-c     Evaluate rmax as a function of rbar
-      CCTK_REAL function nov_rmax(rbar)
-      implicit none
-      CCTK_REAL rbar
-      nov_rmax = (1.0D0 + 2.0D0*rbar)**2 / (4.0D0 * rbar)
-      return
-      end
diff --git a/src/rob-wal.F b/src/rob-wal.F
deleted file mode 100644
index b14bb03..0000000
--- a/src/rob-wal.F
+++ /dev/null
@@ -1,71 +0,0 @@
-C Robertson-Walker universe pure radiation case
-C near t=0, presure = 1/3*density , and k=0
-C We use same parameters as desitter case !
-C See: J.N. Islam, An Introduction to 
-C Mathematical Cosmology, Cambridge, 1992
-C Author : D. Vulcanov (Timisoara, Romania)
-
-#include "cctk.h"
-#include "cctk_Parameters.h"
-
-      subroutine  rob_wal(
-     $     x, y, z, t,
-     $     gdtt, gdtx, gdty, gdtz, 
-     $     gdxx, gdyy, gdzz, gdxy, gdyz, gdzx,
-     $     gutt, gutx, guty, gutz, 
-     $     guxx, guyy, guzz, guxy, guyz, guzx)
-
-      implicit none
-
-      DECLARE_CCTK_PARAMETERS
-
-      CCTK_REAL x, y, z, t
-      CCTK_REAL gdtt, gdtx, gdty, gdtz, 
-     $       gdxx, gdyy, gdzz, gdxy, gdyz, gdzx,
-     $       gutt, gutx, guty, gutz, 
-     $       guxx, guyy, guzz, guxy, guyz, guzx
-
-      logical firstcall
-    
-      CCTK_REAL arad,r2,x2,y2,z2,am
-
-      data firstcall /.true./
-      save firstcall
-      if (firstcall) then
-          arad = desitt_a
-          PRINT *, firstcall, arad, desitt_a
-          firstcall = .false.
-      end if
- 
-      x2=x*x
-      y2=y*y
-      z2=z*z
-      r2 = x2+y2+z2
-      am=arad*t     
-
-      gdtt = -1.d0 
-      gdtx = 0.d0 
-      gdty = 0.d0
-      gdtz = 0.d0
-      gdxx = 1.d0+(-1.d0 +am)*x2/r2 
-      gdyy = 1.d0+(-1.d0 +am)*y2/r2
-      gdzz = 1.d0+(-1.d0 +am)*z2/r2
-      gdxy = (-1.d0   +am)*x*y/r2
-      gdyz = (-1.d0   +am)*y*z/r2
-      gdzx = (-1.d0   +am)*x*z/r2
-
-      gutt = -1.d0
-      gutx = 0.d0
-      guty = 0.d0
-      gutz = 0.d0
-      guxx = 1.d0+(-1.d0 +1.d0/am)*x2/r2 
-      guyy = 1.d0+(-1.d0 +1.d0/am)*y2/r2
-      guzz = 1.d0+(-1.d0 +1.d0/am)*z2/r2
-      guxy = (-1.d0  +1.d0/am)*x*y/r2
-      guyz = (-1.d0  +1.d0/am)*y*z/r2
-      guzx = (-1.d0  +1.d0/am)*x*z/r2
- 
-
-
-      return
-      end
diff --git a/src/starschwarz.F b/src/starschwarz.F
deleted file mode 100644
index 1c44941..0000000
--- a/src/starschwarz.F
+++ /dev/null
@@ -1,126 +0,0 @@
-#include "cctk.h"
-#include "cctk_Parameters.h"
-
-      subroutine starschwarz(
-     $     x, y, z, t,
-     $     gdtt, gdtx, gdty, gdtz, 
-     $     gdxx, gdyy, gdzz, gdxy, gdyz, gdzx,
-     $     gutt, gutx, guty, gutz, 
-     $     guxx, guyy, guzz, guxy, guyz, guzx)
-
-c     The metric given here corresponds to a constant
-c     density star, also known as a "Schwarzschild" star.
-c     It is obviously not a vacuum solution, so if someone
-c     wants to evolve it, we will need to add here the
-c     matter variables.  For the moment, only the metric
-c     is given.
-
-C     OK, I added the matter variables 
-C     (see include/Scalar_CalcTmunu.inc) - Mitica Vulcanov 
-
-c
-c     The metric is given as a conformally flat metric.
-c     Turns out that in the original areal radius, the
-c     metric variables have a kink at the surface of the
-c     star, but they are smooth in the conformal form.
-c
-c     Thanks to Philippos Papadopoulos for suggesting
-c     the use of this metric.
-
-      implicit none
-
-      DECLARE_CCTK_PARAMETERS
-
-c     Input.
-
-      CCTK_REAL x, y, z, t
-
-c     Output.
-
-      CCTK_REAL gdtt, gdtx, gdty, gdtz, 
-     $       gdxx, gdyy, gdzz, gdxy, gdyz, gdzx,
-     $       gutt, gutx, guty, gutz, 
-     $       guxx, guyy, guzz, guxy, guyz, guzx
-
-c     Internal.
-
-      logical firstcall
-
-      CCTK_REAL mass,radius
-      CCTK_REAL r,c,psi4
-      CCTK_REAL zero,one,two
-
-      data firstcall /.true./
-      save firstcall
-
-c     Find {zero,one,two}.
-
-      zero = 0.0D0
-      one  = 1.0D0
-      two  = 2.0D0
-
-c     Get parameters of the metric.
-
-      if (firstcall) then
-
-         mass   = starschwarz_m
-         radius = starschwarz_r
-
-         firstcall = .false.
-
-      end if
-
-c     Find r.
-
-      r = dsqrt(x**2 + y**2 + z**2)
-
-c     Find conformal factor.
-
-      if (r.le.radius) then
-
-         c = mass/(two*radius)
-
-         psi4 = (one + c)**6/(one + c*(r/radius)**2)**2
-
-      else
-
-         c = mass/(two*r)
-
-         psi4 = (one + c)**4
-
-      end if
-
-c     Find metric components.
-
-      gdtt = - psi4
-
-      gdtx = zero
-      gdty = zero
-      gdtz = zero
-
-      gdxx = psi4
-      gdyy = psi4
-      gdzz = psi4
-
-      gdxy = zero
-      gdyz = zero
-      gdzx = zero
-
-c     Find inverse metric.
-
-      gutt = -one/psi4
-
-      gutx = zero
-      guty = zero
-      gutz = zero
-
-      guxx = one/psi4
-      guyy = one/psi4
-      guzz = one/psi4
-
-      guxy = zero
-      guyz = zero
-      guzx = zero
-
-      return
-      end