#include "cctk.h"

c     ========================================================================

      SUBROUTINE D3_to_D2(cctkGH,do_ADMmass,do_momentum,do_spin,
     &     Psi_power,origin,myproc,Dx,Dy,Dz,Psi,
     &     g00,gxx,gxy,gxz,gyy,gyz,gzz,hxx,hxy,hxz,hyy,hyz,hzz,
     &     x,y,z,eta,Nt,Np,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     ------------------------------------------------------------------------
c
c     Project the 3-metric and its 1st radial derivatives onto a 
c     2-sphere.
c
c     ------------------------------------------------------------------------

      USE met_rad_der_int
      USE ADMmass_integrand3D_int
      USE momentum_integrand3D_int
      USE spin_integrand3D_int

      IMPLICIT NONE

c     Input variables

      CCTK_POINTER :: cctkGH

      INTEGER,INTENT(IN) :: 
     &     myproc,Psi_power
      CCTK_INT, INTENT(IN) ::
     &     Nt,Np,do_momentum,do_spin
      CCTK_REAL,INTENT(IN) :: 
     &     origin(3),Dx,Dy,Dz,eta
      CCTK_REAL,INTENT(IN),DIMENSION(:) :: 
     &     theta,phi,x,y,z
      CCTK_REAL,INTENT(IN),DIMENSION(:,:,:) ::  
     &     Psi,g00,gxx,gxy,gxz,gyy,gyz,gzz,
     &     hxx,hxy,hxz,hyy,hyz,hzz
      CCTK_REAL,INTENT(INOUT),DIMENSION(:,:,:) ::
     &     Extract_temp3d
      LOGICAL ::
     &     do_ADMmass(2)

c     Output variables

      CCTK_REAL,INTENT(OUT),DIMENSION(:,:) ::     
     &     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


c     Local variables, passed on

      LOGICAL :: 
     &    err_flag
      INTEGER :: 
     &    iorder,npoints,Nx,Ny,Nz
      INTEGER,DIMENSION(Nt,Np) :: 
     &    ib,jb,kb
      CCTK_REAL,DIMENSION(Nt,Np) :: 
     &    xs,ys,zs,ux,uy,uz,xb,yb,zb

c     Local variables, here only

      INTEGER :: 
     &    i,j,handle,ierror
      CCTK_REAL ::
     &    xmin,xmax,ymin,ymax,zmin,zmax

c     ------------------------------------------------------------------------
      
c     Initial Stuff
c     -------------
      Nx = SIZE(x)
      Ny = SIZE(y)
      Nz = SIZE(z)
      

c     Compute Cartesian coordinates on the surface
c     --------------------------------------------
      DO j = 1, Np
         DO i = 1, Nt

            xs(i,j) = origin(1)+eta*SIN(theta(i))*COS(phi(j))
            ys(i,j) = origin(2)+eta*SIN(theta(i))*SIN(phi(j))
            zs(i,j) = origin(3)+eta*COS(theta(i))

         ENDDO
      ENDDO


c     Only do interpolation on one processor
c     --------------------------------------
      SELECT CASE (myproc) 

      CASE (0)
         npoints = Np*Nt

      CASE DEFAULT
         npoints = 0

      END SELECT 


c     Choose the interpolator
c     -----------------------    
      call CCTK_InterpHandle (handle, "simple")


c     Find the local origin of the spatial coordinates
c     ------------------------------------------------
      call CCTK_CoordLocalRange(ierror,cctkGH,xmin,xmax,-1,"x","cart3d")
      call CCTK_CoordLocalRange(ierror,cctkGH,ymin,ymax,-1,"y","cart3d")
      call CCTK_CoordLocalRange(ierror,cctkGH,zmin,zmax,-1,"z","cart3d")


c     Project un-physical metric and conformal factor  onto sphere
c     ------------------------------------------------------------
      call CCTK_Interp (ierror, cctkGH, handle, npoints, 3, 8, 8,
     $          nx, ny, nz,
     $          xs, ys, zs,
     $          CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
     $          CCTK_VARIABLE_REAL,
     $          xmin, ymin, zmin,
     $          dx, dy, dz,
     $          Psi,g00,gxx,gxy,gxz,gyy,gyz,gzz,
     $          CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,
     $          CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,
     $          CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,
     $          Psis,g00s,gxxs,gxys,gxzs,gyys,gyzs,gzzs,
     $          CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,
     $          CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,
     $          CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL
     $     )


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(cctkGH,"extract::temps")   
      call CCTK_Interp (ierror, cctkGH, handle, npoints, 3, 1, 1,
     $          nx, ny, nz,
     $          xs, ys, zs,
     $          CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,
     $          xmin, ymin, zmin,
     $          dx, dy, dz,
     $          Extract_temp3d,
     $          CCTK_VARIABLE_REAL,
     $          dPsis,
     $          CCTK_VARIABLE_REAL
     $     )

      CALL met_rad_der(origin,Dx,Dy,Dz,x,y,z,gxx,Extract_temp3d)
      CALL CCTK_SyncGroup(cctkGH,"extract::temps")
      call CCTK_Interp (ierror, cctkGH, handle, npoints, 3, 1, 1,
     $          nx, ny, nz,
     $          xs, ys, zs,
     $          CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,
     $          xmin, ymin, zmin,
     $          dx, dy, dz,
     $          Extract_temp3d,
     $          CCTK_VARIABLE_REAL,
     $          dgxxs,
     $          CCTK_VARIABLE_REAL
     $     )

      CALL met_rad_der(origin,Dx,Dy,Dz,x,y,z,gxy,Extract_temp3d)
      CALL CCTK_SyncGroup(cctkGH,"extract::temps")
      call CCTK_Interp (ierror, cctkGH, handle, npoints, 3, 1, 1,
     $          nx, ny, nz,
     $          xs, ys, zs,
     $          CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,
     $          xmin, ymin, zmin,
     $          dx, dy, dz,
     $          Extract_temp3d,
     $          CCTK_VARIABLE_REAL,
     $          dgxys,
     $          CCTK_VARIABLE_REAL
     $     )

      CALL met_rad_der(origin,Dx,Dy,Dz,x,y,z,gxz,Extract_temp3d)
      CALL CCTK_SyncGroup(cctkGH,"extract::temps")
      call CCTK_Interp (ierror, cctkGH, handle, npoints, 3, 1, 1,
     $          nx, ny, nz,
     $          xs, ys, zs,
     $          CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,
     $          xmin, ymin, zmin,
     $          dx, dy, dz,
     $          Extract_temp3d,
     $          CCTK_VARIABLE_REAL,
     $          dgxzs,
     $          CCTK_VARIABLE_REAL
     $     )

      CALL met_rad_der(origin,Dx,Dy,Dz,x,y,z,gyy,Extract_temp3d)
      CALL CCTK_SyncGroup(cctkGH,"extract::temps")
      call CCTK_Interp (ierror, cctkGH, handle, npoints, 3, 1, 1,
     $          nx, ny, nz,
     $          xs, ys, zs,
     $          CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,
     $          xmin, ymin, zmin,
     $          dx, dy, dz,
     $          Extract_temp3d,
     $          CCTK_VARIABLE_REAL,
     $          dgyys,
     $          CCTK_VARIABLE_REAL
     $     )

      CALL met_rad_der(origin,Dx,Dy,Dz,x,y,z,gyz,Extract_temp3d)
      CALL CCTK_SyncGroup(cctkGH,"extract::temps")
      call CCTK_Interp (ierror, cctkGH, handle, npoints, 3, 1, 1,
     $          nx, ny, nz,
     $          xs, ys, zs,
     $          CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,
     $          xmin, ymin, zmin,
     $          dx, dy, dz,
     $          Extract_temp3d,
     $          CCTK_VARIABLE_REAL,
     $          dgyzs,
     $          CCTK_VARIABLE_REAL
     $     )


      CALL met_rad_der(origin,Dx,Dy,Dz,x,y,z,gzz,Extract_temp3d)
      CALL CCTK_SyncGroup(cctkGH,"extract::temps")
      call CCTK_Interp (ierror, cctkGH, handle, npoints, 3, 1, 1,
     $          nx, ny, nz,
     $          xs, ys, zs,
     $          CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,
     $          xmin, ymin, zmin,
     $          dx, dy, dz,
     $          Extract_temp3d,
     $          CCTK_VARIABLE_REAL,
     $          dgzzs,
     $          CCTK_VARIABLE_REAL
     $     )


c     Calculate integrands for ADM masses
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)
         CALL CCTK_SyncGroup(cctkGH,"extract::temps")
         call CCTK_Interp (ierror, cctkGH, handle, npoints, 3, 1, 1,
     $        nx, ny, nz,
     $        xs, ys, zs,
     $        CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,
     $        xmin, ymin, zmin,
     $        dx, dy, dz,
     $        Extract_temp3d,
     $        CCTK_VARIABLE_REAL,
     $        ADMmass_int1,
     $        CCTK_VARIABLE_REAL
     $        )

      END IF

c     Conformal equation
      IF (do_ADMmass(2)) THEN
             ADMmass_int2 = -eta**2/2D0/3.1416D0*dPsis
      ENDIF	

c     Calculate integrands for momentum
c     ---------------------------------
      IF (do_momentum==1) THEN

         CALL momentum_integrand3D(origin,Dx,Dy,Dz,x,y,z,
     &        gxx,gxy,gxz,gyy,gyz,gzz,hxx,hxy,hxz,hyy,hyz,hzz,
     &        Extract_temp3d,Psi,Psi_power)
         CALL CCTK_SyncGroup(cctkGH,"extract::temps")
         call CCTK_Interp (ierror, cctkGH, handle, npoints, 3, 1, 1,
     $        nx, ny, nz,
     $        xs, ys, zs,
     $        CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,
     $        xmin, ymin, zmin,
     $        dx, dy, dz,
     $        Extract_temp3d,
     $        CCTK_VARIABLE_REAL,
     $        momentum_int1,
     $        CCTK_VARIABLE_REAL
     $        )

         CALL momentum_integrand3D(origin,Dx,Dy,Dz,x,y,z,
     &        gxx,gxy,gxz,gyy,gyz,gzz,hxx,hxy,hxz,hyy,hyz,hzz,
     &        Extract_temp3d,Psi,Psi_power)
         CALL CCTK_SyncGroup(cctkGH,"extract::temps")
         call CCTK_Interp (ierror, cctkGH, handle, npoints, 3, 1, 1,
     $        nx, ny, nz,
     $        xs, ys, zs,
     $        CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,
     $        xmin, ymin, zmin,
     $        dx, dy, dz,
     $        Extract_temp3d,
     $        CCTK_VARIABLE_REAL,
     $        momentum_int2,
     $        CCTK_VARIABLE_REAL
     $        )


         CALL momentum_integrand3D(origin,Dx,Dy,Dz,x,y,z,
     &        gxx,gxy,gxz,gyy,gyz,gzz,hxx,hxy,hxz,hyy,hyz,hzz,
     &        Extract_temp3d,Psi,Psi_power)
         CALL CCTK_SyncGroup(cctkGH,"extract::temps")
         call CCTK_Interp (ierror, cctkGH, handle, npoints, 3, 1, 1,
     $        nx, ny, nz,
     $        xs, ys, zs,
     $        CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,
     $        xmin, ymin, zmin,
     $        dx, dy, dz,
     $        Extract_temp3d,
     $        CCTK_VARIABLE_REAL,
     $        momentum_int3,
     $        CCTK_VARIABLE_REAL
     $        )

      END IF

c     Calculate integrands for spin
c     -----------------------------
      IF (do_spin==1) THEN

         CALL spin_integrand3D(origin,x,y,z,
     &        gxx,gxy,gxz,gyy,gyz,gzz,hxx,hxy,hxz,hyy,hyz,hzz,
     &        Extract_temp3d,Psi,Psi_power)
         CALL CCTK_SyncGroup(cctkGH,"extract::temps")
         call CCTK_Interp (ierror, cctkGH, handle, npoints, 3, 1, 1,
     $        nx, ny, nz,
     $        xs, ys, zs,
     $        CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,
     $        xmin, ymin, zmin,
     $        dx, dy, dz,
     $        Extract_temp3d,
     $        CCTK_VARIABLE_REAL,
     $        spin_int1,
     $        CCTK_VARIABLE_REAL
     $        )

         CALL spin_integrand3D(origin,x,y,z,
     &        gxx,gxy,gxz,gyy,gyz,gzz,hxx,hxy,hxz,hyy,hyz,hzz,
     &        Extract_temp3d,Psi,Psi_power)
         CALL CCTK_SyncGroup(cctkGH,"extract::temps")
         call CCTK_Interp (ierror, cctkGH, handle, npoints, 3, 1, 1,
     $        nx, ny, nz,
     $        xs, ys, zs,
     $        CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,
     $        xmin, ymin, zmin,
     $        dx, dy, dz,
     $        Extract_temp3d,
     $        CCTK_VARIABLE_REAL,
     $        spin_int2,
     $        CCTK_VARIABLE_REAL
     $        )


         CALL spin_integrand3D(origin,x,y,z,
     &        gxx,gxy,gxz,gyy,gyz,gzz,hxx,hxy,hxz,hyy,hyz,hzz,
     &        Extract_temp3d,Psi,Psi_power)
         CALL CCTK_SyncGroup(cctkGH,"extract::temps")
         call CCTK_Interp (ierror, cctkGH, handle, npoints, 3, 1, 1,
     $        nx, ny, nz,
     $        xs, ys, zs,
     $        CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,CCTK_VARIABLE_REAL,
     $        xmin, ymin, zmin,
     $        dx, dy, dz,
     $        Extract_temp3d,
     $        CCTK_VARIABLE_REAL,
     $        spin_int3,
     $        CCTK_VARIABLE_REAL
     $        )

      END IF

      END SUBROUTINE D3_to_D2