#include "cctk.h"

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

      SUBROUTINE D3_to_D2(cctkGH,conformal_state,do_ADMmass,do_momentum,do_spin,
     &     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,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 2-sphere.
c
c     ------------------------------------------------------------------------

      USE met_rad_der_int
      USE ADMmass_integrand3D_int
      USE momentum_integrand3D_int
      USE spin_integrand3D_int

      IMPLICIT NONE

      DECLARE_CCTK_FUNCTIONS

c     Input variables

      CCTK_POINTER :: cctkGH

      INTEGER,INTENT(IN) ::
     &     conformal_state,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)
      CCTK_STRING,INTENT(IN) :: 
     &     interpolation_operator

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, num_arrays, ierror
      INTEGER :: interp_handle, param_table_handle, coord_system_handle
      CCTK_POINTER, dimension(3) :: interp_coords
      INTEGER,   dimension(8)    :: in_array_indices
      CCTK_POINTER, dimension(8) :: out_arrays
      CCTK_INT, dimension(8)     :: out_array_type_codes

      character(128) options_string
      character(128) operator
      CCTK_INT nchars

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, parameter table, and coordinate system handles
c     --------------------------------------------------------------------
      interp_handle = -1
      param_table_handle = -1
      coord_system_handle = -1


      call CCTK_FortranString(nchars,interpolation_operator,operator)
      call CCTK_InterpHandle (interp_handle,operator)
      if (interp_handle .lt. 0) then
        call CCTK_WARN(0,"Cannot get handle for interpolation ! Forgot to activate an implementation providing interpolation operators ??")
      endif

      options_string = "order = " // char(ichar('0') + interpolation_order)
      call Util_TableCreateFromString (param_table_handle, options_string)
      if (param_table_handle .lt. 0) then
        call CCTK_WARN(0,"Cannot create parameter table for interpolator")
      endif

      call CCTK_CoordSystemHandle (coord_system_handle, "cart3d")
      if (coord_system_handle .lt. 0) then
        call CCTK_WARN(0,"Cannot get handle for cart3d coordinate system ! Forgot to activate an implementation providing coordinates ??")
      endif


c     fill in the input/output arrays for the interpolator
c     ----------------------------------------------------
      interp_coords(1) = CCTK_PointerTo(xs)
      interp_coords(2) = CCTK_PointerTo(ys)
      interp_coords(3) = CCTK_PointerTo(zs)

      call CCTK_VarIndex(in_array_indices(1), "extract::g00")
      call CCTK_VarIndex(in_array_indices(2), "admbase::gxx")
      call CCTK_VarIndex(in_array_indices(3), "admbase::gxy")
      call CCTK_VarIndex(in_array_indices(4), "admbase::gxz")
      call CCTK_VarIndex(in_array_indices(5), "admbase::gyy")
      call CCTK_VarIndex(in_array_indices(6), "admbase::gyz")
      call CCTK_VarIndex(in_array_indices(7), "admbase::gzz")
      call CCTK_VarIndex(in_array_indices(8), "staticconformal::psi")

      out_arrays(1) = CCTK_PointerTo(g00s)
      out_arrays(2) = CCTK_PointerTo(gxxs)
      out_arrays(3) = CCTK_PointerTo(gxys)
      out_arrays(4) = CCTK_PointerTo(gxzs)
      out_arrays(5) = CCTK_PointerTo(gyys)
      out_arrays(6) = CCTK_PointerTo(gyzs)
      out_arrays(7) = CCTK_PointerTo(gzzs)
      out_arrays(8) = CCTK_PointerTo(Psis)

      out_array_type_codes = CCTK_VARIABLE_REAL


c     Project un-physical metric and conformal factor  onto sphere
c     ------------------------------------------------------------
      if (conformal_state > 0) then
        num_arrays = 8
      else
        num_arrays = 7
      end if

      call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
     $                            param_table_handle, coord_system_handle,
     $                            npoints, CCTK_VARIABLE_REAL, interp_coords,
     $                            num_arrays, in_array_indices,
     $                            num_arrays, out_array_type_codes, out_arrays)
      if (ierror < 0) then
        call CCTK_WARN (1, "interpolator call returned an error code");
      endif


c     Calculate radial derivatives and project onto sphere
c     ----------------------------------------------------
      call CCTK_VarIndex(in_array_indices(1), "extract::temp3d")

      if (conformal_state > 0) then
         CALL met_rad_der(origin,Dx,Dy,Dz,x,y,z,Psi,Extract_temp3d)
         CALL CCTK_SyncGroup(ierror,cctkGH,"extract::temps")

         out_arrays(1) = CCTK_PointerTo(dPsis)
         call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
     $                               param_table_handle, coord_system_handle,
     $                               npoints, CCTK_VARIABLE_REAL, interp_coords,
     $                               1, in_array_indices(1),
     $                               1, out_array_type_codes, out_arrays(1))
         if (ierror < 0) then
           call CCTK_WARN (1, "interpolator call returned an error code");
         endif
      end if

      CALL met_rad_der(origin,Dx,Dy,Dz,x,y,z,gxx,Extract_temp3d)
      CALL CCTK_SyncGroup(ierror,cctkGH,"extract::temps")

      out_arrays(1) = CCTK_PointerTo(dgxxs)
      call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
     $                            param_table_handle, coord_system_handle,
     $                            npoints, CCTK_VARIABLE_REAL, interp_coords,
     $                            1, in_array_indices(1),
     $                            1, out_array_type_codes, out_arrays(1))
      if (ierror < 0) then
        call CCTK_WARN (1, "interpolator call returned an error code");
      endif

      CALL met_rad_der(origin,Dx,Dy,Dz,x,y,z,gxy,Extract_temp3d)
      CALL CCTK_SyncGroup(ierror,cctkGH,"extract::temps")

      out_arrays(1) = CCTK_PointerTo(dgxys)
      call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
     $                            param_table_handle, coord_system_handle,
     $                            npoints, CCTK_VARIABLE_REAL, interp_coords,
     $                            1, in_array_indices(1),
     $                            1, out_array_type_codes, out_arrays(1))
      if (ierror < 0) then
        call CCTK_WARN (1, "interpolator call returned an error code");
      endif

      CALL met_rad_der(origin,Dx,Dy,Dz,x,y,z,gxz,Extract_temp3d)
      CALL CCTK_SyncGroup(ierror,cctkGH,"extract::temps")

      out_arrays(1) = CCTK_PointerTo(dgxzs)
      call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
     $                            param_table_handle, coord_system_handle,
     $                            npoints, CCTK_VARIABLE_REAL, interp_coords,
     $                            1, in_array_indices(1),
     $                            1, out_array_type_codes, out_arrays(1))
      if (ierror < 0) then
        call CCTK_WARN (1, "interpolator call returned an error code");
      endif

      CALL met_rad_der(origin,Dx,Dy,Dz,x,y,z,gyy,Extract_temp3d)
      CALL CCTK_SyncGroup(ierror,cctkGH,"extract::temps")

      out_arrays(1) = CCTK_PointerTo(dgyys)
      call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
     $                            param_table_handle, coord_system_handle,
     $                            npoints, CCTK_VARIABLE_REAL, interp_coords,
     $                            1, in_array_indices(1),
     $                            1, out_array_type_codes, out_arrays(1))
      if (ierror < 0) then
        call CCTK_WARN (1, "interpolator call returned an error code");
      endif

      CALL met_rad_der(origin,Dx,Dy,Dz,x,y,z,gyz,Extract_temp3d)
      CALL CCTK_SyncGroup(ierror,cctkGH,"extract::temps")

      out_arrays(1) = CCTK_PointerTo(dgyzs)
      call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
     $                            param_table_handle, coord_system_handle,
     $                            npoints, CCTK_VARIABLE_REAL, interp_coords,
     $                            1, in_array_indices(1),
     $                            1, out_array_type_codes, out_arrays(1))
      if (ierror < 0) then
        call CCTK_WARN (1, "interpolator call returned an error code");
      endif

      CALL met_rad_der(origin,Dx,Dy,Dz,x,y,z,gzz,Extract_temp3d)
      CALL CCTK_SyncGroup(ierror,cctkGH,"extract::temps")

      out_arrays(1) = CCTK_PointerTo(dgzzs)
      call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
     $                            param_table_handle, coord_system_handle,
     $                            npoints, CCTK_VARIABLE_REAL, interp_coords,
     $                            1, in_array_indices(1),
     $                            1, out_array_type_codes, out_arrays(1))
      if (ierror < 0) then
        call CCTK_WARN (1, "interpolator call returned an error code");
      endif


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,conformal_state)
         CALL CCTK_SyncGroup(ierror,cctkGH,"extract::temps")

         out_arrays(1) = CCTK_PointerTo(ADMmass_int1)
         call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
     $                               param_table_handle, coord_system_handle,
     $                               npoints, CCTK_VARIABLE_REAL, interp_coords,
     $                               1, in_array_indices(1),
     $                               1, out_array_type_codes, out_arrays(1))
         if (ierror < 0) then
           call CCTK_WARN (1, "interpolator call returned an error code");
         endif
      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,conformal_state)
         CALL CCTK_SyncGroup(ierror,cctkGH,"extract::temps")

         out_arrays(1) = CCTK_PointerTo(momentum_int1)
         call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
     $                               param_table_handle, coord_system_handle,
     $                               npoints, CCTK_VARIABLE_REAL, interp_coords,
     $                               1, in_array_indices(1),
     $                               1, out_array_type_codes, out_arrays(1))
         if (ierror < 0) then
           call CCTK_WARN (1, "interpolator call returned an error code");
         endif

         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,conformal_state)
         CALL CCTK_SyncGroup(ierror,cctkGH,"extract::temps")

         out_arrays(1) = CCTK_PointerTo(momentum_int2)
         call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
     $                               param_table_handle, coord_system_handle,
     $                               npoints, CCTK_VARIABLE_REAL, interp_coords,
     $                               1, in_array_indices(1),
     $                               1, out_array_type_codes, out_arrays(1))
         if (ierror < 0) then
           call CCTK_WARN (1, "interpolator call returned an error code");
         endif

         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,conformal_state)
         CALL CCTK_SyncGroup(ierror,cctkGH,"extract::temps")

         out_arrays(1) = CCTK_PointerTo(momentum_int3)
         call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
     $                               param_table_handle, coord_system_handle,
     $                               npoints, CCTK_VARIABLE_REAL, interp_coords,
     $                               1, in_array_indices(1),
     $                               1, out_array_type_codes, out_arrays(1))
         if (ierror < 0) then
           call CCTK_WARN (1, "interpolator call returned an error code");
         endif
      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,conformal_state)
         CALL CCTK_SyncGroup(ierror,cctkGH,"extract::temps")

         out_arrays(1) = CCTK_PointerTo(spin_int1)
         call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
     $                               param_table_handle, coord_system_handle,
     $                               npoints, CCTK_VARIABLE_REAL, interp_coords,
     $                               1, in_array_indices(1),
     $                               1, out_array_type_codes, out_arrays(1))
         if (ierror < 0) then
           call CCTK_WARN (1, "interpolator call returned an error code");
         endif

         CALL spin_integrand3D(origin,x,y,z,
     &        gxx,gxy,gxz,gyy,gyz,gzz,hxx,hxy,hxz,hyy,hyz,hzz,
     &        Extract_temp3d,Psi,Psi_power,conformal_state)
         CALL CCTK_SyncGroup(ierror,cctkGH,"extract::temps")

         out_arrays(1) = CCTK_PointerTo(spin_int2)
         call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
     $                               param_table_handle, coord_system_handle,
     $                               npoints, CCTK_VARIABLE_REAL, interp_coords,
     $                               1, in_array_indices(1),
     $                               1, out_array_type_codes, out_arrays(1))
         if (ierror < 0) then
           call CCTK_WARN (1, "interpolator call returned an error code");
         endif

         CALL spin_integrand3D(origin,x,y,z,
     &        gxx,gxy,gxz,gyy,gyz,gzz,hxx,hxy,hxz,hyy,hyz,hzz,
     &        Extract_temp3d,Psi,Psi_power,conformal_state)
         CALL CCTK_SyncGroup(ierror,cctkGH,"extract::temps")

         out_arrays(1) = CCTK_PointerTo(spin_int3)
         call CCTK_InterpGridArrays (ierror, cctkGH, 3, interp_handle,
     $                               param_table_handle, coord_system_handle,
     $                               npoints, CCTK_VARIABLE_REAL, interp_coords,
     $                               1, in_array_indices(1),
     $                               1, out_array_type_codes, out_arrays(1))
         if (ierror < 0) then
           call CCTK_WARN (1, "interpolator call returned an error code");
         endif

      END IF

      END SUBROUTINE D3_to_D2