This commit was generated by cvs2svn to compensate for changes in r2, which

included commits to RCS files with non-trunk default branches. git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/Extract/trunk@3 5301f0c2-dbc4-4cee-b2f5-8d7afba4d129
author: allen <allen@5301f0c2-dbc4-4cee-b2f5-8d7afba4d129> 1999-11-01 14:56:51 +0000
committer: allen <allen@5301f0c2-dbc4-4cee-b2f5-8d7afba4d129> 1999-11-01 14:56:51 +0000
commit: 661c56c847f4334feb53fd8d94551b2ab1efdb9a (patch)
tree: 3dbdeebf2c0b8f0f0fe30218e9a58fc0cd1638a1 /src
parent: d8cf32e113a96b24744d231013e3d0b81bf55c30 (diff)
22 files changed, 3686 insertions, 0 deletions
diff --git a/src/ADMmass_integrand3D.F b/src/ADMmass_integrand3D.F
new file mode 100644
index 0000000..06b1198
--- /dev/null
+++ b/src/ADMmass_integrand3D.F
@@ -0,0 +1,170 @@
+#include "cctk.h"
+
+c     ========================================================================
+
+      SUBROUTINE ADMmass_integrand3D(origin,Dx,Dy,Dz,x,y,z,gxx,gxy,gxz,
+     &     gyy,gyz,gzz,ADMmass_int,Psi,Psi_power)
+
+c     ------------------------------------------------------------------------
+c
+c     Estimates the ADM mass at a given radius using Equation (7) from
+c     O Murchadha and York, Phys Rev D, 10, 1974 page 2345
+c
+c     ------------------------------------------------------------------------
+      
+      IMPLICIT NONE
+
+c     Input variables
+      INTEGER,INTENT(IN) ::
+     &     Psi_power
+      CCTK_REAL,INTENT(IN) :: 
+     &     Dx,Dy,Dz,origin(3)
+      CCTK_REAL,DIMENSION(:),INTENT(IN) :: 
+     &     x,y,z
+      CCTK_REAL,DIMENSION(:,:,:),INTENT(IN) :: 
+     &     gxx,gxy,gxz,gyy,gyz,gzz,Psi
+      
+c     Output variables
+      CCTK_REAL,DIMENSION(:,:,:),INTENT(OUT) :: 
+     &     ADMmass_int
+      
+c     Local variables, here only
+      INTEGER :: 
+     &     i,j,k,ip
+      CCTK_REAL,PARAMETER ::
+     &     half = 0.5D0
+      CCTK_REAL :: 
+     &     rad,ux,uy,uz,det,dxx,dxy,dxz,dyy,dyz,dzz,uxx,uxy,uxz,uyy,
+     &     uyz,uzz,term1,term2,term3,dxx_y,dxx_z,dxy_x,dxy_y,dxy_z,
+     &     dyy_x,dxz_x,dxz_y,dxz_z,dyz_x,dzz_x,dyy_z,dyz_y,dyz_z,dzz_y,
+     &     Pi,idet,p,pip,pim,pjp,pjm,pkp,pkm
+
+c     ------------------------------------------------------------------------
+
+      Pi = ACOS(-1D0)
+
+c     Because other codes evolve Psi**4
+      SELECT CASE (Psi_power)
+
+      CASE (1)
+         ip = 4
+
+      CASE (4)
+         ip = 1
+
+      CASE DEFAULT
+         WRITE(*,*) "This value of Psi_power is not supported"
+
+      END SELECT   
+      
+
+      DO k = 2, SIZE(z)-1
+         DO j = 2, SIZE(y)-1
+            DO i = 2, SIZE(x)-1
+               
+               rad = SQRT((x(i)-origin(1))**2 
+     &              +(y(j)-origin(2))**2
+     &              +(z(k)-origin(3))**2)
+               
+               IF (rad.NE.0) THEN
+
+                  ux = (x(i)-origin(1))/rad
+                  uy = (y(j)-origin(2))/rad
+                  uz = (z(k)-origin(3))/rad
+
+c                 Abbreviations for metric coefficients
+c                 -------------------------------------
+                  p = psi(i,j,k)**ip
+
+                  dxx = p*gxx(i,j,k); dxy = p*gxy(i,j,k)
+                  dxz = p*gxz(i,j,k); dyy = p*gyy(i,j,k)
+                  dyz = p*gyz(i,j,k); dzz = p*gzz(i,j,k)
+
+c                 Determinant of 3-metric
+c                 -----------------------
+                  det = (dxx*dyy*dzz + 2.0D0*dxy*dxz*dyz
+     &                 - (dxx*dyz**2 + dyy*dxz**2 + dzz*dxy**2))
+
+                  idet = 1.0/det
+
+c                 Inverse 3-metric
+c                 ----------------
+                  uxx = idet*(dyy*dzz - dyz**2)
+                  uyy = idet*(dxx*dzz - dxz**2)
+                  uzz = idet*(dxx*dyy - dxy**2)
+                  uxy = idet*(dxz*dyz - dzz*dxy)
+                  uxz = idet*(dxy*dyz - dyy*dxz)
+                  uyz = idet*(dxy*dxz - dxx*dyz)
+
+c                 Derivatives of the 3-metric
+c                 ---------------------------
+                  pip = psi(i+1,j,k)**ip
+                  pim = psi(i-1,j,k)**ip
+                  pjp = psi(i,j+1,k)**ip
+                  pjm = psi(i,j-1,k)**ip
+                  pkp = psi(i,j,k+1)**ip
+                  pkm = psi(i,j,k-1)**ip
+
+                  dxx_y = half/Dy*(pjp*gxx(i,j+1,k)-pjm*gxx(i,j-1,k))
+                  dxx_z = half/Dz*(pkp*gxx(i,j,k+1)-pkm*gxx(i,j,k-1))
+                  dxy_x = half/Dx*(pip*gxy(i+1,j,k)-pim*gxy(i-1,j,k))
+                  dxy_y = half/Dy*(pjp*gxy(i,j+1,k)-pjm*gxy(i,j-1,k))
+                  dxy_z = half/Dz*(pkp*gxy(i,j,k+1)-pkm*gxy(i,j,k-1))
+                  dyy_x = half/Dx*(pip*gyy(i+1,j,k)-pim*gyy(i-1,j,k))
+                  dyy_z = half/Dz*(pkp*gyy(i,j,k+1)-pkm*gyy(i,j,k-1))
+                  dxz_x = half/Dx*(pip*gxz(i+1,j,k)-pim*gxz(i-1,j,k))
+                  dxz_y = half/Dy*(pjp*gxz(i,j+1,k)-pjm*gxz(i,j-1,k))
+                  dxz_z = half/Dz*(pkp*gxz(i,j,k+1)-pkm*gxz(i,j,k-1))
+                  dyz_x = half/Dx*(pip*gyz(i+1,j,k)-pim*gyz(i-1,j,k))
+                  dyz_y = half/Dy*(pjp*gyz(i,j+1,k)-pjm*gyz(i,j-1,k))
+                  dyz_z = half/Dz*(pkp*gyz(i,j,k+1)-pkm*gyz(i,j,k-1))
+                  dzz_x = half/Dx*(pip*gzz(i+1,j,k)-pim*gzz(i-1,j,k))
+                  dzz_y = half/Dy*(pjp*gzz(i,j+1,k)-pjm*gzz(i,j-1,k))
+
+                  term1 = uxy*(dxx_y-dxy_x)+uxz*(dxx_z-dxz_x)
+     &                 +uyy*(dxy_y-dyy_x)+uyz*(dxy_z-dyz_x)
+     &                 +uyz*(dxz_y-dyz_x)+uzz*(dxz_z-dzz_x)
+
+                  term2 = uyz*(dyy_z-dyz_y)+uxy*(dyy_x-dxy_y)
+     &                 +uzz*(dyz_z-dzz_y)+uxz*(dyz_x-dxz_y)
+     &                 +uxz*(dxy_z-dxz_y)+uxx*(dxy_x-dxx_y)
+
+                  term3 = uxz*(dzz_x-dxz_z)+uyz*(dzz_y-dyz_z)
+     &                 +uxx*(dxz_x-dxx_z)+uxy*(dxz_y-dxy_z)
+     &                 +uxy*(dyz_x-dxy_z)+uyy*(dyz_y-dyy_z)
+
+                  ADMmass_int(i,j,k) = 1.0D0/16.0D0/Pi*(ux*term1+
+     &                 uy*term2+uz*term3)*SQRT(det)*rad**2
+
+               ELSE
+                  
+                  ADMmass_int(i,j,k) = 0.0D0
+                  
+               ENDIF
+               
+          ENDDO
+        ENDDO
+      ENDDO
+
+
+c     This is needed when the grid is an octant, but it does not hurt
+c     if it is not
+
+      DO k = 2, size(z)-1
+         DO j = 2, size(y)-1
+            ADMmass_int(1,j,k) = ADMmass_int(2,j,k)
+         ENDDO
+      ENDDO
+      DO k = 2, size(z)-1
+         DO i = 1, size(x)-1
+            ADMmass_int(i,1,k) = ADMmass_int(i,2,k)
+         ENDDO
+      ENDDO
+      DO j = 1, size(y)-1
+         DO i = 1, size(x)-1
+            ADMmass_int(i,j,1) = ADMmass_int(i,j,2)
+         ENDDO
+      ENDDO
+
+
+      END SUBROUTINE ADMmass_integrand3D
diff --git a/src/ADMmass_integrand3D_int.F b/src/ADMmass_integrand3D_int.F
new file mode 100644
index 0000000..c103f55
--- /dev/null
+++ b/src/ADMmass_integrand3D_int.F
@@ -0,0 +1,27 @@
+
+#include "cctk.h"
+ 
+      MODULE ADMmass_integrand3D_int
+
+c     ------------------------------------------------------------------
+
+      INTERFACE
+
+          SUBROUTINE  ADMmass_integrand3D(origin,Dx,Dy,Dz,x,y,z,gxx,gxy,
+     &        gxz,gyy,gyz,gzz,ADMmass_int,Psi,Psi_power)
+          IMPLICIT NONE
+          INTEGER,INTENT(IN) ::
+     &         Psi_power
+          CCTK_REAL,INTENT(IN) :: 
+     &        Dx,Dy,Dz,origin(3)
+          CCTK_REAL,DIMENSION(:),INTENT(IN) :: 
+     &        x,y,z
+          CCTK_REAL,DIMENSION(:,:,:),INTENT(IN) :: 
+     &        gxx,gxy,gxz,gyy,gyz,gzz,Psi
+          CCTK_REAL,DIMENSION(:,:,:),INTENT(OUT) :: 
+     &        ADMmass_int
+          END SUBROUTINE  ADMmass_integrand3D
+
+      END INTERFACE
+
+      END MODULE ADMmass_integrand3D_int
diff --git a/src/D2_extract.F b/src/D2_extract.F
new file mode 100644
index 0000000..3b24ae0
--- /dev/null
+++ b/src/D2_extract.F
@@ -0,0 +1,918 @@
+
+#include "cctk.h"
+
+c     ========================================================================
+
+      SUBROUTINE D2_extract(eta,igrid,Nt,Np,theta,phi,all_modes,l,m,g00,g11,
+     &     g12,g13,g22,g23,g33,dg22,dg23,dg33,
+     &     do_ADMmass,ADMmass_int1,ADMmass_int2,
+     &     do_momentum,momentum_int1,momentum_int2,momentum_int3,
+     &     do_spin,spin_int1,spin_int2,spin_int3,
+     &     mass,rsch,Qodd,Qeven,ADMmass,momentum,spin,dtaudt)
+
+c     ------------------------------------------------------------------------
+c
+c     Extract the odd and even parity gauge invariant variables
+c     at a given isotropic radius (eta) from the center of a 
+c     slightly perturbed Schwarzschild spacetime.
+c     
+c     The polar coordinates and metric variables must be given on 
+c     an equally spaced grid, which is offset from the axis.
+c
+c     Notice that for octant symmetry we only need calculate the real 
+c     parts of the gauge invariant variables for the even-l,m modes.
+c     This subroutine is optimised (!!) for octant symmetry.
+c
+c     ------------------------------------------------------------------------
+
+      IMPLICIT NONE
+
+c     Input variables
+      INTEGER,INTENT(IN) :: 
+     &     igrid,l,m
+      CCTK_INT,INTENT(IN) ::
+     &     Nt,Np,all_modes,do_momentum,do_spin
+      INTEGER,INTENT(IN) ::
+     &     do_ADMmass(2)
+      CCTK_REAL,INTENT(IN) :: 
+     &     eta,theta(:),phi(:)
+      CCTK_REAL,DIMENSION (:,:) :: 
+     &     g00,g11,g12,g13,g22,g23,g33,dg22,dg23,dg33,
+     &     ADMmass_int1,ADMmass_int2,
+     &     momentum_int1,momentum_int2,momentum_int3,
+     &     spin_int1,spin_int2,spin_int3
+
+c     Output variables
+      CCTK_REAL,INTENT(OUT) :: 
+     &     dtaudt,ADMmass(2),mass,rsch,Qodd(:,2:,0:),Qeven(:,2:,0:),
+     &     momentum(3),spin(3)
+
+c     Local Variables
+      LOGICAL, PARAMETER :: 
+     &    debug = .FALSE.
+      INTEGER ::  
+     &    i,j,il,im,lmin,lmax,mmin,mmax,lstep,mstep
+      CCTK_REAL,PARAMETER ::
+     &    zero  = 0.0D0,
+     &    half  = 0.5D0,
+     &    one   = 1.0D0,
+     &    two   = 2.0D0,
+     &    four  = 4.0D0,
+     &    six   = 6.0D0,
+     &    eight = 8.0D0
+      CCTK_REAL :: 
+     &     Dt,Dp,st,ist,drschdri,dridrsch,
+     &     g22comb,g23comb,g33comb,S,rsch2,Lambda,
+     &     Pi,
+     &     fac_h1,fac_H2,fac_K,fac_G,fac_c1,fac_c2,
+     &     fac_dG,fac_dK,fac_dc2,fac1,fac2,fac3,max_g11,
+     &     sph_g11,sph_g22,sph_g33,sph_dg22,sph_dg33,error,
+     &     sphere_int
+      CCTK_REAL,DIMENSION(2) :: 
+     &     Y,Y1,Y2,Y3,Y4,h1,H2,K,G,c1,c2,dG,dK,dc2
+      CCTK_REAL,DIMENSION(Nt+1,Np+1) :: 
+     &     temp,h1i,H2i,Gi,Ki,c1i,c2i,dGi,dKi,dc2i
+
+
+c     ------------------------------------------------------------------------
+c
+c     0. Initial things
+c
+c     ------------------------------------------------------------------------
+
+      Pi    = two*ASIN(one)
+
+      SELECT CASE (igrid)
+      CASE (0)                  ! full grid
+        Dt = Pi/DBLE(Nt)
+        Dp = two*Pi/DBLE(Np)
+        lstep = 1
+        mstep = 1
+      CASE (1)                  ! octant grid
+        Dt = half*Pi/DBLE(Nt)
+        Dp = half*Pi/DBLE(Np)
+        lstep = 2
+        mstep = 2
+      CASE (2)                  ! cartoon
+        Dt = Pi/DBLE(Nt)
+        Dp = 2D0*Pi
+        lstep = 2
+        mstep = 2
+      CASE (3)                  ! bitant
+        Dt = half*Pi/DBLE(Nt)
+        Dp = two*Pi/DBLE(Np)         
+        lstep = 1
+        mstep = 1
+      CASE DEFAULT
+        WRITE (*,*) "Grid incorrectly set in D2_extract"
+        STOP
+      END SELECT
+         
+c     Calculate ADM mass
+c     ------------------
+      IF (do_ADMmass(1) == 1) THEN
+
+c	Calculate ADM mass using standard formula 
+         DO j = 2, Np+1		
+            DO i = 2, Nt+1
+               temp(i,j) = ADMmass_int1(i-1,j-1)*DSIN(theta(i-1))
+            ENDDO
+         ENDDO
+         ADMmass(1) = sphere_int(temp,igrid,Nt+1,Np+1,Dt,Dp)
+
+      ENDIF
+      IF (do_ADMmass(2) == 1) THEN
+ 
+c       Calculate ADM mass using conformal formula
+         DO j = 2, Np+1		
+            DO i = 2, Nt+1
+               temp(i,j) = ADMmass_int2(i-1,j-1)*DSIN(theta(i-1))
+            ENDDO
+         ENDDO
+         ADMmass(2) = sphere_int(temp,igrid,Nt+1,Np+1,Dt,Dp)
+         
+      END IF
+
+c     Calculate momentum
+c     ------------------
+      IF (do_momentum == 1) THEN
+
+         DO j = 2, Np+1		
+            DO i = 2, Nt+1
+               temp(i,j) = momentum_int1(i-1,j-1)*DSIN(theta(i-1))
+            ENDDO
+         ENDDO
+         momentum(1) = sphere_int(temp,igrid,Nt+1,Np+1,Dt,Dp)
+         
+         DO j = 2, Np+1		
+            DO i = 2, Nt+1
+               temp(i,j) = momentum_int2(i-1,j-1)*DSIN(theta(i-1))
+            ENDDO
+         ENDDO
+         momentum(2) = sphere_int(temp,igrid,Nt+1,Np+1,Dt,Dp)
+ 
+         DO j = 2, Np+1		
+            DO i = 2, Nt+1
+               temp(i,j) = momentum_int3(i-1,j-1)*DSIN(theta(i-1))
+            ENDDO
+         ENDDO
+         momentum(3) = sphere_int(temp,igrid,Nt+1,Np+1,Dt,Dp)
+ 
+      ENDIF
+
+c     Calculate spin
+c     --------------
+      IF (do_spin == 1) THEN
+
+         DO j = 2, Np+1		
+            DO i = 2, Nt+1
+               temp(i,j) = spin_int1(i-1,j-1)*DSIN(theta(i-1))
+            ENDDO
+         ENDDO
+         spin(1) = sphere_int(temp,igrid,Nt+1,Np+1,Dt,Dp)
+         
+         DO j = 2, Np+1		
+            DO i = 2, Nt+1
+               temp(i,j) = spin_int2(i-1,j-1)*DSIN(theta(i-1))
+            ENDDO
+         ENDDO
+         spin(2) = sphere_int(temp,igrid,Nt+1,Np+1,Dt,Dp)
+ 
+         DO j = 2, Np+1		
+            DO i = 2, Nt+1
+               temp(i,j) = spin_int3(i-1,j-1)*DSIN(theta(i-1))
+            ENDDO
+         ENDDO
+         spin(3) = sphere_int(temp,igrid,Nt+1,Np+1,Dt,Dp)
+ 
+      ENDIF
+
+c     ------------------------------------------------------------------
+c
+c     1. Find Spherical Parts of Metric Components
+c
+c     ------------------------------------------------------------------
+
+c     ... g00
+      DO j = 2, Np+1
+        DO i = 2, Nt+1
+          temp(i,j) = g00(i-1,j-1)*DSIN(theta(i-1))
+        ENDDO
+      ENDDO
+      dtaudt = sqrt(one/(four*Pi)*sphere_int(temp,igrid,Nt+1,Np+1,Dt,Dp))
+
+c     ... g11
+      DO j = 2, Np+1
+         DO i = 2, Nt+1
+            temp(i,j) = g11(i-1,j-1)*DSIN(theta(i-1))
+         ENDDO
+      ENDDO
+      sph_g11 = one/(four*Pi)*sphere_int(temp,igrid,Nt+1,Np+1,Dt,Dp)
+
+c     ... g22
+      DO j = 2, Np+1
+         DO i = 2, Nt+1
+            temp(i,j) = g22(i-1,j-1)*DSIN(theta(i-1))
+         ENDDO
+      ENDDO
+      sph_g22 = one/(four*Pi)*sphere_int(temp,igrid,Nt+1,Np+1,Dt,Dp)
+
+c     ... dg22
+      DO j = 2, Np+1
+         DO i = 2, Nt+1
+            temp(i,j) = dg22(i-1,j-1)*DSIN(theta(i-1))
+         ENDDO
+      ENDDO
+      sph_dg22 = one/(four*Pi)*sphere_int(temp,igrid,Nt+1,Np+1,Dt,Dp)
+
+
+c     ... g33/sin(theta)**2
+      DO j = 2, Np+1
+         DO i = 2, Nt+1
+            temp(i,j) = g33(i-1,j-1)/DSIN(theta(i-1))
+         ENDDO
+      ENDDO
+      sph_g33 = one/(four*Pi)*sphere_int(temp,igrid,Nt+1,Np+1,Dt,Dp)
+
+
+c     Calculate error in orthonormality of spherical harmonic
+c      DO j = 2, Np+1
+c         DO i = 2, Nt+1
+c            CALL spher_harm_combs(theta(i),phi(j),l,m,Y,Y1,Y2,Y3,Y4)
+c            temp(i,j) = (Y(1)*Y(1)+Y(2)*Y(2))*sin(theta(i-1))
+c         ENDDO
+c      ENDDO
+c      error = sphere_int(temp,igrid,Nt+1,Np+1,Dt,Dp)
+c      print *,"Error is",one-error
+
+
+
+c     ------------------------------------------------------------------
+c
+c     1. Compute the Schwarzschild radius
+c   
+c     ------------------------------------------------------------------
+
+      DO j = 2, Np+1
+         DO i = 2, Nt+1
+c            temp(i,j) = g22(i-1,j-1)*SIN(theta(i-1))
+            temp(i,j) = (g22(i-1,j-1)+g33(i-1,j-1)/SIN(theta(i-1))**2)
+     &           *SIN(theta(i-1))/(two)
+         ENDDO
+      ENDDO
+
+      rsch2= one/(four*Pi)*sphere_int(temp,igrid,Nt+1,Np+1,Dt,Dp)
+
+      rsch = DSQRT(rsch2)
+
+
+c     ------------------------------------------------------------------
+c
+c     2. Compute the derivative of the schwarzschild radius with  
+c        respect to the isotropic radius eta.
+c
+c     ------------------------------------------------------------------
+
+      DO j = 2, Np+1         
+         DO i = 2, Nt+1
+c           temp(i,j) = (dg22(i-1,j-1))*SIN(theta(i-1))
+           temp(i,j) = (dg22(i-1,j-1)+dg33(i-1,j-1)/SIN(theta(i-1))**2)
+     &             *SIN(theta(i-1))/(two)
+         ENDDO
+      ENDDO
+
+      drschdri = one/(eight*Pi*rsch)*
+     &     sphere_int(temp,igrid,Nt+1,Np+1,Dt,Dp)
+
+      dridrsch = one/drschdri
+
+
+c     ------------------------------------------------------------------
+c
+c     3. Calculate the Schwarzschild mass and S factor
+c
+c     ------------------------------------------------------------------
+
+      S = drschdri**2/sph_g11
+      mass = rsch*(one-S)/two
+
+
+
+c     ------------------------------------------------------------------
+c
+c     4. Subtract spherical part from metric
+c        Do not know how important this is
+c
+c     ------------------------------------------------------------------
+
+      DO i=1,Nt
+         DO j=1,Np
+            st  = SIN(theta(i))
+            g11(i,j) = g11(i,j) - sph_g11
+            g22(i,j) = g22(i,j) - sph_g22
+            g33(i,j) = g33(i,j) - sph_g22*st**2
+            dg22(i,j) = dg22(i,j) - sph_dg22
+            dg33(i,j) = dg33(i,j) - sph_dg22*st**2
+         ENDDO
+      ENDDO
+
+
+      IF (all_modes == 0) THEN
+         lmin = l ; lmax = l
+      ELSE
+         lmin = 2 ; lmax = l
+      ENDIF
+      
+      loop_l: DO il = lmin,lmax,lstep
+         
+         IF (all_modes == 0) THEN
+            mmin = m ; mmax = m
+         ELSE
+            mmin = 0 ; mmax = il
+         END IF
+         
+         loop_m: DO im = mmin,mmax,mstep
+
+c     ------------------------------------------------------------------
+c     
+c     5. Calculate all Regge-Wheeler variables
+c
+c     ------------------------------------------------------------------
+
+c           Factors independent of angular coordinates
+c           ------------------------------------------
+            fac_h1  = dridrsch/DBLE(il*(il+1))
+            fac_H2  = S*dridrsch**2
+            fac_G   = one/(rsch**2*DBLE(il*(il+1)*(il-1)*(il+2)))
+            fac_K   = half/rsch**2
+            fac_c1  = dridrsch/DBLE(il*(il+1))
+            fac_c2  = two/DBLE(il*(il+1)*(il-1)*(il+2))
+            fac_dG  = fac_G
+            fac_dK  = fac_K
+            fac_dc2 = fac_c2
+  
+          
+c     ------------------------------------------------------------------
+c   
+c     5a. Real parts of the Regge-Wheeler variables
+c
+c     ------------------------------------------------------------------
+
+c           Calculate integrands
+
+            loop_phi1: DO j = 1, Np
+               loop_theta1: DO i = 1, Nt
+
+                  st  = SIN(theta(i))
+                  ist = one/st
+                  
+                  CALL spher_harm_combs(theta(i),phi(j),il,im,Y,Y1,Y2,Y3
+     &                 ,Y4)
+                  
+                  h1i(i+1,j+1) =  st*g12(i,j)*Y1(1)+ist*g13(i,j)*Y2(1)    
+                  H2i(i+1,j+1) =  st*g11(i,j)*Y(1)
+                  Gi(i+1,j+1)  =  (st*g22(i,j)-ist*g33(i,j))*Y3(1)
+     &                 +four*ist*g23(i,j)*Y4(1)
+                  Ki(i+1,j+1) =  (st*g22(i,j)+ist*g33(i,j))*Y(1)              
+                  c1i(i+1,j+1) =  g13(i,j)*Y1(1)-g12(i,j)*Y2(1)           
+                  c2i(i+1,j+1) =  (g22(i,j)-ist**2*g33(i,j))*Y4(1)
+     &                 -g23(i,j)*Y3(1)
+                  
+                  g22comb = dridrsch*dg22(i,j)-two/rsch*g22(i,j)
+                  g23comb = dridrsch*dg23(i,j)-two/rsch*g23(i,j)
+                  g33comb = dridrsch*dg33(i,j)-two/rsch*g33(i,j)
+                  
+                  dGi(i+1,j+1) =  (st*g22comb-ist*g33comb)*Y3(1)
+     &                 +four*ist*g23comb*Y4(1)                  
+                  dKi(i+1,j+1) =  (st*g22comb+ist*g33comb)*Y(1)
+                  dc2i(i+1,j+1) =  (dg22(i,j)-ist**2*dg33(i,j))*Y4(1)
+     &                 -dg23(i,j)*Y3(1)
+
+               END DO loop_theta1    
+            END DO loop_phi1
+
+c           Integrations over the 2-sphere
+
+            h1(1) = fac_h1*sphere_int(h1i,igrid,Nt+1,Np+1,Dt,Dp)
+            H2(1) = fac_H2*sphere_int(H2i,igrid,Nt+1,Np+1,Dt,Dp)
+            G(1) = fac_G*sphere_int(Gi,igrid,Nt+1,Np+1,Dt,Dp)
+            K(1) = fac_K*sphere_int(Ki,igrid,Nt+1,Np+1,Dt,Dp)
+     &           +DBLE(il*(il+1))*half*G(1)
+            c1(1) = fac_c1*sphere_int(c1i,igrid,Nt+1,Np+1,Dt,Dp)
+            c2(1) = fac_c2*sphere_int(c2i,igrid,Nt+1,Np+1,Dt,Dp)
+            dG(1) = fac_dG*sphere_int(dGi,igrid,Nt+1,Np+1,Dt,Dp)
+            dK(1) = fac_dK*sphere_int(dKi,igrid,Nt+1,Np+1,Dt,Dp)
+     &           +DBLE(il*(il+1))*half*dG(1)
+            dc2(1) = fac_dc2*sphere_int(dc2i,igrid,Nt+1,Np+1,Dt,Dp)
+
+
+c     ------------------------------------------------------------------
+c   
+c     5b. Imaginary parts of the Regge-Wheeler variables
+c
+c     ------------------------------------------------------------------
+
+            complex_parts: IF (igrid /= 1 .AND. igrid /=2) THEN
+
+c              Calculate integrands
+
+               loop_phi2: DO j = 1, Np
+
+                  loop_theta2: DO i = 1, Nt
+                  
+                     st  = SIN(theta(i))
+                     ist = one/st
+                  
+                     CALL spher_harm_combs(theta(i),phi(j),il,im,Y,Y1,Y2
+     &                    ,Y3,Y4)
+                  
+                     h1i(i+1,j+1) =-st*g12(i,j)*Y1(2)-ist*g13(i,j)*Y2(2)
+                     H2i(i+1,j+1) = -st*g11(i,j)*Y(2) 
+                     Gi(i+1,j+1)  = -(st*g22(i,j)-ist*g33(i,j))*Y3(2)
+     &                    -four*ist*g23(i,j)*Y4(2)                  
+                     Ki(i+1,j+1) = -(st*g22(i,j)+ist*g33(i,j))*Y(2)   
+                     c1i(i+1,j+1) = -g13(i,j)*Y1(2)+g12(i,j)*Y2(2) 
+                     c2i(i+1,j+1) = -(g22(i,j)-ist**2*g33(i,j))*Y4(2)
+     &                    +g23(i,j)*Y3(2)
+                     
+                     g22comb = dridrsch*dg22(i,j)-two/rsch*g22(i,j)
+                     g23comb = dridrsch*dg23(i,j)-two/rsch*g23(i,j)
+                     g33comb = dridrsch*dg33(i,j)-two/rsch*g33(i,j)
+                     
+                     dGi(i+1,j+1) = -(st*g22comb-ist*g33comb)*Y3(2)
+     &                    -four*ist*g23comb*Y4(2) 
+                     dKi(i+1,j+1) = -(st*g22comb+ist*g33comb)*Y(2)  
+                     dc2i(i+1,j+1) = -(dg22(i,j)-ist**2*dg33(i,j))*Y4(2)
+     &                    +dg23(i,j)*Y3(2)
+                  
+                  END DO loop_theta2    
+               END DO loop_phi2
+
+
+c              Integrate over 2-sphere
+
+               h1(2) = fac_h1*sphere_int(h1i,igrid,Nt+1,Np+1,Dt,Dp)
+               H2(2) = fac_H2*sphere_int(H2i,igrid,Nt+1,Np+1,Dt,Dp)
+               G(2) = fac_G*sphere_int(Gi,igrid,Nt+1,Np+1,Dt,Dp)
+               K(2) = fac_K*sphere_int(Ki,igrid,Nt+1,Np+1,Dt,Dp)
+     &              +DBLE(il*(il+1))*half*G(2)
+               c1(2) = fac_c1*sphere_int(c1i,igrid,Nt+1,Np+1,Dt,Dp)
+               c2(2) = fac_c2*sphere_int(c2i,igrid,Nt+1,Np+1,Dt,Dp)
+               dG(2) = fac_dG*sphere_int(dGi,igrid,Nt+1,Np+1,Dt,Dp)
+               dK(2) = fac_dK*sphere_int(dKi,igrid,Nt+1,Np+1,Dt,Dp)
+     &              +DBLE(il*(il+1))*half*dG(2)
+               dc2(2) = fac_dc2*sphere_int(dc2i,igrid,Nt+1,Np+1,Dt,Dp)
+
+            ELSE
+
+               h1(2) = zero
+               H2(2) = zero
+               G(2) = zero
+               K(2) = zero
+               c1(2) = zero
+               c2(2) = zero
+               dG(2) = zero
+               dK(2) = zero
+               dc2(2) = zero
+               
+
+            END IF complex_parts
+            
+
+
+c     ------------------------------------------------------------------
+c     
+c     6. Construct gauge invariant variables
+c
+c     ------------------------------------------------------------------
+
+            Qodd(:,il,im) = SQRT(two*DBLE((il+2)*(il+1)*il*(il-1)))*S/
+     &           rsch*(c1+half*(dc2-two/rsch*c2))
+            
+            Lambda = DBLE((il-1)*(il+2))+3.0D0*(one-S)
+            
+            Qeven(:,il,im) = one/Lambda*SQRT(two*DBLE((il-1)*(il+2))/
+     &           DBLE(il*(il+1)))*(DBLE(il*(il+1))*S*(rsch**2*dG-
+     &           two*h1)+two*rsch*S*(H2-rsch*dK)
+     &           +Lambda*rsch*K)
+            
+
+         END DO loop_m
+         
+      END DO loop_l
+
+
+
+c     ------------------------------------------------------------------
+c     
+c     7. Write some diagnostics if required (for last l,m to be used)
+c
+c     ------------------------------------------------------------------
+
+      IF (debug) THEN
+         WRITE(101,*) eta,h1
+         WRITE(102,*) eta,H2
+         WRITE(103,*) eta,G
+         WRITE(104,*) eta,K
+         WRITE(105,*) eta,c1
+         WRITE(106,*) eta,c2
+	 WRITE(107,*) eta,dG
+         WRITE(108,*) eta,dK
+         WRITE(109,*) eta,dc2
+         WRITE(201,*) eta,mass
+	 WRITE(202,*) eta,rsch
+         WRITE(203,*) eta,drschdri
+	 WRITE(301,*) eta,Qeven(:,lmax,mmax)
+         WRITE(302,*) eta,Qodd(:,lmax,mmax)
+      ENDIF
+
+
+      END SUBROUTINE D2_extract
+
+c     ===============================================================   
+c
+c
+c
+c
+c
+c     ===============================================================
+      
+      SUBROUTINE simpsons(n,Dx,f,intf)
+      
+c     Simpsons rule to evaluate 1D integral equation
+      
+c     Variables In:
+c     n     number of points [must be odd]
+c     Dx    step size
+c     f(n)  integrand at each abscissa
+      
+c     Variables Out:
+c     intf  evaluated integral
+      
+c     Variables Local:
+c     i     index
+      
+c     ---------------------------------------------------------------
+      
+      IMPLICIT NONE
+      
+c     Input variables
+      INTEGER,INTENT(IN) :: 
+     &     n        
+      CCTK_REAL,INTENT(IN) :: 
+     &     Dx,f(n)      
+      
+c     Output variables
+      CCTK_REAL,INTENT(OUT) :: 
+     &     intf    
+      
+c     Local variables
+      INTEGER :: 
+     &     i
+      CCTK_REAL,PARAMETER ::
+     &     two   = 2.0D0,
+     &     three = 3.0D0,
+     &     four  = 4.0D0
+      
+c     ---------------------------------------------------------------
+
+      IF (MOD(n,2) == 0) THEN
+        WRITE(*,*) "Call to -simpsons- with even number of points"
+        STOP
+      END IF
+            
+      intf = f(1) + four*f(n-1) + f(n)
+      
+      DO i = 2, n-3,2
+         
+         intf = intf + four*f(i) + two*f(i+1)
+         
+      ENDDO
+      
+      intf = intf*Dx/three
+      
+      END SUBROUTINE simpsons
+      
+c     ===============================================================
+c
+c
+c
+c
+c
+c     ===============================================================
+
+      FUNCTION sphere_int(temp,igrid,Nt,Np,Dt,Dp)
+
+c     ---------------------------------------------------------------
+c
+c     Integration the function held in temp over the sphere. 
+c     If the full grid is being used then Simpsons rule is used, 
+c     if an octant is used then the symmetry means Simpsons rule 
+c     becomes a simple sum.
+c
+c     ---------------------------------------------------------------
+
+      IMPLICIT NONE
+
+c     Input variables
+      INTEGER,INTENT(IN) :: 
+     &   Nt,Np,igrid
+      CCTK_REAL,INTENT(IN) :: 
+     &   temp(Nt,Np),Dt,Dp
+
+c     Output variables
+      CCTK_REAL ::
+     &   sphere_int
+
+c     Local variables
+      INTEGER ::
+     &   i,j,Np_start
+      CCTK_REAL :: 
+     &   ft(Nt),fp(Np)
+      CCTK_REAL,PARAMETER ::
+     &   two  = 2.0D0,
+     &   four = 4.0D0
+   
+c     ---------------------------------------------------------------  
+
+      DO j = 2, Np
+
+c        Integrate with respect to theta
+
+         DO i = 2, Nt
+            ft(i) = temp(i,j)
+         ENDDO
+         ft(1) = ft(2)     ! from symmetry across axis
+         
+         SELECT CASE (igrid)
+         CASE (0)
+            CALL simpsons(Nt,Dt,ft,fp(j))
+         CASE (1)
+            fp(j)=two*Dt*SUM(ft(2:Nt))
+         CASE (2)
+            CALL simpsons(Nt,Dt,ft,fp(j))
+         CASE (3)
+            fp(j)=two*Dt*SUM(ft(2:Nt))
+         END SELECT
+
+      ENDDO
+
+c     Integrate with respect to phi
+
+      fp(1) = fp(2)             ! from symmetry across axis
+
+
+      SELECT CASE (igrid)
+      CASE (0)
+         CALL simpsons(Np,Dp,fp,sphere_int)
+      CASE (1)
+         sphere_int=four*Dp*SUM(fp(2:Np))
+      CASE (2)
+         sphere_int=Dp*fp(1)
+      CASE (3)
+         CALL simpsons(Np,Dp,fp,sphere_int)
+      END SELECT   
+
+
+      END FUNCTION sphere_int     
+
+c     ==================================================================
+c
+c
+c
+c
+c
+c     ==================================================================
+
+      SUBROUTINE spherical_harmonic(l,m,theta,phi,Ylm)
+
+c     ------------------------------------------------------------------
+c
+c     Calculate the (l,m) spherical harmonic at given angular
+c     coordinates. This number is in general complex, and
+c
+c      Ylm = Ylm(1) + i Ylm(2)
+c
+c      where
+c     
+c                a    ( 2 l + 1 (l-|m|)! )                      i m phi
+c      Ylm = (-1) SQRT( ------- -------  ) P_l|m| (cos(theta)) e
+c                     (   4 Pi  (l+|m|)! ) 
+c
+c      and where
+c
+c      a = m/2 (sign(m)+1)
+c
+c     ------------------------------------------------------------------
+
+      IMPLICIT NONE
+
+c     Input variables
+      INTEGER ::
+     &   l,m
+      CCTK_REAL ::  
+     &   theta,phi
+
+c     Output variables
+      CCTK_REAL :: 
+     &   Ylm(2)
+
+c     Local variables
+      INTEGER ::  
+     &    i
+      CCTK_REAL,PARAMETER ::
+     &    one  = 1.0D0,
+     &    two  = 2.0D0,
+     &    four = 4.0D0
+      CCTK_REAL ::
+     &    a,Pi,fac,plgndr
+
+c     ------------------------------------------------------------------
+
+      Pi = ACOS(-one)
+
+      fac = one
+      DO i = l-ABS(m)+1,l+ABS(m)
+        fac = fac*DBLE(i)
+      ENDDO
+      fac = one/fac
+
+c      a = (-one)**((m*ISIGN(m,1)/ABS(m)+m)/2)*SQRT(DBLE(2*l+1)
+c     &    /four/Pi*fac)*plgndr(l,ABS(m),cos(theta))
+
+      a = (-one)**MAX(m,0)*SQRT(DBLE(2*l+1)
+     &    /four/Pi*fac)*plgndr(l,ABS(m),cos(theta))
+      Ylm(1) = a*COS(DBLE(m)*phi)
+      Ylm(2) = a*SIN(DBLE(m)*phi)
+
+      END SUBROUTINE spherical_harmonic
+
+c     ==================================================================
+c
+c
+c
+c
+c
+c     ==================================================================
+
+      FUNCTION plgndr(l,m,x)
+
+c     ------------------------------------------------------------------
+c
+c     From Numerical Recipes
+c
+c     Calculates the associated Legendre polynomial Plm(x).
+c     Here m and l are integers satisfying 0 <= m <= l,
+c     while x lies in the range -1 <= x <= 1
+c
+c     ------------------------------------------------------------------
+
+c     Input variables
+      INTEGER,INTENT(IN) :: 
+     &    l,m
+      CCTK_REAL,INTENT(IN) :: 
+     &    x
+
+c     Output variables
+      CCTK_REAL :: 
+     &    plgndr
+
+c     Local Variables
+      INTEGER :: 
+     &    i,ll
+      CCTK_REAL,PARAMETER ::
+     &    one = 1.0D0,
+     &    two = 2.0D0
+      CCTK_REAL :: 
+     &    pmm,somx2,fact,pmmp1,pll
+
+c     ------------------------------------------------------------------
+
+      pmm = one
+
+      IF (m.GT.0) THEN
+        somx2=SQRT((one-x)*(one+x))
+        fact=one
+        DO i=1,m
+          pmm  = -pmm*fact*somx2
+          fact = fact+two
+        END DO
+      ENDIF
+
+      IF (l.EQ.m) THEN
+         plgndr = pmm
+      ELSE
+         pmmp1 = x*(two*m+one)*pmm
+         IF (l.EQ.m+1) THEN
+            plgndr=pmmp1
+         ELSE
+            DO ll=m+2,l
+               pll   = (x*DBLE(2*ll-1)*pmmp1-
+     &              DBLE(ll+m-1)*pmm)/DBLE(ll-m)
+               pmm   = pmmp1
+               pmmp1 = pll
+            END DO
+            plgndr = pll
+         ENDIF
+      ENDIF
+
+      END FUNCTION plgndr
+
+c     ==================================================================
+c
+c
+c
+c
+c
+c     ==================================================================
+
+      SUBROUTINE spher_harm_combs(theta,phi,l,m,Y,Y1,Y2,Y3,Y4)
+c     ------------------------------------------------------------------
+c
+c     Calculates the various combinations of spherical harmonics needed
+c     for the extraction (all are complex):     
+c
+c       Y  = Ylm
+c       Y1 = Ylm,theta
+c       Y2 = Ylm,phi
+c       Y3 = Ylm,theta,theta-cot theta Ylm,theta-Ylm,phi,phi/sin^2 theta
+c       Y4 = Ylm,theta,phi-cot theta Ylm,phi
+c
+c     The local variables Yplus is the spherical harmonic at (l+1,m)
+c
+c     ------------------------------------------------------------------
+
+      IMPLICIT NONE
+
+c     Input variables
+      INTEGER ::
+     &     l,m
+      CCTK_REAL ::
+     &     theta,phi
+
+c     Output variables
+      CCTK_REAL,DIMENSION(2) ::
+     &     Y,Y1,Y2,Y3,Y4
+
+c     Local variables
+      INTEGER :: 
+     &     i
+      CCTK_REAL ::  
+     &     Yplus(2),rl,rm,cot_theta,Pi
+      CCTK_REAL,PARAMETER ::
+     &     half = 0.5D0,
+     &     one  = 1.0D0,
+     &     two  = 2.0D0
+
+c     ------------------------------------------------------------------
+
+      Pi = DACOS(-one)
+
+      rl = DBLE(l)
+      rm = DBLE(m)
+
+      cot_theta = DCOS(theta)/DSIN(theta)
+
+      CALL spherical_harmonic(l+1,m,theta,phi,Yplus)
+
+
+c     Find Y ...
+
+      CALL spherical_harmonic(l,m,theta,phi,Y)
+
+
+c     Find Y1 ...
+
+      DO i = 1,2
+        Y1(i) = -(rl+one)*cot_theta*Y(i)+Yplus(i)/DSIN(theta)*
+     &          DSQRT(((rl+one)**2-rm**2)*(rl+half)/(rl+one+half))
+      ENDDO
+
+
+c     Find Y2 ...
+
+      Y2(1) = -rm*Y(2)
+      Y2(2) =  rm*Y(1)
+
+
+c     Find Y3 ...
+
+      DO i = 1,2
+        Y3(i) = -two*cot_theta*Y1(i)+(two*rm*rm/(DSIN(theta)**2)
+     &          -rl*(rl+one))*Y(i)
+      ENDDO
+
+
+c     Find Y4 ...
+
+      Y4(1) = rm*(cot_theta*Y(2)-Y1(2))
+      Y4(2) = rm*(Y1(1)-cot_theta*Y(1))
+
+
+      END SUBROUTINE spher_harm_combs
+
+c     ==================================================================
+
+
+
+
diff --git a/src/D2_extract_int.F b/src/D2_extract_int.F
new file mode 100644
index 0000000..43f6906
--- /dev/null
+++ b/src/D2_extract_int.F
@@ -0,0 +1,42 @@
+
+#include "cctk.h"
+
+      MODULE D2_extract_int
+
+c     ------------------------------------------------------------------
+
+      INTERFACE 
+
+      SUBROUTINE D2_extract(eta,igrid,Nt,Np,theta,phi,all_modes,l,m,g00,g11,
+     &        g12,g13,g22,g23,g33,dg22,dg23,dg33,
+     &        do_ADMmass,ADMmass_int1,ADMmass_int2,
+     &        do_momentum,momentum_int1,momentum_int2,momentum_int3,
+     &        do_spin,spin_int1,spin_int2,spin_int3,
+     &        mass,rsch,Qodd,Qeven,ADMmass,momentum,spin,dtaudt)
+
+      IMPLICIT NONE
+      INTEGER,INTENT(IN) :: 
+     &     igrid,l,m
+      CCTK_INT,INTENT(IN) ::
+     &     Nt,Np,all_modes,do_momentum,do_spin
+      INTEGER,INTENT(IN) ::
+     &     do_ADMmass(2)
+      CCTK_REAL,INTENT(IN) :: 
+     &     eta
+      CCTK_REAL,INTENT(IN),DIMENSION (:) :: 
+     &     theta,phi
+      CCTK_REAL,INTENT(IN),DIMENSION (:,:) :: 
+     &     g00,g11,g12,g13,g22,g23,g33,dg22,dg23,dg33,
+     &     ADMmass_int1, ADMmass_int2,
+     &     momentum_int1,momentum_int2,momentum_int3,
+     &     spin_int1,spin_int2,spin_int3
+      CCTK_REAL,INTENT(OUT) :: 
+     &     ADMmass(2),mass,rsch,Qodd(:,:,:),Qeven(:,:,:),dtaudt,
+     &     momentum(3),spin(3)
+      END SUBROUTINE
+      
+      END INTERFACE
+      
+      END MODULE D2_extract_int
+
+
diff --git a/src/D3_extract.F b/src/D3_extract.F
new file mode 100644
index 0000000..f5f6850
--- /dev/null
+++ b/src/D3_extract.F
@@ -0,0 +1,221 @@
+
+#include "cctk.h"
+
+c     ==================================================================
+
+      SUBROUTINE D3_extract(cctkGH,do_ADMmass,do_momentum,do_spin,igrid,
+     &     origin,myproc,Nt,Np,all_modes,
+     &     l,m,x,y,z,Dx,Dy,Dz,Psi_power,Psi,
+     &     g00,gxx,gxy,gxz,gyy,gyz,gzz,hxx,hxy,hxz,hyy,hyz,hzz,
+     &     eta,ADMmass,momentum,spin,mass,rsch,Qodd,Qeven,Extract_temp3d,dtaudt)
+
+c     ------------------------------------------------------------------
+c
+c     Extract odd and even parity gauge invariant variables on a 
+c     2-surface defined by constant isotropic radius ETA from ORIGIN,
+c     also returns the corresponding Schwarzschild radius RSCH.
+c   
+c     Works for a full grid (IGRID=0) or an octant (IGRID=1), given 
+c     non-conformal (PSI=1,Gij) or conformal (PSI^PSI_POWER,Gij).
+c
+c     IMPORTANT : Nt and Np must both be even numbers (so that Simpsons
+c                 rule works later)
+c
+c     Variables in:
+c     ------------
+c     igrid     Grid type (0=full,1=octant)
+c               [Could extract in an octant from a full grid, but why
+c                would you want to]
+c     origin    The position of the origin of extraction in the x,y,z
+c               coordinate system
+c     myproc    Cactus variable, must be zero if not using Cactus
+c     Nt,Np     Number of theta,phi grid DIVISIONS to use [even numbers]
+c     all_modes if true the extract all l,m modes up to l
+c     l,m       Spherical harmonic to extract if ALL_MODES is false. 
+c               otherwise m is ignored and l is the maximum l mode
+c     x,y,z     Cartesian coordinates in the 3-space, hopefully 
+c               isotropic about ORIGIN
+c     Dx,Dy,Dz  Grid spacings
+c     Psi_power The power of Psi which is being passed in, usually one
+c               but may be four
+c     Psi       The conformal factor, set it to one if non-conformal 
+c               data is passed in
+c     gij       The 3-metric components in cartesian coordinates
+c     eta       The coordinate radius from ORIGIN for extraction, must
+c               obviously give a 2-surface which is contained in the
+c               grid
+c
+c     Variables out:
+c     -------------
+c     ADMmass   Estimate of ADM mass 
+c     momentum  Estimate of momentum
+c     spin      Estimate of spin
+c     mass      The extracted mass
+c     rsch      The extracted Schwarzschild radius
+c     Qodd      The extracted odd parity gauge invariant variable
+c     Qeven     The extracted even parity gauge invariant variable     
+c     
+c     ------------------------------------------------------------------
+
+      USE D3_to_D2_int
+      USE cartesian_to_spherical_int
+      USE unphysical_to_physical_int
+      USE D2_extract_int
+
+c     ------------------------------------------------------------------
+
+      IMPLICIT NONE
+
+c     Input variables
+
+      CCTK_POINTER :: cctkGH
+
+      INTEGER,INTENT(IN) :: 
+     &     igrid,l,m,Psi_power,myproc
+      CCTK_INT,INTENT(IN) :: 
+     &     Nt,Np,all_modes,do_momentum,do_spin
+      INTEGER,INTENT(IN) ::
+     &     do_ADMmass(2)
+      CCTK_REAL,INTENT(IN) :: 
+     &     origin(3),Dx,Dy,Dz,eta
+      CCTK_REAL,INTENT(IN),DIMENSION(:,:,:) :: 
+     &     Psi,g00,gxx,gxy,gxz,gyy,gyz,gzz,
+     &     hxx,hxy,hxz,hyy,hyz,hzz
+      CCTK_REAL,INTENT(IN),DIMENSION(:) :: 
+     &     x,y,z
+      CCTK_REAL,INTENT(INOUT),DIMENSION(:,:,:) ::
+     &     Extract_temp3d
+
+c     Output variables
+      CCTK_REAL,INTENT(OUT) :: 
+     &     ADMmass(2),mass,rsch,Qodd(:,2:,0:),Qeven(:,2:,0:),dtaudt,
+     &     momentum(3),spin(3)
+
+c     Local variables passed on
+      CCTK_REAL :: 
+     &    theta(Nt),phi(Np)
+      CCTK_REAL,DIMENSION(Nt,Np) :: 
+     &    Psis,g00s,gxxs,gxys,gxzs,gyys,gyzs,gzzs,dPsis,dgxxs,dgxys,dgxzs,
+     &    dgyys,dgyzs,dgzzs,grr,grt,grp,gtt,gtp,gpp,dgtt,dgtp,dgpp,
+     &    ADMmass_int1,ADMmass_int2,
+     &    momentum_int1,momentum_int2,momentum_int3,
+     &    spin_int1,spin_int2,spin_int3
+
+c     Local variables here only
+      INTEGER :: 
+     &    i,j
+      CCTK_REAL,PARAMETER ::
+     &    half = 0.5D0,
+     &    one  = 1.0D0,
+     &    two  = 2.0D0
+      CCTK_REAL :: 
+     &    Pi,Dt,Dp
+
+c     ------------------------------------------------------------------
+
+
+c     ------------------------------------------------------------------
+c
+c     1. Specify polar coordinates on chosen 2-surface
+c
+c     ------------------------------------------------------------------
+
+      Pi = two*ASIN(one)
+
+c     Set polar gridspacing
+
+      SELECT CASE (igrid)
+      CASE (0)                  ! full grid
+        Dt = Pi/DBLE(Nt)
+        Dp = two*Pi/DBLE(Np)
+      CASE (1)                  ! octant grid
+        Dt = half*Pi/DBLE(Nt)
+        Dp = half*Pi/DBLE(Np)
+      CASE (2)                  ! cartoon
+        Dt = Pi/DBLE(Nt)
+        Dp = 2D0*Pi
+      CASE (3)                  ! bitant
+        Dt = half*Pi/DBLE(Nt)
+        Dp = two*Pi/DBLE(Np)
+      CASE DEFAULT
+        WRITE (*,*) "Grid incorrectly set in D3_extract"
+        STOP
+      END SELECT
+         
+
+c     Set coordinates 
+
+      DO i = 1, Nt
+         theta(i) = (DBLE(i) - half)* Dt
+      ENDDO
+      DO j = 1, Np
+         phi(j) = (DBLE(j) - half)* Dp
+      ENDDO
+
+      IF (igrid == 2) phi(1) = 0D0
+      
+
+c     ------------------------------------------------------------------
+c
+c     2. Project quantities onto the 2-surface
+c
+c     ------------------------------------------------------------------
+
+      CALL D3_to_D2(cctkGH,do_ADMmass,do_momentum,do_spin,
+     &     Psi_power,origin,myproc,Dx,Dy,Dz,Psi,
+     &     g00,gxx,gxy,gxz,gyy,gyz,gzz,hxx,hxy,hxz,hyy,hyz,hzz,
+     &     x,y,z,eta,Nt,Np,theta,phi,Psis,g00s,gxxs,gxys,gxzs,
+     &     gyys,gyzs,gzzs,dPsis,dgxxs,dgxys,dgxzs,dgyys,dgyzs,dgzzs,
+     &     ADMmass_int1,ADMmass_int2,momentum_int1,momentum_int2,
+     &     momentum_int3,spin_int1,spin_int2,spin_int3,Extract_temp3d)
+
+c     ------------------------------------------------------------------
+c     Now this is done the grid on processor 0 has correct values and the 
+c     grid on all other processors have junk. This means that processor
+c     0 needs to calculate the wave forms and the other processors have
+c     to sit and do nothing. So we shall put this in a conditional
+c     to only happen on processor 0.
+c     PW Jan 31 98
+
+      if (myproc .eq. 0) then 
+c     ------------------------------------------------------------------
+c
+c     3. Transform tensor components to polar coordinates on 2-surface
+c
+c     ------------------------------------------------------------------
+
+         CALL cartesian_to_spherical(theta,phi,eta,gxxs,gxys,gxzs,gyys,
+     &        gyzs,gzzs,grr,grt,grp,gtt,gtp,gpp,dgxxs,dgxys,dgxzs,dgyys,
+     &        dgyzs,dgzzs,dgtt,dgtp,dgpp)
+
+
+c     ------------------------------------------------------------------
+c
+c     4. Calculate physical quantities on 2-surface
+c
+c     ------------------------------------------------------------------
+
+         CALL unphysical_to_physical(grr,grt,grp,gtt,gtp,gpp,dgtt,dgtp,
+     &        dgpp,Psis,dPsis,Psi_power)
+
+
+c     ------------------------------------------------------------------
+c
+c     5. Extract gauge invariant quantities on 2-surface
+c
+c     ------------------------------------------------------------------
+
+         CALL D2_extract(eta,igrid,Nt,Np,theta,phi,all_modes,l,m,g00s,grr,
+     &        grt,grp,gtt,gtp,gpp,dgtt,dgtp,dgpp,
+     &        do_ADMmass,ADMmass_int1,ADMmass_int2,
+     &        do_momentum,momentum_int1,momentum_int2,momentum_int3,
+     &        do_spin,spin_int1,spin_int2,spin_int3,
+     &        mass,rsch,Qodd,Qeven,ADMmass,momentum,spin,dtaudt)
+
+c     End of myproc = 0 conditional
+      endif
+
+      END SUBROUTINE D3_extract
+
+
+
diff --git a/src/D3_extract_int.F b/src/D3_extract_int.F
new file mode 100644
index 0000000..3bb617a
--- /dev/null
+++ b/src/D3_extract_int.F
@@ -0,0 +1,44 @@
+
+#include "cctk.h"
+
+      MODULE D3_extract_int
+
+c     ------------------------------------------------------------------
+
+      INTERFACE
+
+      SUBROUTINE D3_extract(cctkGH,do_ADMmass,do_momentum,do_spin,igrid,
+     &        origin,myproc,
+     &        Nt,Np,all_modes,l,m,x,y,z,Dx,Dy,Dz,Psi_power,Psi,
+     &        g00,gxx,gxy,gxz,gyy,gyz,gzz,hxx,hxy,hxz,hyy,hyz,hzz,
+     &        eta,ADMmass,momentum,spin,mass,rsch,Qodd,Qeven,
+     &        Extract_temp3d,dtaudt)
+
+      IMPLICIT NONE
+
+      CCTK_POINTER :: cctkGH
+
+      INTEGER,INTENT(IN) :: 
+     &     igrid,l,m,Psi_power,myproc
+      CCTK_INT,INTENT(IN) ::
+     &     Nt,Np,all_modes,do_momentum,do_spin
+      INTEGER,INTENT(IN) ::
+     &     do_ADMmass(2)
+      CCTK_REAL,INTENT(IN) :: 
+     &     origin(3),Dx,Dy,Dz,eta
+      CCTK_REAL,INTENT(IN),DIMENSION(:) :: 
+     &     x,y,z
+      CCTK_REAL,INTENT(IN),DIMENSION(:,:,:) :: 
+     &     Psi,g00,gxx,gxy,gxz,gyy,gyz,gzz,
+     &     hxx,hxy,hxz,hyy,hyz,hzz
+      CCTK_REAL,INTENT(INOUT),DIMENSION(:,:,:) ::
+     &     Extract_Temp3d
+      CCTK_REAL,INTENT(OUT) :: 
+     &     ADMmass(2),mass,rsch,Qodd(:,:,:),Qeven(:,:,:),dtaudt,
+     &     momentum(3),spin(3)
+
+      END SUBROUTINE
+         
+      END INTERFACE
+
+      END MODULE D3_extract_int
diff --git a/src/D3_to_D2.F b/src/D3_to_D2.F
new file mode 100644
index 0000000..884f8b8
--- /dev/null
+++ b/src/D3_to_D2.F
@@ -0,0 +1,398 @@
+
+#include "cctk.h"
+
+c     ========================================================================
+
+      SUBROUTINE D3_to_D2(cctkGH,do_ADMmass,do_momentum,do_spin,
+     &     Psi_power,origin,myproc,Dx,Dy,Dz,Psi,
+     &     g00,gxx,gxy,gxz,gyy,gyz,gzz,hxx,hxy,hxz,hyy,hyz,hzz,
+     &     x,y,z,eta,Nt,Np,theta,phi,
+     &     Psis,g00s,gxxs,gxys,gxzs,gyys,gyzs,gzzs,dPsis,dgxxs,dgxys,dgxzs,
+     &     dgyys,dgyzs,dgzzs,ADMmass_int1,ADMmass_int2,
+     &     momentum_int1,momentum_int2,momentum_int3,
+     &     spin_int1,spin_int2,spin_int3,Extract_temp3d)
+ 
+c     ------------------------------------------------------------------------
+c
+c     Project the 3-metric and its 1st radial derivatives onto a 
+c     2-sphere.
+c
+c     ------------------------------------------------------------------------
+
+      USE met_rad_der_int
+      USE ADMmass_integrand3D_int
+      USE momentum_integrand3D_int
+      USE spin_integrand3D_int
+
+      IMPLICIT NONE
+
+c     Input variables
+
+      CCTK_POINTER :: cctkGH
+
+      INTEGER,INTENT(IN) :: 
+     &     myproc,Psi_power
+      CCTK_INT, INTENT(IN) ::
+     &     Nt,Np,do_momentum,do_spin
+      CCTK_REAL,INTENT(IN) :: 
+     &     origin(3),Dx,Dy,Dz,eta
+      CCTK_REAL,INTENT(IN),DIMENSION(:) :: 
+     &     theta,phi,x,y,z
+      CCTK_REAL,INTENT(IN),DIMENSION(:,:,:) ::  
+     &     Psi,g00,gxx,gxy,gxz,gyy,gyz,gzz,
+     &     hxx,hxy,hxz,hyy,hyz,hzz
+      CCTK_REAL,INTENT(INOUT),DIMENSION(:,:,:) ::
+     &     Extract_temp3d
+      LOGICAL ::
+     &     do_ADMmass(2)
+
+c     Output variables
+
+      CCTK_REAL,INTENT(OUT),DIMENSION(:,:) ::     
+     &     Psis,g00s,gxxs,gxys,gxzs,gyys,
+     &     gyzs,gzzs,dPsis,dgxxs,dgxys,dgxzs,dgyys,dgyzs,dgzzs,
+     &     ADMmass_int1,ADMmass_int2,
+     &     momentum_int1,momentum_int2,momentum_int3,
+     &     spin_int1,spin_int2,spin_int3
+
+
+c     Local variables, passed on
+
+      LOGICAL :: 
+     &    err_flag
+      INTEGER :: 
+     &    iorder,npoints,Nx,Ny,Nz
+      INTEGER,DIMENSION(Nt,Np) :: 
+     &    ib,jb,kb
+      CCTK_REAL,DIMENSION(Nt,Np) :: 
+     &    xs,ys,zs,ux,uy,uz,xb,yb,zb
+
+c     Local variables, here only
+
+      INTEGER :: 
+     &    i,j,handle,ierror
+      CCTK_REAL ::
+     &    xmin,xmax,ymin,ymax,zmin,zmax
+
+c     ------------------------------------------------------------------------
+      
+c     Initial Stuff
+c     -------------
+      Nx = SIZE(x)
+      Ny = SIZE(y)
+      Nz = SIZE(z)
+      
+
+c     Compute Cartesian coordinates on the surface
+c     --------------------------------------------
+      DO j = 1, Np
+         DO i = 1, Nt
+
+            xs(i,j) = origin(1)+eta*SIN(theta(i))*COS(phi(j))
+            ys(i,j) = origin(2)+eta*SIN(theta(i))*SIN(phi(j))
+            zs(i,j) = origin(3)+eta*COS(theta(i))
+
+         ENDDO
+      ENDDO
+
+
+c     Only do interpolation on one processor
+c     --------------------------------------
+      SELECT CASE (myproc) 
+
+      CASE (0)
+         npoints = Np*Nt
+
+      CASE DEFAULT
+         npoints = 0
+
+      END SELECT 
+
+
+c     Choose the interpolator
+c     -----------------------    
+      call CCTK_GetInterpHandle (handle, "simple")
+
+c     Find the origin of the spatial coordinates
+c     ------------------------------------------
+      call CCTK_CoordRange(ierror,cctkGH,xmin,xmax,"x")
+      call CCTK_CoordRange(ierror,cctkGH,ymin,ymax,"y")
+      call CCTK_CoordRange(ierror,cctkGH,zmin,zmax,"z")
+
+
+c     Project un-physical metric and conformal factor  onto sphere
+c     ------------------------------------------------------------
+      call CCTK_Interp (ierror, cctkGH, handle, npoints, 3, 8, 8,
+     $          nx, ny, nz,
+     $          xs, ys, zs,
+     $          CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
+     $          CCTK_VARIABLE_REAL,
+     $          xmin, ymin, zmin,
+     $          dx, dy, dz,
+     $          Psi,g00,gxx,gxy,gxz,gyy,gyz,gzz,
+     $          CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,
+     $          CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,
+     $          CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,
+     $          Psis,g00s,gxxs,gxys,gxzs,gyys,gyzs,gzzs,
+     $          CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,
+     $          CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,
+     $          CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL
+     $     )
+
+
+c     Calculate radial derivatives and project onto sphere
+c     ----------------------------------------------------
+      CALL met_rad_der(origin,Dx,Dy,Dz,x,y,z,Psi,Extract_temp3d)
+      CALL CCTK_SyncGroup(cctkGH,"extract::temps")   
+      call CCTK_Interp (ierror, cctkGH, handle, npoints, 3, 1, 1,
+     $          nx, ny, nz,
+     $          xs, ys, zs,
+     $          CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,
+     $          xmin, ymin, zmin,
+     $          dx, dy, dz,
+     $          Extract_temp3d,
+     $          CCTK_VARIABLE_REAL,
+     $          dPsis,
+     $          CCTK_VARIABLE_REAL
+     $     )
+
+      CALL met_rad_der(origin,Dx,Dy,Dz,x,y,z,gxx,Extract_temp3d)
+      CALL CCTK_SyncGroup(cctkGH,"extract::temps")
+      call CCTK_Interp (ierror, cctkGH, handle, npoints, 3, 1, 1,
+     $          nx, ny, nz,
+     $          xs, ys, zs,
+     $          CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,
+     $          xmin, ymin, zmin,
+     $          dx, dy, dz,
+     $          Extract_temp3d,
+     $          CCTK_VARIABLE_REAL,
+     $          dgxxs,
+     $          CCTK_VARIABLE_REAL
+     $     )
+
+      CALL met_rad_der(origin,Dx,Dy,Dz,x,y,z,gxy,Extract_temp3d)
+      CALL CCTK_SyncGroup(cctkGH,"extract::temps")
+      call CCTK_Interp (ierror, cctkGH, handle, npoints, 3, 1, 1,
+     $          nx, ny, nz,
+     $          xs, ys, zs,
+     $          CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,
+     $          xmin, ymin, zmin,
+     $          dx, dy, dz,
+     $          Extract_temp3d,
+     $          CCTK_VARIABLE_REAL,
+     $          dgxys,
+     $          CCTK_VARIABLE_REAL
+     $     )
+
+      CALL met_rad_der(origin,Dx,Dy,Dz,x,y,z,gxz,Extract_temp3d)
+      CALL CCTK_SyncGroup(cctkGH,"extract::temps")
+      call CCTK_Interp (ierror, cctkGH, handle, npoints, 3, 1, 1,
+     $          nx, ny, nz,
+     $          xs, ys, zs,
+     $          CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,
+     $          xmin, ymin, zmin,
+     $          dx, dy, dz,
+     $          Extract_temp3d,
+     $          CCTK_VARIABLE_REAL,
+     $          dgxzs,
+     $          CCTK_VARIABLE_REAL
+     $     )
+
+      CALL met_rad_der(origin,Dx,Dy,Dz,x,y,z,gyy,Extract_temp3d)
+      CALL CCTK_SyncGroup(cctkGH,"extract::temps")
+      call CCTK_Interp (ierror, cctkGH, handle, npoints, 3, 1, 1,
+     $          nx, ny, nz,
+     $          xs, ys, zs,
+     $          CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,
+     $          xmin, ymin, zmin,
+     $          dx, dy, dz,
+     $          Extract_temp3d,
+     $          CCTK_VARIABLE_REAL,
+     $          dgyys,
+     $          CCTK_VARIABLE_REAL
+     $     )
+
+      CALL met_rad_der(origin,Dx,Dy,Dz,x,y,z,gyz,Extract_temp3d)
+      CALL CCTK_SyncGroup(cctkGH,"extract::temps")
+      call CCTK_Interp (ierror, cctkGH, handle, npoints, 3, 1, 1,
+     $          nx, ny, nz,
+     $          xs, ys, zs,
+     $          CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,
+     $          xmin, ymin, zmin,
+     $          dx, dy, dz,
+     $          Extract_temp3d,
+     $          CCTK_VARIABLE_REAL,
+     $          dgyzs,
+     $          CCTK_VARIABLE_REAL
+     $     )
+
+
+      CALL met_rad_der(origin,Dx,Dy,Dz,x,y,z,gzz,Extract_temp3d)
+      CALL CCTK_SyncGroup(cctkGH,"extract::temps")
+      call CCTK_Interp (ierror, cctkGH, handle, npoints, 3, 1, 1,
+     $          nx, ny, nz,
+     $          xs, ys, zs,
+     $          CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,
+     $          xmin, ymin, zmin,
+     $          dx, dy, dz,
+     $          Extract_temp3d,
+     $          CCTK_VARIABLE_REAL,
+     $          dgzzs,
+     $          CCTK_VARIABLE_REAL
+     $     )
+
+
+c     Calculate integrands for ADM masses
+c     -----------------------------------
+c     Standard equation
+      IF (do_ADMmass(1)) THEN
+         CALL ADMmass_integrand3D(origin,Dx,Dy,Dz,x,y,z,gxx,gxy,
+     &        gxz,gyy,gyz,gzz,Extract_temp3d,Psi,Psi_power)
+         CALL CCTK_SyncGroup(cctkGH,"extract::temps")
+         call CCTK_Interp (ierror, cctkGH, handle, npoints, 3, 1, 1,
+     $        nx, ny, nz,
+     $        xs, ys, zs,
+     $        CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,
+     $        xmin, ymin, zmin,
+     $        dx, dy, dz,
+     $        Extract_temp3d,
+     $        CCTK_VARIABLE_REAL,
+     $        ADMmass_int1,
+     $        CCTK_VARIABLE_REAL
+     $        )
+
+      END IF
+
+c     Conformal equation
+      IF (do_ADMmass(2)) THEN
+             ADMmass_int2 = -eta**2/2D0/3.1416D0*dPsis
+      ENDIF	
+
+c     Calculate integrands for momentum
+c     ---------------------------------
+      IF (do_momentum==1) THEN
+
+         CALL momentum_integrand3D(origin,Dx,Dy,Dz,x,y,z,
+     &        gxx,gxy,gxz,gyy,gyz,gzz,hxx,hxy,hxz,hyy,hyz,hzz,
+     &        Extract_temp3d,Psi,Psi_power)
+         CALL CCTK_SyncGroup(cctkGH,"extract::temps")
+         call CCTK_Interp (ierror, cctkGH, handle, npoints, 3, 1, 1,
+     $        nx, ny, nz,
+     $        xs, ys, zs,
+     $        CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,
+     $        xmin, ymin, zmin,
+     $        dx, dy, dz,
+     $        Extract_temp3d,
+     $        CCTK_VARIABLE_REAL,
+     $        momentum_int1,
+     $        CCTK_VARIABLE_REAL
+     $        )
+
+         CALL momentum_integrand3D(origin,Dx,Dy,Dz,x,y,z,
+     &        gxx,gxy,gxz,gyy,gyz,gzz,hxx,hxy,hxz,hyy,hyz,hzz,
+     &        Extract_temp3d,Psi,Psi_power)
+         CALL CCTK_SyncGroup(cctkGH,"extract::temps")
+         call CCTK_Interp (ierror, cctkGH, handle, npoints, 3, 1, 1,
+     $        nx, ny, nz,
+     $        xs, ys, zs,
+     $        CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,
+     $        xmin, ymin, zmin,
+     $        dx, dy, dz,
+     $        Extract_temp3d,
+     $        CCTK_VARIABLE_REAL,
+     $        momentum_int2,
+     $        CCTK_VARIABLE_REAL
+     $        )
+
+
+         CALL momentum_integrand3D(origin,Dx,Dy,Dz,x,y,z,
+     &        gxx,gxy,gxz,gyy,gyz,gzz,hxx,hxy,hxz,hyy,hyz,hzz,
+     &        Extract_temp3d,Psi,Psi_power)
+         CALL CCTK_SyncGroup(cctkGH,"extract::temps")
+         call CCTK_Interp (ierror, cctkGH, handle, npoints, 3, 1, 1,
+     $        nx, ny, nz,
+     $        xs, ys, zs,
+     $        CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,
+     $        xmin, ymin, zmin,
+     $        dx, dy, dz,
+     $        Extract_temp3d,
+     $        CCTK_VARIABLE_REAL,
+     $        momentum_int3,
+     $        CCTK_VARIABLE_REAL
+     $        )
+
+      END IF
+
+c     Calculate integrands for spin
+c     -----------------------------
+      IF (do_momentum==1) THEN
+
+         CALL spin_integrand3D(origin,Dx,Dy,Dz,x,y,z,
+     &        gxx,gxy,gxz,gyy,gyz,gzz,hxx,hxy,hxz,hyy,hyz,hzz,
+     &        Extract_temp3d,Psi,Psi_power)
+         CALL CCTK_SyncGroup(cctkGH,"extract::temps")
+         call CCTK_Interp (ierror, cctkGH, handle, npoints, 3, 1, 1,
+     $        nx, ny, nz,
+     $        xs, ys, zs,
+     $        CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,
+     $        xmin, ymin, zmin,
+     $        dx, dy, dz,
+     $        Extract_temp3d,
+     $        CCTK_VARIABLE_REAL,
+     $        spin_int1,
+     $        CCTK_VARIABLE_REAL
+     $        )
+
+         CALL spin_integrand3D(origin,Dx,Dy,Dz,x,y,z,
+     &        gxx,gxy,gxz,gyy,gyz,gzz,hxx,hxy,hxz,hyy,hyz,hzz,
+     &        Extract_temp3d,Psi,Psi_power)
+         CALL CCTK_SyncGroup(cctkGH,"extract::temps")
+         call CCTK_Interp (ierror, cctkGH, handle, npoints, 3, 1, 1,
+     $        nx, ny, nz,
+     $        xs, ys, zs,
+     $        CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,
+     $        xmin, ymin, zmin,
+     $        dx, dy, dz,
+     $        Extract_temp3d,
+     $        CCTK_VARIABLE_REAL,
+     $        spin_int2,
+     $        CCTK_VARIABLE_REAL
+     $        )
+
+
+         CALL spin_integrand3D(origin,Dx,Dy,Dz,x,y,z,
+     &        gxx,gxy,gxz,gyy,gyz,gzz,hxx,hxy,hxz,hyy,hyz,hzz,
+     &        Extract_temp3d,Psi,Psi_power)
+         CALL CCTK_SyncGroup(cctkGH,"extract::temps")
+         call CCTK_Interp (ierror, cctkGH, handle, npoints, 3, 1, 1,
+     $        nx, ny, nz,
+     $        xs, ys, zs,
+     $        CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,
+     $        xmin, ymin, zmin,
+     $        dx, dy, dz,
+     $        Extract_temp3d,
+     $        CCTK_VARIABLE_REAL,
+     $        spin_int3,
+     $        CCTK_VARIABLE_REAL
+     $        )
+
+      END IF
+
+      END SUBROUTINE D3_to_D2
+
+
+      
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/src/D3_to_D2_int.F b/src/D3_to_D2_int.F
new file mode 100644
index 0000000..297e22e
--- /dev/null
+++ b/src/D3_to_D2_int.F
@@ -0,0 +1,45 @@
+
+#include "cctk.h"
+
+      MODULE D3_to_D2_int
+
+c     ------------------------------------------------------------------
+
+      INTERFACE
+
+      SUBROUTINE D3_to_D2(cctkGH,do_ADMmass,do_momentum,do_spin,
+     &        Psi_power,origin,myproc,Dx,Dy,Dz,Psi,
+     &        g00,gxx,gxy,gxz,gyy,gyz,gzz,hxx,hxy,hxz,hyy,hyz,hzz,
+     &        x,y,z,eta,Nt,Np,theta,phi,Psis,g00s,gxxs,gxys,
+     &        gxzs,gyys,gyzs,gzzs,dPsis,dgxxs,dgxys,dgxzs,dgyys,dgyzs,
+     &        dgzzs,ADMmass_int1,ADMmass_int2,
+     &        momentum_int1,momentum_int2,momentum_int3,
+     &        spin_int1,spin_int2,spin_int3,Extract_temp3d)
+
+      IMPLICIT NONE
+      CCTK_POINTER :: cctkGH
+      INTEGER,INTENT(IN) :: 
+     &     myproc,Psi_power
+      CCTK_INT, INTENT(IN) ::
+     &     Nt,Np,do_momentum,do_spin
+      INTEGER ::
+     &     do_ADMmass(2)
+      CCTK_REAL,INTENT(IN) :: 
+     &     origin(3),Dx,Dy,Dz,eta
+      CCTK_REAL,INTENT(IN),DIMENSION(:,:,:) :: 
+     &     Psi,g00,gxx,gxy,gxz,gyy,gyz,gzz,
+     &     hxx,hxy,hxz,hyy,hyz,hzz
+      CCTK_REAL,INTENT(IN),DIMENSION(:) :: 
+     &     x,y,z,theta,phi
+      CCTK_REAL,INTENT(OUT),DIMENSION(:,:) :: 
+     &     Psis,g00s,gxxs,gxys,gxzs,gyys,gyzs,gzzs,dPsis,dgxxs,dgxys,dgxzs,
+     &     dgyys,dgyzs,dgzzs,ADMmass_int1,ADMmass_int2,
+     &     momentum_int1,momentum_int2,momentum_int3,
+     &     spin_int1,spin_int2,spin_int3
+      CCTK_REAL,INTENT(INOUT),DIMENSION(:,:,:) ::
+     &     Extract_temp3d
+      END SUBROUTINE
+
+      END INTERFACE
+
+      END MODULE D3_to_D2_int
diff --git a/src/Extract.F b/src/Extract.F
new file mode 100644
index 0000000..1aa8534
--- /dev/null
+++ b/src/Extract.F
@@ -0,0 +1,814 @@
+c/*@@
+c  @file      Extract.F
+c  @date      28th October 1997   
+c  @author    Gab Allen
+c  @desc      Waveform extraction
+c  
+c  @enddesc 
+c@@*/
+
+#include "cctk.h"
+#include "cctk_parameters.h"
+#include "cctk_arguments.h"
+
+c/*@@
+c  @routine    Extract
+c  @date       28th October 1997
+c  @author     Gab Allen
+c  @desc       Entry point for waveform extraction routines
+c  
+c  @enddesc 
+c  @calls      D3_extract
+c  @calledby   No idea
+c  @history    
+c
+c  @endhistory 
+c@@*/
+
+
+      SUBROUTINE Extract(CCTK_FARGUMENTS)
+
+      USE D3_extract_int
+
+      IMPLICIT NONE
+
+      DECLARE_CCTK_FARGUMENTS
+      DECLARE_CCTK_PARAMETERS
+      DECLARE_CCTK_FUNCTIONS
+
+c     Non-Cactus input variables for D3_extract
+
+      INTEGER :: 
+     &     igrid,lmode,mmode,Psi_power=1
+      INTEGER :: 
+     &     do_ADMmass(2)
+      CCTK_REAL :: 
+     &     orig(3),radius 
+      CCTK_REAL :: 
+     &     x_1d(cctk_lsh(1)),y_1d(cctk_lsh(2)),z_1d(cctk_lsh(3))
+
+      
+c     Output variables from D3_extract
+
+      CCTK_REAL :: 
+     &     dtaudt,ADMmass(2),mass,rsch,momentum(3),spin(3)
+      CCTK_REAL,ALLOCATABLE ::
+     &     Qodd(:,:,:), Qeven(:,:,:)
+
+ 
+c     Local variables
+
+      INTEGER :: 
+     &     ix,iy,iz,fn1,fn2,lmin,lmax,mmin,mmax,lstep,mstep,
+     &     il,im,it,ndet,idet,nx_parfile,nchar,iconv,ioutput
+      INTEGER,SAVE :: openfiles = 1
+      INTEGER,PARAMETER ::
+     &     max_detectors = 9,number_timecoords = 2
+      CCTK_REAL,PARAMETER ::
+     &     two = 2.0D0
+      CCTK_REAL ::
+     &     time,r1,r2,dr,r_max,xmin,xmax,ymin,ymax,zmin,zmax,gf_min,gf_max
+      CCTK_REAL ::
+     &     detector(max_detectors)
+      CHARACTER*100 ::
+     &     out1,out2,massfile,rschfile,ADMmassfile1,ADMmassfile2,
+     &     momentumfile1,momentumfile2,momentumfile3,spinfile1,spinfile2,spinfile3
+      CHARACTER*3 ::
+     &     timestring
+
+      LOGICAL :: use_proper,use_coordinate,do_output
+      LOGICAL,SAVE :: firsttime
+      CCTK_REAL,DIMENSION(10),SAVE :: proper_time
+      CCTK_REAL,SAVE :: last_time
+      INTEGER, save :: open_file_level(10) = 0
+
+      INTEGER myproc,ierr
+      CCTK_REAL dx,dy,dz
+      character*200 filestr
+      character*80 infoline
+
+c     ------------------------------------------------------------------
+
+c     Get the output directory sorted out
+      call CCTK_FortranString(nchar,outdir,filestr)
+      if (firsttime) then
+         call CCTK_mkdir(ierr,filestr);
+         firsttime = .FALSE.
+      end if
+
+c     ------------------------------------------------------------------
+
+      myproc = CCTK_MyProc(cctkGH)
+      dx = cctk_delta_space(1)
+      dy = cctk_delta_space(2)
+      dz = cctk_delta_space(3)
+
+c     ------------------------------------------------------------------
+c     
+c     0. Check to see if Extract should have been called
+c
+c     ------------------------------------------------------------------
+
+      IF (MODULO(int(cctk_iteration),int(itout)) .NE. 0) RETURN
+
+      if (verbose == 1) then
+         write(infoline,'(A24,G12.7)') 
+     &        'Calling Extract at time ',cctk_time
+         call CCTK_INFO(infoline)
+      end if
+
+c     ------------------------------------------------------------------
+c     
+c     1. Initial Stuff
+c
+c     ------------------------------------------------------------------
+
+
+c     Get the number of polar divisions
+
+      IF (MOD(int(Nt),2) == 1 .OR. MOD(int(Np),2) == 1 .OR. 
+     &     Nt < 0 .OR. Np < 0) THEN
+        call CCTK_WARN(0,"Error in Nt or Np in Extract")
+      END IF
+
+
+c     Get the origin of spherical symmetry
+
+      orig(1) = origin_x 
+      orig(2) = origin_y 
+      orig(3) = origin_z
+      
+
+c     Set the value of igrid
+
+      IF (CCTK_EQUALS(domain,"octant")) THEN
+         igrid = 1
+      ELSEIF (CCTK_EQUALS(domain,"bitant")) THEN
+         igrid = 3
+      ELSE
+         igrid = 0
+      END IF  
+
+#ifdef THORN_CARTOON_2D
+      IF (contains("cartoon_active","yes") .NE. 0) THEN
+         igrid = 2
+         Np = 1
+      END IF
+#endif
+   
+c     Create 1D coordinate arrays
+
+      do ix = 1, cctk_lsh(1)
+        x_1d(ix) = x(ix,1,1)
+      enddo
+      do iy = 1, cctk_lsh(2)
+        y_1d(iy) = y(1,iy,1)
+      enddo
+      do iz = 1, cctk_lsh(3)
+        z_1d(iz) = z(1,1,iz)
+      enddo
+
+c     See if the ADM mass should be calculated and output
+      IF (doADMmass == 1) THEN
+         do_ADMmass(1) = 1
+	 IF (use_conformal == 1) THEN
+            do_ADMmass(2) = 1
+         ELSE
+            do_ADMmass(2) = 0
+	 ENDIF
+      ELSE
+         do_ADMmass(:) = 0
+      END IF
+
+c     Set array containing gtt ... for the moment I am lazy and
+c     just make it the lapse and do not include the shift!
+      g00 = alp
+
+c     What kind of timecoordinate do I want to use
+      if (CCTK_EQUALS(timecoord,"both")) then
+         use_proper = .true.
+         use_coordinate = .true.
+         if (cctk_iteration == 0) then
+            proper_time = 0
+            last_time = cctk_time
+         end if
+      elseif (CCTK_EQUALS(timecoord,"proper")) then
+         use_proper = .true.
+         use_coordinate = .false.
+         if (cctk_iteration == 0) then
+            proper_time = 0
+            last_time = cctk_time
+         end if
+      elseif (CCTK_EQUALS(timecoord,"coord")) then
+         use_proper = .false.
+         use_coordinate = .true.
+      else
+         use_proper = .false.
+         use_coordinate = .true.
+      endif
+
+c     ------------------------------------------------------------------
+c     
+c     2. Sort out which modes to extract
+c
+c     ------------------------------------------------------------------
+
+      lmode  = l_mode
+
+      IF (all_modes == 1) THEN
+         lmin = 2
+	 lmax = lmode
+c     Get m_mode in this case too, or the T3E has a junk value (PW)
+         mmode = 0
+         ALLOCATE(Qodd(2,2:lmode,0:lmode),Qeven(2,2:lmode,0:lmode))
+      ELSE
+         lmin = lmode
+         lmax = lmode
+         mmode = m_mode
+         ALLOCATE(Qodd(2,2:lmode,0:mmode),Qeven(2,2:lmode,0:mmode))
+      ENDIF
+
+
+c     Check do not have bad values for the modes
+
+      IF (lmode < 2 .OR. mmode > lmode .OR. mmode < 0) THEN
+         call CCTK_WARN(0,"Error in lmode,mmode in Extract")
+      END IF
+      IF (lmode > 9) THEN
+         call CCTK_WARN(4,"Filenames will probably be crazy in Extract")
+      END IF
+
+c     If using an octant, do not do odd modes
+
+      IF (igrid == 1 .OR. igrid == 2) THEN
+         lstep = 2 ; mstep = 2
+      ELSE
+         lstep = 1 ; mstep = 1
+      END IF   
+
+c     If we have cartoon, then jump past all the m-modes in loop
+      IF (igrid == 2) mstep = lmode+1
+
+
+c     ------------------------------------------------------------------
+c     
+c     3. Find maximum radius of extraction on grid
+c
+c     ------------------------------------------------------------------
+    
+      call CCTK_CoordRange(ierr,cctkGH,xmin,xmax,"x")
+      call CCTK_CoordRange(ierr,cctkGH,ymin,ymax,"y")
+      call CCTK_CoordRange(ierr,cctkGH,zmin,zmax,"z")
+      
+      IF (igrid == 2) THEN 
+         r_max = MIN(xmax-two*Dx-orig(1),
+     &        zmax-two*Dz-orig(3))
+      ELSE IF (igrid == 1) THEN 
+         r_max = MIN(zmax-two*Dx-orig(1),
+     &        ymax-two*Dy-orig(2),
+     &        zmax-two*Dz-orig(3))
+      ELSE IF (igrid == 3) THEN 
+         r_max = MIN(xmax-two*Dx-orig(1),
+     &        ymax-two*Dy-orig(2),
+     &        zmax-two*Dz-orig(3),
+     &        DABS(xmin+two*Dx-orig(1)),
+     &        DABS(ymin+two*Dy-orig(2)))
+      ELSE
+         r_max = MIN(xmax-two*Dx-orig(1),
+     &        ymax-two*Dy-orig(2),
+     &        zmax-two*Dz-orig(3),
+     &        DABS(xmin+two*Dx-orig(1)),
+     &        DABS(ymin+two*Dy-orig(2)),
+     &        DABS(zmin+two*Dz-orig(3)))
+      END IF
+
+
+
+
+c     ------------------------------------------------------------------
+c     
+c     4. Extract Cauchy initial data for a linear wave equation
+c
+c     ------------------------------------------------------------------
+
+      test_Cauchy: IF (Cauchy == 1) THEN
+
+      it = Cauchy_timestep
+
+      
+c     check_timestep: IF (cctk_iteration == it) THEN
+
+      check_timestep: IF (open_file_level(1) .eq. 0) THEN
+
+c           Open output files
+            write (*,*) "Open output file on level ", 1
+            open_file_level(1)=1;
+
+            test_myproc1: IF (CCTK_MyProc(cctkGH) == 0 ) THEN
+
+               out1 = filestr(1:nchar)//"/rsch_ini.rl"
+               OPEN(UNIT= 77,file = out1,STATUS="unknown")
+               out1 = filestr(1:nchar)//"/mass_ini.rl"
+               OPEN(UNIT= 78,file = out1,STATUS="unknown")
+               IF (do_ADMmass(1) == 1) THEN
+                  out1 = filestr(1:nchar)//"/ADMmass_ini.rl"
+                  OPEN(UNIT= 79,file = out1,STATUS="unknown")
+               ENDIF
+               IF (do_ADMmass(2) == 1) THEN
+                  out1 = filestr(1:nchar)//"/ADMmassc_ini.rl"
+                  OPEN(UNIT= 80,file = out1,STATUS="unknown")
+               ENDIF
+
+               IF (do_momentum == 1) THEN
+                  out1 = filestr(1:nchar)//"/momentum_x_ini.rl"
+                  OPEN(UNIT=781,file = out1,STATUS="unknown")
+                  out1 = filestr(1:nchar)//"/momentum_y_ini.rl"
+                  OPEN(UNIT=782,file = out1,STATUS="unknown")
+                  out1 = filestr(1:nchar)//"/momentum_z_ini.rl"
+                  OPEN(UNIT=783,file = out1,STATUS="unknown")
+               ENDIF
+               IF (do_spin == 1) THEN
+                  out1 = filestr(1:nchar)//"/spin_x_ini.rl"
+                  OPEN(UNIT=784,file = out1,STATUS="unknown")
+                  out1 = filestr(1:nchar)//"/spin_y_ini.rl"
+                  OPEN(UNIT=785,file = out1,STATUS="unknown")
+                  out1 = filestr(1:nchar)//"/spin_z_ini.rl"
+                  OPEN(UNIT=786,file = out1,STATUS="unknown")
+               ENDIF
+
+               loop_l1: DO il = lmin,lmax,lstep
+
+                  IF (all_modes == 0) THEN
+                     mmin = mmode ; mmax = mmode
+                  ELSE
+                     mmin = 0 ; mmax = il
+                  END IF
+
+                  loop_m1: DO im = mmin,mmax,lstep
+
+
+                     fn1 = il*10+im
+                     out1 = filestr(1:nchar)//"/Qeven_ini_"//
+     &                    CHAR(il+48)//CHAR(im+48)//".rl"
+                     OPEN(UNIT=fn1,FILE=out1,STATUS="unknown")
+
+c                    Only print odd-parity if full grid
+                     IF (igrid == 0) THEN
+                        fn2 = 100+il*10+im
+                        out2 = filestr(1:nchar)//"/Qodd_ini_"//
+     &                       CHAR(il+48)//CHAR(im+48)//".rl"
+                        OPEN(UNIT=fn2,FILE=out2,STATUS="unknown")
+                     END IF
+
+                  END DO loop_m1
+
+               END DO loop_l1
+
+            END IF test_myproc1
+            
+            
+c           Find range of extraction radii
+            
+            r1 = Cauchy_r1
+            r2 = r_max 
+            dr = Cauchy_dr
+            
+            
+            
+c           Do extraction at each radius
+            
+            IF (verbose == 1) THEN
+               WRITE(*,*)
+               WRITE(*,*) "Extracting Cauchy initial data"
+               WRITE(*,*) "------------------------------"
+               WRITE(*,*) "From r =",r1
+               WRITE(*,*) "To r =",r2
+               WRITE(*,*) "In steps of ",dr
+               WRITE(*,*)
+            END IF
+            
+            radius = r1
+             
+            extract_at_each_radius: DO WHILE (radius < r2)
+         
+               CALL D3_extract(cctkGH,do_ADMmass,do_momentum,do_spin,
+     &           igrid,orig,myproc,Nt,Np,all_modes,lmode,
+     &           mmode,x_1d,y_1d,z_1d,Dx,Dy,Dz,Psi_power,Psi,g00,
+     &           gxx,gxy,gxz,gyy,gyz,gzz,kxx,kxy,kxz,kyy,kyz,kzz,
+     &           radius,ADMmass,momentum,spin,mass,rsch,
+     &           Qodd,Qeven,temp3d,dtaudt)
+
+               IF (verbose == 1) THEN
+                  WRITE(*,*) "Extracted at r =",radius
+                  WRITE(*,*) "Schwarzschild radius =",rsch
+                  WRITE(*,*) "Schwarzschild mass =",mass
+                  if (do_ADMmass(1) == 1) then
+                     WRITE(*,*) "ADM mass =",ADMmass(1)
+                  end if
+                  if (do_ADMmass(2) == 1) then
+                     WRITE(*,*) "ADM mass (conformal) =",ADMmass(2)
+                  end if
+                  if (do_momentum == 1) then
+                     WRITE(*,*) "momentum = (",momentum(1),",",
+     &                momentum(2),",",momentum(3),")" 
+                  end if
+                  if (do_spin == 1) then
+                     WRITE(*,*) "spin = (",spin(1),",",
+     &                spin(2),",",spin(3),")" 
+                  end if
+               ENDIF
+               
+               
+c              Write to file at each radius
+               
+               test_myproc2: IF (CCTK_MyProc(cctkGH) == 0) THEN
+
+                  WRITE( 77,*) radius,rsch
+                  WRITE( 78,*) radius,mass
+                  IF (do_ADMmass(1) == 1) WRITE( 79,*) radius,ADMmass(1)
+                  IF (do_ADMmass(2) == 1) WRITE( 80,*) radius,ADMmass(2)
+                  IF (do_momentum == 1) then
+                     WRITE(781,*) radius,momentum(1)
+                     WRITE(782,*) radius,momentum(2)
+                     WRITE(783,*) radius,momentum(3)
+                  end if
+
+                  IF (do_spin == 1) THEN 
+                     WRITE(784,*) radius,spin(1)
+                     WRITE(785,*) radius,spin(2)
+                     WRITE(786,*) radius,spin(3)
+                  END IF
+
+                  loop_l2: DO il = lmin,lmax,lstep
+
+                     IF (all_modes == 0) THEN
+                        mmin = mmode ; mmax = mmode
+                     ELSE
+                        mmin = 0 ; mmax = il
+                     END IF
+
+                     loop_m2: DO im = mmin,mmax,mstep
+
+                        fn1 = il*10+im
+                        WRITE(fn1,*) rsch,Qeven(:,il,im)
+
+c                       Only print odd-parity if full grid
+                        IF (igrid == 0) THEN
+                           fn2 = 100+il*10+im
+                           WRITE(fn2,*) rsch,Qeven(:,il,im)
+                        END IF
+
+                     END DO loop_m2
+
+                  END DO loop_l2    
+
+               END IF test_myproc2 
+            
+               radius = radius+dr
+            
+            END DO extract_at_each_radius
+         
+         
+c           Now close all the files
+
+            test_myproc3: IF (myproc == 0) THEN
+
+               CLOSE( 77)   ! Schwarzschild radius 
+               CLOSE( 78)   ! Schwarzschild mass
+               IF (do_ADMmass(1) == 1) CLOSE( 79)   ! ADM mass
+               IF (do_ADMmass(2) == 1) CLOSE( 80)   ! ADM mass
+
+               IF (do_momentum == 1) THEN
+                  CLOSE(781)    ! momentum
+                  CLOSE(782)    ! momentum
+                  CLOSE(783)    ! momentum
+               end if
+
+               IF (do_spin == 1) THEN
+                  CLOSE(784)    ! spin
+                  CLOSE(785)    ! spin
+                  CLOSE(786)    ! spin
+               END IF
+
+               loop_l3: DO il = lmin,lmax,lstep
+
+                  IF (all_modes == 0) THEN
+                     mmin = mmode ; mmax = mmode
+                  ELSE
+                     mmin = 0 ; mmax = il
+                  END IF
+
+                  loop_m3: DO im = mmin,mmax,mstep
+
+                     fn1 = il*10+im
+                     CLOSE(fn1) ! Qeven_ini_lm.dat
+
+                     IF (igrid == 0) THEN
+                        fn2 = 100+il*10+im
+                        CLOSE(fn2) ! Qodd_ini_lm.dat
+                     END IF
+
+                  END DO loop_m3
+
+               END DO loop_l3    
+               
+            END IF test_myproc3
+            
+         END IF check_timestep
+         
+      END IF test_Cauchy
+     
+
+c     ------------------------------------------------------------------
+c     
+c     5. Extract waveforms at detectors
+c
+c     ------------------------------------------------------------------
+
+c     Cannot use the conformal equation for ADM mass now
+      do_ADMmass(2) = 0
+
+      ndet = num_detectors
+
+      test_detectors: IF (ndet > 0) THEN
+
+         IF (ndet > 9) THEN
+            call CCTK_WARN(0,"Too many detectors in Extract")
+         END IF
+      
+         detector(1) = detector1
+         detector(2) = detector2
+         detector(3) = detector3
+         detector(4) = detector4
+         detector(5) = detector5
+         detector(6) = detector6
+         detector(7) = detector7
+         detector(8) = detector8
+         detector(9) = detector9
+
+         DO idet = ndet+1, max_detectors
+            detector(idet) = 0.0D0
+         END DO
+
+         DO idet = 1, ndet
+            IF (detector(idet) > r_max) THEN
+               IF (openfiles == 1 .OR. cctk_iteration == it) THEN
+                  call CCTK_WARN(8,"Removing detectors outside coordinate range")
+c                  IF (iconv == 0) STOP
+               END IF
+               ndet =idet-1
+               GOTO 404
+            ENDIF
+         END DO
+ 404     CONTINUE
+
+         IF (verbose == 1) THEN
+            IF (openfiles == 1) THEN
+               WRITE(*,*)
+	       WRITE(*,*) "Extracting at detectors"
+               WRITE(*,*) "-----------------------"
+               WRITE(*,*) "Using ",ndet," detectors at"
+               DO idet=1,ndet
+                  WRITE(*,*) "r = ",detector(idet)
+               ENDDO
+	    END IF
+         END IF
+
+         detector_loop: DO idet = 1,ndet
+
+            radius = detector(idet)
+
+            IF (verbose == 1) THEN
+              WRITE(*,*) "Calling extract at radius ... ",radius
+            END IF
+
+            CALL D3_extract(cctkGH,do_ADMmass,do_momentum,do_spin,
+     &           igrid,orig,myproc,Nt,
+     &           Np,all_modes,lmode,mmode,x_1d,y_1d,z_1d,
+     &           Dx,Dy,Dz,Psi_power,Psi,g00,
+     &           gxx,gxy,gxz,gyy,gyz,gzz,kxx,kxy,kxz,kyy,kyz,kzz,
+     &           radius,ADMmass,momentum,spin,mass,rsch,
+     &           Qodd,Qeven,temp3d,dtaudt)
+
+            test_verbose: 
+     &           IF (verbose == 1) THEN
+
+                   WRITE(*,*)
+                   WRITE(*,*) "Detector ",idet," ..."
+                   WRITE(*,*) "Schwarzschild radius =",rsch
+                   WRITE(*,*) "Schwarzschild mass =",mass
+
+               loop_l4: DO il = lmin,lmax,lstep
+
+                  IF (all_modes == 0) THEN
+                     mmin = mmode ; mmax = mmode
+                  ELSE
+                     mmin = 0 ; mmax = il
+                  END IF
+
+                  loop_m4: DO im = mmin,mmax,mstep
+
+                     WRITE(*,*) filestr(1:nchar)//"/Qeven_"//
+     &                 CHAR(il+48)//CHAR(im+48),Qeven(:,il,im)
+                     WRITE(*,*) filestr(1:nchar)//"/Qodd_"//
+     &                 CHAR(il+48)//CHAR(im+48),Qodd(:,il,im)
+
+                  END DO loop_m4
+
+               END DO loop_l4
+
+            END IF test_verbose
+                  
+c           Output to files
+
+            test_myproc4: IF (myproc == 0) THEN
+
+               DO ioutput=1,number_timecoords
+
+                  do_output = .false.
+                  timestring = "ERR"
+                  time = 0D0
+                  IF (ioutput == 1 .AND. use_coordinate) THEN
+                     timestring = ".tl"
+                     do_output = .true.
+                     time = cctk_time
+                  ELSEIF (ioutput == 2 .AND. use_proper) THEN
+                     timestring = ".ul"
+                     do_output = .true.
+                     time = proper_time(idet) + (cctk_time-last_time)*dtaudt
+                     proper_time(idet) = time
+                     if (idet == ndet) last_time = cctk_time
+                  ENDIF
+
+                  IF (do_output) THEN
+
+c                   Output extracted radius and mass
+                     rschfile = filestr(1:nchar)//"/rsch_"//"R"//
+     &                    CHAR(idet+48)//timestring
+                     massfile = filestr(1:nchar)//"/mass_"//"R"//
+     &                    CHAR(idet+48)//timestring
+                     ADMmassfile1 = filestr(1:nchar)//"/ADMmass_"//"R"//
+     &                    CHAR(idet+48)//timestring
+                     ADMmassfile2 = filestr(1:nchar)//"/ADMmassc_"//"R"//
+     &                    CHAR(idet+48)//timestring
+                     momentumfile1 = filestr(1:nchar)//"/momentum_x_"//"R"//
+     &                    CHAR(idet+48)//timestring
+                     momentumfile2 = filestr(1:nchar)//"/momentum_y_"//"R"//
+     &                    CHAR(idet+48)//timestring
+                     momentumfile3 = filestr(1:nchar)//"/momentum_z_"//"R"//
+     &                    CHAR(idet+48)//timestring
+                     spinfile1 = filestr(1:nchar)//"/spin_x_"//"R"//
+     &                    CHAR(idet+48)//timestring
+                     spinfile2 = filestr(1:nchar)//"/spin_y_"//"R"//
+     &                    CHAR(idet+48)//timestring
+                     spinfile3 = filestr(1:nchar)//"/spin_z_"//"R"//
+     &                    CHAR(idet+48)//timestring
+                  
+                     IF (openfiles == 1) THEN
+                        OPEN(UNIT= 77,FILE=rschfile,STATUS="unknown")
+                        OPEN(UNIT= 78,FILE=massfile,STATUS="unknown")
+                        IF (do_ADMmass(1) == 1) 
+     &                       OPEN(UNIT= 79,FILE=ADMmassfile1,STATUS="unknown")
+                        IF (do_ADMmass(2) == 1) 
+     &                       OPEN(UNIT= 80,FILE=ADMmassfile2,STATUS="unknown")
+                        IF (do_momentum == 1) THEN
+                           OPEN(UNIT=781,FILE=momentumfile1,STATUS="unknown")
+                           OPEN(UNIT=782,FILE=momentumfile2,STATUS="unknown")
+                           OPEN(UNIT=783,FILE=momentumfile3,STATUS="unknown")
+                        END IF
+                        
+                        IF (do_spin == 1) THEN
+                           OPEN(UNIT=784,FILE=spinfile1,STATUS="unknown")
+                           OPEN(UNIT=785,FILE=spinfile2,STATUS="unknown")
+                           OPEN(UNIT=786,FILE=spinfile3,STATUS="unknown")
+                        END IF
+                     ELSE
+                        OPEN(UNIT= 77,FILE=rschfile,STATUS="old",
+     &                       POSITION="append")
+                        OPEN(UNIT= 78,FILE=massfile,STATUS="old",
+     &                       POSITION="append")
+                        IF (do_ADMmass(1) == 1)
+     &                       OPEN(UNIT= 79,FILE=ADMmassfile1,STATUS="old",
+     &                       POSITION="append")
+                        IF (do_ADMmass(2) == 1)
+     &                       OPEN(UNIT= 80,FILE=ADMmassfile2,STATUS="old",
+     &                       POSITION="append")
+                        IF (do_momentum == 1) THEN
+                           OPEN(UNIT=781,FILE=momentumfile1,STATUS="old",
+     &                          POSITION="append")
+                           OPEN(UNIT=782,FILE=momentumfile2,STATUS="old",
+     &                          POSITION="append")
+                           OPEN(UNIT=783,FILE=momentumfile3,STATUS="old",
+     &                          POSITION="append")
+                        END IF
+                        
+                        IF (do_spin == 1)THEN
+                           OPEN(UNIT=784,FILE=spinfile1,STATUS="old",
+     &                          POSITION="append")
+                           OPEN(UNIT=785,FILE=spinfile2,STATUS="old",
+     &                          POSITION="append")
+                           OPEN(UNIT=786,FILE=spinfile3,STATUS="old",
+     &                          POSITION="append")
+                        END IF
+                        
+                     ENDIF
+                     WRITE( 77,21) time,rsch
+                     WRITE( 78,21) time,mass
+                     IF (do_ADMmass(1) == 1) WRITE( 79,21) time,ADMmass(1)
+                     IF (do_ADMmass(2) == 1) WRITE( 80,21) time,ADMmass(2)
+                     IF (do_momentum == 1) THEN
+                        WRITE(781,*) cctk_time,momentum(1)
+                        WRITE(782,*) cctk_time,momentum(2)
+                        WRITE(783,*) cctk_time,momentum(3)
+                     END IF
+                     IF (do_spin == 1) THEN
+                        WRITE(784,*) cctk_time,spin(1)
+                        WRITE(785,*) cctk_time,spin(2)
+                        WRITE(786,*) cctk_time,spin(3)
+                     END IF
+
+                     CLOSE( 77)
+                     CLOSE( 78)
+                     IF (do_ADMmass(1) == 1) CLOSE( 79)
+                     IF (do_ADMmass(2) == 1) CLOSE( 80)                     
+
+                     IF (do_momentum == 1) THEN
+                        CLOSE(781)
+                        CLOSE(782)
+                        CLOSE(783)
+                     END IF
+                     IF (do_spin == 1) THEN
+                        CLOSE(782)
+                        CLOSE(783)
+                        CLOSE(784)
+                     END IF
+
+
+c                    Output gauge invariant variables
+
+                     loop_l5: DO il = lmin,lmax,lstep
+                     
+                       IF (all_modes == 0) THEN
+                          mmin = mmode ; mmax = mmode
+                       ELSE
+                          mmin = 0 ; mmax = il
+                       END IF
+                       
+                       loop_m5: DO im = mmin,mmax,mstep
+
+                         fn1 = 11
+                         fn2 = 12
+                         
+                         out1 = filestr(1:nchar)//"/Qeven_"//"R"
+     &                        //CHAR(idet+48)//"_"//CHAR(il+48)//
+     &                        CHAR(im+48)//timestring
+                         out2 = filestr(1:nchar)//"/Qodd_"//"R"
+     &                        //CHAR(idet+48)//"_"//CHAR(il+48)//
+     &                        CHAR(im+48)//timestring
+
+c                        Write even parity waveforms
+                         IF (openfiles == 1) THEN
+                            OPEN(UNIT=fn1,FILE=out1,STATUS="unknown")
+                         ELSE
+                            OPEN(UNIT=fn1,FILE=out1,STATUS="old",
+     &                           POSITION="append")
+                         ENDIF
+                         
+                         WRITE(fn1,20) time,Qeven(:,il,im)
+
+                         CLOSE(fn1)
+
+c                        Only write odd parity waveforms if full grid
+                         IF (igrid == 0) THEN
+                            IF (openfiles == 1) THEN
+                               OPEN(UNIT=fn2,FILE=out2,STATUS="unknown")
+                            ELSE
+                               OPEN(UNIT=fn2,FILE=out2,STATUS="old",
+     &                              POSITION="append")
+                            ENDIF
+                        
+                            WRITE(fn2,20) time,Qodd(:,il,im)
+                        
+                            CLOSE(fn2)
+                         END IF
+                         
+                      END DO loop_m5
+
+                   END DO loop_l5
+         
+                ENDIF
+             END DO
+          END IF test_myproc4
+
+       END DO detector_loop
+       
+      END IF test_detectors
+      
+      DEALLOCATE(Qodd,Qeven)
+
+      openfiles = 0
+
+ 20   format(e24.15,e24.15,e24.15)
+ 21   format(e24.15,e24.15)
+
+      END SUBROUTINE Extract
diff --git a/src/cartesian_to_spherical.F b/src/cartesian_to_spherical.F
new file mode 100644
index 0000000..d55f84f
--- /dev/null
+++ b/src/cartesian_to_spherical.F
@@ -0,0 +1,116 @@
+
+#include "cctk.h"
+
+c     ==================================================================
+
+      SUBROUTINE cartesian_to_spherical(theta,phi,r,gxxs,gxys,gxzs,
+     &     gyys,gyzs, gzzs,grrs,grts,grps,gtts,gtps,gpps,dgxxs,dgxys,
+     &     dgxzs,dgyys,dgyzs,dgzzs,dgtts,dgtps,dgpps)
+
+c     ------------------------------------------------------------------
+c
+c     Convert components of a symmetric 2nd rank 3D tensor on a 2-sphere 
+c     from Cartesian to spherical polar coordinates
+c
+c                (x,y,z)        ---->         (r,t,p)
+c
+c           ( gxx  gxy gxz )             ( grr  grt grp )
+c           (  .   gyy gyz )    ---->    (  .   gpp gtp )   
+c           (  .    .  gzz )             (  .    .  gpp ) 
+c
+c     Also convert radial derivatives of the tensor components
+c
+c     ------------------------------------------------------------------
+
+      IMPLICIT NONE
+
+c     Input variables
+      CCTK_REAL,INTENT(IN) ::  
+     &    r
+      CCTK_REAL,INTENT(IN),DIMENSION (:) ::  
+     &    theta,phi
+      CCTK_REAL,INTENT(IN),DIMENSION (:,:) :: 
+     &    gxxs,gxys,gxzs,gyys,gyzs,gzzs,dgxxs,dgxys,dgxzs,dgyys,dgyzs,
+     &    dgzzs
+
+c     Output variables
+      CCTK_REAL,INTENT(OUT),DIMENSION (:,:) :: 
+     &    grrs,grts,grps,gtts,gtps,
+     &    gpps,dgtts,dgtps,dgpps
+
+c     Local variables
+      INTEGER ::
+     &    it,ip,Nt,Np
+      CCTK_REAL :: 
+     &    ct,st,cp,sp,ct2,st2,cp2,sp2,r2,gxx,gxy,gxz,gyy,gyz,gzz,dgxx,
+     &    dgxy,dgxz,dgyy,dgyz,dgzz
+      CCTK_REAL,PARAMETER ::
+     &    two = 2.0D0
+
+c     ------------------------------------------------------------------
+
+      Nt = SIZE(theta)
+      Np = SIZE(phi)
+
+      r2 = r*r
+
+      DO it = 1, Nt
+ 
+	ct = COS(theta(it)) ; ct2 = ct*ct
+        st = SIN(theta(it)) ; st2 = st*st
+
+          DO ip = 1, Np
+
+          cp = COS(phi(ip)) ; cp2 = cp*cp
+          sp = SIN(phi(ip)) ; sp2 = sp*sp
+
+	  gxx = gxxs(it,ip) ; dgxx = dgxxs(it,ip)
+	  gxy = gxys(it,ip) ; dgxy = dgxys(it,ip)
+	  gxz = gxzs(it,ip) ; dgxz = dgxzs(it,ip)
+          gyy = gyys(it,ip) ; dgyy = dgyys(it,ip)
+          gyz = gyzs(it,ip) ; dgyz = dgyzs(it,ip)
+          gzz = gzzs(it,ip) ; dgzz = dgzzs(it,ip)
+
+          grrs(it,ip) = st2*cp2*gxx+st2*sp2*gyy+ct2*gzz
+     &       +two*st2*cp*sp*gxy+two*st*cp*ct*gxz+two*st*ct*sp*gyz
+
+          grts(it,ip) = r*(st*cp2*ct*gxx+two*st*ct*sp*cp*gxy
+     &       +cp*(ct2-st2)*gxz+st*sp2*ct*gyy+sp*(ct2-st2)*gyz-ct*st*gzz)
+
+          grps(it,ip) = r*st*(-st*sp*cp*gxx-st*(sp2-cp2)*gxy-sp*ct*gxz
+     &       +st*sp*cp*gyy+ct*cp*gyz)
+     
+          gtts(it,ip) = r2*(ct2*cp2*gxx+two*ct2*sp*cp*gxy
+     &       -two*st*ct*cp*gxz+ct2*sp2*gyy-two*st*sp*ct*gyz+st2*gzz)
+     
+          gtps(it,ip) = r2*st*(-cp*sp*ct*gxx-ct*(sp2-cp2)*gxy
+     &       +st*sp*gxz+cp*sp*ct*gyy-st*cp*gyz)
+     
+          gpps(it,ip) = r2*st2*(sp2*gxx-two*cp*sp*gxy+cp2*gyy)
+     
+          dgtts(it,ip) = two/r*gtts(it,ip)+r2*(ct2*cp2*dgxx
+     &       +two*ct2*sp*cp*dgxy-two*st*ct*cp*dgxz+ct2*sp2*dgyy
+     &       -two*st*sp*ct*dgyz+st2*dgzz)
+
+          dgtps(it,ip) = two/r*gtps(it,ip)+r2*st*(-cp*sp*ct*dgxx
+     &       -ct*(sp2-cp2)*dgxy+st*sp*dgxz+cp*sp*ct*dgyy-st*cp*dgyz)
+     
+          dgpps(it,ip) = two/r*gpps(it,ip) +r2*st2*(sp2*dgxx
+     &       -two*cp*sp*dgxy+cp2*dgyy)
+
+        END DO
+	  				 
+      END DO
+
+
+      END SUBROUTINE cartesian_to_spherical
+
+
+
+
+
+
+
+
+
+
diff --git a/src/cartesian_to_spherical_int.F b/src/cartesian_to_spherical_int.F
new file mode 100644
index 0000000..a4db3b2
--- /dev/null
+++ b/src/cartesian_to_spherical_int.F
@@ -0,0 +1,24 @@
+
+#include "cctk.h"
+
+      MODULE cartesian_to_spherical_int
+
+c     ------------------------------------------------------------------
+
+      INTERFACE
+
+      SUBROUTINE cartesian_to_spherical(theta,phi,rad,gxxs,gxys,gxzs,
+     &     gyys,gyzs,gzzs,grrs,grts,grps,gtts,gtps,gpps,dgxxs,dgxys,
+     &     dgxzs,dgyys,dgyzs,dgzzs,dgtts,dgtps,dgpps)
+      IMPLICIT NONE
+      CCTK_REAL  :: 
+     &     rad,theta(:),phi(:)
+      CCTK_REAL, DIMENSION (:,:) :: 
+     &     gxxs,gxys,gxzs,gyys,gyzs,gzzs,dgxxs,dgxys,dgxzs,dgyys,
+     &     dgyzs,dgzzs,grrs,grts,grps,gtts,gtps,gpps,dgtts,dgtps,
+     &     dgpps
+      END SUBROUTINE
+
+      END INTERFACE
+
+      END MODULE cartesian_to_spherical_int
diff --git a/src/energy.f b/src/energy.f
new file mode 100644
index 0000000..19ef5e2
--- /dev/null
+++ b/src/energy.f
@@ -0,0 +1,205 @@
+      PROGRAM energy
+
+c     -----------------------------------------------------------------
+c
+c     Calculate energy flux from even parity variable (only real part)
+c
+c     Now copes with having unequally spaced timesteps, it works out
+c     the smallest timestep, and uses this to interpolate to a new
+c     array with equally spaced timesteps.
+c
+c     -----------------------------------------------------------------
+
+      IMPLICIT NONE
+
+      INTEGER,PARAMETER :: nmax=5000
+      INTEGER :: i,itotal,it1,it2,Nt
+      CHARACTER*50 filename
+      REAL*8 t(nmax),Qplus(nmax),dEdt(nmax),E,ddQ(nmax),Qnew(nmax)
+      REAL*8 Dtnew,fac,time,t1,t2,dQdt_l,dQdt_r,dtmin,dt
+     
+      WRITE(6,*)
+      WRITE(6,*) "Starting PROGRAM energy ..."
+      WRITE(6,*) "---------------------------"
+
+      WRITE(6,131) 
+ 131  FORMAT("Filename containing Q ? ",$)
+      READ(5,*) filename
+      WRITE(6,*) "Filename is ... ",filename
+
+      OPEN(UNIT=12,FILE=filename,STATUS="old")
+
+      DO i=1,nmax
+         READ(12,*,ERR=101) t(i),Qplus(i)
+         itotal = i
+      END DO
+      WRITE(6,*) "Didn't read all of data file"
+      STOP
+
+ 101  CLOSE(12) 
+
+c     Give information about the data
+      WRITE(6,*) "We have data from t=",t(1)," to ",t(itotal)
+
+c     Find minimum time step
+      dtmin = 100000.0
+      do i=2,itotal
+         dt = t(i)-t(i-1)
+         if (dt<dtmin) dtmin = dt
+      end do
+      WRITE(6,*) "Minimum timestep is ",dtmin
+
+c     Initial time for energy flux
+ 201  WRITE(6,132) 
+ 132  FORMAT(" Initial time ? ",$)
+      READ(5,*) t1
+      IF (t1 < t(1)) GOTO 201
+     
+c     Calculate initial timeslice
+      do i=1,itotal
+         if (t(i)<t1) it1=i
+      end do
+      WRITE(6,*) "Initial timeslice is ...",it1
+
+c     Final time for energy flux
+      WRITE(6,133) 
+ 133  FORMAT(" Final time ? ",$)
+ 301  READ(5,*) t2
+      IF (t2 > t(itotal)) GOTO 301
+
+c     Calculate final timeslice
+      do i=1,itotal
+         if (t(i)<t2) it2=i
+      end do
+      WRITE(6,*) "Final timeslice is ...",it2
+
+
+      Nt = (t2-t1)/dtmin
+      WRITE(6,*) "Interpolating to ",Nt," timesteps"
+      
+
+c     Do a spline
+      Dtnew = Dtmin
+      dQdt_l = ( Qplus(2)             - Qplus(1)                ) / 
+     &           ( t(2)       - t(1)          )
+      dQdt_r = (Qplus(itotal)       - Qplus(itotal-1)       ) /
+     &           (t(itotal) - t(itotal-1)) 
+      CALL SPLINE(t,Qplus,itotal,dQdt_l,dQdt_r,ddQ)  
+      DO i = 1,Nt 
+        time = t1+DBLE(i-1)*Dtnew
+        CALL SPLINT(t,Qplus,ddQ,itotal,time,Qnew(i))
+      ENDDO
+
+c     This is the normalisation used
+      fac = 1.0D0/32.0D0/ACOS(-1.0D0)
+
+c     Calculate energy flux
+      dEdt(1) = fac*((Qnew(2)-Qnew(1))/Dtnew)**2
+      DO i = 2,Nt-1
+         dEdt(i) = fac*((Qnew(i+1)-Qnew(i-1))*0.5D0/Dtnew)**2
+      END DO
+      dEdt(Nt) = fac*((Qnew(Nt)-Qnew(Nt-1))/Dtnew)**2
+
+c     Write results to file
+      OPEN(UNIT=101,FILE="dEdt",STATUS="unknown")
+      DO i = 1,Nt
+         time = t1+DBLE(i-1)*Dtnew
+         WRITE(101,*) time,dEdt(i)
+      END DO
+      CLOSE(101)
+
+c     Use Simpsons rule to calculate integral
+      IF ( MOD(Nt,2) /= 0) THEN
+         WRITE(6,*) "Using ",Nt," timeslices. Removing one" 
+         WRITE(6,*) "to give Simpsons rule an odd number ..."
+      END IF
+
+      E = dEdt(1) + 4.0D0*dEdt(Nt-1)+dEdt(Nt)
+      DO i=2,Nt-3,2
+         E = E + 4.0D0*dEdt(i) + 2.0D0*dEdt(i+1)
+      ENDDO
+      E = E*Dtnew/3.0D0
+
+c     Write out total energy
+      WRITE(6,*)
+      WRITE(6,*) "Total energy in this mode is : ",E
+      WRITE(6,*) "****************************"
+
+      END PROGRAM energy
+
+
+
+
+      SUBROUTINE spline(x,y,n,yp1,ypn,y2)
+
+      INTEGER, PARAMETER :: nmax=5000
+
+      INTEGER   n,i,k
+      REAL*8    yp1,ypn,x(nmax),y(nmax),y2(nmax)
+      REAL*8    p,qn,sig,un,u(nmax)
+
+      IF (yp1.GT..99d30) THEN
+        y2(1)=0.d0
+        u(1)=0.d0
+      ELSE
+        y2(1)=-0.5d0
+        u(1)=(3.d0/(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)
+      ENDIF
+
+      DO i=2,n-1
+        sig   = (x(i)-x(i-1))/(x(i+1)-x(i-1))
+        p     = sig*y2(i-1)+2.d0
+        y2(i) = (sig-1.d0)/p
+        u(i)  = (6.d0*((y(i+1)-y(i))/(x(i+1)-x(i))-(y(i)-y(i-1))
+     &          /(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*u(i-1))/p
+      ENDDO
+
+      IF (ypn.GT..99d30) then
+        qn=0.d0
+        un=0.d0
+      ELSE
+        qn=0.5d0
+        un=(3.d0/(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))
+      ENDIF
+
+      y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.d0)
+
+      DO k=n-1,1,-1
+        y2(k)=y2(k)*y2(k+1)+u(k)
+      ENDDO
+
+      RETURN
+      END
+
+
+
+      SUBROUTINE splint(xa,ya,y2a,n,x,y)
+ 
+      INTEGER n 
+      REAL*8  x,y,xa(n),y2a(n),ya(n) 
+      INTEGER k,khi,klo 
+      REAL*8  a,b,h 
+
+      klo=1 
+      khi=n 
+1     if (khi-klo.gt.1) then 
+        k=(khi+klo)/2 
+        if(xa(k).gt.x)then 
+          khi=k 
+        else 
+          klo=k 
+        endif 
+      goto 1 
+      endif 
+      h=xa(khi)-xa(klo) 
+      if (h.eq.0.d0) pause 'bad xa input in splint' 
+      a=(xa(khi)-x)/h 
+      b=(x-xa(klo))/h 
+      y=a*ya(klo)+b*ya(khi)+((a**3-a)*y2a(klo)+(b**3-b)*y2a(khi))*(h** 
+     *2)/6.d0 
+
+      RETURN 
+      END 
+
+
+
diff --git a/src/make.code.defn b/src/make.code.defn
new file mode 100644
index 0000000..a0b7189
--- /dev/null
+++ b/src/make.code.defn
@@ -0,0 +1,27 @@
+# Main make.code.defn file for thorn Extract
+# $Header$
+
+# Source files in this directory
+SRCS =  D2_extract_int.F \
+	D3_extract_int.F \
+	D3_to_D2_int.F \
+	cartesian_to_spherical_int.F \
+	unphysical_to_physical_int.F \
+	ADMmass_integrand3D_int.F \
+	momentum_integrand3D_int.F\
+	spin_integrand3D_int.F\
+	met_rad_der_int.F \
+	Extract.F \
+	D3_extract.F \
+	cartesian_to_spherical.F \
+	unphysical_to_physical.F \
+	D2_extract.F \
+	met_rad_der.F \
+	D3_to_D2.F \
+	ADMmass_integrand3D.F \
+	momentum_integrand3D.F \
+	spin_integrand3D.F 
+
+# Subdirectories containing source files
+SUBDIRS = 
+
diff --git a/src/make.code.deps b/src/make.code.deps
new file mode 100644
index 0000000..126c2e9
--- /dev/null
+++ b/src/make.code.deps
@@ -0,0 +1,2 @@
+Extract.o:	D3_extract_int.o
+ADMmass_integrand3D.o: ADMmass_integrand_int3D.o
diff --git a/src/met_rad_der.F b/src/met_rad_der.F
new file mode 100644
index 0000000..31821f0
--- /dev/null
+++ b/src/met_rad_der.F
@@ -0,0 +1,91 @@
+
+#include "cctk.h"
+
+c     ========================================================================
+
+      SUBROUTINE met_rad_der(origin,Dx,Dy,Dz,x,y,z,g,dg)
+
+c     ------------------------------------------------------------------------
+c
+c     Calculate isotropic radial derivative of unphysical metric,
+c     component or conformal factor at each point on the 3D grid, 
+c     using derivatives in the Cartesian directions.
+c
+c     ------------------------------------------------------------------------
+      
+      IMPLICIT NONE
+
+c     Input variables
+      CCTK_REAL,INTENT(IN) :: 
+     &     Dx,Dy,Dz,origin(3)
+      CCTK_REAL,DIMENSION(:),INTENT(IN) :: 
+     &     x,y,z
+      CCTK_REAL,DIMENSION(:,:,:),INTENT(IN) :: 
+     &     g
+      
+c     Output variables
+      CCTK_REAL,DIMENSION(:,:,:),INTENT(OUT) :: 
+     &     dg
+      
+c     Local variables, here only
+      INTEGER :: 
+     &     i,j,k
+      CCTK_REAL,PARAMETER ::
+     &     zero = 0.0D0,
+     &     half = 0.5D0
+      CCTK_REAL :: 
+     &     rad
+
+c     ------------------------------------------------------------------------
+
+      DO k = 2, SIZE(z)-1
+         DO j = 2, SIZE(y)-1
+            DO i = 2, SIZE(x)-1
+               
+               rad = SQRT((x(i)-origin(1))**2 
+     &              +(y(j)-origin(2))**2
+     &              +(z(k)-origin(3))**2)
+               
+               IF (rad.NE.0) THEN
+                  
+                  dg(i,j,k) = half/rad*
+     &                 ((x(i)-origin(1))/Dx*(g(i+1,j,k)-g(i-1,j,k))
+     &                 +(y(j)-origin(2))/Dy*(g(i,j+1,k)-g(i,j-1,k))
+     &                 +(z(k)-origin(3))/Dz*(g(i,j,k+1)-g(i,j,k-1)))
+                  
+               ELSE
+                  
+                  dg(i,j,k) = zero
+                  
+               ENDIF
+
+          ENDDO
+        ENDDO
+      ENDDO
+
+c     This is needed when the grid is an octant, but it does not hurt
+c     if it is not
+
+      DO k = 2, size(z)-1
+         DO j = 2, size(y)-1
+            dg(1,j,k) = dg(2,j,k)
+         ENDDO
+      ENDDO
+      DO k = 2, size(z)-1
+         DO i = 1, size(x)-1
+            dg(i,1,k) = dg(i,2,k)
+         ENDDO
+      ENDDO
+      DO j = 1, size(y)-1
+         DO i = 1, size(x)-1
+            dg(i,j,1) = dg(i,j,2)
+         ENDDO
+      ENDDO
+
+      END SUBROUTINE met_rad_der
+
+
+
+
+
+
diff --git a/src/met_rad_der_int.F b/src/met_rad_der_int.F
new file mode 100644
index 0000000..d7b2f32
--- /dev/null
+++ b/src/met_rad_der_int.F
@@ -0,0 +1,24 @@
+
+#include "cctk.h"
+ 
+      MODULE met_rad_der_int
+
+c     ------------------------------------------------------------------
+
+      INTERFACE
+
+          SUBROUTINE met_rad_der(origin,Dx,Dy,Dz,x,y,z,g,dg)
+          IMPLICIT NONE
+          CCTK_REAL,INTENT(IN) :: 
+     &        Dx,Dy,Dz,origin(3)
+          CCTK_REAL,DIMENSION(:),INTENT(IN) :: 
+     &        x,y,z
+          CCTK_REAL,DIMENSION(:,:,:),INTENT(IN) :: 
+     &        g
+          CCTK_REAL,DIMENSION(:,:,:),INTENT(OUT) :: 
+     &        dg
+          END SUBROUTINE met_rad_der
+
+      END INTERFACE
+
+      END MODULE met_rad_der_int
diff --git a/src/momentum_integrand3D.F b/src/momentum_integrand3D.F
new file mode 100644
index 0000000..011fa68
--- /dev/null
+++ b/src/momentum_integrand3D.F
@@ -0,0 +1,188 @@
+
+#include "cctk.h"
+
+c     ========================================================================
+
+      SUBROUTINE momentum_integrand3D(origin,Dx,Dy,Dz,x,y,z,gxx,gxy,gxz,
+     &     gyy,gyz,gzz,hxx,hxy,hxz,hyy,hyz,hzz,
+     &     momentum_int,Psi,Psi_power)
+
+c     ------------------------------------------------------------------------
+c
+c     Estimates the momentum at a given radius using Equation (11.2.14) from
+c     Walds General Relativity
+c
+c     -----------------------------------------------------------------------
+      
+      IMPLICIT NONE
+
+c     Input variables
+      INTEGER,INTENT(IN) ::
+     &     Psi_power
+      CCTK_REAL,INTENT(IN) :: 
+     &     Dx,Dy,Dz,origin(3)
+      CCTK_REAL,DIMENSION(:),INTENT(IN) :: 
+     &     x,y,z
+      CCTK_REAL,DIMENSION(:,:,:),INTENT(IN) :: 
+     &     gxx,gxy,gxz,gyy,gyz,gzz,Psi,
+     &     hxx,hxy,hxz,hyy,hyz,hzz
+
+      
+c     Output variables
+      CCTK_REAL,DIMENSION(:,:,:),INTENT(OUT) :: 
+     &     momentum_int
+      
+c     Local variables, here only
+      INTEGER :: 
+     &     i,j,k,ip,count
+
+      CCTK_REAL,PARAMETER ::
+     &     half = 0.5D0
+      CCTK_REAL :: 
+     &     rad,ux,uy,uz,det,dxx,dxy,dxz,dyy,dyz,dzz,uxx,uxy,uxz,uyy,
+     &     uyz,uzz,
+     &     Pi,idet,p,term1,
+     &     tracek
+
+      data count / 1 /
+      save count
+c     ------------------------------------------------------------------------
+
+      Pi = ACOS(-1D0)
+
+c     Because other codes evolve Psi**4
+      SELECT CASE (Psi_power)
+
+      CASE (1)
+         ip = 4
+
+      CASE (4)
+         ip = 1
+
+      CASE DEFAULT
+         WRITE(*,*) "This value of Psi_power is not supported"
+
+      END SELECT   
+      
+
+      DO k = 2, SIZE(z)-1
+         DO j = 2, SIZE(y)-1
+            DO i = 2, SIZE(x)-1
+               
+               rad = SQRT((x(i)-origin(1))**2 
+     &              +(y(j)-origin(2))**2
+     &              +(z(k)-origin(3))**2)
+               
+               IF (rad.NE.0) THEN
+
+                  ux = (x(i)-origin(1))/rad
+                  uy = (y(j)-origin(2))/rad
+                  uz = (z(k)-origin(3))/rad
+
+c                 Abbreviations for metric coefficients
+c                 -------------------------------------
+                  p = psi(i,j,k)**ip
+
+                  dxx = p*gxx(i,j,k); dxy = p*gxy(i,j,k)
+                  dxz = p*gxz(i,j,k); dyy = p*gyy(i,j,k)
+                  dyz = p*gyz(i,j,k); dzz = p*gzz(i,j,k)
+
+c                 Determinant of 3-metric
+c                 -----------------------
+                  det = (dxx*dyy*dzz + 2.0D0*dxy*dxz*dyz
+     &                 - (dxx*dyz**2 + dyy*dxz**2 + dzz*dxy**2))
+
+                  idet = 1.0/det
+
+c                 Inverse 3-metric
+c                 ----------------
+                  uxx = idet*(dyy*dzz - dyz**2)
+                  uyy = idet*(dxx*dzz - dxz**2)
+                  uzz = idet*(dxx*dyy - dxy**2)
+                  uxy = idet*(dxz*dyz - dzz*dxy)
+                  uxz = idet*(dxy*dyz - dyy*dxz)
+                  uyz = idet*(dxy*dxz - dxx*dyz)
+
+
+c                 Trace of extrinsic curvature
+c                 ----------------------------
+
+                  tracek = uxx*hxx(i,j,k) + uyy*hyy(i,j,k) +
+     &                 uzz*hzz(i,j,k)
+
+                  
+
+c                 Integrands
+c                 ----------
+
+                  if (count.eq.1) then
+                     term1 = ux*hxx(i,j,k)+uy*hxy(i,j,k)+uz*hxz(i,j,k)
+
+c                     momentum_int(i,j,k) = 1.0D0/8.0D0/Pi*
+c     &                    (term1 - tracek*ux)*rad**2
+                     momentum_int(i,j,k) = 1.0D0/8.0D0/Pi*
+     &                    (term1 - tracek*(dxx*ux+dxy*uy+dxz*uz))*rad**2
+
+                  else if (count.eq.2) then
+                     term1 = ux*hxy(i,j,k)+uy*hyy(i,j,k)+uz*hyz(i,j,k)
+
+c                     momentum_int(i,j,k) = 1.0D0/8.0D0/Pi*
+c     &                    (term1 - tracek*uy)*rad**2
+                     momentum_int(i,j,k) = 1.0D0/8.0D0/Pi*
+     &                    (term1 - tracek*(dxy*ux+dyy*uy+dyz*uz))*rad**2
+
+                  else if (count.eq.3) then
+                     term1 = ux*hxz(i,j,k)+uy*hyz(i,j,k)+uz*hzz(i,j,k)
+
+c                     momentum_int(i,j,k) = 1.0D0/8.0D0/Pi*
+c     &                    (term1 - tracek*uz)*rad**2
+                      momentum_int(i,j,k) = 1.0D0/8.0D0/Pi*
+     &                    (term1 - tracek*(dxz*ux+dyz*uy+dzz*uz))*rad**2
+                   end if
+
+
+               ELSE
+                  
+                  momentum_int(i,j,k) = 0.0D0
+                  
+               ENDIF
+               
+          ENDDO
+        ENDDO
+      ENDDO
+
+
+c     This is needed when the grid is an octant, but it does not hurt
+c     if it is not
+
+      DO k = 2, size(z)-1
+         DO j = 2, size(y)-1
+            momentum_int(1,j,k) = momentum_int(2,j,k)
+         ENDDO
+      ENDDO
+      DO k = 2, size(z)-1
+         DO i = 1, size(x)-1
+            momentum_int(i,1,k) = momentum_int(i,2,k)
+         ENDDO
+      ENDDO
+      DO j = 1, size(y)-1
+         DO i = 1, size(x)-1
+            momentum_int(i,j,1) = momentum_int(i,j,2)
+         ENDDO
+      ENDDO
+
+
+c     setting counter
+c     ---------------
+
+      if (count.eq.3) then
+         count = 1
+      else
+         count = count + 1
+      end if
+
+
+c     end
+c     ---
+
+      END SUBROUTINE momentum_integrand3D
diff --git a/src/momentum_integrand3D_int.F b/src/momentum_integrand3D_int.F
new file mode 100644
index 0000000..e97174f
--- /dev/null
+++ b/src/momentum_integrand3D_int.F
@@ -0,0 +1,30 @@
+
+#include "cctk.h"
+ 
+      MODULE momentum_integrand3D_int
+
+c     ------------------------------------------------------------------
+
+      INTERFACE
+
+          SUBROUTINE  momentum_integrand3D(origin,Dx,Dy,Dz,x,y,z,gxx,
+     &        gxy,gxz,gyy,gyz,gzz,hxx,hxy,hxz,hyy,hyz,hzz,momentum_int,
+     &        Psi,Psi_power)
+
+          IMPLICIT NONE
+          INTEGER,INTENT(IN) ::
+     &         Psi_power
+          CCTK_REAL,INTENT(IN) :: 
+     &        Dx,Dy,Dz,origin(3)
+          CCTK_REAL,DIMENSION(:),INTENT(IN) :: 
+     &        x,y,z
+          CCTK_REAL,DIMENSION(:,:,:),INTENT(IN) :: 
+     &         gxx,gxy,gxz,gyy,gyz,gzz,Psi,
+     &         hxx,hxy,hxz,hyy,hyz,hzz
+          CCTK_REAL,DIMENSION(:,:,:),INTENT(OUT) :: 
+     &        momentum_int
+          END SUBROUTINE  momentum_integrand3D
+
+      END INTERFACE
+
+      END MODULE momentum_integrand3D_int
diff --git a/src/spin_integrand3D.F b/src/spin_integrand3D.F
new file mode 100644
index 0000000..3f1afc3
--- /dev/null
+++ b/src/spin_integrand3D.F
@@ -0,0 +1,193 @@
+
+#include "cctk.h"
+
+c     ========================================================================
+
+      SUBROUTINE spin_integrand3D(origin,Dx,Dy,Dz,x,y,z,gxx,gxy,gxz,
+     &     gyy,gyz,gzz,hxx,hxy,hxz,hyy,hyz,hzz,
+     &     spin_int,Psi,Psi_power)
+
+c     ------------------------------------------------------------------------
+c
+c     Estimates the spin at a given radius using Equation (12) from
+c     <Time-asymmetric initial data for black holes and black-hole 
+c     collisions>, Bowen and York, Phys. Rev. D, 21 2047(1980)
+c
+c     -----------------------------------------------------------------------
+      
+      IMPLICIT NONE
+
+c     Input variables
+      INTEGER,INTENT(IN) ::
+     &     Psi_power
+      CCTK_REAL,INTENT(IN) :: 
+     &     Dx,Dy,Dz,origin(3)
+      CCTK_REAL,DIMENSION(:),INTENT(IN) :: 
+     &     x,y,z
+      CCTK_REAL,DIMENSION(:,:,:),INTENT(IN) :: 
+     &     gxx,gxy,gxz,gyy,gyz,gzz,Psi,
+     &     hxx,hxy,hxz,hyy,hyz,hzz
+
+      
+c     Output variables
+      CCTK_REAL,DIMENSION(:,:,:),INTENT(OUT) :: 
+     &     spin_int
+      
+c     Local variables, here only
+      INTEGER :: 
+     &     i,j,k,ip,count
+
+      CCTK_REAL,PARAMETER ::
+     &     half = 0.5D0
+      CCTK_REAL :: 
+     &     rad,ux,uy,uz,det,dxx,dxy,dxz,dyy,dyz,dzz,uxx,uxy,uxz,uyy,
+     &     uyz,uzz,
+     &     Pi,idet,p,
+     &     term1,term2
+
+      data count / 1 /
+      save count
+c     ------------------------------------------------------------------------
+
+      Pi = ACOS(-1D0)
+
+      SELECT CASE (Psi_power)
+
+      CASE (1)
+         ip = 4
+
+      CASE (4)
+         ip = 1
+
+      CASE DEFAULT
+         WRITE(*,*) "This value of Psi_power is not supported"
+
+      END SELECT   
+      
+
+      DO k = 2, SIZE(z)-1
+         DO j = 2, SIZE(y)-1
+            DO i = 2, SIZE(x)-1
+               
+               rad = SQRT((x(i)-origin(1))**2 
+     &              +(y(j)-origin(2))**2
+     &              +(z(k)-origin(3))**2)
+               
+               IF (rad.NE.0) THEN
+
+                  ux = (x(i)-origin(1))/rad
+                  uy = (y(j)-origin(2))/rad
+                  uz = (z(k)-origin(3))/rad
+
+c                 Abbreviations for metric coefficients
+c                 -------------------------------------
+                  p = psi(i,j,k)**ip
+
+                  dxx = p*gxx(i,j,k); dxy = p*gxy(i,j,k)
+                  dxz = p*gxz(i,j,k); dyy = p*gyy(i,j,k)
+                  dyz = p*gyz(i,j,k); dzz = p*gzz(i,j,k)
+
+c                 Determinant of 3-metric
+c                 -----------------------
+                  det = (dxx*dyy*dzz + 2.0D0*dxy*dxz*dyz
+     &                 - (dxx*dyz**2 + dyy*dxz**2 + dzz*dxy**2))
+
+                  idet = 1.0/det
+
+c                 Inverse 3-metric
+c                 ----------------
+                  uxx = idet*(dyy*dzz - dyz**2)
+                  uyy = idet*(dxx*dzz - dxz**2)
+                  uzz = idet*(dxx*dyy - dxy**2)
+                  uxy = idet*(dxz*dyz - dzz*dxy)
+                  uxz = idet*(dxy*dyz - dyy*dxz)
+                  uyz = idet*(dxy*dxz - dxx*dyz)
+
+
+c                 Integrands
+c                 ----------
+
+                  if (count.eq.1) then
+              term1 = uxz*(ux*hxx(i,j,k)+uy*hxy(i,j,k)+uz*hxz(i,j,k))
+     &                +uyz*(ux*hxy(i,j,k)+uy*hyy(i,j,k))
+     &                +uzz*(ux*hxz(i,j,k)+uy*hyz(i,j,k)+uz*hzz(i,j,k))
+
+              term2 = uxy*(ux*hxx(i,j,k)+uy*hxy(i,j,k)+uz*hxz(i,j,k))
+     &                +uyy*(ux*hxy(i,j,k)+uy*hyy(i,j,k)+uz*hyz(i,j,k))
+     &                +uyz*(ux*hxz(i,j,k)+uz*hzz(i,j,k))
+
+                     spin_int(i,j,k) = 1.0D0/8.0D0/Pi*
+     &                    (uy*term1-uz*term2)*rad**3
+
+                  else if (count.eq.2) then
+              term1 = uxx*(ux*hxx(i,j,k)+uy*hxy(i,j,k)+uz*hxz(i,j,k))
+     &                +uxy*(ux*hxy(i,j,k)+uy*hyy(i,j,k)+uz*hyz(i,j,k))
+     &                +uxz*(uy*hyz(i,j,k)+uz*hzz(i,j,k))
+
+              term2 = uxz*(ux*hxx(i,j,k)+uy*hxy(i,j,k))
+     &                +uyz*(ux*hxy(i,j,k)+uy*hyy(i,j,k)+uz*hyz(i,j,k))
+     &                +uzz*(ux*hxz(i,j,k)+uy*hyz(i,j,k)+uz*hzz(i,j,k))
+
+                     spin_int(i,j,k) = 1.0D0/8.0D0/Pi*
+     &                    (uz*term1-ux*term2)*rad**3
+
+                  else if (count.eq.3) then
+              term1 = uxy*(ux*hxx(i,j,k)+uz*hxz(i,j,k))
+     &                +uyy*(ux*hxy(i,j,k)+uy*hyy(i,j,k)+uz*hyz(i,j,k))
+     &                +uyz*(ux*hxz(i,j,k)+uy*hyz(i,j,k)+uz*hzz(i,j,k))
+
+              term2 = uxx*(ux*hxx(i,j,k)+uy*hxy(i,j,k)+uz*hxz(i,j,k))
+     &                +uxy*(uy*hyy(i,j,k)+uz*hyz(i,j,k))
+     &                +uxz*(ux*hxz(i,j,k)+uy*hyz(i,j,k)+uz*hzz(i,j,k))
+
+                     spin_int(i,j,k) = 1.0D0/8.0D0/Pi*
+     &                    (ux*term1-uy*term2)*rad**3
+
+                   end if
+
+
+               ELSE
+                  
+                  spin_int(i,j,k) = 0.0D0
+                  
+               ENDIF
+               
+          ENDDO
+        ENDDO
+      ENDDO
+
+
+c     This is needed when the grid is an octant, but it does not hurt
+c     if it is not
+
+      DO k = 2, size(z)-1
+         DO j = 2, size(y)-1
+            spin_int(1,j,k) = spin_int(2,j,k)
+         ENDDO
+      ENDDO
+      DO k = 2, size(z)-1
+         DO i = 1, size(x)-1
+            spin_int(i,1,k) = spin_int(i,2,k)
+         ENDDO
+      ENDDO
+      DO j = 1, size(y)-1
+         DO i = 1, size(x)-1
+            spin_int(i,j,1) = spin_int(i,j,2)
+         ENDDO
+      ENDDO
+
+
+c     setting counter
+c     ---------------
+
+      if (count.eq.3) then
+         count = 1
+      else
+         count = count + 1
+      end if
+
+
+c     end
+c     ---
+
+      END SUBROUTINE spin_integrand3D
diff --git a/src/spin_integrand3D_int.F b/src/spin_integrand3D_int.F
new file mode 100644
index 0000000..cb2a1ca
--- /dev/null
+++ b/src/spin_integrand3D_int.F
@@ -0,0 +1,30 @@
+
+#include "cctk.h"
+ 
+      MODULE spin_integrand3D_int
+
+c     ------------------------------------------------------------------
+
+      INTERFACE
+
+          SUBROUTINE  spin_integrand3D(origin,Dx,Dy,Dz,x,y,z,gxx,
+     &        gxy,gxz,gyy,gyz,gzz,hxx,hxy,hxz,hyy,hyz,hzz,spin_int,
+     &        Psi,Psi_power)
+
+          IMPLICIT NONE
+          INTEGER,INTENT(IN) ::
+     &         Psi_power
+          CCTK_REAL,INTENT(IN) :: 
+     &        Dx,Dy,Dz,origin(3)
+          CCTK_REAL,DIMENSION(:),INTENT(IN) :: 
+     &        x,y,z
+          CCTK_REAL,DIMENSION(:,:,:),INTENT(IN) :: 
+     &         gxx,gxy,gxz,gyy,gyz,gzz,Psi,
+     &         hxx,hxy,hxz,hyy,hyz,hzz
+          CCTK_REAL,DIMENSION(:,:,:),INTENT(OUT) :: 
+     &        spin_int
+          END SUBROUTINE  spin_integrand3D
+
+      END INTERFACE
+
+      END MODULE spin_integrand3D_int
diff --git a/src/unphysical_to_physical.F b/src/unphysical_to_physical.F
new file mode 100644
index 0000000..74b4ad8
--- /dev/null
+++ b/src/unphysical_to_physical.F
@@ -0,0 +1,56 @@
+
+#include "cctk.h"
+
+c     ==================================================================
+
+      SUBROUTINE unphysical_to_physical(grr,grt,grp,gtt,gtp,gpp,dgtt,
+     &     dgtp,dgpp,Psis,dPsis,Psi_power)
+
+c     ------------------------------------------------------------------
+c
+c     Convert unphysical metric components and radial (eta) derivatives 
+c     on the 2-sphere into physical quantities, using the conformal factor 
+c     on the sphere
+c
+c     ------------------------------------------------------------------
+
+      IMPLICIT NONE
+
+c     Input variables
+      INTEGER :: 
+     &    Psi_power
+      CCTK_REAL,INTENT(INOUT),DIMENSION (:,:) :: 
+     &    grr,grt,grp,gtt,gtp,gpp,dgtt,dgtp,dgpp
+      CCTK_REAL,INTENT(IN),DIMENSION (:,:) :: 
+     &    Psis,dPsis
+
+c     Output variables
+c     WE ARE CHANGING THE INPUT VARIABLES !!!!
+
+c     ------------------------------------------------------------------
+
+      dgtt = Psis**4*dgtt + 4.0d0*Psis**3*dPsis*gtt
+      dgtp = Psis**4*dgtp + 4.0d0*Psis**3*dPsis*gtp
+      dgpp = Psis**4*dgpp + 4.0d0*Psis**3*dPsis*gpp
+
+      grr  = Psis**4*grr
+      grt  = Psis**4*grt
+      grp  = Psis**4*grp
+      gtt  = Psis**4*gtt
+      gtp  = Psis**4*gtp
+      gpp  = Psis**4*gpp
+
+
+      END SUBROUTINE unphysical_to_physical
+
+c     ==================================================================
+
+      
+
+
+
+
+
+
+
+
diff --git a/src/unphysical_to_physical_int.F b/src/unphysical_to_physical_int.F
new file mode 100644
index 0000000..b53fa0b
--- /dev/null
+++ b/src/unphysical_to_physical_int.F
@@ -0,0 +1,21 @@
+
+#include "cctk.h"
+
+      MODULE unphysical_to_physical_int
+
+c     ------------------------------------------------------------------
+
+      INTERFACE
+
+      SUBROUTINE unphysical_to_physical(grr,grt,grp,gtt,gtp,gpp,dgtt,
+     &     dgtp,dgpp,Psis,dPsis,Psi_power)
+      IMPLICIT NONE
+      INTEGER :: 
+     &     Psi_power
+      CCTK_REAL, DIMENSION (:,:) :: 
+     &     grr,grt,grp,gtt,gtp,gpp,dgtt,dgtp,dgpp,Psis,dPsis
+      END SUBROUTINE
+      
+      END INTERFACE
+      
+      END MODULE unphysical_to_physical_int
author	allen <allen@5301f0c2-dbc4-4cee-b2f5-8d7afba4d129>	1999-11-01 14:56:51 +0000
committer	allen <allen@5301f0c2-dbc4-4cee-b2f5-8d7afba4d129>	1999-11-01 14:56:51 +0000
commit	661c56c847f4334feb53fd8d94551b2ab1efdb9a (patch)
tree	3dbdeebf2c0b8f0f0fe30218e9a58fc0cd1638a1 /src
parent	d8cf32e113a96b24744d231013e3d0b81bf55c30 (diff)