remove from here.

git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinInitialData/DistortedBHIVP/trunk@11 971fb155-194f-0410-9daf-e2eca44e59f5

2 files changed, 0 insertions, 945 deletions

diff --git a/src/DistortedBHIVP.F.sav b/src/DistortedBHIVP.F.sav
deleted file mode 100644
index 815d25d..0000000
--- a/src/DistortedBHIVP.F.sav
+++ /dev/null
@@ -1,681 +0,0 @@
-c     /*@@
-c     @file      DistortedBHIVP.F
-c     @date      
-c     @author    
-c     @desc 
-c     
-c     @enddesc 
-c@@   */
-      
-#include "cctk.h" 
-#include "cctk_Parameters.h"
-#include "cctk_Arguments.h"
-      
-c     /*@@
-c     @routine    DistortedBHIVP
-c     @date       
-c     @author     
-c     @desc 
-c     
-c     @enddesc 
-c     @calls     
-c     @calledby   
-c     @history 
-c     
-c     @endhistory 
-c@@   */
-
-c     Need include file from Einstein
-#include "CactusEinstein/Einstein/src/Einstein.h"
-
-      subroutine DistortedBHIVP(CCTK_FARGUMENTS)
-
-      implicit none
-
-      DECLARE_CCTK_FARGUMENTS
-      DECLARE_CCTK_PARAMETERS
-
-c     Perhaps this and others should go into cctk.h
-      integer CCTK_Equals
-      
-      real*8 :: amp,eta0,sigma,c,etamax,deta,dq,dphi
-      integer :: ne,nq,np,n
-      real*8, allocatable :: ac(:,:,:),ae(:,:,:),aw(:,:,:),an(:,:,:),
-     $     as(:,:,:),aq(:,:,:),ab(:,:,:),rhs(:,:,:),
-     $     qf(:,:,:),qfetaeta(:,:,:),qfqq(:,:,:),qfphi(:,:,:),
-     $     qfphiphi(:,:,:),
-     $     psisph(:,:,:),psiprim(:,:,:),detapsisph(:,:,:),
-     $     dqpsisph(:,:,:),dphipsisph(:,:,:),detaetapsisph(:,:,:),
-     $     detaqpsisph(:,:,:),detaphipsisph(:,:,:),dqqpsisph(:,:,:),
-     $     dqphipsisph(:,:,:),dphiphipsisph(:,:,:),
-     $     psi3d(:,:,:),detapsi3d(:,:,:),dqpsi3d(:,:,:),
-     $     dphipsi3d(:,:,:),detaetapsi3d(:,:,:),detaqpsi3d(:,:,:),
-     $     detaphipsi3d(:,:,:),dqqpsi3d(:,:,:),dqphipsi3d(:,:,:),
-     $     dphiphipsi3d(:,:,:)
-      real*8, allocatable :: etagrd(:),qgrd(:),phigrd(:)
-      real*8, allocatable :: eta(:,:,:),abseta(:,:,:),sign_eta(:,:,:),
-     $     q(:,:,:),phi(:,:,:)
-      real*8 o1,o2,o3,o4,o5,o6,o7,o8,o9,o10,o11,o12,o13,o14,o15,o16,o17,
-     $     o18,o19,o20,o21,o22,o23,o24,o25,o26,o27,o28,o29,o30,o31,o32,
-     $     o33,o34,o35,o36,o37,o38,o39,o40,o41,o42,o43,o44,o45,o46,o47,
-     $     o48,o49,o50,o51,o52,o53,o54,o55,o56,o57,o58,o59,o60,o61,o62,
-     $     o63,o64,o65,o66,o67
-      real*8 rmax,adm
-      real*8,parameter :: dbh_eps = 1.0d-9
-      real*8 r,pi
-      integer :: nx,ny,nz
-      integer i,j,k,ier,nquads,nocts,order
-      integer npoints,handle,ierror
-      
-      conformal_state = CONFORMAL_METRIC     
-      
-      pi = 4.0d0*atan(1.0d0)
-      
-c     Set up the grid spacings
-      nx = cctk_lsh(1)
-      ny = cctk_lsh(2)
-      nz = cctk_lsh(3)      
-      
-c     Distorted Schwarzchild BH parameters
-      amp = distortedbh_amp
-      eta0 = distortedbh_eta0
-      c = distortedbh_c
-      sigma = distortedbh_sigma
-      n = distortedbh_n
-      etamax = distortedbh_etamax
-      
-      print *,"Brill wave + Distorted BH solve"
-      write(*,123)amp,eta0,c,sigma,n
-      print*,'etamax=',etamax
- 123  format(1x, 'Pars: amp',f8.5,' eta0',f8.5,' c',f8.5,' sigma',f8.5,' n',i3)
-      
-c     Sovle on this sized cartesian grid
-c     3D grid size NE x NT x NP
-c     Add 2 zones for eta coordinate and 4 for theta
-c     and phi coordenate.
-      ne = distortedbh_ne + 2
-      nq = distortedbh_nq + 2
-      np = distortedbh_np + 2
-c     
-      allocate(ac(ne,nq,np),ae(ne,nq,np),aw(ne,nq,np),an(ne,nq,np),
-     $     as(ne,nq,np),aq(ne,nq,np),ab(ne,nq,np),rhs(ne,nq,np),
-     $     qf(ne,nq,np),qfetaeta(ne,nq,np),qfqq(ne,nq,np),
-     $     qfphi(ne,nq,np),qfphiphi(ne,nq,np),
-     $     psisph(ne,nq,np),psiprim(ne,nq,np),detapsisph(ne,nq,np),
-     $     dqpsisph(ne,nq,np),dphipsisph(ne,nq,np),
-     $     detaetapsisph(ne,nq,np),detaqpsisph(ne,nq,np),
-     $     detaphipsisph(ne,nq,np),dqqpsisph(ne,nq,np),
-     $     dqphipsisph(ne,nq,np),dphiphipsisph(ne,nq,np),
-     $     psi3d(nx,ny,nz),detapsi3d(nx,ny,nz),dqpsi3d(nx,ny,nz),
-     $     dphipsi3d(nx,ny,nz),detaetapsi3d(nx,ny,nz),
-     $     detaqpsi3d(nx,ny,nz),detaphipsi3d(nx,ny,nz),
-     $     dqqpsi3d(nx,ny,nz),dqphipsi3d(nx,ny,nz),
-     $     dphiphipsi3d(nx,ny,nz))
-      allocate(etagrd(ne),qgrd(nq),phigrd(np))
-      allocate(eta(nx,ny,nz),abseta(nx,ny,nz),sign_eta(nx,ny,nz),
-     $     q(nx,ny,nz),phi(nx,ny,nz))
-c     
-c     Initialize some array
-c     
-      nocts = 4
-      nquads = 2      
-      dphi = nocts*0.5*pi/(np-2)
-      dq = nquads*0.5*pi/(nq-2)
-      deta = etamax/(ne-3)
-      
-      do k = 1,np
-         phigrd(k) = (k-1.5)*dphi
-         do j = 1,nq
-            qgrd(j) = (j-1.5)*dq
-            do i=1,ne
-               etagrd(i) = (i-2)*deta
-            enddo
-         enddo
-      enddo
-c     
-c Initialize q-function and its derivatives: should be generalized
-c     
-      do k = 1,np
-         do j = 1,nq
-            do i = 1,ne
-#include "CactusEinstein/DistortedBHIVP/src/qfunc.x"
-            enddo
-         enddo
-      enddo
-      
-c     
-c     Initialize psi to the Schwarzschild solution:
-c     
-      do k = 1,np
-         do j = 1,nq
-            do i = 1,ne
-               psisph(i,j,k) = 2.*cosh(0.5*etagrd(i))
-               psiprim(i,j,k) = 0.
-            enddo
-         enddo
-      enddo
-c     
-c     Initialize stencil coefficients:
-c     
-      do k = 1,np
-         do j = 1,nq
-            do i = 1,ne
-               ac(i,j,k) = 0.
-               ae(i,j,k) = 0.
-               aw(i,j,k) = 0.
-               an(i,j,k) = 0.
-               as(i,j,k) = 0.
-               aq(i,j,k) = 0.
-               ab(i,j,k) = 0.
-               rhs(i,j,k) = 0.
-            enddo
-         enddo
-      enddo
-c     
-      do k = 2,np-1
-         do j = 2,nq-1
-            do i = 2,ne-1
-               ac(i,j,k) = -2./deta**2-2./dq**2-2.*exp(2.*
-     &              qf(i,j,k))/(dphi**2*
-     &              sin(qgrd(j))**2)+0.25*
-     &              (qfetaeta(i,j,k)+qfqq(i,j,k)+2.*
-     &              exp(2.*qf(i,j,k))*qfphiphi(i,j,k)/
-     &              sin(qgrd(j))**2+3.*exp(2.*
-     &              qf(i,j,k))*qfphi(i,j,k)**2/
-     &              sin(qgrd(j))**2-1.)
-               
-               ae(i,j,k) = 1./deta**2
-               
-               aw(i,j,k) = 1./deta**2
-               
-               an(i,j,k) = 1./dq**2+0.5/(dq*tan(qgrd(j)))
-               
-               as(i,j,k) = 1./dq**2-0.5/(dq*tan(qgrd(j)))
-               
-               aq(i,j,k) = exp(2.*qf(i,j,k))/(dphi**2*
-     &              sin(qgrd(j))**2)+exp(2.*
-     &              qf(i,j,k))*qfphi(i,j,k)/(dphi*
-     &              sin(qgrd(j))**2)
-               
-               ab(i,j,k) = exp(2.*qf(i,j,k))/(dphi**2*
-     &              sin(qgrd(j))**2)-exp(2.*
-     &              qf(i,j,k))*qfphi(i,j,k)/(dphi*
-     &              sin(qgrd(j))**2)
-               
-               rhs(i,j,k) = -0.25*(qfetaeta(i,j,k)+
-     &              qfqq(i,j,k)+2.*exp(2.*qf(i,j,k))*
-     &              qfphiphi(i,j,k)/sin(qgrd(j))**2+3.*
-     &              exp(2.*qf(i,j,k))*qfphi(i,j,k)**2/
-     &              sin(qgrd(j))**2)*psisph(i,j,k)
-            enddo
-         enddo
-      enddo
-c
-c Apply boundary conditions to the faces of the cube:
-c
-c    i=2:
-      do k = 3,np-2
-         do j = 3,nq-2
-            ae(2,j,k) = ae(2,j,k) + aw(2,j,k)
-            aw(2,j,k) = 0.
-c
-c    i=ne-1:
-            ac(ne-1,j,k) = ac(ne-1,j,k)+4.*ae(ne-1,j,k)/(3.+deta)
-            aw(ne-1,j,k) = aw(ne-1,j,k) - ae(ne-1,j,k)/(3.+deta)
-            ae(ne-1,j,k) = 0.
-         enddo
-      enddo
-c
-c    j=2:
-      do k = 3,np-2
-         do i = 3,ne-2
-            ac(i,2,k) = ac(i,2,k) + as(i,2,k)
-            as(i,2,k) = 0.
-c
-c    j=nq-1:
-            ac(i,nq-1,k) = ac(i,nq-1,k) + an(i,nq-1,k)
-            an(i,nq-1,k) = 0.
-         enddo
-      enddo
-c
-c    k=2:
-      do j = 3,nq-2
-         do i = 3,ne-2
-            ac(i,j,2) = ac(i,j,2) + ab(i,j,2)
-            ab(i,j,2) = 0.
-c
-c    k=np-1:
-            ac(i,j,np-1) = ac(i,j,np-1) + aq(i,j,np-1)
-            aq(i,j,np-1) = 0.
-         enddo
-      enddo
-c
-c Apply boundary conditions to the edges of the cube:
-c
-c    i=2, j=2:
-      do k = 3,np-2
-         ae(2,2,k) = ae(2,2,k) + aw(2,2,k)
-         ac(2,2,k) = ac(2,2,k) + as(2,2,k)
-         aw(2,2,k) = 0.
-         as(2,2,k) = 0.
-c
-c    i=ne-1, j=2:
-         aw(ne-1,2,k) = aw(ne-1,2,k) - ae(ne-1,2,k)/(3.+deta)
-         ac(ne-1,2,k) = ac(ne-1,2,k) + as(ne-1,2,k) +
-     &        4.*ae(ne-1,2,k)/(3.+deta)
-         ae(ne-1,2,k) = 0.
-         as(ne-1,2,k) = 0.
-c
-c    i=2, j=nq-1:
-         ae(2,nq-1,k) = ae(2,nq-1,k) + aw(2,nq-1,k)
-         ac(2,nq-1,k) = ac(2,nq-1,k) + an(2,nq-1,k)
-         aw(2,nq-1,k) = 0.
-         an(2,nq-1,k) = 0.
-c
-c    i=ne-1, j=nq-1:
-         aw(ne-1,nq-1,k) = aw(ne-1,nq-1,k) - ae(ne-1,nq-1,k)/
-     &        (3.+deta)
-         ac(ne-1,nq-1,k) = ac(ne-1,nq-1,k) + an(ne-1,nq-1,k) +
-     &        4.*ae(ne-1,nq-1,k)/(3.+deta)
-         ae(ne-1,nq-1,k) = 0.
-         an(ne-1,nq-1,k) = 0.
-      enddo
-c
-c    i=2, k=2:
-      do j = 3,nq-2
-         ae(2,j,2) = ae(2,j,2) + aw(2,j,2)
-         ac(2,j,2) = ac(2,j,2) + ab(2,j,2)
-         aw(2,j,2) = 0.
-         ab(2,j,2) = 0.
-c
-c    i=ne-1, k=2:
-         aw(ne-1,j,2) = aw(ne-1,j,2) - ae(ne-1,j,2)/(3.+deta)
-         ac(ne-1,j,2) = ac(ne-1,j,2) + ab(ne-1,j,2) +
-     &        4.*ae(ne-1,j,2)/(3.+deta)
-         ae(ne-1,j,2) = 0.
-         ab(ne-1,j,2) = 0.
-c
-c    i=2, k=np-1:
-         ae(2,j,np-1) = ae(2,j,np-1) + aw(2,j,np-1)
-         ac(2,j,np-1) = ac(2,j,np-1) + aq(2,j,np-1)
-         aw(2,j,np-1) = 0.
-         aq(2,j,np-1) = 0.
-c
-c    i=ne-1, k=np-1:
-         aw(ne-1,j,np-1) = aw(ne-1,j,np-1) - ae(ne-1,j,np-1)/
-     &        (3.+deta)
-         ac(ne-1,j,np-1) = ac(ne-1,j,np-1) + aq(ne-1,j,np-1) +
-     &        4.*ae(ne-1,j,np-1)/(3.+deta)
-         ae(ne-1,j,np-1) = 0.
-         aq(ne-1,j,np-1) = 0.
-      enddo
-c
-c    j=2, k=2:
-      do i = 3,ne-2
-         ac(i,2,2) = ac(i,2,2) + as(i,2,2) + ab(i,2,2)
-         as(i,2,2) = 0.
-         ab(i,2,2) = 0.
-c
-c    j=nq-1, k=2:
-         ac(i,nq-1,2) = ac(i,nq-1,2) + an(i,nq-1,2) +
-     &        ab(i,nq-1,2)
-         an(i,nq-1,2) = 0.
-         ab(i,nq-1,2) = 0.
-c
-c    j=2, k=np-1:
-         ac(i,2,np-1) = ac(i,2,np-1) + as(i,2,np-1) +
-     &        aq(i,2,np-1)
-         as(i,2,np-1) = 0.
-         aq(i,2,np-1) = 0.
-c
-c    j=nq-1, k=np-1:
-         ac(i,nq-1,np-1) = ac(i,nq-1,np-1) + an(i,nq-1,np-1) +
-     &        aq(i,nq-1,np-1)
-         an(i,nq-1,np-1) = 0.
-         aq(i,nq-1,np-1) = 0.
-      enddo
-c
-c Apply boundary conditions to the corners of the cube:
-c
-c    i=2, j=2, k=2:
-      ae(2,2,2) = ae(2,2,2) + aw(2,2,2)
-      ac(2,2,2) = ac(2,2,2) + as(2,2,2) + ab(2,2,2)
-      aw(2,2,2) = 0.
-      as(2,2,2) = 0.
-      ab(2,2,2) = 0.
-c
-c    i=ne-1, j=2, k=2:
-      aw(ne-1,2,2) = aw(ne-1,2,2) - ae(ne-1,2,2)/(3.+deta)
-      ac(ne-1,2,2) = ac(ne-1,2,2) + as(ne-1,2,2) + ab(ne-1,2,2) +
-     &     4.*ae(ne-1,2,2)/(3.+deta)
-      ae(ne-1,2,2) = 0.
-      as(ne-1,2,2) = 0.
-      ab(ne-1,2,2) = 0.
-c
-c    i=2, j=nq-1, k=2:
-      ae(2,nq-1,2) = ae(2,nq-1,2) + aw(2,nq-1,2)
-      ac(2,nq-1,2) = ac(2,nq-1,2) + an(2,nq-1,2) + ab(2,nq-1,2)
-      aw(2,nq-1,2) = 0.
-      an(2,nq-1,2) = 0.
-      ab(2,nq-1,2) = 0.
-c
-c    i=2, j=2, k=np-1:
-      ae(2,2,np-1) = ae(2,2,np-1) + aw(2,2,np-1)
-      ac(2,2,np-1) = ac(2,2,np-1) + as(2,2,np-1) + aq(2,2,np-1)
-      aw(2,2,np-1) = 0.
-      as(2,2,np-1) = 0.
-      aq(2,2,np-1) = 0.
-c
-c    i=ne-1, j=nq-1, k=2:
-      aw(ne-1,nq-1,2) = aw(ne-1,nq-1,2) - ae(ne-1,nq-1,2)/(3.+deta)
-      ac(ne-1,nq-1,2) = ac(ne-1,nq-1,2) + an(ne-1,nq-1,2) + ab(ne-1,nq-1,2) +
-     &     4.*ae(ne-1,nq-1,2)/(3.+deta)
-      ae(ne-1,nq-1,2) = 0.
-      an(ne-1,nq-1,2) = 0.
-      ab(ne-1,nq-1,2) = 0.
-c
-c    i=ne-1, j=2, k=np-1:
-      aw(ne-1,2,np-1) = aw(ne-1,2,np-1) - ae(ne-1,2,np-1)/(3.+deta)
-      ac(ne-1,2,np-1) = ac(ne-1,2,np-1) + as(ne-1,2,np-1) + aq(ne-1,2,np-1) +
-     &     4.*ae(ne-1,2,np-1)/(3.+deta)
-      ae(ne-1,2,np-1) = 0.
-      as(ne-1,2,np-1) = 0.
-      aq(ne-1,2,np-1) = 0.
-c
-c    i=2, j=nq-1, k=np-1:
-      ae(2,nq-1,np-1) = ae(2,nq-1,np-1) + aw(2,nq-1,np-1)
-      ac(2,nq-1,np-1) = ac(2,nq-1,np-1) + an(2,nq-1,np-1) + aq(2,nq-1,np-1)
-      aw(2,nq-1,np-1) = 0.
-      an(2,nq-1,np-1) = 0.
-      aq(2,nq-1,np-1) = 0.
-c
-c    i=ne-1, j=nq-1, k=np-1:
-      aw(ne-1,nq-1,np-1) = aw(ne-1,nq-1,np-1) - ae(ne-1,nq-1,np-1)/(3.+deta)
-      ac(ne-1,nq-1,np-1) = ac(ne-1,nq-1,np-1) + an(ne-1,nq-1,np-1) +
-     &     aq(ne-1,nq-1,np-1) + 4.*ae(ne-1,nq-1,np-1)/(3.+deta)
-      ae(ne-1,nq-1,np-1) = 0.
-      an(ne-1,nq-1,np-1) = 0.
-      aq(ne-1,nq-1,np-1) = 0.
-c
-c Solve for psi:
-c
-      call bicgst3d(ac,ae,aw,an,as,aq,ab,psiprim,rhs,dbh_eps,rmax,ier,
-     $     ne,nq,np)
-c
-      if (rmax.gt.dbh_eps) then
-         write(*,*) '***WARNING: bicgst3d did not converge.'
-      endif
-      if (ier.eq.-1) then
-         write(*,*) '***WARNING: ier=-1'
-      endif
-      print *,'psiprim = ',maxval(psiprim),' ',minval(psiprim)
-c
-c Now, apply boundary conditions to psiprim:
-c
-      do k = 1,np
-         do j = 1,nq
-            psiprim(1,j,k) = psiprim(3,j,k)
-            psiprim(ne,j,k) = (4.*psiprim(ne-1,j,k)-psiprim(ne-2,j,k))/
-     $           (3.+deta)
-         enddo
-      enddo
-      do k = 1,np
-         do i = 1,ne
-            psiprim(i,1,k) = psiprim(i,2,k)
-            psiprim(i,nq,k) = psiprim(i,nq-1,k)
-         enddo
-      enddo
-      do j = 1,nq
-         do i = 1,ne
-            psiprim(i,j,1) = psiprim(i,j,2)
-            psiprim(i,j,np) = psiprim(i,j,np-1)
-         enddo
-      enddo
-c     
-c     Here, compute the derivatives of the spherical conformal factor
-c
-c      goto 110 
-
-      do k = 1,np
-         do j = 1,nq
-            do i = 2,ne-1
-               detapsisph(i,j,k)=0.5*(psiprim(i+1,j,k)-psiprim(i-1,j,k))
-     $              /deta + sinh(0.5*etagrd(i))
-            enddo
-            detapsisph(1,j,k) = -detapsisph(3,j,k)
-         enddo
-      enddo
-c     
-      do k = 1,np
-         do j = 2,nq-1
-            do i = 1,ne
-               dqpsisph(i,j,k)=0.5*(psiprim(i,j+1,k)-psiprim(i,j-1,k))/
-     $              dq
-            enddo
-         enddo
-         do i = 1,ne
-            dqpsisph(i,1,k) = -dqpsisph(i,2,k)
-            dqpsisph(i,nq,k) = -dqpsisph(i,nq-1,k)
-         enddo
-      enddo
-c
-      do k = 2,np-1
-         do j = 1,nq
-            do i = 1,ne
-               dphipsisph(i,j,k)=0.5*(psiprim(i,j,k+1)-psiprim(i,j,k-1))
-     $              /dphi
-            enddo
-         enddo
-      enddo
-      do j = 1,nq
-         do i = 1,ne
-            dphipsisph(i,j,1) = -dphipsisph(i,j,2)
-            dphipsisph(i,j,np) = -dphipsisph(i,j,np-1)
-         enddo
-      enddo
-c     
-      do k = 1,np
-         do j = 1,nq
-            do i = 2,ne-1
-               detaetapsisph(i,j,k)=(psiprim(i+1,j,k)-2.*psiprim(i,j,k)+
-     &              psiprim(i-1,j,k))/deta**2+sqrt(0.25)*
-     &              cosh(0.5*etagrd(i))
-            enddo
-            detaetapsisph(1,j,k) = detaetapsisph(3,j,k)
-         enddo
-      enddo
-c     
-      do k = 1,np
-         do j = 2,nq-1
-            do i = 1,ne
-               detaqpsisph(i,j,k)=0.5*(detapsisph(i,j+1,k)-detapsisph(i,
-     $              j-1,k))/dq
-            enddo
-         enddo
-         do i = 1,ne
-            detaqpsisph(i,1,k) = -detaqpsisph(i,2,k)
-            detaqpsisph(i,nq,k) = -detaqpsisph(i,nq-1,k)
-         enddo
-      enddo
-c
-      do k = 2,np-1
-         do j = 1,nq
-            do i = 1,ne
-               detaphipsisph(i,j,k)=0.5*(detapsisph(i,j,k+1)-detapsisph(
-     $              i,j,k-1))/dphi
-            enddo
-         enddo
-      enddo
-      do j = 1,nq
-         do i = 1,ne
-            detaphipsisph(i,j,1) = -detaphipsisph(i,j,2)
-            detaphipsisph(i,j,np) = -detaphipsisph(i,j,np-1)
-         enddo
-      enddo
-c
-      do k = 1,np
-         do j = 2,nq-1
-            do i = 1,ne
-               dqqpsisph(i,j,k)=0.5*(dqpsisph(i,j+1,k)-dqpsisph(i,j-1,k))/
-     $              dq
-            enddo
-         enddo
-         do i = 1,ne
-            dqqpsisph(i,1,k) = dqqpsisph(i,2,k)
-            dqqpsisph(i,nq,k) = dqqpsisph(i,nq-1,k)
-         enddo
-      enddo
-c
-      do k = 2,np-1
-         do j = 1,nq
-            do i = 1,ne
-               dqphipsisph(i,j,k)=0.5*(dqpsisph(i,j,k+1)-dqpsisph(i,j,k-1))/
-     $              dphi
-            enddo
-         enddo
-      enddo
-      do j = 1,nq
-         do i = 1,ne
-            dqphipsisph(i,j,1) = -dqphipsisph(i,j,2)
-            dqphipsisph(i,j,np) = -dqphipsisph(i,j,np-1)
-         enddo
-      enddo
-c     
-      do k = 2,np-1
-         do j = 1,nq
-            do i = 1,ne
-               dphiphipsisph(i,j,k)=0.5*(dphipsisph(i,j,k+1)-
-     $              dphipsisph(i,j,k-1))/dphi
-            enddo
-         enddo
-      enddo
-      do j = 1,nq
-         do i = 1,ne
-            dphiphipsisph(i,j,1) = dphiphipsisph(i,j,2)
-            dphiphipsisph(i,j,np) = dphiphipsisph(i,j,np-1)
-         enddo
-      enddo
-c     
-      do k = 1,np
-         do j = 1,nq
-            psisph(:,j,k)=psiprim(:,j,k)+2.0*cosh(0.5*etagrd)
-         enddo
-      enddo
-      
-c     
-c     Now compute on the Cartesian coordinate.
-c     
-c     Compute eta,q,phi at the each points of cartesian grid
-      do k = 1,nz
-         do j = 1,ny
-            do i = 1,nx
-               eta(i,j,k) = 0.5d0*dlog(x(i,j,k)**2+y(i,j,k)**2+z(i,j,k)**2)
-               abseta(i,j,k) = abs(eta(i,j,k))
-               if(eta(i,j,k) .lt. 0)then
-                  sign_eta(i,j,k) = -1.0d0
-               else
-                  sign_eta(i,j,k) = 1.0d0
-               endif
-               q(i,j,k) = datan2(sqrt(x(i,j,k)**2+y(i,j,k)**2),z(i,j,k))
-               phi(i,j,k) = datan2(y(i,j,k),x(i,j,k))
-            enddo
-         enddo
-      enddo
-      
-      call CCTK_GetInterpHandle (handle, "simple_local")
-      
-      npoints = nx*ny*nz
-
-      call CCTK_Interp (cctkGH,ierror,handle,npoints,3,10,10,
-     $     ne,nq,np,abseta,q,phi,
-     $     CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,
-     $     etagrd(1),qgrd(1),phigrd(1)-pi,deta,dq,dphi,
-     $     psisph,detapsisph,dqpsisph,dphipsisph,detaetapsisph,
-     $     detaqpsisph,detaphipsisph,dqqpsisph,dqphipsisph,dphiphipsisph,
-     $     CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,
-     $     CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,
-     $     CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,
-     $     CCTK_VARIABLE_REAL,
-     $     psi3d,detapsi3d,dqpsi3d,dphipsi3d,detaetapsi3d,detaqpsi3d,
-     $     detaphipsi3d,dqqpsi3d,dqphipsi3d,dphiphipsi3d,
-     $     CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,
-     $     CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,
-     $     CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,
-     $     CCTK_VARIABLE_REAL)
-
-      do k=1,nz
-         do j=1,ny
-            do i=1,nx
-               psi(i,j,k) = psi3d(i,j,k)*exp(-0.5*eta(i,j,k))
-               
-               detapsi3d(i,j,k) = sign_eta(i,j,k)*detapsi3d(i,j,k)
-               
-c     psix = \partial psi / \partial x / psi
-#include "CactusEinstein/DistortedBHIVP/src/psi_1st_deriv.x"
-               
-               detaqpsi3d(i,j,k) = sign_eta(i,j,k)*detaqpsi3d(i,j,k)
-               detaphipsi3d(i,j,k) = sign_eta(i,j,k)*detaphipsi3d(i,j,k)
-
-c     psixx = \partial^2\psi / \partial x^2 / psi
-#include "CactusEinstein/DistortedBHIVP/src/psi_2nd_deriv.x"
-            enddo
-         enddo
-      enddo
-      
-c     Conformal metric
-c     gxx = ...
-      
-c     Derivatives of the metric
-c     dxgxx = 1/2 \partial gxx / \partial x
-c     
-      do k=1,nz
-         do j=1,ny
-            do i=1,nx
-               eta(i,j,k) = 0.5d0*dlog(x(i,j,k)**2+y(i,j,k)**2+z(i,j,k)**2)
-               q(i,j,k) = datan2(sqrt(x(i,j,k)**2+y(i,j,k)**2),z(i,j,k))
-               phi(i,j,k) = datan2(y(i,j,k),x(i,j,k))
-#include "CactusEinstein/DistortedBHIVP/src/gij.x"
-            enddo
-         enddo
-      enddo
-      
-c     Courvature
-      kxx = 0.0d0
-      kxy = 0.0d0
-      kxz = 0.0d0
-      kyy = 0.0d0
-      kyz = 0.0d0
-      kzz = 0.0d0
-      
- 110  continue
-c     Set ADM mass
-      i = ne-15
-      adm = 0.0
-      do k=2,np-1
-         do j=2,nq-1
-            adm=adm+(psisph(i,j,k)-(psisph(i+1,j,k)-psisph(i-1,j,k))/
-     $           deta)*exp(0.5*etagrd(i))
-         enddo
-      enddo
-      adm=adm/(nq-2)/(np-2)
-      print *,'ADM mass: ',adm
-      if (CCTK_Equals(initial_lapse,"schwarz")==1) then 
-         write (*,*)"Initial with schwarzschild-like lapse"
-         write (*,*)"using alp = (2.*r - adm)/(2.*r+adm)."
-         alp = (2.*r - adm)/(2.*r+adm)
-      endif
-      
-      deallocate(ac,ae,aw,an,as,aq,ab,rhs,
-     $     qf,qfetaeta,qfqq,qfphi,qfphiphi,
-     $     psisph,psiprim,detapsisph,dqpsisph,dphipsisph,detaetapsisph,
-     $     detaqpsisph,detaphipsisph,dqqpsisph,dqphipsisph,
-     $     dphiphipsisph,
-     $     psi3d,detapsi3d,dqpsi3d,
-     $     dphipsi3d,detaetapsi3d,detaqpsi3d,
-     $     detaphipsi3d,dqqpsi3d,dqphipsi3d,
-     $     dphiphipsi3d)
-      deallocate(etagrd,qgrd,phigrd)
-      deallocate(eta,abseta,sign_eta,q,phi)
-      
-      return
-      end
diff --git a/src/interp3.F b/src/interp3.F
deleted file mode 100644
index 0bbdf90..0000000
--- a/src/interp3.F
+++ /dev/null
@@ -1,264 +0,0 @@
-c /*@@
-c   @routine    interp3
-c   @date       Fri Feb 14 08:46:53 1997
-c   @author     Ryoji Takahashi
-c   @desc 
-c      Interpolates from 3D data var with coordinates x, y, and z and
-c      sizes nx , ny, and nz onto 1D data out with position outx, outy
-c      ,and outz nout points.
-c      <p>
-c      This has linear, quadratic and cubic interpolators in it.
-c      Or will one day very soon.
-c   @enddesc 
-c   @calls     
-c   @calledby   numerical_nonaxi
-c@@*/
-
-      subroutine interp3(dat,x,y,z,nx,ny,nz,out,outx,outy,outz,nout
-     $     ,order)
-      implicit none 
-      integer nx,ny,nz,nout
-      real*8 dat(nx,ny,nz), x(nx), y(ny), z(nz)
-      real*8 out(nout),outx(nout),outy(nout),outz(nout)
-      integer order
-c Interpolation goes from ibelow to ibelow+[1,2,3] depending on order
-      integer i,j,k,l,m,ibelow,jbelow,kbelow,pt
-      real*8 xsym,ysym,zsym,findx,findy,findz,frac
-      real*8 ydir(order+1)
-      real*8 zdir(order+1)
-      real*8 f(4), fp(4,4), fl(4)
-      real*8 dbh_linear, dbh_quad, dbh_cubic
-      real*8 dx, dy, dz, PI
-      integer twobhjsad
-
-      pi = 4.0d0*atan(1.0d0)
-      
-      dx = x(2) - x(1)
-      dy = y(2) - y(1)
-      dz = z(2) - z(1)
-      
-c Loop over all out points
-      do pt=1,nout
-         zsym = 1.0D0
-         ysym = 1.0D0
-         xsym = 1.0D0
-c Check bounds
-         findx = outx(pt)
-         if (findx .lt. x(1)) then 
-            write (*,*) "Below inner bound at ",pt,outx(pt)
-            STOP
-         endif
-         if (findx .gt. x(nx)) then
-            write (*,*) "Above x bounds at ",pt,outx(pt),x(nx)
-            STOP
-         endif
-         findy = outy(pt)
-         if (findy .lt. y(1)) then 
-            write (*,*) "Below  y inner bound at ",pt,outy(pt)
-            STOP
-         endif
-         if (findy .gt. y(ny)) then 
-            write (*,*) "Below  y inner bound at ",pt,outy(pt)
-            STOP
-         endif
-         findz = outz(pt)+pi
-         if (findz .lt. z(1)) then 
-            write (*,*) "Below  z inner bound at ",pt,outz(pt)
-            STOP
-         endif
-         if (findz .gt. z(nz)) then 
-            write (*,*) "Below  z inner bound at ",pt,outz(pt)
-            STOP
-         endif
-
-c Locate ourselves in i,j space
-c         do i=1,nx
-c            if (x(i) .lt. findx) then
-c               ibelow = i
-c            endif
-c         enddo
-c     Assume a regular grid
-         ibelow = (findx-x(1))/dx+1
-
-         if (order .eq. 3 .and. ibelow .gt. 1) then
-            ibelow = ibelow - 1
-         endif
-         if (ibelow + order .gt. nx) then
-            ibelow = nx - order
-         endif
-
-c         do i=1,ny
-c            if (y(i) .lt. findy) then
-c     jbelow = i
-c     endif
-c         enddo
-c     Assume a regular grid
-         jbelow = (findy-y(1))/dy+1
-         if (order .eq. 3 .and. jbelow .gt. 1) then
-            jbelow = jbelow - 1
-         endif
-         if (jbelow + order .gt. ny) then
-            jbelow = ny - order
-         endif
-
-c         do i=1,nz
-c            if (z(i) .lt. findz) then
-c               kbelow = i
-c            endif
-c         enddo
-c     Assume a regular grid
-         kbelow = (findz-z(1))/dz+1
-         if (order .eq. 3 .and. kbelow .gt. 1) then
-            kbelow = kbelow - 1
-         endif
-         if (kbelow + order .gt. nz) then
-            kbelow = nz - order
-         endif
-         
-c     write (*,*) "PT :",findx,findy
-c     write (*,*) "SYM:",sym
-c     write (*,*) "BOUND X ",ibelow,x(ibelow),x(ibelow+1)
-c     write (*,*) "BOUND Y ",jbelow,y(jbelow),y(jbelow+1)
-c     write (*,*) "BOUND z ",kbelow,z(kbelow),z(kbelow+1)
-         
-c     So do the interpolation
-         if (order .eq. 1) then 
-c     Interp in the x direction
-!            frac = (findx-x(ibelow))/(x(ibelow+1)-x(ibelow))
-            do l = 1,2
-               do m = 1,2
-                  f(1) = dat(ibelow,jbelow+l-1,kbelow+m-1)
-                  f(2) = dat(ibelow+1,jbelow+l-1,kbelow+m-1)
-                  fp(l,m) = dbh_linear(f,x(ibelow),dx,findx)
-               enddo
-            enddo
-c     Now take our 2x2 plane and interp to the center of both
-c     in the y direction
-!            frac = (findy-y(jbelow))/(y(jbelow+1)-y(jbelow))
-            do m = 1,2
-               f(1) = fp(1,m)
-               f(2) = fp(2,m)
-               fl(m) = dbh_linear(f,y(jbelow),dy,findy)
-            enddo
-c     And finally, interp in the z direction
-            out(pt) = xsym * ysym * zsym * 
-     $           dbh_linear(fl,z(kbelow),dz,findz)
-         else if (order .eq. 2) then
-c     Load up for calls to poly2inter
-            do l=1,3
-               do m=1,3
-                  f(1) = dat(ibelow,jbelow+l-1,kbelow+m-1)
-                  f(2) = dat(ibelow+1,jbelow+l-1,kbelow+m-1)
-                  f(3) = dat(ibelow+2,jbelow+l-1,kbelow+m-1)
-                  fp(l,m) = dbh_quad(f,x(ibelow),dx,findx)
-               enddo
-            enddo
-c     Now take our 2x2 plane and interp to the center of both
-c     in the y direction
-            do m=1,3
-               f(1) = fp(1,m)
-               f(2) = fp(2,m)
-               f(3) = fp(3,m)
-               fl(m) = dbh_quad(f,y(jbelow),dy,findy)
-            enddo
-c     And finally, interp in the z direction
-            out(pt) = xsym * ysym * zsym *
-     $           dbh_quad(fl,z(kbelow),dz,findz)
-         else if (order .eq. 3) then
-c     Load up for calls to cubic
-            do l=1,4
-               do m=1,4
-                  f(1) = dat(ibelow,jbelow+l-1,kbelow+m-1)
-                  f(2) = dat(ibelow+1,jbelow+l-1,kbelow+m-1)
-                  f(3) = dat(ibelow+2,jbelow+l-1,kbelow+m-1)
-                  f(4) = dat(ibelow+3,jbelow+l-1,kbelow+m-1)
-                  fp(l,m) = dbh_cubic(f,x(ibelow),dx,findx)
-               enddo
-            enddo
-c     Now take our 2x2 plane and interp to the center of both
-c     in the y direction
-            do m=1,4
-               f(1) = fp(1,m)
-               f(2) = fp(2,m)
-               f(3) = fp(3,m)
-               f(4) = fp(4,m)
-               fl(m) = dbh_cubic(f,y(jbelow),dy,findz)
-            enddo
-c     And finally, interp in the z direction
-            out(pt) = xsym * ysym * zsym * 
-     $           dbh_cubic(fl,z(kbelow),dz,findz)
-         else
-            write (*,*) "ORDER set wrong in interp3d",order
-            stop
-         endif
-
-      enddo
-      return
-      end
-
-      real*8 function dbh_linear(f, x0, dx, findx)
-      implicit none
-      real*8 f(2),x0,dx,findx
-      real*8 frac
-
-      frac = (findx-x0)/dx
-      dbh_linear = (frac)*f(2) + (1.0-frac)*f(1)
-
-      return
-      end
-
-      real*8 function dbh_quad(f, x0, dx, findx)
-      implicit none
-      real*8 f(3),x0, dx, findx, dbh_quad
-      real*8 f0,f1,f2
-      real*8 a,b,c, dx2, x02, o2dx2
-c Mathematica tells us
-c     -  List(List(Rule(c,(2*dx**2*f0 + 3*dx*f0*x0 - 4*dx*f1*x0 
-c     -        + dx*f2*x0 +  f0*x0**2 - 2*f1*x0**2 + f2*x0**2)
-c     -        /(2*dx**2)), Rule(b,(-3*dx*f0 + 4*dx*f1 - dx*f2 - 2*f0*x0 
-c     -        + 4*f1*x0 - 2*f2*x0)/(2*dx**2)),Rule(a,(f0 - 2*f1 + 
-c     -        f2)/(2*dx**2))))
-
-	f0 = f(1)
-	f1 = f(2)
-	f2 = f(3)
-        dx2 = dx**2
-        x02 = x0**2
-        o2dx2 = 1.0D0/(2.0D0*dx2)
-
-      c = (2.0D0*dx2*f0 + dx*x0*(3.0D0*f0 - 4.0D0*f1 + f2) + 
-     $       x02*(f0 - 2.0D0*f1 + f2))*o2dx2
-      b = (dx * (-3.0D0*f0 + 4.0D0*f1 - f2) + x0 * (- 2.0D0*f0 + 
-     $     4.0D0*f1 - 2.0D0*f2))*o2dx2
-
-      a = (f0 - 2.0D0*f1 + f2)*o2dx2
-
-      dbh_quad = a*findx**2 + b*findx + c
-
-      end
-
-
-      real*8 function dbh_cubic(f, x0, dx, findx)
-      implicit none
-      real*8 a,b,c,d,dbh_cubic
-      real*8 f(4),x0,dx,findx
-
-      a = -(f(1)-3.0*f(2)+3.0*f(3)-f(4)) / (6.0*(dx**3))
-      
-      b = (f(1)-2.0*f(2)+f(3))/(2.0*(dx**2)) +
-     $     (f(1)-3.0*f(2)+3.0*f(3)-f(4))*(dx+x0)/(2.0*(dx**3))
-      
-      c = ((dx**2)*(-11.0*f(1) + 18.0*f(2) -  9.0*f(3) +  2.0*f(4)) +
-     $     dx*x0*  (-12.0*f(1) + 30.0*f(2) - 24.0*f(3) +  6.0*f(4)) +
-     $     (x0**2)*( -3.0*f(1) +  9.0*f(2) -  9.0*f(3) +  3.0*f(4))) /
-     $     (6.0*(dx**3))
-      
-      d = ((dx**3)*   (  6.0*f(1) ) +
-     $     (dx**2)*x0*( 11.0*f(1) - 18.0*f(2) +  9.0*f(3) -  2.0*f(4)) +
-     $     (x0**2)*dx*(  6.0*f(1) - 15.0*f(2) + 12.0*f(3) -  3.0*f(4)) +
-     $     (x0**3)*   (  1.0*f(1) -  3.0*f(2) +  3.0*f(3) -  1.0*f(4)))/
-     $     (6.0*(dx**3))
-         
-      dbh_cubic = ((a*findx + b)*findx + c)*findx + d
-      return 
-      end