#include "cctk.h"

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

      SUBROUTINE D3_to_D2(cctkGH,do_ADMmass,do_momentum,do_spin,
     &     Psi_power,origin,myproc,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,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,interpolation_order
      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,interp_handle,coord_system_handle,ierror
      INTEGER,   DIMENSION(8)   :: in_array_indices


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     Get the interpolator and coordinate system handle
c     -------------------------------------------------
      interp_handle = -1
      coord_system_handle = -1

      if (interpolation_order .eq. 1) then
        call CCTK_InterpHandle (interp_handle, "first-order uniform cartesian")
      else if (interpolation_order .eq. 2) then
        call CCTK_InterpHandle (interp_handle, "second-order uniform cartesian")
      else if (interpolation_order .eq. 3) then
        call CCTK_InterpHandle (interp_handle, "third-order uniform cartesian")
      endif

      call CCTK_CoordSystemHandle (coord_system_handle, "cart3d")

      if (interp_handle .lt. 0 .or. coord_system_handle .lt. 0) then
        call CCTK_WARN (0, "Couldn't get handles for interpolation operator and/or coordinate system")
      endif

c     Get indices of GFs to interpolate
c     ---------------------------------
      call CCTK_VarIndex(in_array_indices(1), "einstein::psi")
      call CCTK_VarIndex(in_array_indices(2), "extract::g00")
      call CCTK_VarIndex(in_array_indices(3), "einstein::gxx")
      call CCTK_VarIndex(in_array_indices(4), "einstein::gxy")
      call CCTK_VarIndex(in_array_indices(5), "einstein::gxz")
      call CCTK_VarIndex(in_array_indices(6), "einstein::gyy")
      call CCTK_VarIndex(in_array_indices(7), "einstein::gyz")
      call CCTK_VarIndex(in_array_indices(8), "einstein::gzz")


c     Project un-physical metric and conformal factor  onto sphere
c     ------------------------------------------------------------
      call CCTK_InterpGV (ierror, cctkGH, interp_handle, coord_system_handle,
     $                    npoints, 8, 8,
     $                    xs, ys, zs,
     $                    CCTK_VARIABLE_REAL, CCTK_VARIABLE_REAL,
     $                    CCTK_VARIABLE_REAL,
     $                    in_array_indices(1), in_array_indices(2),
     $                    in_array_indices(3), in_array_indices(4),
     $                    in_array_indices(5), in_array_indices(6),
     $                    in_array_indices(7), in_array_indices(8),
     $                    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(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)

      CALL met_rad_der(origin,Dx,Dy,Dz,x,y,z,gxx,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),
     $                    dgxxs,
     $                    CCTK_VARIABLE_REAL)

      CALL met_rad_der(origin,Dx,Dy,Dz,x,y,z,gxy,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),
     $                    dgxys,
     $                    CCTK_VARIABLE_REAL)

      CALL met_rad_der(origin,Dx,Dy,Dz,x,y,z,gxz,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),
     $                    dgxzs,
     $                    CCTK_VARIABLE_REAL)

      CALL met_rad_der(origin,Dx,Dy,Dz,x,y,z,gyy,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),
     $                    dgyys,
     $                    CCTK_VARIABLE_REAL)

      CALL met_rad_der(origin,Dx,Dy,Dz,x,y,z,gyz,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),
     $                    dgyzs,
     $                    CCTK_VARIABLE_REAL)

      CALL met_rad_der(origin,Dx,Dy,Dz,x,y,z,gzz,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),
     $                    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(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),
     $                       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(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),
     $                       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(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),
     $                       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(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),
     $                       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(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),
     $                       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(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),
     $                       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(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),
     $                       spin_int3,
     $                       CCTK_VARIABLE_REAL)

      END IF

      END SUBROUTINE D3_to_D2