


















































































































 


































































































































































 






























































































































































































 



















































 










 








































































































































































































































 










                                
                                

                                
                                

                                
                                

                                





































c

      SUBROUTINE D3_extract(cctkGH,conformal_state,do_ADMmass,do_momentu
     &m,do_spin,igrid,
     &     origin,myproc,interpolation_operator,interpolation_order,Nt,N
     &p,
     &     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)

c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c
c

      USE D3_to_D2_int
      USE cartesian_to_spherical_int
      USE unphysical_to_physical_int
      USE D2_extract_int

c

      IMPLICIT NONE

c

      integer*8 :: cctkGH

      INTEGER*4,INTENT(IN) :: 
     &     conformal_state
      INTEGER,INTENT(IN) :: 
     &     igrid,l,m,Psi_power,myproc
      INTEGER*4,INTENT(IN) :: 
     &     Nt,Np,nx,ny,nz,all_modes,do_momentum,do_spin,interpolation_or
     &der
      INTEGER,INTENT(IN) ::
     &     do_ADMmass(2)
      REAL*8,INTENT(IN) :: 
     &     origin(3),Dx,Dy,Dz,eta
      REAL*8,INTENT(IN),DIMENSION(nx,ny,nz) :: 
     &     Psi,g00,gxx,gxy,gxz,gyy,gyz,gzz,
     &     hxx,hxy,hxz,hyy,hyz,hzz
      REAL*8,INTENT(IN),DIMENSION(nx) :: x
      REAL*8,INTENT(IN),DIMENSION(ny) :: y
      REAL*8,INTENT(IN),DIMENSION(nz) :: z
      REAL*8,INTENT(INOUT),DIMENSION(nx,ny,nz) ::
     &     Extract_temp3d
      integer*8,INTENT(IN) ::
     &     interpolation_operator

c
      REAL*8,INTENT(OUT) :: 
     &     ADMmass(2),mass,rsch,Qodd(:,2:,0:),Qeven(:,2:,0:),dtaudt,
     &     momentum(3),spin(3)

c
      REAL*8 :: 
     &    theta(Nt),phi(Np)
      REAL*8,DIMENSION(Nt,Np) :: 
     &    Psis,g00s,gxxs,gxys,gxzs,gyys,gyzs,gzzs,dPsis,dgxxs,dgxys,dgxz
     &s,
     &    dgyys,dgyzs,dgzzs,grr,grt,grp,gtt,gtp,gpp,dgtt,dgtp,dgpp,
     &    ADMmass_int1,ADMmass_int2,
     &    momentum_int1,momentum_int2,momentum_int3,
     &    spin_int1,spin_int2,spin_int3

c
      INTEGER :: 
     &    i,j
      REAL*8,PARAMETER ::
     &    half = 0.5D0,
     &    one  = 1.0D0,
     &    two  = 2.0D0
      REAL*8 :: 
     &    Pi,Dt,Dp

c


c
c
c
c
c

c
c
c
c
c
      Pi = two*ASIN(one)

c

      SELECT CASE (igrid)
      CASE (0)                  ! full grid
        Dt = Pi/DBLE(Nt)
        Dp = two*Pi/DBLE(Np)
      CASE (1)                  ! octant grid
        Dt = half*Pi/DBLE(Nt)
        Dp = half*Pi/DBLE(Np)
      CASE (2)                  ! cartoon
        Dt = Pi/DBLE(Nt)
        Dp = 2D0*Pi
      CASE (3)                  ! bitant
        Dt = half*Pi/DBLE(Nt)
        Dp = two*Pi/DBLE(Np)
      CASE DEFAULT
        WRITE (*,*) "Grid incorrectly set in D3_extract"
        STOP
      END SELECT
         

c

      DO i = 1, Nt
         theta(i) = (DBLE(i) - half)* Dt
      ENDDO
      DO j = 1, Np
         phi(j) = (DBLE(j) - half)* Dp
      ENDDO

      IF (igrid == 2) phi(1) = 0D0
      

c
c
c
c
c

c
c
      CALL D3_to_D2(cctkGH,conformal_state,do_ADMmass,do_momentum,do_spi
     &n,
     &     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,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


      if (myproc .eq. 0) then 

c
c
c
c
c

         CALL cartesian_to_spherical(theta,phi,eta,gxxs,gxys,gxzs,gyys,
     &        gyzs,gzzs,grr,grt,grp,gtt,gtp,gpp,dgxxs,dgxys,dgxzs,dgyys,
     &
     &        dgyzs,dgzzs,dgtt,dgtp,dgpp)


c
c
c
c
c

         if (conformal_state > 0) then
           CALL unphysical_to_physical(grr,grt,grp,gtt,gtp,gpp,dgtt,dgtp
     &,
     &        dgpp,Psis,dPsis,Psi_power)
         end if

c
c
c
c
c

c
         CALL D2_extract(eta,igrid,Nt,Np,theta,phi,all_modes,l,m,g00s,gr
     &r,
     &        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
c
c
      

c
      endif

      END SUBROUTINE D3_extract



