Fix extract for physical metrics

Closes Cactus/1196


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

3 files changed, 60 insertions, 26 deletions

diff --git a/src/ADMmass_integrand3D.F b/src/ADMmass_integrand3D.F
index 06b1198..a501ef1 100644
--- a/src/ADMmass_integrand3D.F
+++ b/src/ADMmass_integrand3D.F
@@ -1,9 +1,10 @@
 #include "cctk.h"
+#include "cctk_Parameters.h"
 
 c     ========================================================================
 
       SUBROUTINE ADMmass_integrand3D(origin,Dx,Dy,Dz,x,y,z,gxx,gxy,gxz,
-     &     gyy,gyz,gzz,ADMmass_int,Psi,Psi_power)
+     &     gyy,gyz,gzz,ADMmass_int,Psi,Psi_power,conformal_state)
 
 c     ------------------------------------------------------------------------
 c
@@ -16,7 +17,7 @@ c     ------------------------------------------------------------------------
 
 c     Input variables
       INTEGER,INTENT(IN) ::
-     &     Psi_power
+     &     Psi_power,conformal_state
       CCTK_REAL,INTENT(IN) :: 
      &     Dx,Dy,Dz,origin(3)
       CCTK_REAL,DIMENSION(:),INTENT(IN) :: 
@@ -39,6 +40,9 @@ c     Local variables, here only
      &     dyy_x,dxz_x,dxz_y,dxz_z,dyz_x,dzz_x,dyy_z,dyz_y,dyz_z,dzz_y,
      &     Pi,idet,p,pip,pim,pjp,pjm,pkp,pkm
 
+      DECLARE_CCTK_PARAMETERS
+      DECLARE_CCTK_FUNCTIONS
+
 c     ------------------------------------------------------------------------
 
       Pi = ACOS(-1D0)
@@ -74,11 +78,24 @@ c     Because other codes evolve Psi**4
 
 c                 Abbreviations for metric coefficients
 c                 -------------------------------------
-                  p = psi(i,j,k)**ip
 
-                  dxx = p*gxx(i,j,k); dxy = p*gxy(i,j,k)
-                  dxz = p*gxz(i,j,k); dyy = p*gyy(i,j,k)
-                  dyz = p*gyz(i,j,k); dzz = p*gzz(i,j,k)
+                  if (conformal_state > 0) then
+  
+                    p = psi(i,j,k)**ip
+
+                    dxx = p*gxx(i,j,k); dxy = p*gxy(i,j,k)
+                    dxz = p*gxz(i,j,k); dyy = p*gyy(i,j,k)
+                    dyz = p*gyz(i,j,k); dzz = p*gzz(i,j,k)
+
+                 else 
+
+                    p = 1
+
+                    dxx = gxx(i,j,k); dxy = gxy(i,j,k)
+                    dxz = gxz(i,j,k); dyy = gyy(i,j,k)
+                    dyz = gyz(i,j,k); dzz = gzz(i,j,k)
+
+                 end if
 
 c                 Determinant of 3-metric
 c                 -----------------------
@@ -98,12 +115,26 @@ c                 ----------------
 
 c                 Derivatives of the 3-metric
 c                 ---------------------------
-                  pip = psi(i+1,j,k)**ip
-                  pim = psi(i-1,j,k)**ip
-                  pjp = psi(i,j+1,k)**ip
-                  pjm = psi(i,j-1,k)**ip
-                  pkp = psi(i,j,k+1)**ip
-                  pkm = psi(i,j,k-1)**ip
+
+                  IF (conformal_state > 0) THEN
+
+                    pip = psi(i+1,j,k)**ip
+                    pim = psi(i-1,j,k)**ip
+                    pjp = psi(i,j+1,k)**ip
+                    pjm = psi(i,j-1,k)**ip
+                    pkp = psi(i,j,k+1)**ip
+                    pkm = psi(i,j,k-1)**ip
+
+                 else
+
+                    pip = 1.0d0
+                    pim = 1.0d0
+                    pjp = 1.0d0
+                    pjm = 1.0d0
+                    pkp = 1.0d0
+                    pkm = 1.0d0
+
+                 end if
 
                   dxx_y = half/Dy*(pjp*gxx(i,j+1,k)-pjm*gxx(i,j-1,k))
                   dxx_z = half/Dz*(pkp*gxx(i,j,k+1)-pkm*gxx(i,j,k-1))
@@ -141,7 +172,7 @@ c                 ---------------------------
                   ADMmass_int(i,j,k) = 0.0D0
                   
                ENDIF
-               
+
           ENDDO
         ENDDO
       ENDDO
diff --git a/src/ADMmass_integrand3D_int.F b/src/ADMmass_integrand3D_int.F
index c103f55..b9f98d3 100644
--- a/src/ADMmass_integrand3D_int.F
+++ b/src/ADMmass_integrand3D_int.F
@@ -8,10 +8,10 @@ c     ------------------------------------------------------------------
       INTERFACE
 
           SUBROUTINE  ADMmass_integrand3D(origin,Dx,Dy,Dz,x,y,z,gxx,gxy,
-     &        gxz,gyy,gyz,gzz,ADMmass_int,Psi,Psi_power)
+     &        gxz,gyy,gyz,gzz,ADMmass_int,Psi,Psi_power,conformal_state)
           IMPLICIT NONE
           INTEGER,INTENT(IN) ::
-     &         Psi_power
+     &         Psi_power,conformal_state
           CCTK_REAL,INTENT(IN) :: 
      &        Dx,Dy,Dz,origin(3)
           CCTK_REAL,DIMENSION(:),INTENT(IN) :: 
diff --git a/src/D3_to_D2.F b/src/D3_to_D2.F
index bc684dc..2948517 100644
--- a/src/D3_to_D2.F
+++ b/src/D3_to_D2.F
@@ -177,17 +177,20 @@ c     ------------------------------------------------------------
 
 c     Calculate radial derivatives and project onto sphere
 c     ----------------------------------------------------
-      CALL met_rad_der(origin,Dx,Dy,Dz,x,y,z,Psi,Extract_temp3d)
-      CALL CCTK_SyncGroup(ierror,cctkGH,"extract::temps")   
       call CCTK_VarIndex(in_array_indices(1), "extract::temp3d")
-      call CCTK_InterpGV (ierror, cctkGH, interp_handle, coord_system_handle,
-     $                    npoints, 1, 1,
-     $                    xs, ys, zs,
-     $                    CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
-     $                    CCTK_VARIABLE_REAL,
-     $                    in_array_indices(1),
-     $                    dPsis,
-     $                    CCTK_VARIABLE_REAL)
+
+      if (conformal_state > 0) then
+         CALL met_rad_der(origin,Dx,Dy,Dz,x,y,z,Psi,Extract_temp3d)
+         CALL CCTK_SyncGroup(ierror,cctkGH,"extract::temps")   
+         call CCTK_InterpGV (ierror, cctkGH, interp_handle, coord_system_handle,
+     $        npoints, 1, 1,
+     $        xs, ys, zs,
+     $        CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
+     $        CCTK_VARIABLE_REAL,
+     $        in_array_indices(1),
+     $        dPsis,
+     $        CCTK_VARIABLE_REAL)
+      end if
 
       CALL met_rad_der(origin,Dx,Dy,Dz,x,y,z,gxx,Extract_temp3d)
       CALL CCTK_SyncGroup(ierror,cctkGH,"extract::temps")
@@ -261,7 +264,7 @@ c     -----------------------------------
 c     Standard equation
       IF (do_ADMmass(1)) THEN
          CALL ADMmass_integrand3D(origin,Dx,Dy,Dz,x,y,z,gxx,gxy,
-     &        gxz,gyy,gyz,gzz,Extract_temp3d,Psi,Psi_power)
+     &        gxz,gyy,gyz,gzz,Extract_temp3d,Psi,Psi_power,conformal_state)
          CALL CCTK_SyncGroup(ierror,cctkGH,"extract::temps")
          call CCTK_InterpGV (ierror, cctkGH, interp_handle, coord_system_handle,
      $                       npoints, 1, 1,