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

      subroutine ivp_init_data(CCTK_FARGUMENTS)

      implicit none

      DECLARE_CCTK_FARGUMENTS
      DECLARE_CCTK_PARAMETERS
      INTEGER CCTK_Equals

/* #define DEBUG_IVP */

#ifdef DEBUG_IVP
      write(*,*)'initializing data in ivp_init_data'
#endif

c     conformal stuff
c      conformal_state = NOCONFORMAL_METRIC 
      conformal_state = CONFORMAL_METRIC 

      psi = 1.0d0
      psix = 0.0d0
      psiy = 0.0d0
      psiz = 0.0d0    
      psixy = 0.0d0
      psixz = 0.0d0
      psiyz = 0.0d0  
      psixx = 0.0d0
      psiyy = 0.0d0
      psizz = 0.0d0  

      if(CCTK_Equals(ivp_init_method,"flat").eq.1) then

         write(*,*)'IVP_TEST_METRIC: flat space'
         
         gxx = 1.0d0
         gxy = 0.0d0
         gxz = 0.0d0
         gyy = 1.0d0
         gyz = 0.0d0
         gzz = 1.0d0
         
         kxx = 0.0d0
         kxy = 0.0d0
         kxz = 0.0d0
         kyy = 0.0d0
         kyz = 0.0d0
         kzz = 0.0d0
         
      elseif(CCTK_Equals(ivp_init_method,"pflat_conf").eq.1) then
         
         write(*,*)'IVP_TEST_METRIC: perturbed flat space'
         write(*,*)'conforamlly flat'
         write(*,*)'and perturbed ext. curv.'
         
         gxx = 1.0d0 + 0.1d0*exp(-r**2/5.0d0)
         gxy = 0.0d0
         gxz = 0.0d0
         gyy = 1.0d0 + 0.1d0*exp(-r**2/5.0d0)
         gyz = 0.0d0
         gzz = 1.0d0 + 0.1d0*exp(-r**2/5.0d0)
         
         kxx = 0.1d0*exp(-r**2/5.0d0)
         kxy = 0.0d0
         kxz = 0.0d0
         kyy = 0.1d0*exp(-r**2/5.0d0)
         kyz = 0.0d0
         kzz = 0.1d0*exp(-r**2/5.0d0)
         
      elseif(CCTK_Equals(ivp_init_method,"pflat_nconf").eq.1) then

         write(*,*)'IVP_TEST_METRIC: perturbed flat space'
         write(*,*)'not conforamlly flat'
         write(*,*)'and perturbed ext. curv.'
         
         gxx = 1.0d0 + 0.1d0*exp(-r**2/5.0d0)
         gxy = 0.0d0
         gxz = 0.0d0
         gyy = 1.0d0 - 0.1d0*exp(-r**2/5.0d0)
         gyz = 0.0d0
         gzz = 1.0d0 

         kxx = 0.1d0*exp(-r**2/5.0d0)
         kxy = 0.0d0
         kxz = 0.0d0
         kyy = -0.1d0*exp(-r**2/5.0d0)
         kyz = 0.0d0
         kzz = 0.0d0

      elseif(CCTK_Equals(ivp_init_method,"pflat_nconf_0K").eq.1) then

         write(*,*)'IVP_TEST_METRIC: perturbed flat space'
         write(*,*)'not conforamlly flat'
         write(*,*)'and zero ext. curv.'
         
         gxx = 1.0d0 + 0.1d0*exp(-r**2/5.0d0)
         gxy = 0.0d0
         gxz = 0.0d0
         gyy = 1.0d0 - 0.1d0*exp(-r**2/5.0d0)
         gyz = 0.0d0
         gzz = 1.0d0 

         kxx = 0.0d0
         kxy = 0.0d0
         kxz = 0.0d0
         kyy = 0.0d0
         kyz = 0.0d0
         kzz = 0.0d0

      elseif(CCTK_Equals(ivp_init_method,"pflat_pK").eq.1) then

         write(*,*)'IVP_TEST_METRIC: flat space'
         write(*,*)'and perturbed ext. curv.'
         
         gxx = 1.0d0 
         gxy = 0.0d0
         gxz = 0.0d0
         gyy = 1.0d0 
         gyz = 0.0d0
         gzz = 1.0d0 

         kxx = 0.1d0*exp(-r**2/5.0d0)
         kxy = 0.0d0
         kxz = 0.0d0
         kyy =-0.1d0*exp(-r**2/5.0d0)
         kyz = 0.0d0
         kzz = 0.0d0

      endif

      return
      end