diff options
Diffstat (limited to 'src/ADMmass_integrand3D.F')
| -rw-r--r-- | src/ADMmass_integrand3D.F | 57 |
1 files changed, 44 insertions, 13 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 |