Correct nans

git-svn-id: https://svn.cct.lsu.edu/repos/numrel/LSUThorns/QuasiLocalMeasures/trunk@5 ef6f4158-a8ee-47d1-ba14-cb73256398e0

12 files changed, 618 insertions, 445 deletions

diff --git a/interface.ccl b/interface.ccl
index fffdbf8..0c01900 100644
--- a/interface.ccl
+++ b/interface.ccl
@@ -2,7 +2,7 @@
 
 IMPLEMENTS: QuasiLocalMeasures
 
-INHERITS: ADMBase SphericalSurface
+INHERITS: ADMBase SphericalSurface TmunuBase
 
 
 
@@ -86,24 +86,6 @@ COMPLEX qlm_tetrad_m[num_surfaces] TYPE=array DIM=2 SIZE=SphericalSurface::maxnt
   qlm_m0 qlm_m1 qlm_m2 qlm_m3       # null normal within surface
 } "Tetrad vector m^mu"
 
-REAL qlm_tetrad_l_p[num_surfaces] TYPE=array DIM=2 SIZE=SphericalSurface::maxntheta,SphericalSurface::maxnphi DISTRIB=constant TAGS='convergence_power=1'
-{
-  qlm_l0_p qlm_l1_p qlm_l2_p qlm_l3_p
-  qlm_l0_p_p qlm_l1_p_p qlm_l2_p_p qlm_l3_p_p
-} "Previous tetrad vectors l^mu"
-
-REAL qlm_tetrad_n_p[num_surfaces] TYPE=array DIM=2 SIZE=SphericalSurface::maxntheta,SphericalSurface::maxnphi DISTRIB=constant TAGS='convergence_power=1'
-{
-  qlm_n0_p qlm_n1_p qlm_n2_p qlm_n3_p
-  qlm_n0_p_p qlm_n1_p_p qlm_n2_p_p qlm_n3_p_p
-} "Previous tetrad vectors n^mu"
-
-COMPLEX qlm_tetrad_m_p[num_surfaces] TYPE=array DIM=2 SIZE=SphericalSurface::maxntheta,SphericalSurface::maxnphi DISTRIB=constant TAGS='convergence_power=1'
-{
-  qlm_m0_p qlm_m1_p qlm_m2_p qlm_m3_p
-  qlm_m0_p_p qlm_m1_p_p qlm_m2_p_p qlm_m3_p_p
-} "Previous tetrad vectors m^mu"
-
 COMPLEX qlm_newman_penrose[num_surfaces] TYPE=array DIM=2 SIZE=SphericalSurface::maxntheta,SphericalSurface::maxnphi DISTRIB=constant TAGS='Checkpoint="no"' TAGS='convergence_power=1'
 {
   qlm_npkappa qlm_nptau qlm_npsigma qlm_nprho
@@ -156,12 +138,6 @@ REAL qlm_killing_vector[num_surfaces] TYPE=array DIM=2 SIZE=SphericalSurface::ma
   qlm_chi
 } "Killing vector field"
 
-REAL qlm_killing_vector_p[num_surfaces] TYPE=array DIM=2 SIZE=SphericalSurface::maxntheta,SphericalSurface::maxnphi DISTRIB=constant TAGS='convergence_power=1'
-{
-  qlm_xi_t_p qlm_xi_p_p
-  qlm_xi_t_p_p qlm_xi_p_p_p
-} "Previous Killing vector field"
-
 REAL qlm_killed_twometric[num_surfaces] TYPE=array DIM=2 SIZE=SphericalSurface::maxntheta,SphericalSurface::maxnphi DISTRIB=constant TAGS='Checkpoint="no"' TAGS='convergence_power=1'
 {
   qlm_lqtt qlm_lqtp qlm_lqpp
diff --git a/schedule.ccl b/schedule.ccl
index 92ca28f..630aa04 100644
--- a/schedule.ccl
+++ b/schedule.ccl
@@ -4,8 +4,8 @@ STORAGE: qlm_state qlm_scalars
 STORAGE: qlm_state_p qlm_scalars_p
 STORAGE: qlm_grid_int qlm_grid_real qlm_grid_real_p
 STORAGE: qlm_shapes qlm_tetrad_l qlm_tetrad_n qlm_tetrad_m
-STORAGE: qlm_shapes_p qlm_tetrad_l_p qlm_tetrad_n_p qlm_tetrad_m_p
-STORAGE: qlm_killing_vector qlm_killing_vector_p
+STORAGE: qlm_shapes_p
+STORAGE: qlm_killing_vector
 
 # Ensure that there is enough storage for the ADM variables
 STORAGE: ADMBase::metric[3] ADMBase::curv[3] ADMBase::lapse[3] ADMBase::shift[3]
@@ -23,15 +23,15 @@ SCHEDULE qlm_paramcheck AT paramcheck
 {
   OPTIONS: global
   LANG: Fortran
-} "Check Isolated and Dynamical Horizon parameter settings"
+} "Check quasi-local parameter settings"
 
 SCHEDULE qlm_init AT initial
 {
   OPTIONS: global
   LANG: Fortran
-} "Initialise Isolated and Dynamical Horizon calculations"
+} "Initialise quasi-local calculations"
 
-SCHEDULE qlm_calculate AT analysis AFTER SphericalSurface_HasBeenSet
+SCHEDULE qlm_calculate AT analysis AFTER (SphericalSurface_HasBeenSet SetTmunu)
 {
   OPTIONS: global
   LANG: Fortran
@@ -47,4 +47,4 @@ SCHEDULE qlm_calculate AT analysis AFTER SphericalSurface_HasBeenSet
   TRIGGERS: qlm_ricci_scalars qlm_twometric qlm_killing_vector qlm_killed_twometric
   TRIGGERS: qlm_invariant_coordinates qlm_multipole_moments qlm_3determinant
   TRIGGERS: qlm_scalars
-} "Calculate Isolated and Dynamical Horizon quantities"
+} "Calculate quasi-local quantities"
diff --git a/src/qlm_3determinant.F90 b/src/qlm_3determinant.F90
index d263009..1b575c1 100644
--- a/src/qlm_3determinant.F90
+++ b/src/qlm_3determinant.F90
@@ -37,7 +37,7 @@ SUBROUTINE qlm_calc_3determinant (CCTK_ARGUMENTS, hn)
   CCTK_REAL    :: theta, phi
   
   IF (veryverbose/=0) THEN
-     CALL CCTK_INFO ("Computing 3-Volume Element of Horizon")
+     CALL CCTK_INFO ("Computing 3-Volume Element of Surface")
   END IF
   
   t0 = qlm_time(hn)
diff --git a/src/qlm_broadcast.c b/src/qlm_broadcast.c
index a491fb5..3266177 100644
--- a/src/qlm_broadcast.c
+++ b/src/qlm_broadcast.c
@@ -147,9 +147,6 @@ CCTK_FNAME(qlm_broadcast) (CCTK_POINTER_TO_CONST * restrict const cctkGH_)
     bcast (cctkGH, "QuasiLocalMeasures::qlm_tetrad_l"             , hn, root);
     bcast (cctkGH, "QuasiLocalMeasures::qlm_tetrad_n"             , hn, root);
     bcast (cctkGH, "QuasiLocalMeasures::qlm_tetrad_m"             , hn, root);
-    bcast (cctkGH, "QuasiLocalMeasures::qlm_tetrad_l_p"           , hn, root);
-    bcast (cctkGH, "QuasiLocalMeasures::qlm_tetrad_n_p"           , hn, root);
-    bcast (cctkGH, "QuasiLocalMeasures::qlm_tetrad_m_p"           , hn, root);
     bcast (cctkGH, "QuasiLocalMeasures::qlm_newman_penrose"       , hn, root);
     bcast (cctkGH, "QuasiLocalMeasures::qlm_weyl_scalars"         , hn, root);
     bcast (cctkGH, "QuasiLocalMeasures::qlm_ricci_scalars"        , hn, root);
diff --git a/src/qlm_calculate.F90 b/src/qlm_calculate.F90
index eb4c4ba..19b0ffb 100644
--- a/src/qlm_calculate.F90
+++ b/src/qlm_calculate.F90
@@ -106,6 +106,9 @@ subroutine qlm_calculate (CCTK_ARGUMENTS)
         end if
         call qlm_calc_3determinant (CCTK_PASS_FTOF, hn)
         call qlm_analyse (CCTK_PASS_FTOF, hn)
+        if (qlm_have_killing_vector(hn) /= 0) then
+           call qlm_multipoles_normalise (CCTK_PASS_FTOF, hn)
+        end if
         
 9999    continue
         
diff --git a/src/qlm_coordinates.F90 b/src/qlm_coordinates.F90
index 6807bcd..17b7b0c 100644
--- a/src/qlm_coordinates.F90
+++ b/src/qlm_coordinates.F90
@@ -9,6 +9,7 @@ subroutine qlm_calc_coordinates (CCTK_ARGUMENTS, hn)
   use cctk
   use constants
   use qlm_boundary
+  use tensor2
   implicit none
   DECLARE_CCTK_ARGUMENTS
   DECLARE_CCTK_FUNCTIONS
@@ -43,7 +44,8 @@ subroutine qlm_calc_coordinates (CCTK_ARGUMENTS, hn)
         qq(1,2) = qlm_qtp(i,j,hn)
         qq(2,2) = qlm_qpp(i,j,hn)
         qq(2,1) = qq(1,2)
-        dtq = qq(1,1) * qq(2,2) - qq(1,2) * qq(2,1)
+        
+        call calc_2det (qq, dtq)
         
         area = area + sqrt(dtq) * qlm_delta_theta(hn) * qlm_delta_phi(hn)
         
@@ -98,7 +100,8 @@ subroutine qlm_calc_coordinates (CCTK_ARGUMENTS, hn)
         qq(1,2) = qlm_qtp(i,j,hn)
         qq(2,2) = qlm_qpp(i,j,hn)
         qq(2,1) = qq(1,2)
-        dtq = qq(1,1) * qq(2,2) - qq(1,2) * qq(2,1)
+        
+        call calc_2det (qq, dtq)
         
         integral_z = integral_z + qlm_inv_z(i,j,hn) * sqrt(dtq) * qlm_delta_theta(hn) * qlm_delta_phi(hn)
         
@@ -172,7 +175,8 @@ contains
     qq(1,2) = qlm_qtp(i,j,hn)
     qq(2,2) = qlm_qpp(i,j,hn)
     qq(2,1) = qq(1,2)
-    dtq = qq(1,1) * qq(2,2) - qq(1,2) * qq(2,1)
+    
+    call calc_2det (qq, dtq)
     
     zdot = vv(1) * (- (1/radius**2) * qlm_xi_p(i,j,hn) * sqrt(dtq)) &
          + vv(2) * (+ (1/radius**2) * qlm_xi_t(i,j,hn) * sqrt(dtq))
diff --git a/src/qlm_import_surface.F90 b/src/qlm_import_surface.F90
index 8eec453..d10b736 100644
--- a/src/qlm_import_surface.F90
+++ b/src/qlm_import_surface.F90
@@ -57,38 +57,6 @@ subroutine qlm_import_surface (CCTK_ARGUMENTS, hn)
      
      qlm_shape_p_p(:,:,hn) = qlm_shape_p(:,:,hn)
      qlm_shape_p  (:,:,hn) = qlm_shape  (:,:,hn)
-       
-     qlm_l0_p_p(:,:,hn) = qlm_l0_p(:,:,hn)
-     qlm_l0_p  (:,:,hn) = qlm_l0  (:,:,hn)
-     qlm_l1_p_p(:,:,hn) = qlm_l1_p(:,:,hn)
-     qlm_l1_p  (:,:,hn) = qlm_l1  (:,:,hn)
-     qlm_l2_p_p(:,:,hn) = qlm_l2_p(:,:,hn)
-     qlm_l2_p  (:,:,hn) = qlm_l2  (:,:,hn)
-     qlm_l3_p_p(:,:,hn) = qlm_l3_p(:,:,hn)
-     qlm_l3_p  (:,:,hn) = qlm_l3  (:,:,hn)
-     
-     qlm_n0_p_p(:,:,hn) = qlm_n0_p(:,:,hn)
-     qlm_n0_p  (:,:,hn) = qlm_n0  (:,:,hn)
-     qlm_n1_p_p(:,:,hn) = qlm_n1_p(:,:,hn)
-     qlm_n1_p  (:,:,hn) = qlm_n1  (:,:,hn)
-     qlm_n2_p_p(:,:,hn) = qlm_n2_p(:,:,hn)
-     qlm_n2_p  (:,:,hn) = qlm_n2  (:,:,hn)
-     qlm_n3_p_p(:,:,hn) = qlm_n3_p(:,:,hn)
-     qlm_n3_p  (:,:,hn) = qlm_n3  (:,:,hn)
-     
-     qlm_m0_p_p(:,:,hn) = qlm_m0_p(:,:,hn)
-     qlm_m0_p  (:,:,hn) = qlm_m0  (:,:,hn)
-     qlm_m1_p_p(:,:,hn) = qlm_m1_p(:,:,hn)
-     qlm_m1_p  (:,:,hn) = qlm_m1  (:,:,hn)
-     qlm_m2_p_p(:,:,hn) = qlm_m2_p(:,:,hn)
-     qlm_m2_p  (:,:,hn) = qlm_m2  (:,:,hn)
-     qlm_m3_p_p(:,:,hn) = qlm_m3_p(:,:,hn)
-     qlm_m3_p  (:,:,hn) = qlm_m3  (:,:,hn)
-     
-     qlm_xi_t_p_p(:,:,hn) = qlm_xi_t_p(:,:,hn)
-     qlm_xi_t_p  (:,:,hn) = qlm_xi_t  (:,:,hn)
-     qlm_xi_p_p_p(:,:,hn) = qlm_xi_p_p(:,:,hn)
-     qlm_xi_p_p  (:,:,hn) = qlm_xi_p  (:,:,hn)
      
   end if
   
diff --git a/src/qlm_interpolate.F90 b/src/qlm_interpolate.F90
index 82ec27e..8317de7 100644
--- a/src/qlm_interpolate.F90
+++ b/src/qlm_interpolate.F90
@@ -29,6 +29,8 @@ subroutine qlm_interpolate (CCTK_ARGUMENTS, hn)
   integer,   parameter :: sk = kind(interpolator)
   CCTK_REAL, parameter :: one = 1
   
+  CCTK_REAL, parameter :: poison_value = -42
+  
   integer      :: len_coordsystem
   integer      :: len_interpolator
   integer      :: len_interpolator_options
@@ -54,14 +56,17 @@ subroutine qlm_interpolate (CCTK_ARGUMENTS, hn)
   integer      :: ind_kxx, ind_kxy, ind_kxz, ind_kyy, ind_kyz, ind_kzz
   integer      :: ind_alpha
   integer      :: ind_betax, ind_betay, ind_betaz
+  integer      :: ind_ttt
+  integer      :: ind_ttx, ind_tty, ind_ttz
+  integer      :: ind_txx, ind_txy, ind_txz, ind_tyy, ind_tyz, ind_tzz
   
   integer      :: coord_type
   CCTK_POINTER :: coords(3)
-  CCTK_INT     :: inputs(16)
-  CCTK_INT     :: output_types(88)
-  CCTK_POINTER :: outputs(88)
-  CCTK_INT     :: operand_indices(88)
-  CCTK_INT     :: operation_codes(88)
+  CCTK_INT     :: inputs(26)
+  CCTK_INT     :: output_types(98)
+  CCTK_POINTER :: outputs(98)
+  CCTK_INT     :: operand_indices(98)
+  CCTK_INT     :: operation_codes(98)
   integer      :: npoints
   
   character    :: msg*1000
@@ -82,6 +87,10 @@ subroutine qlm_interpolate (CCTK_ARGUMENTS, hn)
      call CCTK_WARN (0, "The shift must have storage")
   end if
   
+!!$  if (stress_energy_state==0) then
+!!$     call CCTK_WARN (0, "The stress-energy tensor must have storage")
+!!$  end if
+  
   
   
   ! Get coordinate system
@@ -137,22 +146,45 @@ subroutine qlm_interpolate (CCTK_ARGUMENTS, hn)
   ! TODO: check the excision mask
   
   ! Get variable indices
-  call CCTK_VarIndex (ind_gxx  , "ADMBase::gxx"  )
-  call CCTK_VarIndex (ind_gxy  , "ADMBase::gxy"  )
-  call CCTK_VarIndex (ind_gxz  , "ADMBase::gxz"  )
-  call CCTK_VarIndex (ind_gyy  , "ADMBase::gyy"  )
-  call CCTK_VarIndex (ind_gyz  , "ADMBase::gyz"  )
-  call CCTK_VarIndex (ind_gzz  , "ADMBase::gzz"  )
-  call CCTK_VarIndex (ind_kxx  , "ADMBase::kxx"  )
-  call CCTK_VarIndex (ind_kxy  , "ADMBase::kxy"  )
-  call CCTK_VarIndex (ind_kxz  , "ADMBase::kxz"  )
-  call CCTK_VarIndex (ind_kyy  , "ADMBase::kyy"  )
-  call CCTK_VarIndex (ind_kyz  , "ADMBase::kyz"  )
-  call CCTK_VarIndex (ind_kzz  , "ADMBase::kzz"  )
-  call CCTK_VarIndex (ind_alpha, "ADMBase::alp"  )
-  call CCTK_VarIndex (ind_betax, "ADMBase::betax")
-  call CCTK_VarIndex (ind_betay, "ADMBase::betay")
-  call CCTK_VarIndex (ind_betaz, "ADMBase::betaz")
+  call CCTK_VarIndex (ind_gxx  , "ADMBase::gxx"   )
+  call CCTK_VarIndex (ind_gxy  , "ADMBase::gxy"   )
+  call CCTK_VarIndex (ind_gxz  , "ADMBase::gxz"   )
+  call CCTK_VarIndex (ind_gyy  , "ADMBase::gyy"   )
+  call CCTK_VarIndex (ind_gyz  , "ADMBase::gyz"   )
+  call CCTK_VarIndex (ind_gzz  , "ADMBase::gzz"   )
+  call CCTK_VarIndex (ind_kxx  , "ADMBase::kxx"   )
+  call CCTK_VarIndex (ind_kxy  , "ADMBase::kxy"   )
+  call CCTK_VarIndex (ind_kxz  , "ADMBase::kxz"   )
+  call CCTK_VarIndex (ind_kyy  , "ADMBase::kyy"   )
+  call CCTK_VarIndex (ind_kyz  , "ADMBase::kyz"   )
+  call CCTK_VarIndex (ind_kzz  , "ADMBase::kzz"   )
+  call CCTK_VarIndex (ind_alpha, "ADMBase::alp"   )
+  call CCTK_VarIndex (ind_betax, "ADMBase::betax" )
+  call CCTK_VarIndex (ind_betay, "ADMBase::betay" )
+  call CCTK_VarIndex (ind_betaz, "ADMBase::betaz" )
+  if (stress_energy_state /= 0) then
+     call CCTK_VarIndex (ind_ttt  , "TmunuBase::eTtt")
+     call CCTK_VarIndex (ind_ttx  , "TmunuBase::eTtx")
+     call CCTK_VarIndex (ind_tty  , "TmunuBase::eTty")
+     call CCTK_VarIndex (ind_ttz  , "TmunuBase::eTtz")
+     call CCTK_VarIndex (ind_txx  , "TmunuBase::eTxx")
+     call CCTK_VarIndex (ind_txy  , "TmunuBase::eTxy")
+     call CCTK_VarIndex (ind_txz  , "TmunuBase::eTxz")
+     call CCTK_VarIndex (ind_tyy  , "TmunuBase::eTyy")
+     call CCTK_VarIndex (ind_tyz  , "TmunuBase::eTyz")
+     call CCTK_VarIndex (ind_tzz  , "TmunuBase::eTzz")
+  else
+     ind_ttt = -1
+     ind_ttx = -1
+     ind_tty = -1
+     ind_ttz = -1
+     ind_txx = -1
+     ind_txy = -1
+     ind_txz = -1
+     ind_tyy = -1
+     ind_tyz = -1
+     ind_tzz = -1
+  end if
   
   
   
@@ -170,7 +202,11 @@ subroutine qlm_interpolate (CCTK_ARGUMENTS, hn)
        ind_gxx, ind_gxy, ind_gxz, ind_gyy, ind_gyz, ind_gzz, &
        ind_kxx, ind_kxy, ind_kxz, ind_kyy, ind_kyz, ind_kzz, &
        ind_alpha, &
-       ind_betax, ind_betay, ind_betaz /)
+       ind_betax, ind_betay, ind_betaz, &
+       ind_ttt, &
+       ind_ttx, ind_tty, ind_ttz, &
+       ind_txx, ind_txy, ind_txz, ind_tyy, ind_tyz, ind_tzz /)
+  
   call CCTK_NumVars (nvars)
   if (nvars < 0) call CCTK_WARN (0, "internal error")
   if (any(inputs /= -1 .and. (inputs < 0 .or. inputs >= nvars))) then
@@ -193,7 +229,10 @@ subroutine qlm_interpolate (CCTK_ARGUMENTS, hn)
        06, 07, 08, 09, 10, 11, &
        06, 07, 08, 09, 10, 11, &
        12, &                    ! alp
-       13, 14, 15 /)            ! beta^i
+       13, 14, 15, &            ! beta^i
+       16, &                    ! T_tt
+       17, 18, 19, &            ! T_ti
+       20, 21, 22, 23, 24, 25 /) ! T_ij
   
   operation_codes = (/ &
        0, 0, 0, 0, 0, 0, &      ! g_ij
@@ -211,7 +250,10 @@ subroutine qlm_interpolate (CCTK_ARGUMENTS, hn)
        2, 2, 2, 2, 2, 2, &
        3, 3, 3, 3, 3, 3, &
        0, &                     ! alp
-       0, 0, 0 /)               ! beta^i
+       0, 0, 0, &               ! beta^i
+       0, &                     ! T_tt
+       0, 0, 0, &               ! T_ti
+       0, 0, 0, 0, 0, 0 /)      ! T_ij
   
   output_types(:) = CCTK_VARIABLE_REAL
   if (hn > 0) then
@@ -231,7 +273,10 @@ subroutine qlm_interpolate (CCTK_ARGUMENTS, hn)
           P(qlm_dkxxy), P(qlm_dkxyy), P(qlm_dkxzy), P(qlm_dkyyy), P(qlm_dkyzy), P(qlm_dkzzy), &
           P(qlm_dkxxz), P(qlm_dkxyz), P(qlm_dkxzz), P(qlm_dkyyz), P(qlm_dkyzz), P(qlm_dkzzz), &
           P(qlm_alpha), &
-          P(qlm_betax), P(qlm_betay), P(qlm_betaz) /)
+          P(qlm_betax), P(qlm_betay), P(qlm_betaz), &
+          P(qlm_ttt), &
+          P(qlm_ttx), P(qlm_tty), P(qlm_ttz), &
+          P(qlm_txx), P(qlm_txy), P(qlm_txz), P(qlm_tyy), P(qlm_tyz), P(qlm_tzz) /)
   else
      outputs(:) = CCTK_NullPointer()
   end if
@@ -245,95 +290,104 @@ subroutine qlm_interpolate (CCTK_ARGUMENTS, hn)
   
 #if 0
   ! Poison the output variables
-#define poison -42
-  qlm_gxx = poison
-  qlm_gxy = poison
-  qlm_gxz = poison
-  qlm_gyy = poison
-  qlm_gyz = poison
-  qlm_gzz = poison
-  qlm_dgxxx = poison
-  qlm_dgxyx = poison
-  qlm_dgxzx = poison
-  qlm_dgyyx = poison
-  qlm_dgyzx = poison
-  qlm_dgzzx = poison
-  qlm_dgxxy = poison
-  qlm_dgxyy = poison
-  qlm_dgxzy = poison
-  qlm_dgyyy = poison
-  qlm_dgyzy = poison
-  qlm_dgzzy = poison
-  qlm_dgxxz = poison
-  qlm_dgxyz = poison
-  qlm_dgxzz = poison
-  qlm_dgyyz = poison
-  qlm_dgyzz = poison
-  qlm_dgzzz = poison
-  qlm_ddgxxxx = poison
-  qlm_ddgxyxx = poison
-  qlm_ddgxzxx = poison
-  qlm_ddgyyxx = poison
-  qlm_ddgyzxx = poison
-  qlm_ddgzzxx = poison
-  qlm_ddgxxxy = poison
-  qlm_ddgxyxy = poison
-  qlm_ddgxzxy = poison
-  qlm_ddgyyxy = poison
-  qlm_ddgyzxy = poison
-  qlm_ddgzzxy = poison
-  qlm_ddgxxxz = poison
-  qlm_ddgxyxz = poison
-  qlm_ddgxzxz = poison
-  qlm_ddgyyxz = poison
-  qlm_ddgyzxz = poison
-  qlm_ddgzzxz = poison
-  qlm_ddgxxyy = poison
-  qlm_ddgxyyy = poison
-  qlm_ddgxzyy = poison
-  qlm_ddgyyyy = poison
-  qlm_ddgyzyy = poison
-  qlm_ddgzzyy = poison
-  qlm_ddgxxyz = poison
-  qlm_ddgxyyz = poison
-  qlm_ddgxzyz = poison
-  qlm_ddgyyyz = poison
-  qlm_ddgyzyz = poison
-  qlm_ddgzzyz = poison
-  qlm_ddgxxzz = poison
-  qlm_ddgxyzz = poison
-  qlm_ddgxzzz = poison
-  qlm_ddgyyzz = poison
-  qlm_ddgyzzz = poison
-  qlm_ddgzzzz = poison
-  qlm_kxx = poison
-  qlm_kxy = poison
-  qlm_kxz = poison
-  qlm_kyy = poison
-  qlm_kyz = poison
-  qlm_kzz = poison
-  qlm_dkxxx = poison
-  qlm_dkxyx = poison
-  qlm_dkxzx = poison
-  qlm_dkyyx = poison
-  qlm_dkyzx = poison
-  qlm_dkzzx = poison
-  qlm_dkxxy = poison
-  qlm_dkxyy = poison
-  qlm_dkxzy = poison
-  qlm_dkyyy = poison
-  qlm_dkyzy = poison
-  qlm_dkzzy = poison
-  qlm_dkxxz = poison
-  qlm_dkxyz = poison
-  qlm_dkxzz = poison
-  qlm_dkyyz = poison
-  qlm_dkyzz = poison
-  qlm_dkzzz = poison
-  qlm_alpha = poison
-  qlm_betax = poison
-  qlm_betay = poison
-  qlm_betaz = poison
+  call poison (qlm_gxx    )
+  call poison (qlm_gxy    )
+  call poison (qlm_gxz    )
+  call poison (qlm_gyy    )
+  call poison (qlm_gyz    )
+  call poison (qlm_gzz    )
+  call poison (qlm_dgxxx  )
+  call poison (qlm_dgxyx  )
+  call poison (qlm_dgxzx  )
+  call poison (qlm_dgyyx  )
+  call poison (qlm_dgyzx  )
+  call poison (qlm_dgzzx  )
+  call poison (qlm_dgxxy  )
+  call poison (qlm_dgxyy  )
+  call poison (qlm_dgxzy  )
+  call poison (qlm_dgyyy  )
+  call poison (qlm_dgyzy  )
+  call poison (qlm_dgzzy  )
+  call poison (qlm_dgxxz  )
+  call poison (qlm_dgxyz  )
+  call poison (qlm_dgxzz  )
+  call poison (qlm_dgyyz  )
+  call poison (qlm_dgyzz  )
+  call poison (qlm_dgzzz  )
+  call poison (qlm_ddgxxxx)
+  call poison (qlm_ddgxyxx)
+  call poison (qlm_ddgxzxx)
+  call poison (qlm_ddgyyxx)
+  call poison (qlm_ddgyzxx)
+  call poison (qlm_ddgzzxx)
+  call poison (qlm_ddgxxxy)
+  call poison (qlm_ddgxyxy)
+  call poison (qlm_ddgxzxy)
+  call poison (qlm_ddgyyxy)
+  call poison (qlm_ddgyzxy)
+  call poison (qlm_ddgzzxy)
+  call poison (qlm_ddgxxxz)
+  call poison (qlm_ddgxyxz)
+  call poison (qlm_ddgxzxz)
+  call poison (qlm_ddgyyxz)
+  call poison (qlm_ddgyzxz)
+  call poison (qlm_ddgzzxz)
+  call poison (qlm_ddgxxyy)
+  call poison (qlm_ddgxyyy)
+  call poison (qlm_ddgxzyy)
+  call poison (qlm_ddgyyyy)
+  call poison (qlm_ddgyzyy)
+  call poison (qlm_ddgzzyy)
+  call poison (qlm_ddgxxyz)
+  call poison (qlm_ddgxyyz)
+  call poison (qlm_ddgxzyz)
+  call poison (qlm_ddgyyyz)
+  call poison (qlm_ddgyzyz)
+  call poison (qlm_ddgzzyz)
+  call poison (qlm_ddgxxzz)
+  call poison (qlm_ddgxyzz)
+  call poison (qlm_ddgxzzz)
+  call poison (qlm_ddgyyzz)
+  call poison (qlm_ddgyzzz)
+  call poison (qlm_ddgzzzz)
+  call poison (qlm_kxx    )
+  call poison (qlm_kxy    )
+  call poison (qlm_kxz    )
+  call poison (qlm_kyy    )
+  call poison (qlm_kyz    )
+  call poison (qlm_kzz    )
+  call poison (qlm_dkxxx  )
+  call poison (qlm_dkxyx  )
+  call poison (qlm_dkxzx  )
+  call poison (qlm_dkyyx  )
+  call poison (qlm_dkyzx  )
+  call poison (qlm_dkzzx  )
+  call poison (qlm_dkxxy  )
+  call poison (qlm_dkxyy  )
+  call poison (qlm_dkxzy  )
+  call poison (qlm_dkyyy  )
+  call poison (qlm_dkyzy  )
+  call poison (qlm_dkzzy  )
+  call poison (qlm_dkxxz  )
+  call poison (qlm_dkxyz  )
+  call poison (qlm_dkxzz  )
+  call poison (qlm_dkyyz  )
+  call poison (qlm_dkyzz  )
+  call poison (qlm_dkzzz  )
+  call poison (qlm_alpha  )
+  call poison (qlm_betax  )
+  call poison (qlm_betay  )
+  call poison (qlm_betaz  )
+  call poison (qlm_ttt    )
+  call poison (qlm_ttx    )
+  call poison (qlm_tty    )
+  call poison (qlm_ttz    )
+  call poison (qlm_txx    )
+  call poison (qlm_txy    )
+  call poison (qlm_txz    )
+  call poison (qlm_tyy    )
+  call poison (qlm_tyz    )
+  call poison (qlm_tzz    )
 #endif
   
   
@@ -368,187 +422,220 @@ subroutine qlm_interpolate (CCTK_ARGUMENTS, hn)
   ! Unpack the variables
   if (hn > 0) then
      
-     call unpack (qlm_gxx         , ni, nj)
-     call unpack (qlm_gxy         , ni, nj)
-     call unpack (qlm_gxz         , ni, nj)
-     call unpack (qlm_gyy         , ni, nj)
-     call unpack (qlm_gyz         , ni, nj)
-     call unpack (qlm_gzz         , ni, nj)
-     call unpack (qlm_dgxxx       , ni, nj)
-     call unpack (qlm_dgxyx       , ni, nj)
-     call unpack (qlm_dgxzx       , ni, nj)
-     call unpack (qlm_dgyyx       , ni, nj)
-     call unpack (qlm_dgyzx       , ni, nj)
-     call unpack (qlm_dgzzx       , ni, nj)
-     call unpack (qlm_dgxxy       , ni, nj)
-     call unpack (qlm_dgxyy       , ni, nj)
-     call unpack (qlm_dgxzy       , ni, nj)
-     call unpack (qlm_dgyyy       , ni, nj)
-     call unpack (qlm_dgyzy       , ni, nj)
-     call unpack (qlm_dgzzy       , ni, nj)
-     call unpack (qlm_dgxxz       , ni, nj)
-     call unpack (qlm_dgxyz       , ni, nj)
-     call unpack (qlm_dgxzz       , ni, nj)
-     call unpack (qlm_dgyyz       , ni, nj)
-     call unpack (qlm_dgyzz       , ni, nj)
-     call unpack (qlm_dgzzz       , ni, nj)
-     call unpack (qlm_ddgxxxx     , ni, nj)
-     call unpack (qlm_ddgxyxx     , ni, nj)
-     call unpack (qlm_ddgxzxx     , ni, nj)
-     call unpack (qlm_ddgyyxx     , ni, nj)
-     call unpack (qlm_ddgyzxx     , ni, nj)
-     call unpack (qlm_ddgzzxx     , ni, nj)
-     call unpack (qlm_ddgxxxy     , ni, nj)
-     call unpack (qlm_ddgxyxy     , ni, nj)
-     call unpack (qlm_ddgxzxy     , ni, nj)
-     call unpack (qlm_ddgyyxy     , ni, nj)
-     call unpack (qlm_ddgyzxy     , ni, nj)
-     call unpack (qlm_ddgzzxy     , ni, nj)
-     call unpack (qlm_ddgxxxz     , ni, nj)
-     call unpack (qlm_ddgxyxz     , ni, nj)
-     call unpack (qlm_ddgxzxz     , ni, nj)
-     call unpack (qlm_ddgyyxz     , ni, nj)
-     call unpack (qlm_ddgyzxz     , ni, nj)
-     call unpack (qlm_ddgzzxz     , ni, nj)
-     call unpack (qlm_ddgxxyy     , ni, nj)
-     call unpack (qlm_ddgxyyy     , ni, nj)
-     call unpack (qlm_ddgxzyy     , ni, nj)
-     call unpack (qlm_ddgyyyy     , ni, nj)
-     call unpack (qlm_ddgyzyy     , ni, nj)
-     call unpack (qlm_ddgzzyy     , ni, nj)
-     call unpack (qlm_ddgxxyz     , ni, nj)
-     call unpack (qlm_ddgxyyz     , ni, nj)
-     call unpack (qlm_ddgxzyz     , ni, nj)
-     call unpack (qlm_ddgyyyz     , ni, nj)
-     call unpack (qlm_ddgyzyz     , ni, nj)
-     call unpack (qlm_ddgzzyz     , ni, nj)
-     call unpack (qlm_ddgxxzz     , ni, nj)
-     call unpack (qlm_ddgxyzz     , ni, nj)
-     call unpack (qlm_ddgxzzz     , ni, nj)
-     call unpack (qlm_ddgyyzz     , ni, nj)
-     call unpack (qlm_ddgyzzz     , ni, nj)
-     call unpack (qlm_ddgzzzz     , ni, nj)
-     call unpack (qlm_kxx         , ni, nj)
-     call unpack (qlm_kxy         , ni, nj)
-     call unpack (qlm_kxz         , ni, nj)
-     call unpack (qlm_kyy         , ni, nj)
-     call unpack (qlm_kyz         , ni, nj)
-     call unpack (qlm_kzz         , ni, nj)
-     call unpack (qlm_dkxxx       , ni, nj)
-     call unpack (qlm_dkxyx       , ni, nj)
-     call unpack (qlm_dkxzx       , ni, nj)
-     call unpack (qlm_dkyyx       , ni, nj)
-     call unpack (qlm_dkyzx       , ni, nj)
-     call unpack (qlm_dkzzx       , ni, nj)
-     call unpack (qlm_dkxxy       , ni, nj)
-     call unpack (qlm_dkxyy       , ni, nj)
-     call unpack (qlm_dkxzy       , ni, nj)
-     call unpack (qlm_dkyyy       , ni, nj)
-     call unpack (qlm_dkyzy       , ni, nj)
-     call unpack (qlm_dkzzy       , ni, nj)
-     call unpack (qlm_dkxxz       , ni, nj)
-     call unpack (qlm_dkxyz       , ni, nj)
-     call unpack (qlm_dkxzz       , ni, nj)
-     call unpack (qlm_dkyyz       , ni, nj)
-     call unpack (qlm_dkyzz       , ni, nj)
-     call unpack (qlm_dkzzz       , ni, nj)
-     call unpack (qlm_alpha       , ni, nj)
-     call unpack (qlm_betax       , ni, nj)
-     call unpack (qlm_betay       , ni, nj)
-     call unpack (qlm_betaz       , ni, nj)
+     call unpack (qlm_gxx    , ni, nj)
+     call unpack (qlm_gxy    , ni, nj)
+     call unpack (qlm_gxz    , ni, nj)
+     call unpack (qlm_gyy    , ni, nj)
+     call unpack (qlm_gyz    , ni, nj)
+     call unpack (qlm_gzz    , ni, nj)
+     call unpack (qlm_dgxxx  , ni, nj)
+     call unpack (qlm_dgxyx  , ni, nj)
+     call unpack (qlm_dgxzx  , ni, nj)
+     call unpack (qlm_dgyyx  , ni, nj)
+     call unpack (qlm_dgyzx  , ni, nj)
+     call unpack (qlm_dgzzx  , ni, nj)
+     call unpack (qlm_dgxxy  , ni, nj)
+     call unpack (qlm_dgxyy  , ni, nj)
+     call unpack (qlm_dgxzy  , ni, nj)
+     call unpack (qlm_dgyyy  , ni, nj)
+     call unpack (qlm_dgyzy  , ni, nj)
+     call unpack (qlm_dgzzy  , ni, nj)
+     call unpack (qlm_dgxxz  , ni, nj)
+     call unpack (qlm_dgxyz  , ni, nj)
+     call unpack (qlm_dgxzz  , ni, nj)
+     call unpack (qlm_dgyyz  , ni, nj)
+     call unpack (qlm_dgyzz  , ni, nj)
+     call unpack (qlm_dgzzz  , ni, nj)
+     call unpack (qlm_ddgxxxx, ni, nj)
+     call unpack (qlm_ddgxyxx, ni, nj)
+     call unpack (qlm_ddgxzxx, ni, nj)
+     call unpack (qlm_ddgyyxx, ni, nj)
+     call unpack (qlm_ddgyzxx, ni, nj)
+     call unpack (qlm_ddgzzxx, ni, nj)
+     call unpack (qlm_ddgxxxy, ni, nj)
+     call unpack (qlm_ddgxyxy, ni, nj)
+     call unpack (qlm_ddgxzxy, ni, nj)
+     call unpack (qlm_ddgyyxy, ni, nj)
+     call unpack (qlm_ddgyzxy, ni, nj)
+     call unpack (qlm_ddgzzxy, ni, nj)
+     call unpack (qlm_ddgxxxz, ni, nj)
+     call unpack (qlm_ddgxyxz, ni, nj)
+     call unpack (qlm_ddgxzxz, ni, nj)
+     call unpack (qlm_ddgyyxz, ni, nj)
+     call unpack (qlm_ddgyzxz, ni, nj)
+     call unpack (qlm_ddgzzxz, ni, nj)
+     call unpack (qlm_ddgxxyy, ni, nj)
+     call unpack (qlm_ddgxyyy, ni, nj)
+     call unpack (qlm_ddgxzyy, ni, nj)
+     call unpack (qlm_ddgyyyy, ni, nj)
+     call unpack (qlm_ddgyzyy, ni, nj)
+     call unpack (qlm_ddgzzyy, ni, nj)
+     call unpack (qlm_ddgxxyz, ni, nj)
+     call unpack (qlm_ddgxyyz, ni, nj)
+     call unpack (qlm_ddgxzyz, ni, nj)
+     call unpack (qlm_ddgyyyz, ni, nj)
+     call unpack (qlm_ddgyzyz, ni, nj)
+     call unpack (qlm_ddgzzyz, ni, nj)
+     call unpack (qlm_ddgxxzz, ni, nj)
+     call unpack (qlm_ddgxyzz, ni, nj)
+     call unpack (qlm_ddgxzzz, ni, nj)
+     call unpack (qlm_ddgyyzz, ni, nj)
+     call unpack (qlm_ddgyzzz, ni, nj)
+     call unpack (qlm_ddgzzzz, ni, nj)
+     call unpack (qlm_kxx    , ni, nj)
+     call unpack (qlm_kxy    , ni, nj)
+     call unpack (qlm_kxz    , ni, nj)
+     call unpack (qlm_kyy    , ni, nj)
+     call unpack (qlm_kyz    , ni, nj)
+     call unpack (qlm_kzz    , ni, nj)
+     call unpack (qlm_dkxxx  , ni, nj)
+     call unpack (qlm_dkxyx  , ni, nj)
+     call unpack (qlm_dkxzx  , ni, nj)
+     call unpack (qlm_dkyyx  , ni, nj)
+     call unpack (qlm_dkyzx  , ni, nj)
+     call unpack (qlm_dkzzx  , ni, nj)
+     call unpack (qlm_dkxxy  , ni, nj)
+     call unpack (qlm_dkxyy  , ni, nj)
+     call unpack (qlm_dkxzy  , ni, nj)
+     call unpack (qlm_dkyyy  , ni, nj)
+     call unpack (qlm_dkyzy  , ni, nj)
+     call unpack (qlm_dkzzy  , ni, nj)
+     call unpack (qlm_dkxxz  , ni, nj)
+     call unpack (qlm_dkxyz  , ni, nj)
+     call unpack (qlm_dkxzz  , ni, nj)
+     call unpack (qlm_dkyyz  , ni, nj)
+     call unpack (qlm_dkyzz  , ni, nj)
+     call unpack (qlm_dkzzz  , ni, nj)
+     call unpack (qlm_alpha  , ni, nj)
+     call unpack (qlm_betax  , ni, nj)
+     call unpack (qlm_betay  , ni, nj)
+     call unpack (qlm_betaz  , ni, nj)
+     if (stress_energy_state /= 0) then
+        call unpack (qlm_ttt    , ni, nj)
+        call unpack (qlm_ttx    , ni, nj)
+        call unpack (qlm_tty    , ni, nj)
+        call unpack (qlm_ttz    , ni, nj)
+        call unpack (qlm_txx    , ni, nj)
+        call unpack (qlm_txy    , ni, nj)
+        call unpack (qlm_txz    , ni, nj)
+        call unpack (qlm_tyy    , ni, nj)
+        call unpack (qlm_tyz    , ni, nj)
+        call unpack (qlm_tzz    , ni, nj)
+     else
+        qlm_ttt = 0
+        qlm_ttx = 0
+        qlm_tty = 0
+        qlm_ttz = 0
+        qlm_txx = 0
+        qlm_txy = 0
+        qlm_txz = 0
+        qlm_tyy = 0
+        qlm_tyz = 0
+        qlm_tzz = 0
+     end if
      
      
      
 #if 0
      ! Check for poison
-     if (any(qlm_gxx == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_gxy == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_gxz == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_gyy == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_gyz == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_gzz == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_dgxxx == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_dgxyx == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_dgxzx == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_dgyyx == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_dgyzx == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_dgzzx == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_dgxxy == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_dgxyy == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_dgxzy == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_dgyyy == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_dgyzy == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_dgzzy == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_dgxxz == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_dgxyz == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_dgxzz == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_dgyyz == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_dgyzz == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_dgzzz == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_ddgxxxx == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_ddgxyxx == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_ddgxzxx == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_ddgyyxx == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_ddgyzxx == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_ddgzzxx == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_ddgxxxy == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_ddgxyxy == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_ddgxzxy == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_ddgyyxy == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_ddgyzxy == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_ddgzzxy == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_ddgxxxz == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_ddgxyxz == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_ddgxzxz == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_ddgyyxz == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_ddgyzxz == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_ddgzzxz == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_ddgxxyy == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_ddgxyyy == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_ddgxzyy == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_ddgyyyy == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_ddgyzyy == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_ddgzzyy == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_ddgxxyz == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_ddgxyyz == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_ddgxzyz == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_ddgyyyz == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_ddgyzyz == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_ddgzzyz == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_ddgxxzz == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_ddgxyzz == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_ddgxzzz == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_ddgyyzz == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_ddgyzzz == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_ddgzzzz == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_kxx == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_kxy == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_kxz == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_kyy == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_kyz == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_kzz == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_dkxxx == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_dkxyx == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_dkxzx == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_dkyyx == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_dkyzx == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_dkzzx == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_dkxxy == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_dkxyy == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_dkxzy == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_dkyyy == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_dkyzy == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_dkzzy == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_dkxxz == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_dkxyz == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_dkxzz == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_dkyyz == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_dkyzz == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_dkzzz == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_alpha == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_betax == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_betay == poison)) call CCTK_WARN (0, "poison found")
-     if (any(qlm_betaz == poison)) call CCTK_WARN (0, "poison found")
+     call poison_check (qlm_gxx    , "qlm_gxx    ")
+     call poison_check (qlm_gxy    , "qlm_gxy    ")
+     call poison_check (qlm_gxz    , "qlm_gxz    ")
+     call poison_check (qlm_gyy    , "qlm_gyy    ")
+     call poison_check (qlm_gyz    , "qlm_gyz    ")
+     call poison_check (qlm_gzz    , "qlm_gzz    ")
+     call poison_check (qlm_dgxxx  , "qlm_dgxxx  ")
+     call poison_check (qlm_dgxyx  , "qlm_dgxyx  ")
+     call poison_check (qlm_dgxzx  , "qlm_dgxzx  ")
+     call poison_check (qlm_dgyyx  , "qlm_dgyyx  ")
+     call poison_check (qlm_dgyzx  , "qlm_dgyzx  ")
+     call poison_check (qlm_dgzzx  , "qlm_dgzzx  ")
+     call poison_check (qlm_dgxxy  , "qlm_dgxxy  ")
+     call poison_check (qlm_dgxyy  , "qlm_dgxyy  ")
+     call poison_check (qlm_dgxzy  , "qlm_dgxzy  ")
+     call poison_check (qlm_dgyyy  , "qlm_dgyyy  ")
+     call poison_check (qlm_dgyzy  , "qlm_dgyzy  ")
+     call poison_check (qlm_dgzzy  , "qlm_dgzzy  ")
+     call poison_check (qlm_dgxxz  , "qlm_dgxxz  ")
+     call poison_check (qlm_dgxyz  , "qlm_dgxyz  ")
+     call poison_check (qlm_dgxzz  , "qlm_dgxzz  ")
+     call poison_check (qlm_dgyyz  , "qlm_dgyyz  ")
+     call poison_check (qlm_dgyzz  , "qlm_dgyzz  ")
+     call poison_check (qlm_dgzzz  , "qlm_dgzzz  ")
+     call poison_check (qlm_ddgxxxx, "qlm_ddgxxxx")
+     call poison_check (qlm_ddgxyxx, "qlm_ddgxyxx")
+     call poison_check (qlm_ddgxzxx, "qlm_ddgxzxx")
+     call poison_check (qlm_ddgyyxx, "qlm_ddgyyxx")
+     call poison_check (qlm_ddgyzxx, "qlm_ddgyzxx")
+     call poison_check (qlm_ddgzzxx, "qlm_ddgzzxx")
+     call poison_check (qlm_ddgxxxy, "qlm_ddgxxxy")
+     call poison_check (qlm_ddgxyxy, "qlm_ddgxyxy")
+     call poison_check (qlm_ddgxzxy, "qlm_ddgxzxy")
+     call poison_check (qlm_ddgyyxy, "qlm_ddgyyxy")
+     call poison_check (qlm_ddgyzxy, "qlm_ddgyzxy")
+     call poison_check (qlm_ddgzzxy, "qlm_ddgzzxy")
+     call poison_check (qlm_ddgxxxz, "qlm_ddgxxxz")
+     call poison_check (qlm_ddgxyxz, "qlm_ddgxyxz")
+     call poison_check (qlm_ddgxzxz, "qlm_ddgxzxz")
+     call poison_check (qlm_ddgyyxz, "qlm_ddgyyxz")
+     call poison_check (qlm_ddgyzxz, "qlm_ddgyzxz")
+     call poison_check (qlm_ddgzzxz, "qlm_ddgzzxz")
+     call poison_check (qlm_ddgxxyy, "qlm_ddgxxyy")
+     call poison_check (qlm_ddgxyyy, "qlm_ddgxyyy")
+     call poison_check (qlm_ddgxzyy, "qlm_ddgxzyy")
+     call poison_check (qlm_ddgyyyy, "qlm_ddgyyyy")
+     call poison_check (qlm_ddgyzyy, "qlm_ddgyzyy")
+     call poison_check (qlm_ddgzzyy, "qlm_ddgzzyy")
+     call poison_check (qlm_ddgxxyz, "qlm_ddgxxyz")
+     call poison_check (qlm_ddgxyyz, "qlm_ddgxyyz")
+     call poison_check (qlm_ddgxzyz, "qlm_ddgxzyz")
+     call poison_check (qlm_ddgyyyz, "qlm_ddgyyyz")
+     call poison_check (qlm_ddgyzyz, "qlm_ddgyzyz")
+     call poison_check (qlm_ddgzzyz, "qlm_ddgzzyz")
+     call poison_check (qlm_ddgxxzz, "qlm_ddgxxzz")
+     call poison_check (qlm_ddgxyzz, "qlm_ddgxyzz")
+     call poison_check (qlm_ddgxzzz, "qlm_ddgxzzz")
+     call poison_check (qlm_ddgyyzz, "qlm_ddgyyzz")
+     call poison_check (qlm_ddgyzzz, "qlm_ddgyzzz")
+     call poison_check (qlm_ddgzzzz, "qlm_ddgzzzz")
+     call poison_check (qlm_kxx    , "qlm_kxx    ")
+     call poison_check (qlm_kxy    , "qlm_kxy    ")
+     call poison_check (qlm_kxz    , "qlm_kxz    ")
+     call poison_check (qlm_kyy    , "qlm_kyy    ")
+     call poison_check (qlm_kyz    , "qlm_kyz    ")
+     call poison_check (qlm_kzz    , "qlm_kzz    ")
+     call poison_check (qlm_dkxxx  , "qlm_dkxxx  ")
+     call poison_check (qlm_dkxyx  , "qlm_dkxyx  ")
+     call poison_check (qlm_dkxzx  , "qlm_dkxzx  ")
+     call poison_check (qlm_dkyyx  , "qlm_dkyyx  ")
+     call poison_check (qlm_dkyzx  , "qlm_dkyzx  ")
+     call poison_check (qlm_dkzzx  , "qlm_dkzzx  ")
+     call poison_check (qlm_dkxxy  , "qlm_dkxxy  ")
+     call poison_check (qlm_dkxyy  , "qlm_dkxyy  ")
+     call poison_check (qlm_dkxzy  , "qlm_dkxzy  ")
+     call poison_check (qlm_dkyyy  , "qlm_dkyyy  ")
+     call poison_check (qlm_dkyzy  , "qlm_dkyzy  ")
+     call poison_check (qlm_dkzzy  , "qlm_dkzzy  ")
+     call poison_check (qlm_dkxxz  , "qlm_dkxxz  ")
+     call poison_check (qlm_dkxyz  , "qlm_dkxyz  ")
+     call poison_check (qlm_dkxzz  , "qlm_dkxzz  ")
+     call poison_check (qlm_dkyyz  , "qlm_dkyyz  ")
+     call poison_check (qlm_dkyzz  , "qlm_dkyzz  ")
+     call poison_check (qlm_dkzzz  , "qlm_dkzzz  ")
+     call poison_check (qlm_alpha  , "qlm_alpha  ")
+     call poison_check (qlm_betax  , "qlm_betax  ")
+     call poison_check (qlm_betay  , "qlm_betay  ")
+     call poison_check (qlm_betaz  , "qlm_betaz  ")
+     call poison_check (qlm_ttt    , "qlm_ttt    ")
+     call poison_check (qlm_ttx    , "qlm_ttx    ")
+     call poison_check (qlm_tty    , "qlm_tty    ")
+     call poison_check (qlm_ttz    , "qlm_ttz    ")
+     call poison_check (qlm_txx    , "qlm_txx    ")
+     call poison_check (qlm_txy    , "qlm_txy    ")
+     call poison_check (qlm_txz    , "qlm_txz    ")
+     call poison_check (qlm_tyy    , "qlm_tyy    ")
+     call poison_check (qlm_tyz    , "qlm_tyz    ")
+     call poison_check (qlm_tzz    , "qlm_tzz    ")
 #endif
      
   end if
@@ -599,4 +686,25 @@ contains
     a = b
   end subroutine copy
   
+  subroutine poison (arr)
+    CCTK_REAL, intent(out) :: arr(:,:)
+    arr = poison_value
+  end subroutine poison
+  
+  subroutine poison_check (arr, name)
+    CCTK_REAL,    intent(in) :: arr(:,:)
+    character(*), intent(in) :: name
+    character*1000 :: msg
+!!$    integer        :: i, j
+    if (any(arr==poison_value)) then
+       write (msg, '("Poison found in ",a)') trim(name)
+       call CCTK_WARN (CCTK_WARN_ALERT, msg)
+!!$       do j=1,size(arr,2)
+!!$          do i=1,size(arr,1)
+!!$             print '(2i6)', i,j
+!!$          end do
+!!$       end do
+    end if
+  end subroutine poison_check
+  
 end subroutine qlm_interpolate
diff --git a/src/qlm_multipoles.F90 b/src/qlm_multipoles.F90
index 79b75d5..3183d90 100644
--- a/src/qlm_multipoles.F90
+++ b/src/qlm_multipoles.F90
@@ -25,9 +25,9 @@ subroutine qlm_multipoles (CCTK_ARGUMENTS, hn)
   CCTK_REAL, parameter :: o128=one/128
   
   CCTK_REAL    :: qq(2,2), dtq, rsc
-  CCTK_REAL    :: ee(3,2)
-  CCTK_REAL    :: kk(3,3)
-  CCTK_REAL    :: ll(0:3), nn(0:3), ss(0:3)
+!!$  CCTK_REAL    :: ee(3,2)
+!!$  CCTK_REAL    :: kk(3,3)
+!!$  CCTK_REAL    :: ll(0:3), nn(0:3), ss(0:3)
   CCTK_COMPLEX :: psi2
   CCTK_REAL    :: zz, dzz(2)
   CCTK_REAL    :: area, mass, spin
@@ -74,31 +74,31 @@ subroutine qlm_multipoles (CCTK_ARGUMENTS, hn)
         
         call calc_2det (qq, dtq)
         
-        ee(1,1:2) = deriv (qlm_x(:,:,hn), i, j, delta_space)
-        ee(2,1:2) = deriv (qlm_y(:,:,hn), i, j, delta_space)
-        ee(3,1:2) = deriv (qlm_z(:,:,hn), i, j, delta_space)
-        
-        kk(1,1) = qlm_kxx(i,j)
-        kk(1,2) = qlm_kxy(i,j)
-        kk(1,3) = qlm_kxz(i,j)
-        kk(2,2) = qlm_kyy(i,j)
-        kk(2,3) = qlm_kyz(i,j)
-        kk(3,3) = qlm_kzz(i,j)
-        kk(2,1) = kk(1,2)
-        kk(3,1) = kk(1,3)
-        kk(3,2) = kk(2,3)
-        
-        ll(0) = qlm_l0(i,j,hn)
-        ll(1) = qlm_l1(i,j,hn)
-        ll(2) = qlm_l2(i,j,hn)
-        ll(3) = qlm_l3(i,j,hn)
-        
-        nn(0) = qlm_n0(i,j,hn)
-        nn(1) = qlm_n1(i,j,hn)
-        nn(2) = qlm_n2(i,j,hn)
-        nn(3) = qlm_n3(i,j,hn)
-        
-        ss = (ll - nn) / sqrt(two)
+!!$        ee(1,1:2) = deriv (qlm_x(:,:,hn), i, j, delta_space)
+!!$        ee(2,1:2) = deriv (qlm_y(:,:,hn), i, j, delta_space)
+!!$        ee(3,1:2) = deriv (qlm_z(:,:,hn), i, j, delta_space)
+!!$        
+!!$        kk(1,1) = qlm_kxx(i,j)
+!!$        kk(1,2) = qlm_kxy(i,j)
+!!$        kk(1,3) = qlm_kxz(i,j)
+!!$        kk(2,2) = qlm_kyy(i,j)
+!!$        kk(2,3) = qlm_kyz(i,j)
+!!$        kk(3,3) = qlm_kzz(i,j)
+!!$        kk(2,1) = kk(1,2)
+!!$        kk(3,1) = kk(1,3)
+!!$        kk(3,2) = kk(2,3)
+!!$        
+!!$        ll(0) = qlm_l0(i,j,hn)
+!!$        ll(1) = qlm_l1(i,j,hn)
+!!$        ll(2) = qlm_l2(i,j,hn)
+!!$        ll(3) = qlm_l3(i,j,hn)
+!!$        
+!!$        nn(0) = qlm_n0(i,j,hn)
+!!$        nn(1) = qlm_n1(i,j,hn)
+!!$        nn(2) = qlm_n2(i,j,hn)
+!!$        nn(3) = qlm_n3(i,j,hn)
+!!$        
+!!$        ss = (ll - nn) / sqrt(two)
         
         rsc = qlm_rsc(i,j,hn)
         psi2 = qlm_psi2(i,j,hn)
@@ -157,50 +157,6 @@ subroutine qlm_multipoles (CCTK_ARGUMENTS, hn)
      end do
   end do
   
-  ! Normalise
-  
-!!$  ! This is the normalisation for I_n and L_n
-!!$  qlm_mp_m0(hn) = qlm_mp_m0(hn) / sqrt(4*pi/ 1)
-!!$  qlm_mp_m1(hn) = qlm_mp_m1(hn) / sqrt(4*pi/ 3)
-!!$  qlm_mp_m2(hn) = qlm_mp_m2(hn) / sqrt(4*pi/ 5)
-!!$  qlm_mp_m3(hn) = qlm_mp_m3(hn) / sqrt(4*pi/ 7)
-!!$  qlm_mp_m4(hn) = qlm_mp_m4(hn) / sqrt(4*pi/ 9)
-!!$  qlm_mp_m5(hn) = qlm_mp_m5(hn) / sqrt(4*pi/11)
-!!$  qlm_mp_m6(hn) = qlm_mp_m6(hn) / sqrt(4*pi/13)
-!!$  qlm_mp_m7(hn) = qlm_mp_m7(hn) / sqrt(4*pi/15)
-!!$  qlm_mp_m8(hn) = qlm_mp_m8(hn) / sqrt(4*pi/17)
-!!$  
-!!$  qlm_mp_j0(hn) = qlm_mp_j0(hn) / sqrt(4*pi/ 1)
-!!$  qlm_mp_j1(hn) = qlm_mp_j1(hn) / sqrt(4*pi/ 3)
-!!$  qlm_mp_j2(hn) = qlm_mp_j2(hn) / sqrt(4*pi/ 5)
-!!$  qlm_mp_j3(hn) = qlm_mp_j3(hn) / sqrt(4*pi/ 7)
-!!$  qlm_mp_j4(hn) = qlm_mp_j4(hn) / sqrt(4*pi/ 9)
-!!$  qlm_mp_j5(hn) = qlm_mp_j5(hn) / sqrt(4*pi/11)
-!!$  qlm_mp_j6(hn) = qlm_mp_j6(hn) / sqrt(4*pi/13)
-!!$  qlm_mp_j7(hn) = qlm_mp_j7(hn) / sqrt(4*pi/15)
-!!$  qlm_mp_j8(hn) = qlm_mp_j8(hn) / sqrt(4*pi/17)
-  
-  ! This is the normalisation for M_n and J_n
-  qlm_mp_m0(hn) = qlm_mp_m0(hn) * qlm_mass(hn) * qlm_radius(hn)**0 / (2*pi)
-  qlm_mp_m1(hn) = qlm_mp_m1(hn) * qlm_mass(hn) * qlm_radius(hn)**1 / (2*pi)
-  qlm_mp_m2(hn) = qlm_mp_m2(hn) * qlm_mass(hn) * qlm_radius(hn)**2 / (2*pi)
-  qlm_mp_m3(hn) = qlm_mp_m3(hn) * qlm_mass(hn) * qlm_radius(hn)**3 / (2*pi)
-  qlm_mp_m4(hn) = qlm_mp_m4(hn) * qlm_mass(hn) * qlm_radius(hn)**4 / (2*pi)
-  qlm_mp_m5(hn) = qlm_mp_m5(hn) * qlm_mass(hn) * qlm_radius(hn)**5 / (2*pi)
-  qlm_mp_m6(hn) = qlm_mp_m6(hn) * qlm_mass(hn) * qlm_radius(hn)**6 / (2*pi)
-  qlm_mp_m7(hn) = qlm_mp_m7(hn) * qlm_mass(hn) * qlm_radius(hn)**7 / (2*pi)
-  qlm_mp_m8(hn) = qlm_mp_m8(hn) * qlm_mass(hn) * qlm_radius(hn)**8 / (2*pi)
-  
-  qlm_mp_j0(hn) = qlm_mp_j0(hn) * qlm_radius(hn)**1 / (4*pi)
-  qlm_mp_j1(hn) = qlm_mp_j1(hn) * qlm_radius(hn)**2 / (4*pi)
-  qlm_mp_j2(hn) = qlm_mp_j2(hn) * qlm_radius(hn)**3 / (4*pi)
-  qlm_mp_j3(hn) = qlm_mp_j3(hn) * qlm_radius(hn)**4 / (4*pi)
-  qlm_mp_j4(hn) = qlm_mp_j4(hn) * qlm_radius(hn)**5 / (4*pi)
-  qlm_mp_j5(hn) = qlm_mp_j5(hn) * qlm_radius(hn)**6 / (4*pi)
-  qlm_mp_j6(hn) = qlm_mp_j6(hn) * qlm_radius(hn)**7 / (4*pi)
-  qlm_mp_j7(hn) = qlm_mp_j7(hn) * qlm_radius(hn)**8 / (4*pi)
-  qlm_mp_j8(hn) = qlm_mp_j8(hn) * qlm_radius(hn)**9 / (4*pi)
-  
 contains
   
   ! Legendre polynomials
@@ -321,3 +277,64 @@ contains
   end function dp8
   
 end subroutine qlm_multipoles
+
+
+
+subroutine qlm_multipoles_normalise (CCTK_ARGUMENTS, hn)
+  use cctk
+  use constants
+  implicit none
+  DECLARE_CCTK_ARGUMENTS
+  DECLARE_CCTK_FUNCTIONS
+  DECLARE_CCTK_PARAMETERS
+  integer :: hn
+  
+  if (veryverbose/=0) then
+     call CCTK_INFO ("Normalising multipole moments")
+  end if
+  
+  ! Normalise
+  
+!!$  ! This is the normalisation for I_n and L_n
+!!$  qlm_mp_m0(hn) = qlm_mp_m0(hn) / sqrt(4*pi/ 1)
+!!$  qlm_mp_m1(hn) = qlm_mp_m1(hn) / sqrt(4*pi/ 3)
+!!$  qlm_mp_m2(hn) = qlm_mp_m2(hn) / sqrt(4*pi/ 5)
+!!$  qlm_mp_m3(hn) = qlm_mp_m3(hn) / sqrt(4*pi/ 7)
+!!$  qlm_mp_m4(hn) = qlm_mp_m4(hn) / sqrt(4*pi/ 9)
+!!$  qlm_mp_m5(hn) = qlm_mp_m5(hn) / sqrt(4*pi/11)
+!!$  qlm_mp_m6(hn) = qlm_mp_m6(hn) / sqrt(4*pi/13)
+!!$  qlm_mp_m7(hn) = qlm_mp_m7(hn) / sqrt(4*pi/15)
+!!$  qlm_mp_m8(hn) = qlm_mp_m8(hn) / sqrt(4*pi/17)
+!!$  
+!!$  qlm_mp_j0(hn) = qlm_mp_j0(hn) / sqrt(4*pi/ 1)
+!!$  qlm_mp_j1(hn) = qlm_mp_j1(hn) / sqrt(4*pi/ 3)
+!!$  qlm_mp_j2(hn) = qlm_mp_j2(hn) / sqrt(4*pi/ 5)
+!!$  qlm_mp_j3(hn) = qlm_mp_j3(hn) / sqrt(4*pi/ 7)
+!!$  qlm_mp_j4(hn) = qlm_mp_j4(hn) / sqrt(4*pi/ 9)
+!!$  qlm_mp_j5(hn) = qlm_mp_j5(hn) / sqrt(4*pi/11)
+!!$  qlm_mp_j6(hn) = qlm_mp_j6(hn) / sqrt(4*pi/13)
+!!$  qlm_mp_j7(hn) = qlm_mp_j7(hn) / sqrt(4*pi/15)
+!!$  qlm_mp_j8(hn) = qlm_mp_j8(hn) / sqrt(4*pi/17)
+  
+  ! This is the normalisation for M_n and J_n
+  qlm_mp_m0(hn) = qlm_mp_m0(hn) * qlm_mass(hn) * qlm_radius(hn)**0 / (2*pi)
+  qlm_mp_m1(hn) = qlm_mp_m1(hn) * qlm_mass(hn) * qlm_radius(hn)**1 / (2*pi)
+  qlm_mp_m2(hn) = qlm_mp_m2(hn) * qlm_mass(hn) * qlm_radius(hn)**2 / (2*pi)
+  qlm_mp_m3(hn) = qlm_mp_m3(hn) * qlm_mass(hn) * qlm_radius(hn)**3 / (2*pi)
+  qlm_mp_m4(hn) = qlm_mp_m4(hn) * qlm_mass(hn) * qlm_radius(hn)**4 / (2*pi)
+  qlm_mp_m5(hn) = qlm_mp_m5(hn) * qlm_mass(hn) * qlm_radius(hn)**5 / (2*pi)
+  qlm_mp_m6(hn) = qlm_mp_m6(hn) * qlm_mass(hn) * qlm_radius(hn)**6 / (2*pi)
+  qlm_mp_m7(hn) = qlm_mp_m7(hn) * qlm_mass(hn) * qlm_radius(hn)**7 / (2*pi)
+  qlm_mp_m8(hn) = qlm_mp_m8(hn) * qlm_mass(hn) * qlm_radius(hn)**8 / (2*pi)
+  
+  qlm_mp_j0(hn) = qlm_mp_j0(hn) * qlm_radius(hn)**1 / (4*pi)
+  qlm_mp_j1(hn) = qlm_mp_j1(hn) * qlm_radius(hn)**2 / (4*pi)
+  qlm_mp_j2(hn) = qlm_mp_j2(hn) * qlm_radius(hn)**3 / (4*pi)
+  qlm_mp_j3(hn) = qlm_mp_j3(hn) * qlm_radius(hn)**4 / (4*pi)
+  qlm_mp_j4(hn) = qlm_mp_j4(hn) * qlm_radius(hn)**5 / (4*pi)
+  qlm_mp_j5(hn) = qlm_mp_j5(hn) * qlm_radius(hn)**6 / (4*pi)
+  qlm_mp_j6(hn) = qlm_mp_j6(hn) * qlm_radius(hn)**7 / (4*pi)
+  qlm_mp_j7(hn) = qlm_mp_j7(hn) * qlm_radius(hn)**8 / (4*pi)
+  qlm_mp_j8(hn) = qlm_mp_j8(hn) * qlm_radius(hn)**9 / (4*pi)
+  
+end subroutine qlm_multipoles_normalise
diff --git a/src/qlm_tetrad.F90 b/src/qlm_tetrad.F90
index 85837dc..d309906 100644
--- a/src/qlm_tetrad.F90
+++ b/src/qlm_tetrad.F90
@@ -46,8 +46,8 @@ subroutine qlm_calc_tetrad (CCTK_ARGUMENTS, hn)
   CCTK_REAL    :: nabla_ll(0:3,0:3), nabla_nn(0:3,0:3)
   CCTK_COMPLEX :: nabla_mm(0:3,0:3)
   
-  CCTK_REAL    :: t0, t1, t2
-  logical      :: ce0, ce1, ce2
+  !CCTK_REAL    :: t0, t1, t2
+  !logical      :: ce0, ce1, ce2
   CCTK_REAL    :: delta_space(2)
   
   CCTK_REAL    :: count, accuracy
@@ -153,16 +153,17 @@ subroutine qlm_calc_tetrad (CCTK_ARGUMENTS, hn)
         ee(1,2) = qlm_y(i,j,hn) - qlm_origin_y(hn)
         ee(1,3) = qlm_z(i,j,hn) - qlm_origin_z(hn)
         
-        ee_p(1,1) = qlm_x_p(i,j,hn) - qlm_origin_x_p(hn)
-        ee_p(1,2) = qlm_y_p(i,j,hn) - qlm_origin_y_p(hn)
-        ee_p(1,3) = qlm_z_p(i,j,hn) - qlm_origin_z_p(hn)
+        !ee_p(1,1) = qlm_x_p(i,j,hn) - qlm_origin_x_p(hn)
+        !ee_p(1,2) = qlm_y_p(i,j,hn) - qlm_origin_y_p(hn)
+        !ee_p(1,3) = qlm_z_p(i,j,hn) - qlm_origin_z_p(hn)
         
-        ee_p_p(1,1) = qlm_x_p_p(i,j,hn) - qlm_origin_x_p_p(hn)
-        ee_p_p(1,2) = qlm_y_p_p(i,j,hn) - qlm_origin_y_p_p(hn)
-        ee_p_p(1,3) = qlm_z_p_p(i,j,hn) - qlm_origin_z_p_p(hn)
+        !ee_p_p(1,1) = qlm_x_p_p(i,j,hn) - qlm_origin_x_p_p(hn)
+        !ee_p_p(1,2) = qlm_y_p_p(i,j,hn) - qlm_origin_y_p_p(hn)
+        !ee_p_p(1,3) = qlm_z_p_p(i,j,hn) - qlm_origin_z_p_p(hn)
         
         dee(1,0,:) = 0
-        dee(1,1:3,0) = timederiv (ee(1,1:3), ee_p(1,1:3), ee_p_p(1,1:3), t0,t1,t2, ce0,ce1,ce2)
+        !dee(1,1:3,0) = timederiv (ee(1,1:3), ee_p(1,1:3), ee_p_p(1,1:3), t0,t1,t2, ce0,ce1,ce2)
+        dee(1,1:3,0) = 0
         dee_spher(1,:,1) = 0    ! this is a choice
         dee_spher(1,1,2:3) = deriv (qlm_x(:,:,hn), i,j, delta_space)
         dee_spher(1,2,2:3) = deriv (qlm_y(:,:,hn), i,j, delta_space)
@@ -184,7 +185,8 @@ subroutine qlm_calc_tetrad (CCTK_ARGUMENTS, hn)
         ee_p_p(2:3,3) = deriv (qlm_z_p_p(:,:,hn), i,j, delta_space)
         
         dee(2:3,0,:) = 0
-        dee(2:3,1:3,0) = timederiv (ee(2:3,1:3), ee_p(2:3,1:3), ee_p_p(2:3,1:3), t0,t1,t2, ce0,ce1,ce2)
+        !dee(2:3,1:3,0) = timederiv (ee(2:3,1:3), ee_p(2:3,1:3), ee_p_p(2:3,1:3), t0,t1,t2, ce0,ce1,ce2)
+        dee(2:3,1:3,0) = 0
         dee_spher(2:3,:,1) = 0  ! this is a choice
         dee_spher(2:3,1,2:3) = deriv2 (qlm_x(:,:,hn), i,j, delta_space)
         dee_spher(2:3,2,2:3) = deriv2 (qlm_y(:,:,hn), i,j, delta_space)
diff --git a/src/qlm_variables.F90 b/src/qlm_variables.F90
index d21842c..9343666 100644
--- a/src/qlm_variables.F90
+++ b/src/qlm_variables.F90
@@ -21,7 +21,10 @@ module qlm_variables
              qlm_dkxxy, qlm_dkxyy, qlm_dkxzy, qlm_dkyyy, qlm_dkyzy, qlm_dkzzy, &
              qlm_dkxxz, qlm_dkxyz, qlm_dkxzz, qlm_dkyyz, qlm_dkyzz, qlm_dkzzz, &
              qlm_alpha, &
-             qlm_betax, qlm_betay, qlm_betaz
+             qlm_betax, qlm_betay, qlm_betaz, &
+             qlm_ttt, &
+             qlm_ttx, qlm_tty, qlm_ttz, &
+             qlm_txx, qlm_txy, qlm_txz, qlm_tyy, qlm_tyz, qlm_tzz
   
   type tetrad_derivs
      ! nabla_ll(a,b) = D_b l_a
@@ -34,8 +37,11 @@ module qlm_variables
 contains
   
   subroutine allocate_variables(ntheta, nphi)
-    integer, intent(in) :: ntheta, nphi
-
+    integer,   intent(in) :: ntheta, nphi
+    CCTK_REAL, parameter  :: zero = 0
+    integer,   parameter  :: rk = kind(zero)
+    type(tetrad_derivs)   :: tetrad_derivs_nan
+    
     allocate(qlm_gxx(ntheta,nphi))
     allocate(qlm_gxy(ntheta,nphi))
     allocate(qlm_gxz(ntheta,nphi))
@@ -124,6 +130,16 @@ contains
     allocate(qlm_betax(ntheta,nphi))
     allocate(qlm_betay(ntheta,nphi))
     allocate(qlm_betaz(ntheta,nphi))
+    allocate(qlm_ttt(ntheta,nphi))
+    allocate(qlm_ttx(ntheta,nphi))
+    allocate(qlm_tty(ntheta,nphi))
+    allocate(qlm_ttz(ntheta,nphi))
+    allocate(qlm_txx(ntheta,nphi))
+    allocate(qlm_txy(ntheta,nphi))
+    allocate(qlm_txz(ntheta,nphi))
+    allocate(qlm_tyy(ntheta,nphi))
+    allocate(qlm_tyz(ntheta,nphi))
+    allocate(qlm_tzz(ntheta,nphi))
     
     allocate(qlm_tetrad_derivs(ntheta,nphi))
     
@@ -215,6 +231,21 @@ contains
     qlm_betax = TAT_nan()
     qlm_betay = TAT_nan()
     qlm_betaz = TAT_nan()
+    qlm_ttt = TAT_nan()
+    qlm_ttx = TAT_nan()
+    qlm_tty = TAT_nan()
+    qlm_ttz = TAT_nan()
+    qlm_txx = TAT_nan()
+    qlm_txy = TAT_nan()
+    qlm_txz = TAT_nan()
+    qlm_tyy = TAT_nan()
+    qlm_tyz = TAT_nan()
+    qlm_tzz = TAT_nan()
+    
+    tetrad_derivs_nan%nabla_ll = TAT_nan()
+    tetrad_derivs_nan%nabla_nn = TAT_nan()
+    tetrad_derivs_nan%nabla_mm = cmplx(TAT_nan(),TAT_nan(),rk)
+    qlm_tetrad_derivs = tetrad_derivs_nan
   end subroutine allocate_variables
 
   subroutine deallocate_variables
@@ -306,6 +337,16 @@ contains
     deallocate(qlm_betax)
     deallocate(qlm_betay)
     deallocate(qlm_betaz)
+    deallocate(qlm_ttt)
+    deallocate(qlm_ttx)
+    deallocate(qlm_tty)
+    deallocate(qlm_ttz)
+    deallocate(qlm_txx)
+    deallocate(qlm_txy)
+    deallocate(qlm_txz)
+    deallocate(qlm_tyy)
+    deallocate(qlm_tyz)
+    deallocate(qlm_tzz)
     
     deallocate(qlm_tetrad_derivs)
   end subroutine deallocate_variables
diff --git a/src/qlm_weyl_scalars.F90 b/src/qlm_weyl_scalars.F90
index d79751c..ade2e31 100644
--- a/src/qlm_weyl_scalars.F90
+++ b/src/qlm_weyl_scalars.F90
@@ -24,6 +24,9 @@ subroutine qlm_calc_weyl_scalars (CCTK_ARGUMENTS, hn)
   
   CCTK_REAL, parameter :: two=2, four=4
   CCTK_REAL    :: gg(3,3), dgg(3,3,3), ddgg(3,3,3,3), gg_dot(3,3), gg_dot2(3,3), dgg_dot(3,3,3)
+  CCTK_REAL    :: kk(3,3), dkk(3,3,3), kk_dot(3,3)
+  CCTK_REAL    :: tt(3,3)
+  CCTK_REAL    :: dtg, gu(3,3), dgu(3,3,3), gamma(3,3,3), dgamma(3,3,3,3), ri(3,3)
   CCTK_REAL    :: g4(0:3,0:3), dg4(0:3,0:3,0:3), ddg4(0:3,0:3,0:3,0:3)
   CCTK_REAL    :: gu4(0:3,0:3), dgu4(0:3,0:3,0:3)
   CCTK_REAL    :: gamma4(0:3,0:3,0:3), dgamma4(0:3,0:3,0:3,0:3)
@@ -122,6 +125,48 @@ subroutine qlm_calc_weyl_scalars (CCTK_ARGUMENTS, hn)
         ddgg(:,:,3,1) = ddgg(:,:,1,3)
         ddgg(:,:,3,2) = ddgg(:,:,2,3)
         
+        kk(1,1) = qlm_kxx(i,j)
+        kk(1,2) = qlm_kxy(i,j)
+        kk(1,3) = qlm_kxz(i,j)
+        kk(2,2) = qlm_kyy(i,j)
+        kk(2,3) = qlm_kyz(i,j)
+        kk(3,3) = qlm_kzz(i,j)
+        kk(2,1) = kk(1,2)
+        kk(3,1) = kk(1,3)
+        kk(3,2) = kk(2,3)
+        
+        dkk(1,1,1) = qlm_dkxxx(i,j)
+        dkk(1,2,1) = qlm_dkxyx(i,j)
+        dkk(1,3,1) = qlm_dkxzx(i,j)
+        dkk(2,2,1) = qlm_dkyyx(i,j)
+        dkk(2,3,1) = qlm_dkyzx(i,j)
+        dkk(3,3,1) = qlm_dkzzx(i,j)
+        dkk(1,1,2) = qlm_dkxxy(i,j)
+        dkk(1,2,2) = qlm_dkxyy(i,j)
+        dkk(1,3,2) = qlm_dkxzy(i,j)
+        dkk(2,2,2) = qlm_dkyyy(i,j)
+        dkk(2,3,2) = qlm_dkyzy(i,j)
+        dkk(3,3,2) = qlm_dkzzy(i,j)
+        dkk(1,1,3) = qlm_dkxxz(i,j)
+        dkk(1,2,3) = qlm_dkxyz(i,j)
+        dkk(1,3,3) = qlm_dkxzz(i,j)
+        dkk(2,2,3) = qlm_dkyyz(i,j)
+        dkk(2,3,3) = qlm_dkyzz(i,j)
+        dkk(3,3,3) = qlm_dkzzz(i,j)
+        dkk(2,1,:) = dkk(1,2,:)
+        dkk(3,1,:) = dkk(1,3,:)
+        dkk(3,2,:) = dkk(2,3,:)
+        
+        tt(1,1) = qlm_txx(i,j)
+        tt(1,2) = qlm_txy(i,j)
+        tt(1,3) = qlm_txz(i,j)
+        tt(2,2) = qlm_tyy(i,j)
+        tt(2,3) = qlm_tyz(i,j)
+        tt(3,3) = qlm_tzz(i,j)
+        tt(2,1) = tt(1,2)
+        tt(3,1) = tt(1,3)
+        tt(3,2) = tt(2,3)
+        
         ll(0) = qlm_l0(i,j,hn)
         ll(1) = qlm_l1(i,j,hn)
         ll(2) = qlm_l2(i,j,hn)
@@ -140,6 +185,18 @@ subroutine qlm_calc_weyl_scalars (CCTK_ARGUMENTS, hn)
         
         
         ! Calculate 4-metric
+        call calc_det (gg, dtg)
+        call calc_inv (gg, dtg, gu)
+        call calc_invderiv (gu, dgg, dgu)
+        call calc_connections (gu, dgg, gamma)
+        call calc_connectionderivs (gu, dgg, dgu, ddgg, dgamma)
+        call calc_ricci (gamma, dgamma, ri)
+        
+        call calc_3metricdot_simple (kk, gg_dot)
+        call calc_3metricderivdot_simple (dkk, dgg_dot)
+        call calc_extcurvdot_simple (gg,gu,ri, kk, tt, kk_dot)
+        call calc_3metricdot2_simple (kk_dot, gg_dot2)
+        
         call calc_4metricderivs2_simple (gg, dgg, &
              ddgg, gg_dot, gg_dot2, dgg_dot, g4,dg4,ddg4)
         call calc_4inv (g4, gu4)