Changed assumed shape arrays into explicitly declared arrays to make

the code run on seaborg (it was segmentation faulting before). This
included passing the nx, ny, and nz of gridfunctions into the subroutines.
Also made sure that the interfaces for the subroutine calls agree with the
actual subroutines.


git-svn-id: http://svn.einsteintoolkit.org/cactus/EinsteinAnalysis/Extract/trunk@79 5301f0c2-dbc4-4cee-b2f5-8d7afba4d129

7 files changed, 274 insertions, 76 deletions

diff --git a/src/D2_extract.F b/src/D2_extract.F
index ca832f5..69f65f1 100644
--- a/src/D2_extract.F
+++ b/src/D2_extract.F
@@ -3,8 +3,8 @@
 
 c     ========================================================================
 
-      SUBROUTINE D2_extract(eta,igrid,Nt,Np,theta,phi,all_modes,l,m,g00,g11,
-     &     g12,g13,g22,g23,g33,dg22,dg23,dg33,
+      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,
@@ -32,19 +32,22 @@ c     Input variables
      &     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,
+      INTEGER,INTENT(IN) :: do_ADMmass(2)
+      CCTK_REAL,INTENT(IN) :: eta
+      CCTK_REAL,INTENT(IN),DIMENSION (Nt) :: theta
+      CCTK_REAL,INTENT(IN),DIMENSION (Np) :: phi
+      CCTK_REAL,INTENT(IN),DIMENSION (Nt,Np) :: 
+     &     g00,g12,g13,g23,dg23,
      &     ADMmass_int1,ADMmass_int2,
      &     momentum_int1,momentum_int2,momentum_int3,
      &     spin_int1,spin_int2,spin_int3
 
+      CCTK_REAL,INTENT(INOUT),DIMENSION (Nt,Np) ::
+     &     g11,g22,g33,dg22,dg33
+
 c     Output variables
       CCTK_REAL,INTENT(OUT) :: 
-     &     dtaudt,ADMmass(2),mass,rsch,Qodd(:,2:,0:),Qeven(:,2:,0:),
+     &     ADMmass(2),mass,rsch,Qodd(:,2:,0:),Qeven(:,2:,0:),dtaudt,
      &     momentum(3),spin(3)
 
 c     Local Variables
@@ -80,6 +83,9 @@ c     0. Initial things
 c
 c     ------------------------------------------------------------------------
 
+c      print*,'D2_extract Qodd'
+c      print*,lbound(Qodd,1),lbound(Qodd,2),lbound(Qodd,3)
+c      print*,size(Qodd,1),size(Qodd,2),size(Qodd,3)
       Pi    = two*ASIN(one)
 
       SELECT CASE (igrid)
@@ -112,6 +118,7 @@ c     Calculate ADM mass
 c     ------------------
       IF (do_ADMmass(1) == 1) THEN
 
+c      print*,'1'
 c	Calculate ADM mass using standard formula 
          DO j = 2, Np+1		
             DO i = 2, Nt+1
@@ -123,6 +130,7 @@ c	Calculate ADM mass using standard formula
       ENDIF
       IF (do_ADMmass(2) == 1) THEN
  
+c      print*,'2'
 c       Calculate ADM mass using conformal formula
          DO j = 2, Np+1		
             DO i = 2, Nt+1
@@ -137,6 +145,7 @@ c     Calculate momentum
 c     ------------------
       IF (do_momentum == 1) THEN
 
+c      print*,'3'
          DO j = 2, Np+1		
             DO i = 2, Nt+1
                temp(i,j) = momentum_int1(i-1,j-1)*SIN(theta(i-1))
@@ -144,6 +153,7 @@ c     ------------------
          ENDDO
          momentum(1) = sphere_int(temp,igrid,Nt+1,Np+1,Dt,Dp)
          
+c      print*,'4'
          DO j = 2, Np+1		
             DO i = 2, Nt+1
                temp(i,j) = momentum_int2(i-1,j-1)*SIN(theta(i-1))
@@ -151,6 +161,7 @@ c     ------------------
          ENDDO
          momentum(2) = sphere_int(temp,igrid,Nt+1,Np+1,Dt,Dp)
  
+c      print*,'5'
          DO j = 2, Np+1		
             DO i = 2, Nt+1
                temp(i,j) = momentum_int3(i-1,j-1)*SIN(theta(i-1))
@@ -164,6 +175,7 @@ c     Calculate spin
 c     --------------
       IF (do_spin == 1) THEN
 
+c      print*,'6'
          DO j = 2, Np+1		
             DO i = 2, Nt+1
                temp(i,j) = spin_int1(i-1,j-1)*SIN(theta(i-1))
@@ -171,6 +183,7 @@ c     --------------
          ENDDO
          spin(1) = sphere_int(temp,igrid,Nt+1,Np+1,Dt,Dp)
          
+c      print*,'7'
          DO j = 2, Np+1		
             DO i = 2, Nt+1
                temp(i,j) = spin_int2(i-1,j-1)*SIN(theta(i-1))
@@ -178,6 +191,7 @@ c     --------------
          ENDDO
          spin(2) = sphere_int(temp,igrid,Nt+1,Np+1,Dt,Dp)
  
+c      print*,'8'
          DO j = 2, Np+1		
             DO i = 2, Nt+1
                temp(i,j) = spin_int3(i-1,j-1)*SIN(theta(i-1))
@@ -194,6 +208,7 @@ c
 c     ------------------------------------------------------------------
 
 c     ... g00
+c      print*,'9'
       DO j = 2, Np+1
         DO i = 2, Nt+1
           temp(i,j) = g00(i-1,j-1)*SIN(theta(i-1))
@@ -202,6 +217,7 @@ c     ... g00
       dtaudt = sqrt(one/(four*Pi)*sphere_int(temp,igrid,Nt+1,Np+1,Dt,Dp))
 
 c     ... g11
+c      print*,'10'
       DO j = 2, Np+1
          DO i = 2, Nt+1
             temp(i,j) = g11(i-1,j-1)*SIN(theta(i-1))
@@ -210,6 +226,7 @@ c     ... g11
       sph_g11 = one/(four*Pi)*sphere_int(temp,igrid,Nt+1,Np+1,Dt,Dp)
 
 c     ... g22
+c      print*,'11'
       DO j = 2, Np+1
          DO i = 2, Nt+1
             temp(i,j) = g22(i-1,j-1)*SIN(theta(i-1))
@@ -218,6 +235,7 @@ c     ... g22
       sph_g22 = one/(four*Pi)*sphere_int(temp,igrid,Nt+1,Np+1,Dt,Dp)
 
 c     ... dg22
+c      print*,'12'
       DO j = 2, Np+1
          DO i = 2, Nt+1
             temp(i,j) = dg22(i-1,j-1)*SIN(theta(i-1))
@@ -227,6 +245,7 @@ c     ... dg22
 
 
 c     ... g33/sin(theta)**2
+c      print*,'13'
       DO j = 2, Np+1
          DO i = 2, Nt+1
             temp(i,j) = g33(i-1,j-1)/SIN(theta(i-1))
@@ -253,6 +272,7 @@ c     1. Compute the Schwarzschild radius
 c   
 c     ------------------------------------------------------------------
 
+c      print*,'14'
       DO j = 2, Np+1
          DO i = 2, Nt+1
 c            temp(i,j) = g22(i-1,j-1)*SIN(theta(i-1))
@@ -273,6 +293,7 @@ c        respect to the isotropic radius eta.
 c
 c     ------------------------------------------------------------------
 
+c      print*,'15'
       DO j = 2, Np+1         
          DO i = 2, Nt+1
 c           temp(i,j) = (dg22(i-1,j-1))*SIN(theta(i-1))
@@ -305,6 +326,7 @@ c        Do not know how important this is
 c
 c     ------------------------------------------------------------------
 
+c      print*,'16'
       DO i=1,Nt
          DO j=1,Np
             st  = SIN(theta(i))
@@ -360,12 +382,14 @@ c     ------------------------------------------------------------------
 
 c           Calculate integrands
 
+c      print*,'17'
             loop_phi1: DO j = 1, Np
                loop_theta1: DO i = 1, Nt
 
                   st  = SIN(theta(i))
                   ist = one/st
                   
+c                  print*,'Calling spher_harm_coms',i,j
                   CALL spher_harm_combs(theta(i),phi(j),il,im,Y,Y1,Y2,Y3
      &                 ,Y4)
                   
@@ -390,20 +414,30 @@ c           Calculate integrands
 
                END DO loop_theta1    
             END DO loop_phi1
+c            print*,'End loop'
 
 c           Integrations over the 2-sphere
 
             h1(1) = fac_h1*sphere_int(h1i,igrid,Nt+1,Np+1,Dt,Dp)
+c            print*,'h1 done'
             H2(1) = fac_H2*sphere_int(H2i,igrid,Nt+1,Np+1,Dt,Dp)
+c            print*,'H2 done'
             G(1) = fac_G*sphere_int(Gi,igrid,Nt+1,Np+1,Dt,Dp)
+c            print*,'G done'
             K(1) = fac_K*sphere_int(Ki,igrid,Nt+1,Np+1,Dt,Dp)
      &           +DBLE(il*(il+1))*half*G(1)
+c            print*,'K done'
             c1(1) = fac_c1*sphere_int(c1i,igrid,Nt+1,Np+1,Dt,Dp)
+c            print*,'c1 done'
             c2(1) = fac_c2*sphere_int(c2i,igrid,Nt+1,Np+1,Dt,Dp)
+c            print*,'c2 done'
             dG(1) = fac_dG*sphere_int(dGi,igrid,Nt+1,Np+1,Dt,Dp)
+c            print*,'dG done'
             dK(1) = fac_dK*sphere_int(dKi,igrid,Nt+1,Np+1,Dt,Dp)
      &           +DBLE(il*(il+1))*half*dG(1)
+c            print*,'dK done'
             dc2(1) = fac_dc2*sphere_int(dc2i,igrid,Nt+1,Np+1,Dt,Dp)
+c            print*,'dc2 done'
 
 
 c     ------------------------------------------------------------------
@@ -412,10 +446,12 @@ c     5b. Imaginary parts of the Regge-Wheeler variables
 c
 c     ------------------------------------------------------------------
 
+c      print*,'17 1/2'
             complex_parts: IF (igrid /= 1 .AND. igrid /=2) THEN
 
 c              Calculate integrands
 
+c      print*,'18'
                loop_phi2: DO j = 1, Np
 
                   loop_theta2: DO i = 1, Nt
@@ -486,11 +522,16 @@ c     6. Construct gauge invariant variables
 c
 c     ------------------------------------------------------------------
 
+c            print*,'Qodd',il,im
+c            print*,lbound(Qodd,1),lbound(Qodd,2),lbound(Qodd,3)
+c            print*,size(Qodd,1),size(Qodd,2),size(Qodd,3)
             Qodd(:,il,im) = SQRT(two*DBLE((il+2)*(il+1)*il*(il-1)))*S/
      &           rsch*(c1+half*(dc2-two/rsch*c2))
             
+c            print*,'Lambda'
             Lambda = DBLE((il-1)*(il+2))+3.0D0*(one-S)
             
+c            print*,'Qeven'
             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)
@@ -526,6 +567,7 @@ c     ------------------------------------------------------------------
          WRITE(302,*) eta,Qodd(:,lmax,mmax)
       ENDIF
 
+c      print*,'Exiting D2_extract'
 
       END SUBROUTINE D2_extract
 
diff --git a/src/D2_extract_int.F b/src/D2_extract_int.F
index 43f6906..58a264d 100644
--- a/src/D2_extract_int.F
+++ b/src/D2_extract_int.F
@@ -7,8 +7,8 @@ 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,
+      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,
@@ -19,19 +19,21 @@ c     ------------------------------------------------------------------
      &     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,
+      INTEGER,INTENT(IN) :: do_ADMmass(2)
+      CCTK_REAL,INTENT(IN) :: eta
+      CCTK_REAL,INTENT(IN),DIMENSION (Nt) :: theta
+      CCTK_REAL,INTENT(IN),DIMENSION (Np) :: phi
+      CCTK_REAL,INTENT(IN),DIMENSION (Nt,Np) :: 
+     &     g00,g12,g13,g23,dg23,
      &     ADMmass_int1, ADMmass_int2,
      &     momentum_int1,momentum_int2,momentum_int3,
      &     spin_int1,spin_int2,spin_int3
+
+      CCTK_REAL,INTENT(INOUT),DIMENSION (Nt,Np) ::
+     &     g11,g22,g33,dg22,dg33
+
       CCTK_REAL,INTENT(OUT) :: 
-     &     ADMmass(2),mass,rsch,Qodd(:,:,:),Qeven(:,:,:),dtaudt,
+     &     ADMmass(2),mass,rsch,Qodd(:,2:,0:),Qeven(:,2:,0:),dtaudt,
      &     momentum(3),spin(3)
       END SUBROUTINE
       
diff --git a/src/D3_extract.F b/src/D3_extract.F
index 47a7fb3..7ef0f0f 100644
--- a/src/D3_extract.F
+++ b/src/D3_extract.F
@@ -4,7 +4,8 @@
 c     ==================================================================
 
       SUBROUTINE D3_extract(cctkGH,conformal_state,do_ADMmass,do_momentum,do_spin,igrid,
-     &     origin,myproc,interpolation_operator,interpolation_order,Nt,Np,all_modes,
+     &     origin,myproc,interpolation_operator,interpolation_order,Nt,Np,
+     &     nx,ny,nz,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)
@@ -73,17 +74,18 @@ c     Input variables
       INTEGER,INTENT(IN) :: 
      &     conformal_state,igrid,l,m,Psi_power,myproc
       CCTK_INT,INTENT(IN) :: 
-     &     Nt,Np,all_modes,do_momentum,do_spin,interpolation_order
+     &     Nt,Np,nx,ny,nz,all_modes,do_momentum,do_spin,interpolation_order
       INTEGER,INTENT(IN) ::
      &     do_ADMmass(2)
       CCTK_REAL,INTENT(IN) :: 
      &     origin(3),Dx,Dy,Dz,eta
-      CCTK_REAL,INTENT(IN),DIMENSION(:,:,:) :: 
+      CCTK_REAL,INTENT(IN),DIMENSION(nx,ny,nz) :: 
      &     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(:,:,:) ::
+      CCTK_REAL,INTENT(IN),DIMENSION(nx) :: x
+      CCTK_REAL,INTENT(IN),DIMENSION(ny) :: y
+      CCTK_REAL,INTENT(IN),DIMENSION(nz) :: z
+      CCTK_REAL,INTENT(INOUT),DIMENSION(nx,ny,nz) ::
      &     Extract_temp3d
       CCTK_STRING,INTENT(IN) ::
      &     interpolation_operator
@@ -122,6 +124,11 @@ c     1. Specify polar coordinates on chosen 2-surface
 c
 c     ------------------------------------------------------------------
 
+c      print*,'Entering D3_extract'
+c      print*,nx,ny,nz
+c      print*,'D3_extract Qodd'
+c      print*,lbound(Qodd,1),lbound(Qodd,2),lbound(Qodd,3)
+c      print*,size(Qodd,1),size(Qodd,2),size(Qodd,3)
       Pi = two*ASIN(one)
 
 c     Set polar gridspacing
@@ -163,14 +170,17 @@ c     2. Project quantities onto the 2-surface
 c
 c     ------------------------------------------------------------------
 
+c      print*,nx,ny,nz,Nt,Np
+c      print*,'Calling D3_to_D2'
       CALL D3_to_D2(cctkGH,conformal_state,do_ADMmass,do_momentum,do_spin,
      &     Psi_power,origin,myproc,interpolation_operator,interpolation_order,
      &     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,
+     &     x,y,z,eta,Nt,Np,nx,ny,nz,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      print*,'Exited D3_to_D2'
 
 
       if (myproc .eq. 0) then 
@@ -203,12 +213,17 @@ c     5. Extract gauge invariant quantities on 2-surface
 c
 c     ------------------------------------------------------------------
 
+c      print*,'Calling D2_extract'
          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      print*,Qodd
+c      print*,Qeven
+c      print*,'Exiting D2_extract'
+      
 
 c     End of myproc = 0 conditional
       endif
diff --git a/src/D3_extract_int.F b/src/D3_extract_int.F
index ce14b37..e00b321 100644
--- a/src/D3_extract_int.F
+++ b/src/D3_extract_int.F
@@ -9,7 +9,7 @@ c     ------------------------------------------------------------------
 
       SUBROUTINE D3_extract(cctkGH,conformal_state,do_ADMmass,do_momentum,do_spin,igrid,
      &        origin,myproc,interpolation_operator,interpolation_order,
-     &        Nt,Np,all_modes,l,m,x,y,z,Dx,Dy,Dz,Psi_power,Psi,
+     &        Nt,Np,nx,ny,nz,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)
@@ -19,26 +19,28 @@ c     ------------------------------------------------------------------
       CCTK_POINTER :: cctkGH
 
       INTEGER,INTENT(IN) :: 
-     &     igrid,l,m,Psi_power,myproc
+     &     conformal_state,igrid,l,m,Psi_power,myproc
       CCTK_INT,INTENT(IN) ::
-     &     conformal_state,Nt,Np,all_modes,do_momentum,do_spin,interpolation_order
+     &     Nt,Np,nx,ny,nz,all_modes,do_momentum,do_spin,interpolation_order
       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(:,:,:) :: 
+      CCTK_REAL,INTENT(IN),DIMENSION(nx,ny,nz) :: 
      &     Psi,g00,gxx,gxy,gxz,gyy,gyz,gzz,
      &     hxx,hxy,hxz,hyy,hyz,hzz
-      CCTK_REAL,INTENT(INOUT),DIMENSION(:,:,:) ::
+      CCTK_REAL,INTENT(IN),DIMENSION(nx) :: x
+      CCTK_REAL,INTENT(IN),DIMENSION(ny) :: y
+      CCTK_REAL,INTENT(IN),DIMENSION(nz) :: z
+      CCTK_REAL,INTENT(INOUT),DIMENSION(nx,ny,nz) ::
      &     Extract_Temp3d
-      CCTK_REAL,INTENT(OUT) :: 
-     &     ADMmass(2),mass,rsch,Qodd(:,:,:),Qeven(:,:,:),dtaudt,
-     &     momentum(3),spin(3)
       CCTK_STRING,INTENT(IN) ::
      &     interpolation_operator
 
+      CCTK_REAL,INTENT(OUT) :: 
+     &     ADMmass(2),mass,rsch,Qodd(:,2:,0:),Qeven(:,2:,0:),dtaudt,
+     &     momentum(3),spin(3)
+
 
       END SUBROUTINE
          
diff --git a/src/D3_to_D2.F b/src/D3_to_D2.F
index da2294e..1e71ff4 100644
--- a/src/D3_to_D2.F
+++ b/src/D3_to_D2.F
@@ -7,7 +7,7 @@ c     ========================================================================
      &     Psi_power,origin,myproc,interpolation_operator,interpolation_order,
      &     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,
+     &     x,y,z,eta,Nt,Np,nx,ny,nz,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,
@@ -35,24 +35,26 @@ c     Input variables
       INTEGER,INTENT(IN) ::
      &     conformal_state,myproc,Psi_power
       CCTK_INT, INTENT(IN) ::
-     &     Nt,Np,do_momentum,do_spin,interpolation_order
+     &     Nt,Np,nx,ny,nz,do_momentum,do_spin,interpolation_order
       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(:,:,:) ::
+      CCTK_REAL,INTENT(IN),DIMENSION(Nt) :: theta
+      CCTK_REAL,INTENT(IN),DIMENSION(Np) :: phi
+      CCTK_REAL,INTENT(IN),DIMENSION(nx) :: x
+      CCTK_REAL,INTENT(IN),DIMENSION(ny) :: y
+      CCTK_REAL,INTENT(IN),DIMENSION(nz) :: z
+      CCTK_REAL,INTENT(IN),DIMENSION(nx,ny,nz) ::
      &     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)
+      INTEGER,INTENT(IN) :: do_ADMmass(2)
       CCTK_STRING,INTENT(IN) :: 
      &     interpolation_operator
+      CCTK_REAL,INTENT(INOUT),DIMENSION(nx,ny,nz) ::
+     &     Extract_temp3d
 
 c     Output variables
 
-      CCTK_REAL,INTENT(OUT),DIMENSION(:,:) ::
+      CCTK_REAL,INTENT(OUT),DIMENSION(Nt,Np) ::
      &     Psis,g00s,gxxs,gxys,gxzs,gyys,
      &     gyzs,gzzs,dPsis,dgxxs,dgxys,dgxzs,dgyys,dgyzs,dgzzs,
      &     ADMmass_int1,ADMmass_int2,
@@ -65,7 +67,7 @@ c     Local variables, passed on
       LOGICAL ::
      &    err_flag
       INTEGER ::
-     &    iorder,npoints,Nx,Ny,Nz
+     &    iorder,npoints
       INTEGER,DIMENSION(Nt,Np) ::
      &    ib,jb,kb
       CCTK_REAL,DIMENSION(Nt,Np) ::
@@ -83,18 +85,14 @@ c     Local variables, here only
       character(128) options_string
       character(128) operator
       CCTK_INT nchars
+      character(len=20) :: error_code
 
 c     ------------------------------------------------------------------------
 
-c     Initial Stuff
-c     -------------
-      Nx = SIZE(x)
-      Ny = SIZE(y)
-      Nz = SIZE(z)
-
-
 c     Compute Cartesian coordinates on the surface
 c     --------------------------------------------
+c      print*,'Entering D3_to_D2'
+c      print*,nx,ny,nz,Np,Nt
       DO j = 1, Np
          DO i = 1, Nt
 
@@ -105,6 +103,7 @@ c     --------------------------------------------
          ENDDO
       ENDDO
 
+c      print*,'Coordinates set'
 
 c     Only do interpolation on one processor
 c     --------------------------------------
@@ -118,6 +117,7 @@ c     --------------------------------------
 
       END SELECT
 
+c      print*,'Number of points set'
 
 c     Get the interpolator, parameter table, and coordinate system handles
 c     --------------------------------------------------------------------
@@ -143,6 +143,7 @@ c     --------------------------------------------------------------------
         call CCTK_WARN(0,"Cannot get handle for cart3d coordinate system ! Forgot to activate an implementation providing coordinates ??")
       endif
 
+c      print*,'Handles obtained'
 
 c     fill in the input/output arrays for the interpolator
 c     ----------------------------------------------------
@@ -179,13 +180,19 @@ c     ------------------------------------------------------------
         num_arrays = 7
       end if
 
+c      print*,'1'
+c      print*,in_array_indices
+c      print*,out_arrays
       call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
      $                            param_table_handle, coord_system_handle,
      $                            npoints, CCTK_VARIABLE_REAL, interp_coords,
      $                            num_arrays, in_array_indices,
      $                            num_arrays, out_array_type_codes, out_arrays)
+c      print*,'1 done'
       if (ierror < 0) then
         call CCTK_WARN (1, "interpolator call returned an error code");
+        write(error_code,'(a7,i13)') 'code = ',ierror
+        call CCTK_WARN (1, error_code );
       endif
 
 
@@ -194,57 +201,87 @@ c     ----------------------------------------------------
       call CCTK_VarIndex(in_array_indices(1), "extract::temp3d")
 
       if (conformal_state > 0) then
+c         print*,'1 1/2'
          CALL met_rad_der(origin,Dx,Dy,Dz,x,y,z,Psi,Extract_temp3d)
+c         print*,'1 1/2 done'
          CALL CCTK_SyncGroup(ierror,cctkGH,"extract::temps")
 
          out_arrays(1) = CCTK_PointerTo(dPsis)
+c         print*,'2'
+c         print*,in_array_indices(1)
+c         print*,out_arrays(1)
          call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
      $                               param_table_handle, coord_system_handle,
      $                               npoints, CCTK_VARIABLE_REAL, interp_coords,
      $                               1, in_array_indices(1),
      $                               1, out_array_type_codes, out_arrays(1))
+c         print*,'2 done'
          if (ierror < 0) then
            call CCTK_WARN (1, "interpolator call returned an error code");
+           write(error_code,'(a7,i13)') 'code = ',ierror
+           call CCTK_WARN (1, error_code );
          endif
       end if
 
+c      print*,'3'
       CALL met_rad_der(origin,Dx,Dy,Dz,x,y,z,gxx,Extract_temp3d)
+c      print*,'3 done'
+c      print*,'4'
       CALL CCTK_SyncGroup(ierror,cctkGH,"extract::temps")
+c      print*,'4 done'
 
       out_arrays(1) = CCTK_PointerTo(dgxxs)
+c      print*,'5'
+c      print*,in_array_indices(1)
+c      print*,out_arrays(1)
       call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
      $                            param_table_handle, coord_system_handle,
      $                            npoints, CCTK_VARIABLE_REAL, interp_coords,
      $                            1, in_array_indices(1),
      $                            1, out_array_type_codes, out_arrays(1))
+c      print*,'5 done'
       if (ierror < 0) then
         call CCTK_WARN (1, "interpolator call returned an error code");
+        write(error_code,'(a7,i13)') 'code = ',ierror
+        call CCTK_WARN (1, error_code );
       endif
 
       CALL met_rad_der(origin,Dx,Dy,Dz,x,y,z,gxy,Extract_temp3d)
       CALL CCTK_SyncGroup(ierror,cctkGH,"extract::temps")
 
       out_arrays(1) = CCTK_PointerTo(dgxys)
+c      print*,'6'
+c      print*,in_array_indices(1)
+c      print*,out_arrays(1)
       call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
      $                            param_table_handle, coord_system_handle,
      $                            npoints, CCTK_VARIABLE_REAL, interp_coords,
      $                            1, in_array_indices(1),
      $                            1, out_array_type_codes, out_arrays(1))
+c      print*,'6 done'
       if (ierror < 0) then
         call CCTK_WARN (1, "interpolator call returned an error code");
+        write(error_code,'(a7,i13)') 'code = ',ierror
+        call CCTK_WARN (1, error_code );
       endif
 
       CALL met_rad_der(origin,Dx,Dy,Dz,x,y,z,gxz,Extract_temp3d)
       CALL CCTK_SyncGroup(ierror,cctkGH,"extract::temps")
 
       out_arrays(1) = CCTK_PointerTo(dgxzs)
+c      print*,'7'
+c      print*,in_array_indices(1)
+c      print*,out_arrays(1)
       call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
      $                            param_table_handle, coord_system_handle,
      $                            npoints, CCTK_VARIABLE_REAL, interp_coords,
      $                            1, in_array_indices(1),
      $                            1, out_array_type_codes, out_arrays(1))
+c      print*,'7 done'
       if (ierror < 0) then
         call CCTK_WARN (1, "interpolator call returned an error code");
+        write(error_code,'(a7,i13)') 'code = ',ierror
+        call CCTK_WARN (1, error_code );
       endif
 
       CALL met_rad_der(origin,Dx,Dy,Dz,x,y,z,gyy,Extract_temp3d)
@@ -258,56 +295,80 @@ c     ----------------------------------------------------
      $                            1, out_array_type_codes, out_arrays(1))
       if (ierror < 0) then
         call CCTK_WARN (1, "interpolator call returned an error code");
+        write(error_code,'(a7,i13)') 'code = ',ierror
+        call CCTK_WARN (1, error_code );
       endif
 
       CALL met_rad_der(origin,Dx,Dy,Dz,x,y,z,gyz,Extract_temp3d)
       CALL CCTK_SyncGroup(ierror,cctkGH,"extract::temps")
 
       out_arrays(1) = CCTK_PointerTo(dgyzs)
+c      print*,'8'
+c      print*,in_array_indices(1)
+c      print*,out_arrays(1)
       call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
      $                            param_table_handle, coord_system_handle,
      $                            npoints, CCTK_VARIABLE_REAL, interp_coords,
      $                            1, in_array_indices(1),
      $                            1, out_array_type_codes, out_arrays(1))
+c      print*,'8 done'
       if (ierror < 0) then
         call CCTK_WARN (1, "interpolator call returned an error code");
+        write(error_code,'(a7,i13)') 'code = ',ierror
+        call CCTK_WARN (1, error_code );
       endif
 
       CALL met_rad_der(origin,Dx,Dy,Dz,x,y,z,gzz,Extract_temp3d)
       CALL CCTK_SyncGroup(ierror,cctkGH,"extract::temps")
 
       out_arrays(1) = CCTK_PointerTo(dgzzs)
+c      print*,'9'
+c      print*,in_array_indices(1)
+c      print*,out_arrays(1)
       call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
      $                            param_table_handle, coord_system_handle,
      $                            npoints, CCTK_VARIABLE_REAL, interp_coords,
      $                            1, in_array_indices(1),
      $                            1, out_array_type_codes, out_arrays(1))
+c      print*,'9 done'
       if (ierror < 0) then
         call CCTK_WARN (1, "interpolator call returned an error code");
+        write(error_code,'(a7,i13)') 'code = ',ierror
+        call CCTK_WARN (1, error_code );
       endif
 
 
 c     Calculate integrands for ADM masses
 c     -----------------------------------
 c     Standard equation
-      IF (do_ADMmass(1)) THEN
+      IF (do_ADMmass(1) == 1) THEN
+c         print*,'10'
+c         print*,in_array_indices(1)
+c         print*,out_arrays(1)
          CALL ADMmass_integrand3D(origin,Dx,Dy,Dz,x,y,z,gxx,gxy,
      &        gxz,gyy,gyz,gzz,Extract_temp3d,Psi,Psi_power,conformal_state)
          CALL CCTK_SyncGroup(ierror,cctkGH,"extract::temps")
+c         print*,'10 done'
 
          out_arrays(1) = CCTK_PointerTo(ADMmass_int1)
+c         print*,'11'
+c         print*,in_array_indices(1)
+c         print*,out_arrays(1)
          call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
      $                               param_table_handle, coord_system_handle,
      $                               npoints, CCTK_VARIABLE_REAL, interp_coords,
      $                               1, in_array_indices(1),
      $                               1, out_array_type_codes, out_arrays(1))
+c         print*,'11 done'
          if (ierror < 0) then
            call CCTK_WARN (1, "interpolator call returned an error code");
+           write(error_code,'(a7,i13)') 'code = ',ierror
+           call CCTK_WARN (1, error_code );
          endif
       END IF
 
 c     Conformal equation
-      IF (do_ADMmass(2)) THEN
+      IF (do_ADMmass(2) == 1) THEN
              ADMmass_int2 = -eta**2/2D0/3.1416D0*dPsis
       ENDIF
 
@@ -315,49 +376,79 @@ c     Calculate integrands for momentum
 c     ---------------------------------
       IF (do_momentum==1) THEN
 
+c         print*,'12'
+c         print*,in_array_indices(1)
+c         print*,out_arrays(1)
          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,conformal_state)
          CALL CCTK_SyncGroup(ierror,cctkGH,"extract::temps")
+c         print*,'12 done'
 
          out_arrays(1) = CCTK_PointerTo(momentum_int1)
+c         print*,'13'
+c         print*,in_array_indices(1)
+c         print*,out_arrays(1)
          call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
      $                               param_table_handle, coord_system_handle,
      $                               npoints, CCTK_VARIABLE_REAL, interp_coords,
      $                               1, in_array_indices(1),
      $                               1, out_array_type_codes, out_arrays(1))
+c         print*,'13 done'
          if (ierror < 0) then
            call CCTK_WARN (1, "interpolator call returned an error code");
+           write(error_code,'(a7,i13)') 'code = ',ierror
+           call CCTK_WARN (1, error_code );
          endif
 
+c         print*,'14'
+c         print*,in_array_indices(1)
+c         print*,out_arrays(1)
          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,conformal_state)
          CALL CCTK_SyncGroup(ierror,cctkGH,"extract::temps")
+c         print*,'14 done'
 
          out_arrays(1) = CCTK_PointerTo(momentum_int2)
+c         print*,'15'
+c         print*,in_array_indices(1)
+c         print*,out_arrays(1)
          call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
      $                               param_table_handle, coord_system_handle,
      $                               npoints, CCTK_VARIABLE_REAL, interp_coords,
      $                               1, in_array_indices(1),
      $                               1, out_array_type_codes, out_arrays(1))
+c         print*,'15 done'
          if (ierror < 0) then
            call CCTK_WARN (1, "interpolator call returned an error code");
+           write(error_code,'(a7,i13)') 'code = ',ierror
+           call CCTK_WARN (1, error_code );
          endif
 
+c         print*,'16'
+c         print*,in_array_indices(1)
+c         print*,out_arrays(1)
          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,conformal_state)
          CALL CCTK_SyncGroup(ierror,cctkGH,"extract::temps")
+c         print*,'16 done'
 
          out_arrays(1) = CCTK_PointerTo(momentum_int3)
+c         print*,'17'
+c         print*,in_array_indices(1)
+c         print*,out_arrays(1)
          call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
      $                               param_table_handle, coord_system_handle,
      $                               npoints, CCTK_VARIABLE_REAL, interp_coords,
      $                               1, in_array_indices(1),
      $                               1, out_array_type_codes, out_arrays(1))
+c         print*,'17 done'
          if (ierror < 0) then
            call CCTK_WARN (1, "interpolator call returned an error code");
+           write(error_code,'(a7,i13)') 'code = ',ierror
+           call CCTK_WARN (1, error_code );
          endif
       END IF
 
@@ -365,49 +456,79 @@ c     Calculate integrands for spin
 c     -----------------------------
       IF (do_spin==1) THEN
 
+c         print*,'18'
+c         print*,in_array_indices(1)
+c         print*,out_arrays(1)
          CALL spin_integrand3D(origin,x,y,z,
      &        gxx,gxy,gxz,gyy,gyz,gzz,hxx,hxy,hxz,hyy,hyz,hzz,
      &        Extract_temp3d,Psi,Psi_power,conformal_state)
          CALL CCTK_SyncGroup(ierror,cctkGH,"extract::temps")
+c         print*,'18 done'
 
          out_arrays(1) = CCTK_PointerTo(spin_int1)
+c         print*,'19'
+c         print*,in_array_indices(1)
+c         print*,out_arrays(1)
          call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
      $                               param_table_handle, coord_system_handle,
      $                               npoints, CCTK_VARIABLE_REAL, interp_coords,
      $                               1, in_array_indices(1),
      $                               1, out_array_type_codes, out_arrays(1))
+c         print*,'19 done'
          if (ierror < 0) then
            call CCTK_WARN (1, "interpolator call returned an error code");
+           write(error_code,'(a7,i13)') 'code = ',ierror
+           call CCTK_WARN (1, error_code );
          endif
 
+c         print*,'20'
+c         print*,in_array_indices(1)
+c         print*,out_arrays(1)
          CALL spin_integrand3D(origin,x,y,z,
      &        gxx,gxy,gxz,gyy,gyz,gzz,hxx,hxy,hxz,hyy,hyz,hzz,
      &        Extract_temp3d,Psi,Psi_power,conformal_state)
          CALL CCTK_SyncGroup(ierror,cctkGH,"extract::temps")
+c         print*,'20 done'
 
          out_arrays(1) = CCTK_PointerTo(spin_int2)
+c         print*,'21'
+c         print*,in_array_indices(1)
+c         print*,out_arrays(1)
          call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
      $                               param_table_handle, coord_system_handle,
      $                               npoints, CCTK_VARIABLE_REAL, interp_coords,
      $                               1, in_array_indices(1),
      $                               1, out_array_type_codes, out_arrays(1))
+c         print*,'21 done'
          if (ierror < 0) then
            call CCTK_WARN (1, "interpolator call returned an error code");
+           write(error_code,'(a7,i13)') 'code = ',ierror
+           call CCTK_WARN (1, error_code );
          endif
 
+c         print*,'22'
+c         print*,in_array_indices(1)
+c         print*,out_arrays(1)
          CALL spin_integrand3D(origin,x,y,z,
      &        gxx,gxy,gxz,gyy,gyz,gzz,hxx,hxy,hxz,hyy,hyz,hzz,
      &        Extract_temp3d,Psi,Psi_power,conformal_state)
          CALL CCTK_SyncGroup(ierror,cctkGH,"extract::temps")
+c         print*,'22 done'
 
          out_arrays(1) = CCTK_PointerTo(spin_int3)
+c         print*,'23'
+c         print*,in_array_indices(1)
+c         print*,out_arrays(1)
          call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
      $                               param_table_handle, coord_system_handle,
      $                               npoints, CCTK_VARIABLE_REAL, interp_coords,
      $                               1, in_array_indices(1),
      $                               1, out_array_type_codes, out_arrays(1))
+c         print*,'23 done'
          if (ierror < 0) then
            call CCTK_WARN (1, "interpolator call returned an error code");
+           write(error_code,'(a7,i13)') 'code = ',ierror
+           call CCTK_WARN (1, error_code );
          endif
 
       END IF
diff --git a/src/D3_to_D2_int.F b/src/D3_to_D2_int.F
index db013a2..4aeef0a 100644
--- a/src/D3_to_D2_int.F
+++ b/src/D3_to_D2_int.F
@@ -11,7 +11,7 @@ c     ------------------------------------------------------------------
      &        Psi_power,origin,myproc,interpolation_operator,interpolation_order,
      &        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,
+     &        x,y,z,eta,Nt,Np,nx,ny,nz,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,
@@ -22,26 +22,30 @@ c     ------------------------------------------------------------------
       INTEGER,INTENT(IN) :: 
      &     conformal_state,myproc,Psi_power
       CCTK_INT, INTENT(IN) ::
-     &     Nt,Np,do_momentum,do_spin,interpolation_order
-      INTEGER ::
-     &     do_ADMmass(2)
+     &     Nt,Np,nx,ny,nz,do_momentum,do_spin,interpolation_order
       CCTK_REAL,INTENT(IN) :: 
      &     origin(3),Dx,Dy,Dz,eta
-      CCTK_REAL,INTENT(IN),DIMENSION(:,:,:) :: 
+      CCTK_REAL,INTENT(IN),DIMENSION(Nt) :: theta
+      CCTK_REAL,INTENT(IN),DIMENSION(Np) :: phi
+      CCTK_REAL,INTENT(IN),DIMENSION(nx) :: x
+      CCTK_REAL,INTENT(IN),DIMENSION(ny) :: y
+      CCTK_REAL,INTENT(IN),DIMENSION(nz) :: z
+      CCTK_REAL,INTENT(IN),DIMENSION(nx,ny,nz) ::
      &     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
+      INTEGER,INTENT(IN) :: do_ADMmass(2)
       CCTK_STRING,INTENT(IN) ::
      &     interpolation_operator
 
+      CCTK_REAL,INTENT(INOUT),DIMENSION(nx,ny,nz) ::
+     &     Extract_temp3d
+
+      CCTK_REAL,INTENT(OUT),DIMENSION(Nt,Np) :: 
+     &     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
       END SUBROUTINE
 
       END INTERFACE
diff --git a/src/Extract.F b/src/Extract.F
index ace99fb..b4ba753 100644
--- a/src/Extract.F
+++ b/src/Extract.F
@@ -83,6 +83,7 @@ c     Local variables
 
       INTEGER myproc,ierr
       CCTK_REAL dx,dy,dz
+      CCTK_INT :: nx, ny, nz
 
       character*80 infoline
 
@@ -93,6 +94,10 @@ c     ------------------------------------------------------------------
       dy = cctk_delta_space(2)
       dz = cctk_delta_space(3)
 
+      nx = cctk_lsh(1)
+      ny = cctk_lsh(2)
+      nz = cctk_lsh(3)
+
 c     ------------------------------------------------------------------
 c
 c     0. Check to see if Extract should have been called
@@ -149,13 +154,13 @@ c     Set the value of igrid
 
 c     Create 1D coordinate arrays
 
-      do ix = 1, cctk_lsh(1)
+      do ix = 1, nx
         x_1d(ix) = x(ix,1,1)
       enddo
-      do iy = 1, cctk_lsh(2)
+      do iy = 1, ny
         y_1d(iy) = y(1,iy,1)
       enddo
-      do iz = 1, cctk_lsh(3)
+      do iz = 1, nz
         z_1d(iz) = z(1,1,iz)
       enddo
 
@@ -218,6 +223,9 @@ c     Get m_mode in this case too, or the T3E has a junk value (PW)
          mmode = m_mode
          ALLOCATE(Qodd(2,2:lmode,0:mmode),Qeven(2,2:lmode,0:mmode))
       ENDIF
+c      print*,'Extract Qodd'
+c      print*,lbound(Qodd,1),lbound(Qodd,2),lbound(Qodd,3)
+c      print*,size(Qodd,1),size(Qodd,2),size(Qodd,3)
 
 
 c     Check do not have bad values for the modes
@@ -377,13 +385,17 @@ c           Do extraction at each radius
 
             extract_at_each_radius: DO WHILE (radius < r2)
 
+c               print*,'Calling D3_extract'
+c               print*,nx,ny,nz,Nt,Np
                CALL D3_extract(cctkGH,conformal_state,do_ADMmass,do_momentum,do_spin,
      &           igrid,orig,myproc,interpolation_operator,interpolation_order,
-     &           Nt,Np,all_modes,lmode,
+     &           Nt,Np,nx,ny,nz,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)
+c               print*,Qodd,Qeven
+c               print*,'Exited D3_extract'
 
                IF (verbose == 1) THEN
                   WRITE(*,*) "Extracted at r =",radius
@@ -551,7 +563,7 @@ c     Cannot use the conformal equation for ADM mass now
 
             CALL D3_extract(cctkGH,conformal_state,do_ADMmass,do_momentum,do_spin,
      &           igrid,orig,myproc,interpolation_operator,interpolation_order,
-     &           Nt,Np,all_modes,
+     &           Nt,Np,nx,ny,nz,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,