aboutsummaryrefslogtreecommitdiff
path: root/src/D2_extract.F
diff options
context:
space:
mode:
Diffstat (limited to 'src/D2_extract.F')
-rw-r--r--src/D2_extract.F60
1 files changed, 51 insertions, 9 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
generated by cgit v1.2.3 (git 2.39.1) at 2024-05-19 13:06:11 +0000