#include "cctk.h"
#include "cctk_Parameters.h"
#include "cctk_Arguments.h"
#include "CactusEinstein/Einstein/src/Einstein.h"

      subroutine setupIVP(CCTK_FARGUMENTS)
      implicit none
      DECLARE_CCTK_FARGUMENTS
c     local vars
      integer i,j,k
      CCTK_REAL det,uxx,uxy,uxz,uyy,uyz,uzz
      CCTK_REAL p4, dp4,two, dx,dy,dz,delta_x,delta_y,delta_z
#include "INCLUDE/ivp_deriv_temps_con.inc"
#include "INCLUDE/ivp_ricci_temps.inc"

      two = 2.0d0
c     dx settings 
      dx=cctk_delta_space(1)
      dy=cctk_delta_space(2)
      dz=cctk_delta_space(3)

c     first, compute the inverse metric
      do k=1,cctk_lsh(3)
      do j=1,cctk_lsh(2)
      do i=1,cctk_lsh(1)
         if (conformal_state .eq. CONFORMAL_METRIC) then
            p4 = psi(i,j,k) ** 4
            dp4 = 1.0D0 / p4
         else
            p4 = 1.0D0
            dp4 = 1.0D0
         endif
         det= -(gxz(i,j,k)**2*gyy(i,j,k)) +
     .        2.0d0*gxy(i,j,k)*gxz(i,j,k)*gyz(i,j,k) -
     .        gxx(i,j,k)*gyz(i,j,k)**2  -
     .        gxy(i,j,k)**2*gzz(i,j,k) +
     .        gxx(i,j,k)*gyy(i,j,k)*gzz(i,j,k)
         uxx=(-gyz(i,j,k)**2 + gyy(i,j,k)*gzz(i,j,k))/det/p4
         uxy=(gxz(i,j,k)*gyz(i,j,k) - gxy(i,j,k)*gzz(i,j,k))/det/p4
         uyy=(-gxz(i,j,k)**2 + gxx(i,j,k)*gzz(i,j,k))/det/p4
         uxz=(-gxz(i,j,k)*gyy(i,j,k) + gxy(i,j,k)*gyz(i,j,k))/det/p4
         uyz=(gxy(i,j,k)*gxz(i,j,k) - gxx(i,j,k)*gyz(i,j,k))/det/p4
         uzz=(-gxy(i,j,k)**2 + gxx(i,j,k)*gyy(i,j,k))/det/p4
         ivp_det(i,j,k) = det
         ivp_uxx(i,j,k) = uxx
         ivp_uxy(i,j,k) = uxy
         ivp_uxz(i,j,k) = uxz
         ivp_uyy(i,j,k) = uyy
         ivp_uyz(i,j,k) = uyz
         ivp_uzz(i,j,k) = uzz
      enddo
      enddo
      enddo

c     now, compute scalar curvature
      do k=2,cctk_lsh(3)-1
      do j=2,cctk_lsh(2)-1
      do i=2,cctk_lsh(1)-1

c        first, compute the 3-scalar curvature
         if (conformal_state .eq. CONFORMAL_METRIC) then
            p4 = psi(i,j,k) ** 4
            dp4 = 1.0D0 / p4
         else
            p4 = 1.0D0
            dp4 = 1.0D0
         endif
         det = ivp_det(i,j,k)
         uxx = ivp_uxx(i,j,k)
         uxy = ivp_uxy(i,j,k)
         uxz = ivp_uxz(i,j,k)
         uyy = ivp_uyy(i,j,k)
         uyz = ivp_uyz(i,j,k)
         uzz = ivp_uzz(i,j,k)
#include "INCLUDE/ivp_setup_derivs_con.inc"
#include "INCLUDE/ivp_ricci_calc.inc"
         ivp_scalarc(i,j,k) = riccixx*uxx + ricciyy*uyy + riccizz*uzz +
     .      two*riccixy*uxy + two*riccixz*uxz + two*ricciyz*uyz
      enddo
      enddo
      enddo

      call CCTK_SyncGroup(cctkGH, "ivp::scalarc")
c     sync the ivp_scalarc GF.

      return
      end