














































c
c
c
c
c

c




































































 


































































































































































 


















































































































































































































 



















































 










 








































































































































































































































 










                                
                                

                                
                                

                                
                                

                                















































      subroutine  Exact__Bona_Masso_data(
     $     decoded_exact_model,
     $     x, y, z, t,
     $     gxx, gyy, gzz, gxy, gyz, gxz,
     $     hxx, hyy, hzz, hxy, hyz, hxz,
     $     psi, psix, psiy, psiz,
     $     psixx, psiyy, psizz, psixy, psiyz, psixz,
     $     dxgxx, dxgyy, dxgzz, dxgxy, dxgyz, dxgxz,
     $     dygxx, dygyy, dygzz, dygxy, dygyz, dygxz,
     $     dzgxx, dzgyy, dzgzz, dzgxy, dzgyz, dzgxz,
     $     alp, dtalp, ax, ay, az,
     $     betax, betay, betaz, dtbetax, dtbetay, dtbetaz,
     $     bxx, bxy, bxz, byx, byy, byz, bzx, bzy, bzz)

      implicit none
      INTEGER*4, parameter :: izero = 0
      integer, parameter :: ik = kind(izero)
      INTEGER*4 decoded_exact_model
      logical psi_on
      REAL*8 x, y, z, t,
     $       gxx, gyy, gzz, gxy, gyz, gxz,
     $       hxx, hyy, hzz, hxy, hyz, hxz,
     $       psi, psix, psiy, psiz,
     $       psixx, psiyy, psizz, psixy, psiyz, psixz,
     $       dxgxx, dxgyy, dxgzz, dxgxy, dxgyz, dxgxz,
     $       dygxx, dygyy, dygzz, dygxy, dygyz, dygxz,
     $       dzgxx, dzgyy, dzgzz, dzgxy, dzgyz, dzgxz,
     $       alp, dtalp, ax, ay, az,
     $       betax, betay, betaz, dtbetax, dtbetay, dtbetaz,
     $       bxx, bxy, bxz, byx, byy, byz, bzx, bzy, bzz

C     gxx is g_xx etc.
C     hxx is K_xx etc.
C     dxgyy is (/2) dg_yy / dx (sic!)
C     alp is N, betax is N^x etc.
C     bxy is (/2) dN^y / dx (sic and sic!)
C     ax is dN / dx / N (sic!)

      REAL*8
     $       gdtt, gdtx, gdty, gdtz, 
     $       gutt, gutx, guty, gutz, 
     $       guxx, guyy, guzz, guxy, guyz, guxz,
     $       dgdtt, dgdtx, dgdty, dgdtz,
     $       dgdxx, dgdyy, dgdzz, dgdxy, dgdyz, dgdxz,
     $       dgutt, dgutx, dguty, dgutz,
     $       dguxx, dguyy, dguzz, dguxy, dguyz, dguxz,
     $       dpsi

C     Save, if we have to provide the conformal factor

      psi_on = psi .gt. 0.0d0

C     Get the spacetime metric and its inverse at the base point.

      call Exact__metric(
     $     decoded_exact_model,
     $     x, y, z, t,
     $     gdtt, gdtx, gdty, gdtz, 
     $     gxx, gyy, gzz, gxy, gyz, gxz,
     $     gutt, gutx, guty, gutz, 
     $     guxx, guyy, guzz, guxy, guyz, guxz,
     $     psi)

C     Calculate lapse and shift from the upper metric.

      alp = 1.d0 / sqrt(-gutt)

      betax = - gutx / gutt
      betay = - guty / gutt
      betaz = - gutz / gutt

C     In order to calculate the derivatives of the lapse and shift from
C     the contravariant metric, use g^tt = -1/N^2 and g^i0 = N^i / N^2
C     Note that ax is dN/dx / N and that bxy is 1/2 dN^y / dx.

C     Calculate x-derivatives. 

      call Exact__metric_deriv(
     $     decoded_exact_model,
     $     1_ik,
     $     x, y, z, t,
     $     dgdtt, dgdtx, dgdty, dgdtz,
     $     dgdxx, dgdyy, dgdzz, dgdxy, dgdyz, dgdxz,
     $     dgutt, dgutx, dguty, dgutz,
     $     dguxx, dguyy, dguzz, dguxy, dguyz, dguxz,
     $     dpsi)

      dxgxx = 0.5d0 * dgdxx
      dxgyy = 0.5d0 * dgdyy
      dxgzz = 0.5d0 * dgdzz
      dxgxy = 0.5d0 * dgdxy
      dxgyz = 0.5d0 * dgdyz
      dxgxz = 0.5d0 * dgdxz

      ax = - 0.5d0 * dgutt / gutt

      bxx = 0.5d0 * (- dgutx * gutt + gutx * dgutt) / gutt**2
      bxy = 0.5d0 * (- dguty * gutt + guty * dgutt) / gutt**2
      bxz = 0.5d0 * (- dgutz * gutt + gutz * dgutt) / gutt**2

      if (psi_on) then
         psix = dpsi
      end if

C     Calculate y-derivatives.

      call Exact__metric_deriv(
     $     decoded_exact_model,
     $     2_ik,
     $     x, y, z, t,
     $     dgdtt, dgdtx, dgdty, dgdtz,
     $     dgdxx, dgdyy, dgdzz, dgdxy, dgdyz, dgdxz,
     $     dgutt, dgutx, dguty, dgutz,
     $     dguxx, dguyy, dguzz, dguxy, dguyz, dguxz,
     $     dpsi)

      dygxx = 0.5d0 * dgdxx
      dygyy = 0.5d0 * dgdyy
      dygzz = 0.5d0 * dgdzz
      dygxy = 0.5d0 * dgdxy
      dygyz = 0.5d0 * dgdyz
      dygxz = 0.5d0 * dgdxz

      ay = - 0.5d0 * dgutt / gutt

      byx = 0.5d0 * (- dgutx * gutt + gutx * dgutt) / gutt**2
      byy = 0.5d0 * (- dguty * gutt + guty * dgutt) / gutt**2
      byz = 0.5d0 * (- dgutz * gutt + gutz * dgutt) / gutt**2

      if (psi_on) then
         psiy = dpsi
      end if

C     Calculate z-derivatives.

      call Exact__metric_deriv(
     $     decoded_exact_model,
     $     3_ik,
     $     x, y, z, t,
     $     dgdtt, dgdtx, dgdty, dgdtz,
     $     dgdxx, dgdyy, dgdzz, dgdxy, dgdyz, dgdxz,
     $     dgutt, dgutx, dguty, dgutz,
     $     dguxx, dguyy, dguzz, dguxy, dguyz, dguxz,
     $     dpsi)

      dzgxx = 0.5d0 * dgdxx
      dzgyy = 0.5d0 * dgdyy
      dzgzz = 0.5d0 * dgdzz
      dzgxy = 0.5d0 * dgdxy
      dzgyz = 0.5d0 * dgdyz
      dzgxz = 0.5d0 * dgdxz

      az = - 0.5d0 * dgutt / gutt

      bzx = 0.5d0 * (- dgutx * gutt + gutx * dgutt) / gutt**2
      bzy = 0.5d0 * (- dguty * gutt + guty * dgutt) / gutt**2
      bzz = 0.5d0 * (- dgutz * gutt + gutz * dgutt) / gutt**2

      if (psi_on) then
         psiz = dpsi
      end if

C     Calculate t-derivatives, and extrinsic curvature.

      call Exact__metric_deriv(
     $     decoded_exact_model,
     $     0_ik,
     $     x, y, z, t,
     $     dgdtt, dgdtx, dgdty, dgdtz,
     $     dgdxx, dgdyy, dgdzz, dgdxy, dgdyz, dgdxz,
     $     dgutt, dgutx, dguty, dgutz,
     $     dguxx, dguyy, dguzz, dguxy, dguyz, dguxz,
     $     dpsi)

      hxx = - 0.5d0 * dgdxx / alp
     $     + (dxgxx * betax + dygxx * betay + dzgxx * betaz
     $     + 2.d0 * (bxx * gxx + bxy * gxy + bxz * gxz)) / alp

      hyy = - 0.5d0 * dgdyy / alp
     $     + (dxgyy * betax + dygyy * betay + dzgyy * betaz
     $     + 2.d0 * (byx * gxy + byy * gyy + byz * gyz)) / alp

      hzz = - 0.5d0 * dgdzz / alp
     $     + (dxgzz * betax + dygzz * betay + dzgzz * betaz
     $     + 2.d0 * (bzx * gxz + bzy * gyz + bzz * gzz)) / alp

      hxy = - 0.5d0 * dgdxy / alp
     $     + (dxgxy * betax + dygxy * betay + dzgxy * betaz
     $     + bxx * gxy + bxy * gyy + bxz * gyz
     $     + byx * gxx + byy * gxy + byz * gxz) / alp

      hyz = - 0.5d0 * dgdyz / alp
     $     + (dxgyz * betax + dygyz * betay + dzgyz * betaz
     $     + byx * gxz + byy * gyz + byz * gzz
     $     + bzx * gxy + bzy * gyy + bzz * gyz) / alp

      hxz = - 0.5d0 * dgdxz / alp
     $     + (dxgxz * betax + dygxz * betay + dzgxz * betaz
     $     + bxx * gxz + bxy * gyz + bxz * gzz
     $     + bzx * gxx + bzy * gxy + bzz * gxz) / alp

C     Calculate time derivatives of  lapse and shift

C     alp = 1.d0 / sqrt(-gutt)
      dtalp = 0.5d0 / sqrt(-gutt)**3 * dgutt

C     betax = - gutx / gutt
C     betay = - guty / gutt
C     betaz = - gutz / gutt
      dtbetax = (- dgutx * gutt + gutx * dgutt) / gutt**2
      dtbetay = (- dguty * gutt + guty * dgutt) / gutt**2
      dtbetaz = (- dgutz * gutt + gutz * dgutt) / gutt**2

C     Calculate second derivatives of the conformal factor

      if (psi_on) then
        call Exact__metric_deriv2(
     $        decoded_exact_model,
     $        1_ik, 1_ik,
     $        x, y, z, t,
     $        dgdtt, dgdtx, dgdty, dgdtz,
     $        dgdxx, dgdyy, dgdzz, dgdxy, dgdyz, dgdxz,
     $        dgutt, dgutx, dguty, dgutz,
     $        dguxx, dguyy, dguzz, dguxy, dguyz, dguxz,
     $        psixx)
        call Exact__metric_deriv2(
     $        decoded_exact_model,
     $        1_ik, 2_ik,
     $        x, y, z, t,
     $        dgdtt, dgdtx, dgdty, dgdtz,
     $        dgdxx, dgdyy, dgdzz, dgdxy, dgdyz, dgdxz,
     $        dgutt, dgutx, dguty, dgutz,
     $        dguxx, dguyy, dguzz, dguxy, dguyz, dguxz,
     $        psixy)
        call Exact__metric_deriv2(
     $        decoded_exact_model,
     $        1_ik, 3_ik,
     $        x, y, z, t,
     $        dgdtt, dgdtx, dgdty, dgdtz,
     $        dgdxx, dgdyy, dgdzz, dgdxy, dgdyz, dgdxz,
     $        dgutt, dgutx, dguty, dgutz,
     $        dguxx, dguyy, dguzz, dguxy, dguyz, dguxz,
     $        psixz)
        call Exact__metric_deriv2(
     $        decoded_exact_model,
     $        2_ik, 2_ik,
     $        x, y, z, t,
     $        dgdtt, dgdtx, dgdty, dgdtz,
     $        dgdxx, dgdyy, dgdzz, dgdxy, dgdyz, dgdxz,
     $        dgutt, dgutx, dguty, dgutz,
     $        dguxx, dguyy, dguzz, dguxy, dguyz, dguxz,
     $        psiyy)
        call Exact__metric_deriv2(
     $        decoded_exact_model,
     $        2_ik, 3_ik,
     $        x, y, z, t,
     $        dgdtt, dgdtx, dgdty, dgdtz,
     $        dgdxx, dgdyy, dgdzz, dgdxy, dgdyz, dgdxz,
     $        dgutt, dgutx, dguty, dgutz,
     $        dguxx, dguyy, dguzz, dguxy, dguyz, dguxz,
     $        psiyz)
        call Exact__metric_deriv2(
     $        decoded_exact_model,
     $        3_ik, 3_ik,
     $        x, y, z, t,
     $        dgdtt, dgdtx, dgdty, dgdtz,
     $        dgdxx, dgdyy, dgdzz, dgdxy, dgdyz, dgdxz,
     $        dgutt, dgutx, dguty, dgutz,
     $        dguxx, dguyy, dguzz, dguxy, dguyz, dguxz,
     $        psizz)
        
      end if
      return
      end



      subroutine Exact__metric_deriv(
     $     decoded_exact_model,
     $     dir,
     $     x, y, z, t,
     $     dgdtt, dgdtx, dgdty, dgdtz,
     $     dgdxx, dgdyy, dgdzz, dgdxy, dgdyz, dgdxz,
     $     dgutt, dgutx, dguty, dgutz,
     $     dguxx, dguyy, dguzz, dguxy, dguyz, dguxz,
     $     dpsi)
      
      implicit none
      REAL*8 Alvi__mass1
      integer, parameter :: cctki_use_Alvi__mass1 = kind(Alvi__mass1)
      REAL*8 Alvi__mass2
      integer, parameter :: cctki_use_Alvi__mass2 = kind(Alvi__mass2)
      REAL*8 Alvi__separation
      integer, parameter :: cctki_use_Alvi__separation = kind(Alvi__sepa
     &ration)
      REAL*8 Bertotti__Lambda
      integer, parameter :: cctki_use_Bertotti__Lambda = kind(Bertotti__
     &Lambda)
      REAL*8 Bianchi_I__scale
      integer, parameter :: cctki_use_Bianchi_I__scale = kind(Bianchi_I_
     &_scale)
      REAL*8 Goedel__scale
      integer, parameter :: cctki_use_Goedel__scale = kind(Goedel__scale
     &)
      REAL*8 Gowdy_wave__amplitude
      integer, parameter :: cctki_use_Gowdy_wave__amplitude = kind(Gowdy
     &_wave__amplitude)
      REAL*8 Kasner_generalized__p1
      integer, parameter :: cctki_use_Kasner_generalized__p1 = kind(Kasn
     &er_generalized__p1)
      REAL*8 Kasner_generalized__p2
      integer, parameter :: cctki_use_Kasner_generalized__p2 = kind(Kasn
     &er_generalized__p2)
      REAL*8 Kasner_like__q
      integer, parameter :: cctki_use_Kasner_like__q = kind(Kasner_like_
     &_q)
      REAL*8 Kerr_BoyerLindquist__mass
      integer, parameter :: cctki_use_Kerr_BoyerLindquist__mass = kind(K
     &err_BoyerLindquist__mass)
      REAL*8 Kerr_BoyerLindquist__spin
      integer, parameter :: cctki_use_Kerr_BoyerLindquist__spin = kind(K
     &err_BoyerLindquist__spin)
      REAL*8 Kerr_KerrSchild__boost_v
      integer, parameter :: cctki_use_Kerr_KerrSchild__boost_v = kind(Ke
     &rr_KerrSchild__boost_v)
      REAL*8 Kerr_KerrSchild__epsilon
      integer, parameter :: cctki_use_Kerr_KerrSchild__epsilon = kind(Ke
     &rr_KerrSchild__epsilon)
      REAL*8 Kerr_KerrSchild__mass
      integer, parameter :: cctki_use_Kerr_KerrSchild__mass = kind(Kerr_
     &KerrSchild__mass)
      REAL*8 Kerr_KerrSchild__spin
      integer, parameter :: cctki_use_Kerr_KerrSchild__spin = kind(Kerr_
     &KerrSchild__spin)
      REAL*8 Kerr_KerrSchild__t
      integer, parameter :: cctki_use_Kerr_KerrSchild__t = kind(Kerr_Ker
     &rSchild__t)
      REAL*8 Kerr_KerrSchild__x
      integer, parameter :: cctki_use_Kerr_KerrSchild__x = kind(Kerr_Ker
     &rSchild__x)
      REAL*8 Kerr_KerrSchild__y
      integer, parameter :: cctki_use_Kerr_KerrSchild__y = kind(Kerr_Ker
     &rSchild__y)
      REAL*8 Kerr_KerrSchild__z
      integer, parameter :: cctki_use_Kerr_KerrSchild__z = kind(Kerr_Ker
     &rSchild__z)
      REAL*8 Lemaitre__Lambda
      integer, parameter :: cctki_use_Lemaitre__Lambda = kind(Lemaitre__
     &Lambda)
      REAL*8 Lemaitre__R0
      integer, parameter :: cctki_use_Lemaitre__R0 = kind(Lemaitre__R0)
      REAL*8 Lemaitre__epsilon0
      integer, parameter :: cctki_use_Lemaitre__epsilon0 = kind(Lemaitre
     &__epsilon0)
      REAL*8 Lemaitre__kappa
      integer, parameter :: cctki_use_Lemaitre__kappa = kind(Lemaitre__k
     &appa)
      REAL*8 Minkowski_conf_wave__amplitude
      integer, parameter :: cctki_use_Minkowski_conf_wave__amplitude = k
     &ind(Minkowski_conf_wave__amplitude)
      REAL*8 Minkowski_conf_wave__wavelength
      integer, parameter :: cctki_use_Minkowski_conf_wave__wavelength = 
     &kind(Minkowski_conf_wave__wavelength)
      REAL*8 Minkowski_funny__amplitude
      integer, parameter :: cctki_use_Minkowski_funny__amplitude = kind(
     &Minkowski_funny__amplitude)
      REAL*8 Minkowski_funny__sigma
      integer, parameter :: cctki_use_Minkowski_funny__sigma = kind(Mink
     &owski_funny__sigma)
      REAL*8 Minkowski_gauge_wave__amplitude
      integer, parameter :: cctki_use_Minkowski_gauge_wave__amplitude = 
     &kind(Minkowski_gauge_wave__amplitude)
      REAL*8 Minkowski_gauge_wave__lambda
      integer, parameter :: cctki_use_Minkowski_gauge_wave__lambda = kin
     &d(Minkowski_gauge_wave__lambda)
      REAL*8 Minkowski_gauge_wave__omega
      integer, parameter :: cctki_use_Minkowski_gauge_wave__omega = kind
     &(Minkowski_gauge_wave__omega)
      REAL*8 Minkowski_gauge_wave__phase
      integer, parameter :: cctki_use_Minkowski_gauge_wave__phase = kind
     &(Minkowski_gauge_wave__phase)
      REAL*8 Minkowski_shift__amplitude
      integer, parameter :: cctki_use_Minkowski_shift__amplitude = kind(
     &Minkowski_shift__amplitude)
      REAL*8 Minkowski_shift__sigma
      integer, parameter :: cctki_use_Minkowski_shift__sigma = kind(Mink
     &owski_shift__sigma)
      REAL*8 Schwarzschild_BL__epsilon
      integer, parameter :: cctki_use_Schwarzschild_BL__epsilon = kind(S
     &chwarzschild_BL__epsilon)
      REAL*8 Schwarzschild_BL__mass
      integer, parameter :: cctki_use_Schwarzschild_BL__mass = kind(Schw
     &arzschild_BL__mass)
      REAL*8 Schwarzschild_EF__epsilon
      integer, parameter :: cctki_use_Schwarzschild_EF__epsilon = kind(S
     &chwarzschild_EF__epsilon)
      REAL*8 Schwarzschild_EF__mass
      integer, parameter :: cctki_use_Schwarzschild_EF__mass = kind(Schw
     &arzschild_EF__mass)
      REAL*8 Schwarzschild_Lemaitre__Lambda
      integer, parameter :: cctki_use_Schwarzschild_Lemaitre__Lambda = k
     &ind(Schwarzschild_Lemaitre__Lambda)
      REAL*8 Schwarzschild_Lemaitre__mass
      integer, parameter :: cctki_use_Schwarzschild_Lemaitre__mass = kin
     &d(Schwarzschild_Lemaitre__mass)
      REAL*8 Schwarzschild_Novikov__epsilon
      integer, parameter :: cctki_use_Schwarzschild_Novikov__epsilon = k
     &ind(Schwarzschild_Novikov__epsilon)
      REAL*8 Schwarzschild_Novikov__mass
      integer, parameter :: cctki_use_Schwarzschild_Novikov__mass = kind
     &(Schwarzschild_Novikov__mass)
      REAL*8 Schwarzschild_PG__epsilon
      integer, parameter :: cctki_use_Schwarzschild_PG__epsilon = kind(S
     &chwarzschild_PG__epsilon)
      REAL*8 Schwarzschild_PG__mass
      integer, parameter :: cctki_use_Schwarzschild_PG__mass = kind(Schw
     &arzschild_PG__mass)
      REAL*8 Thorne_fakebinary__Omega0
      integer, parameter :: cctki_use_Thorne_fakebinary__Omega0 = kind(T
     &horne_fakebinary__Omega0)
      REAL*8 Thorne_fakebinary__epsilon
      integer, parameter :: cctki_use_Thorne_fakebinary__epsilon = kind(
     &Thorne_fakebinary__epsilon)
      REAL*8 Thorne_fakebinary__mass
      integer, parameter :: cctki_use_Thorne_fakebinary__mass = kind(Tho
     &rne_fakebinary__mass)
      REAL*8 Thorne_fakebinary__separation
      integer, parameter :: cctki_use_Thorne_fakebinary__separation = ki
     &nd(Thorne_fakebinary__separation)
      REAL*8 Thorne_fakebinary__smoothing
      integer, parameter :: cctki_use_Thorne_fakebinary__smoothing = kin
     &d(Thorne_fakebinary__smoothing)
      REAL*8 anti_de_Sitter_Lambda__scale
      integer, parameter :: cctki_use_anti_de_Sitter_Lambda__scale = kin
     &d(anti_de_Sitter_Lambda__scale)
      REAL*8 boost_rotation_symmetric__amp
      integer, parameter :: cctki_use_boost_rotation_symmetric__amp = ki
     &nd(boost_rotation_symmetric__amp)
      REAL*8 boost_rotation_symmetric__min_d
      integer, parameter :: cctki_use_boost_rotation_symmetric__min_d = 
     &kind(boost_rotation_symmetric__min_d)
      REAL*8 boost_rotation_symmetric__scale
      integer, parameter :: cctki_use_boost_rotation_symmetric__scale = 
     &kind(boost_rotation_symmetric__scale)
      REAL*8 bowl__center
      integer, parameter :: cctki_use_bowl__center = kind(bowl__center)
      REAL*8 bowl__sigma
      integer, parameter :: cctki_use_bowl__sigma = kind(bowl__sigma)
      REAL*8 bowl__sigma_t
      integer, parameter :: cctki_use_bowl__sigma_t = kind(bowl__sigma_t
     &)
      REAL*8 bowl__strength
      integer, parameter :: cctki_use_bowl__strength = kind(bowl__streng
     &th)
      REAL*8 bowl__t0
      integer, parameter :: cctki_use_bowl__t0 = kind(bowl__t0)
      REAL*8 bowl__x_scale
      integer, parameter :: cctki_use_bowl__x_scale = kind(bowl__x_scale
     &)
      REAL*8 bowl__y_scale
      integer, parameter :: cctki_use_bowl__y_scale = kind(bowl__y_scale
     &)
      REAL*8 bowl__z_scale
      integer, parameter :: cctki_use_bowl__z_scale = kind(bowl__z_scale
     &)
      REAL*8 constant_density_star__mass
      integer, parameter :: cctki_use_constant_density_star__mass = kind
     &(constant_density_star__mass)
      REAL*8 constant_density_star__radius
      integer, parameter :: cctki_use_constant_density_star__radius = ki
     &nd(constant_density_star__radius)
      REAL*8 de_Sitter_Lambda__scale
      integer, parameter :: cctki_use_de_Sitter_Lambda__scale = kind(de_
     &Sitter_Lambda__scale)
      REAL*8 de_Sitter__scale
      integer, parameter :: cctki_use_de_Sitter__scale = kind(de_Sitter_
     &_scale)
      REAL*8 multi_BH__Hubble
      integer, parameter :: cctki_use_multi_BH__Hubble = kind(multi_BH__
     &Hubble)
      REAL*8 multi_BH__mass1
      integer, parameter :: cctki_use_multi_BH__mass1 = kind(multi_BH__m
     &ass1)
      REAL*8 multi_BH__mass2
      integer, parameter :: cctki_use_multi_BH__mass2 = kind(multi_BH__m
     &ass2)
      REAL*8 multi_BH__mass3
      integer, parameter :: cctki_use_multi_BH__mass3 = kind(multi_BH__m
     &ass3)
      REAL*8 multi_BH__mass4
      integer, parameter :: cctki_use_multi_BH__mass4 = kind(multi_BH__m
     &ass4)
      REAL*8 multi_BH__x1
      integer, parameter :: cctki_use_multi_BH__x1 = kind(multi_BH__x1)
      REAL*8 multi_BH__x2
      integer, parameter :: cctki_use_multi_BH__x2 = kind(multi_BH__x2)
      REAL*8 multi_BH__x3
      integer, parameter :: cctki_use_multi_BH__x3 = kind(multi_BH__x3)
      REAL*8 multi_BH__x4
      integer, parameter :: cctki_use_multi_BH__x4 = kind(multi_BH__x4)
      REAL*8 multi_BH__y1
      integer, parameter :: cctki_use_multi_BH__y1 = kind(multi_BH__y1)
      REAL*8 multi_BH__y2
      integer, parameter :: cctki_use_multi_BH__y2 = kind(multi_BH__y2)
      REAL*8 multi_BH__y3
      integer, parameter :: cctki_use_multi_BH__y3 = kind(multi_BH__y3)
      REAL*8 multi_BH__y4
      integer, parameter :: cctki_use_multi_BH__y4 = kind(multi_BH__y4)
      REAL*8 multi_BH__z1
      integer, parameter :: cctki_use_multi_BH__z1 = kind(multi_BH__z1)
      REAL*8 multi_BH__z2
      integer, parameter :: cctki_use_multi_BH__z2 = kind(multi_BH__z2)
      REAL*8 multi_BH__z3
      integer, parameter :: cctki_use_multi_BH__z3 = kind(multi_BH__z3)
      REAL*8 multi_BH__z4
      integer, parameter :: cctki_use_multi_BH__z4 = kind(multi_BH__z4)
      integer*8 Minkowski_gauge_wave__what_fn
      integer, parameter :: cctki_use_Minkowski_gauge_wave__what_fn = ki
     &nd(Minkowski_gauge_wave__what_fn)
      integer*8 Thorne_fakebinary__atype
      integer, parameter :: cctki_use_Thorne_fakebinary__atype = kind(Th
     &orne_fakebinary__atype)
      integer*8 bowl__shape
      integer, parameter :: cctki_use_bowl__shape = kind(bowl__shape)
      integer*8 exact_model
      integer, parameter :: cctki_use_exact_model = kind(exact_model)
      INTEGER*4 Kerr_KerrSchild__parabolic
      integer, parameter :: cctki_use_Kerr_KerrSchild__parabolic = kind(
     &Kerr_KerrSchild__parabolic)
      INTEGER*4 Kerr_KerrSchild__power
      integer, parameter :: cctki_use_Kerr_KerrSchild__power = kind(Kerr
     &_KerrSchild__power)
      INTEGER*4 Minkowski_conf_wave__direction
      integer, parameter :: cctki_use_Minkowski_conf_wave__direction = k
     &ind(Minkowski_conf_wave__direction)
      INTEGER*4 Minkowski_gauge_wave__diagonal
      integer, parameter :: cctki_use_Minkowski_gauge_wave__diagonal = k
     &ind(Minkowski_gauge_wave__diagonal)
      INTEGER*4 Thorne_fakebinary__retarded
      integer, parameter :: cctki_use_Thorne_fakebinary__retarded = kind
     &(Thorne_fakebinary__retarded)
      INTEGER*4 bowl__evolve
      integer, parameter :: cctki_use_bowl__evolve = kind(bowl__evolve)
      INTEGER*4 multi_BH__nBH
      integer, parameter :: cctki_use_multi_BH__nBH = kind(multi_BH__nBH
     &)
      COMMON /exactrest/Alvi__mass1, Alvi__mass2, Alvi__separation, Bert
     &otti__Lambda, Bianchi_I__scale, Goedel__scale, Gowdy_wave__amplitu
     &de, Kasner_generalized__p1, Kasner_generalized__p2, Kasner_like__q
     &, Kerr_BoyerLindquist__mass, Kerr_BoyerLindquist__spin, Kerr_KerrS
     &child__boost_v, Kerr_KerrSchild__epsilon, Kerr_KerrSchild__mass, K
     &err_KerrSchild__spin, Kerr_KerrSchild__t, Kerr_KerrSchild__x, Kerr
     &_KerrSchild__y, Kerr_KerrSchild__z, Lemaitre__Lambda, Lemaitre__R0
     &, Lemaitre__epsilon0, Lemaitre__kappa, Minkowski_conf_wave__amplit
     &ude, Minkowski_conf_wave__wavelength, Minkowski_funny__amplitude, 
     &Minkowski_funny__sigma, Minkowski_gauge_wave__amplitude, Minkowski
     &_gauge_wave__lambda, Minkowski_gauge_wave__omega, Minkowski_gauge_
     &wave__phase, Minkowski_shift__amplitude, Minkowski_shift__sigma, S
     &chwarzschild_BL__epsilon, Schwarzschild_BL__mass, Schwarzschild_EF
     &__epsilon, Schwarzschild_EF__mass, Schwarzschild_Lemaitre__Lambda,
     & Schwarzschild_Lemaitre__mass, Schwarzschild_Novikov__epsilon, Sch
     &warzschild_Novikov__mass, Schwarzschild_PG__epsilon, Schwarzschild
     &_PG__mass, Thorne_fakebinary__Omega0, Thorne_fakebinary__epsilon, 
     &Thorne_fakebinary__mass, Thorne_fakebinary__separation, Thorne_fak
     &ebinary__smoothing, anti_de_Sitter_Lambda__scale, boost_rotation_s
     &ymmetric__amp, boost_rotation_symmetric__min_d, boost_rotation_sym
     &metric__scale, bowl__center, bowl__sigma, bowl__sigma_t, bowl__str
     &ength, bowl__t0, bowl__x_scale, bowl__y_scale, bowl__z_scale, cons
     &tant_density_star__mass, constant_density_star__radius, de_Sitter_
     &Lambda__scale, de_Sitter__scale, multi_BH__Hubble, multi_BH__mass1
     &, multi_BH__mass2, multi_BH__mass3, multi_BH__mass4, multi_BH__x1,
     & multi_BH__x2, multi_BH__x3, multi_BH__x4, multi_BH__y1, multi_BH_
     &_y2, multi_BH__y3, multi_BH__y4, multi_BH__z1, multi_BH__z2, multi
     &_BH__z3, multi_BH__z4, Minkowski_gauge_wave__what_fn, Thorne_fakeb
     &inary__atype, bowl__shape, exact_model, Kerr_KerrSchild__parabolic
     &, Kerr_KerrSchild__power, Minkowski_conf_wave__direction, Minkowsk
     &i_gauge_wave__diagonal, Thorne_fakebinary__retarded, bowl__evolve,
     & multi_BH__nBH
      REAL*8 boost_vx
      integer, parameter :: cctki_use_boost_vx = kind(boost_vx)
      REAL*8 boost_vy
      integer, parameter :: cctki_use_boost_vy = kind(boost_vy)
      REAL*8 boost_vz
      integer, parameter :: cctki_use_boost_vz = kind(boost_vz)
      REAL*8 exact_eps
      integer, parameter :: cctki_use_exact_eps = kind(exact_eps)
      REAL*8 exblend_rout
      integer, parameter :: cctki_use_exblend_rout = kind(exblend_rout)
      REAL*8 exblend_width
      integer, parameter :: cctki_use_exblend_width = kind(exblend_width
     &)
      REAL*8 rotation_euler_phi
      integer, parameter :: cctki_use_rotation_euler_phi = kind(rotation
     &_euler_phi)
      REAL*8 rotation_euler_psi
      integer, parameter :: cctki_use_rotation_euler_psi = kind(rotation
     &_euler_psi)
      REAL*8 rotation_euler_theta
      integer, parameter :: cctki_use_rotation_euler_theta = kind(rotati
     &on_euler_theta)
      REAL*8 shift_add_x
      integer, parameter :: cctki_use_shift_add_x = kind(shift_add_x)
      REAL*8 shift_add_y
      integer, parameter :: cctki_use_shift_add_y = kind(shift_add_y)
      REAL*8 shift_add_z
      integer, parameter :: cctki_use_shift_add_z = kind(shift_add_z)
      REAL*8 slice_Gauss_ampl
      integer, parameter :: cctki_use_slice_Gauss_ampl = kind(slice_Gaus
     &s_ampl)
      REAL*8 slice_Gauss_width
      integer, parameter :: cctki_use_slice_Gauss_width = kind(slice_Gau
     &ss_width)
      integer*8 overwrite_boundary
      integer, parameter :: cctki_use_overwrite_boundary = kind(overwrit
     &e_boundary)
      INTEGER*4 exact_order
      integer, parameter :: cctki_use_exact_order = kind(exact_order)
      INTEGER*4 exblend_Ks
      integer, parameter :: cctki_use_exblend_Ks = kind(exblend_Ks)
      INTEGER*4 exblend_gauge
      integer, parameter :: cctki_use_exblend_gauge = kind(exblend_gauge
     &)
      INTEGER*4 exblend_gs
      integer, parameter :: cctki_use_exblend_gs = kind(exblend_gs)
      COMMON /Exactpriv/boost_vx, boost_vy, boost_vz, exact_eps, exblend
     &_rout, exblend_width, rotation_euler_phi, rotation_euler_psi, rota
     &tion_euler_theta, shift_add_x, shift_add_y, shift_add_z, slice_Gau
     &ss_ampl, slice_Gauss_width, overwrite_boundary, exact_order, exble
     &nd_Ks, exblend_gauge, exblend_gs
      integer*8 CCTKH0
      integer, parameter :: cctki_use_CCTKH0 = kind(CCTKH0)
      integer*8 dtlapse_evolution_method
      integer, parameter :: cctki_use_dtlapse_evolution_method = kind(dt
     &lapse_evolution_method)
      integer*8 dtshift_evolution_method
      integer, parameter :: cctki_use_dtshift_evolution_method = kind(dt
     &shift_evolution_method)
      integer*8 evolution_method
      integer, parameter :: cctki_use_evolution_method = kind(evolution_
     &method)
      integer*8 initial_data
      integer, parameter :: cctki_use_initial_data = kind(initial_data)
      integer*8 initial_dtlapse
      integer, parameter :: cctki_use_initial_dtlapse = kind(initial_dtl
     &apse)
      integer*8 initial_dtshift
      integer, parameter :: cctki_use_initial_dtshift = kind(initial_dts
     &hift)
      integer*8 initial_lapse
      integer, parameter :: cctki_use_initial_lapse = kind(initial_lapse
     &)
      integer*8 initial_shift
      integer, parameter :: cctki_use_initial_shift = kind(initial_shift
     &)
      integer*8 lapse_evolution_method
      integer, parameter :: cctki_use_lapse_evolution_method = kind(laps
     &e_evolution_method)
      integer*8 CCTKH1
      integer, parameter :: cctki_use_CCTKH1 = kind(CCTKH1)
      integer*8 CCTKH3
      integer, parameter :: cctki_use_CCTKH3 = kind(CCTKH3)
      integer*8 metric_type
      integer, parameter :: cctki_use_metric_type = kind(metric_type)
      integer*8 shift_evolution_method
      integer, parameter :: cctki_use_shift_evolution_method = kind(shif
     &t_evolution_method)
      integer*8 CCTKH5
      integer, parameter :: cctki_use_CCTKH5 = kind(CCTKH5)
      INTEGER*4 CCTKH2
      integer, parameter :: cctki_use_CCTKH2 = kind(CCTKH2)
      INTEGER*4 CCTKH4
      integer, parameter :: cctki_use_CCTKH4 = kind(CCTKH4)
      INTEGER*4 CCTKH6
      integer, parameter :: cctki_use_CCTKH6 = kind(CCTKH6)
      COMMON /ADMBASErest/CCTKH0, dtlapse_evolution_method, dtshift_evol
     &ution_method, evolution_method, initial_data, initial_dtlapse, ini
     &tial_dtshift, initial_lapse, initial_shift, lapse_evolution_method
     &, CCTKH1, CCTKH3, metric_type, shift_evolution_method, CCTKH5, CCT
     &KH2, CCTKH4, CCTKH6
      integer*8 conformal_storage
      integer, parameter :: cctki_use_conformal_storage = kind(conformal
     &_storage)
      COMMON /STATICCONFORMALrest/conformal_storage
      
      INTEGER*4
     $     decoded_exact_model,
     $     dir
      REAL*8
     $     x, y, z, t,
     $     dgdtt, dgdtx, dgdty, dgdtz,
     $     dgdxx, dgdyy, dgdzz, dgdxy, dgdyz, dgdxz,
     $     dgutt, dgutx, dguty, dgutz,
     $     dguxx, dguyy, dguzz, dguxy, dguyz, dguxz,
     $     dpsi,
     $     gdtt_p, gdtx_p, gdty_p, gdtz_p, 
     $     gdxx_p, gdyy_p, gdzz_p, gdxy_p, gdyz_p, gdxz_p,
     $     gutt_p, gutx_p, guty_p, gutz_p, 
     $     guxx_p, guyy_p, guzz_p, guxy_p, guyz_p, guxz_p,
     $     psi_p,
     $     gdtt_m, gdtx_m, gdty_m, gdtz_m, 
     $     gdxx_m, gdyy_m, gdzz_m, gdxy_m, gdyz_m, gdxz_m,
     $     gutt_m, gutx_m, guty_m, gutz_m, 
     $     guxx_m, guyy_m, guzz_m, guxy_m, guyz_m, guxz_m,
     $     psi_m,
     $     gdtt_p_p, gdtx_p_p, gdty_p_p, gdtz_p_p, 
     $     gdxx_p_p, gdyy_p_p, gdzz_p_p, gdxy_p_p, gdyz_p_p, gdxz_p_p,
     $     gutt_p_p, gutx_p_p, guty_p_p, gutz_p_p, 
     $     guxx_p_p, guyy_p_p, guzz_p_p, guxy_p_p, guyz_p_p, guxz_p_p,
     $     psi_p_p,
     $     gdtt_m_m, gdtx_m_m, gdty_m_m, gdtz_m_m, 
     $     gdxx_m_m, gdyy_m_m, gdzz_m_m, gdxy_m_m, gdyz_m_m, gdxz_m_m,
     $     gutt_m_m, gutx_m_m, guty_m_m, gutz_m_m, 
     $     guxx_m_m, guyy_m_m, guzz_m_m, guxy_m_m, guyz_m_m, guxz_m_m,
     $     psi_m_m,
     $     eps,
     $     dx, dy, dz, dt
      
      eps = exact_eps
      
      dx = 0
      dy = 0
      dz = 0
      dt = 0
      if (dir.eq.0) dt = eps
      if (dir.eq.1) dx = eps
      if (dir.eq.2) dy = eps
      if (dir.eq.3) dz = eps
      
      if (exact_order .eq. 2) then
         
         call Exact__metric(
     $        decoded_exact_model,
     $        x-dx, y-dy, z-dz, t-dt,
     $        gdtt_m, gdtx_m, gdty_m, gdtz_m,
     $        gdxx_m, gdyy_m, gdzz_m, gdxy_m, gdyz_m, gdxz_m,
     $        gutt_m, gutx_m, guty_m, gutz_m,
     $        guxx_m, guyy_m, guzz_m, guxy_m, guyz_m, guxz_m,
     $        psi_m)
         call Exact__metric(
     $        decoded_exact_model,
     $        x+dx, y+dy, z+dz, t+dt,
     $        gdtt_p, gdtx_p, gdty_p, gdtz_p,
     $        gdxx_p, gdyy_p, gdzz_p, gdxy_p, gdyz_p, gdxz_p,
     $        gutt_p, gutx_p, guty_p, gutz_p,
     $        guxx_p, guyy_p, guzz_p, guxy_p, guyz_p, guxz_p,
     $        psi_p)
         
         dgdtt = (gdtt_p - gdtt_m) / (2*eps)
         dgdtx = (gdtx_p - gdtx_m) / (2*eps)
         dgdty = (gdty_p - gdty_m) / (2*eps)
         dgdtz = (gdtz_p - gdtz_m) / (2*eps)
         dgdxx = (gdxx_p - gdxx_m) / (2*eps)
         dgdyy = (gdyy_p - gdyy_m) / (2*eps)
         dgdzz = (gdzz_p - gdzz_m) / (2*eps)
         dgdxy = (gdxy_p - gdxy_m) / (2*eps)
         dgdyz = (gdyz_p - gdyz_m) / (2*eps)
         dgdxz = (gdxz_p - gdxz_m) / (2*eps)
         dgutt = (gutt_p - gutt_m) / (2*eps)
         dgutx = (gutx_p - gutx_m) / (2*eps)
         dguty = (guty_p - guty_m) / (2*eps)
         dgutz = (gutz_p - gutz_m) / (2*eps)
         dguxx = (guxx_p - guxx_m) / (2*eps)
         dguyy = (guyy_p - guyy_m) / (2*eps)
         dguzz = (guzz_p - guzz_m) / (2*eps)
         dguxy = (guxy_p - guxy_m) / (2*eps)
         dguyz = (guyz_p - guyz_m) / (2*eps)
         dguxz = (guxz_p - guxz_m) / (2*eps)
         dpsi  = (psi_p  - psi_m ) / (2*eps)
         
      else if (exact_order .eq. 4) then
         
         call Exact__metric(
     $        decoded_exact_model,
     $        x-2*dx, y-2*dy, z-2*dz, t-2*dt,
     $        gdtt_m_m, gdtx_m_m, gdty_m_m, gdtz_m_m,
     $        gdxx_m_m, gdyy_m_m, gdzz_m_m, gdxy_m_m, gdyz_m_m, gdxz_m_m
     &,
     $        gutt_m_m, gutx_m_m, guty_m_m, gutz_m_m,
     $        guxx_m_m, guyy_m_m, guzz_m_m, guxy_m_m, guyz_m_m, guxz_m_m
     &,
     $        psi_m_m)
         call Exact__metric(
     $        decoded_exact_model,
     $        x-dx, y-dy, z-dz, t-dt,
     $        gdtt_m, gdtx_m, gdty_m, gdtz_m,
     $        gdxx_m, gdyy_m, gdzz_m, gdxy_m, gdyz_m, gdxz_m,
     $        gutt_m, gutx_m, guty_m, gutz_m,
     $        guxx_m, guyy_m, guzz_m, guxy_m, guyz_m, guxz_m,
     $        psi_m)
         call Exact__metric(
     $        decoded_exact_model,
     $        x+dx, y+dy, z+dz, t+dt,
     $        gdtt_p, gdtx_p, gdty_p, gdtz_p,
     $        gdxx_p, gdyy_p, gdzz_p, gdxy_p, gdyz_p, gdxz_p,
     $        gutt_p, gutx_p, guty_p, gutz_p,
     $        guxx_p, guyy_p, guzz_p, guxy_p, guyz_p, guxz_p,
     $        psi_p)
         call Exact__metric(
     $        decoded_exact_model,
     $        x+2*dx, y+2*dy, z+2*dz, t+2*dt,
     $        gdtt_p_p, gdtx_p_p, gdty_p_p, gdtz_p_p,
     $        gdxx_p_p, gdyy_p_p, gdzz_p_p, gdxy_p_p, gdyz_p_p, gdxz_p_p
     &,
     $        gutt_p_p, gutx_p_p, guty_p_p, gutz_p_p,
     $        guxx_p_p, guyy_p_p, guzz_p_p, guxy_p_p, guyz_p_p, guxz_p_p
     &,
     $        psi_p_p)
         
         dgdtt = (- gdtt_p_p + 8*gdtt_p - 8*gdtt_m + gdtt_m_m) / (12*eps
     &)
         dgdtx = (- gdtx_p_p + 8*gdtx_p - 8*gdtx_m + gdtx_m_m) / (12*eps
     &)
         dgdty = (- gdty_p_p + 8*gdty_p - 8*gdty_m + gdty_m_m) / (12*eps
     &)
         dgdtz = (- gdtz_p_p + 8*gdtz_p - 8*gdtz_m + gdtz_m_m) / (12*eps
     &)
         dgdxx = (- gdxx_p_p + 8*gdxx_p - 8*gdxx_m + gdxx_m_m) / (12*eps
     &)
         dgdyy = (- gdyy_p_p + 8*gdyy_p - 8*gdyy_m + gdyy_m_m) / (12*eps
     &)
         dgdzz = (- gdzz_p_p + 8*gdzz_p - 8*gdzz_m + gdzz_m_m) / (12*eps
     &)
         dgdxy = (- gdxy_p_p + 8*gdxy_p - 8*gdxy_m + gdxy_m_m) / (12*eps
     &)
         dgdyz = (- gdyz_p_p + 8*gdyz_p - 8*gdyz_m + gdyz_m_m) / (12*eps
     &)
         dgdxz = (- gdxz_p_p + 8*gdxz_p - 8*gdxz_m + gdxz_m_m) / (12*eps
     &)
         dgutt = (- gutt_p_p + 8*gutt_p - 8*gutt_m + gutt_m_m) / (12*eps
     &)
         dgutx = (- gutx_p_p + 8*gutx_p - 8*gutx_m + gutx_m_m) / (12*eps
     &)
         dguty = (- guty_p_p + 8*guty_p - 8*guty_m + guty_m_m) / (12*eps
     &)
         dgutz = (- gutz_p_p + 8*gutz_p - 8*gutz_m + gutz_m_m) / (12*eps
     &)
         dguxx = (- guxx_p_p + 8*guxx_p - 8*guxx_m + guxx_m_m) / (12*eps
     &)
         dguyy = (- guyy_p_p + 8*guyy_p - 8*guyy_m + guyy_m_m) / (12*eps
     &)
         dguzz = (- guzz_p_p + 8*guzz_p - 8*guzz_m + guzz_m_m) / (12*eps
     &)
         dguxy = (- guxy_p_p + 8*guxy_p - 8*guxy_m + guxy_m_m) / (12*eps
     &)
         dguyz = (- guyz_p_p + 8*guyz_p - 8*guyz_m + guyz_m_m) / (12*eps
     &)
         dguxz = (- guxz_p_p + 8*guxz_p - 8*guxz_m + guxz_m_m) / (12*eps
     &)
         dpsi  = (- psi_p_p  + 8*psi_p  - 8*psi_m  + psi_m_m ) / (12*eps
     &)
         
      else
         call CCTK_Warn(0,445,"Bona_Masso_data.F","Exact", "internal err
     &or")
      end if
      
      end



      subroutine Exact__metric_deriv2(
     $     decoded_exact_model,
     $     dir1, dir2,
     $     x, y, z, t,
     $     dgdtt, dgdtx, dgdty, dgdtz,
     $     dgdxx, dgdyy, dgdzz, dgdxy, dgdyz, dgdxz,
     $     dgutt, dgutx, dguty, dgutz,
     $     dguxx, dguyy, dguzz, dguxy, dguyz, dguxz,
     $     dpsi)
      
      implicit none
      REAL*8 Alvi__mass1
      integer, parameter :: cctki_use_Alvi__mass1 = kind(Alvi__mass1)
      REAL*8 Alvi__mass2
      integer, parameter :: cctki_use_Alvi__mass2 = kind(Alvi__mass2)
      REAL*8 Alvi__separation
      integer, parameter :: cctki_use_Alvi__separation = kind(Alvi__sepa
     &ration)
      REAL*8 Bertotti__Lambda
      integer, parameter :: cctki_use_Bertotti__Lambda = kind(Bertotti__
     &Lambda)
      REAL*8 Bianchi_I__scale
      integer, parameter :: cctki_use_Bianchi_I__scale = kind(Bianchi_I_
     &_scale)
      REAL*8 Goedel__scale
      integer, parameter :: cctki_use_Goedel__scale = kind(Goedel__scale
     &)
      REAL*8 Gowdy_wave__amplitude
      integer, parameter :: cctki_use_Gowdy_wave__amplitude = kind(Gowdy
     &_wave__amplitude)
      REAL*8 Kasner_generalized__p1
      integer, parameter :: cctki_use_Kasner_generalized__p1 = kind(Kasn
     &er_generalized__p1)
      REAL*8 Kasner_generalized__p2
      integer, parameter :: cctki_use_Kasner_generalized__p2 = kind(Kasn
     &er_generalized__p2)
      REAL*8 Kasner_like__q
      integer, parameter :: cctki_use_Kasner_like__q = kind(Kasner_like_
     &_q)
      REAL*8 Kerr_BoyerLindquist__mass
      integer, parameter :: cctki_use_Kerr_BoyerLindquist__mass = kind(K
     &err_BoyerLindquist__mass)
      REAL*8 Kerr_BoyerLindquist__spin
      integer, parameter :: cctki_use_Kerr_BoyerLindquist__spin = kind(K
     &err_BoyerLindquist__spin)
      REAL*8 Kerr_KerrSchild__boost_v
      integer, parameter :: cctki_use_Kerr_KerrSchild__boost_v = kind(Ke
     &rr_KerrSchild__boost_v)
      REAL*8 Kerr_KerrSchild__epsilon
      integer, parameter :: cctki_use_Kerr_KerrSchild__epsilon = kind(Ke
     &rr_KerrSchild__epsilon)
      REAL*8 Kerr_KerrSchild__mass
      integer, parameter :: cctki_use_Kerr_KerrSchild__mass = kind(Kerr_
     &KerrSchild__mass)
      REAL*8 Kerr_KerrSchild__spin
      integer, parameter :: cctki_use_Kerr_KerrSchild__spin = kind(Kerr_
     &KerrSchild__spin)
      REAL*8 Kerr_KerrSchild__t
      integer, parameter :: cctki_use_Kerr_KerrSchild__t = kind(Kerr_Ker
     &rSchild__t)
      REAL*8 Kerr_KerrSchild__x
      integer, parameter :: cctki_use_Kerr_KerrSchild__x = kind(Kerr_Ker
     &rSchild__x)
      REAL*8 Kerr_KerrSchild__y
      integer, parameter :: cctki_use_Kerr_KerrSchild__y = kind(Kerr_Ker
     &rSchild__y)
      REAL*8 Kerr_KerrSchild__z
      integer, parameter :: cctki_use_Kerr_KerrSchild__z = kind(Kerr_Ker
     &rSchild__z)
      REAL*8 Lemaitre__Lambda
      integer, parameter :: cctki_use_Lemaitre__Lambda = kind(Lemaitre__
     &Lambda)
      REAL*8 Lemaitre__R0
      integer, parameter :: cctki_use_Lemaitre__R0 = kind(Lemaitre__R0)
      REAL*8 Lemaitre__epsilon0
      integer, parameter :: cctki_use_Lemaitre__epsilon0 = kind(Lemaitre
     &__epsilon0)
      REAL*8 Lemaitre__kappa
      integer, parameter :: cctki_use_Lemaitre__kappa = kind(Lemaitre__k
     &appa)
      REAL*8 Minkowski_conf_wave__amplitude
      integer, parameter :: cctki_use_Minkowski_conf_wave__amplitude = k
     &ind(Minkowski_conf_wave__amplitude)
      REAL*8 Minkowski_conf_wave__wavelength
      integer, parameter :: cctki_use_Minkowski_conf_wave__wavelength = 
     &kind(Minkowski_conf_wave__wavelength)
      REAL*8 Minkowski_funny__amplitude
      integer, parameter :: cctki_use_Minkowski_funny__amplitude = kind(
     &Minkowski_funny__amplitude)
      REAL*8 Minkowski_funny__sigma
      integer, parameter :: cctki_use_Minkowski_funny__sigma = kind(Mink
     &owski_funny__sigma)
      REAL*8 Minkowski_gauge_wave__amplitude
      integer, parameter :: cctki_use_Minkowski_gauge_wave__amplitude = 
     &kind(Minkowski_gauge_wave__amplitude)
      REAL*8 Minkowski_gauge_wave__lambda
      integer, parameter :: cctki_use_Minkowski_gauge_wave__lambda = kin
     &d(Minkowski_gauge_wave__lambda)
      REAL*8 Minkowski_gauge_wave__omega
      integer, parameter :: cctki_use_Minkowski_gauge_wave__omega = kind
     &(Minkowski_gauge_wave__omega)
      REAL*8 Minkowski_gauge_wave__phase
      integer, parameter :: cctki_use_Minkowski_gauge_wave__phase = kind
     &(Minkowski_gauge_wave__phase)
      REAL*8 Minkowski_shift__amplitude
      integer, parameter :: cctki_use_Minkowski_shift__amplitude = kind(
     &Minkowski_shift__amplitude)
      REAL*8 Minkowski_shift__sigma
      integer, parameter :: cctki_use_Minkowski_shift__sigma = kind(Mink
     &owski_shift__sigma)
      REAL*8 Schwarzschild_BL__epsilon
      integer, parameter :: cctki_use_Schwarzschild_BL__epsilon = kind(S
     &chwarzschild_BL__epsilon)
      REAL*8 Schwarzschild_BL__mass
      integer, parameter :: cctki_use_Schwarzschild_BL__mass = kind(Schw
     &arzschild_BL__mass)
      REAL*8 Schwarzschild_EF__epsilon
      integer, parameter :: cctki_use_Schwarzschild_EF__epsilon = kind(S
     &chwarzschild_EF__epsilon)
      REAL*8 Schwarzschild_EF__mass
      integer, parameter :: cctki_use_Schwarzschild_EF__mass = kind(Schw
     &arzschild_EF__mass)
      REAL*8 Schwarzschild_Lemaitre__Lambda
      integer, parameter :: cctki_use_Schwarzschild_Lemaitre__Lambda = k
     &ind(Schwarzschild_Lemaitre__Lambda)
      REAL*8 Schwarzschild_Lemaitre__mass
      integer, parameter :: cctki_use_Schwarzschild_Lemaitre__mass = kin
     &d(Schwarzschild_Lemaitre__mass)
      REAL*8 Schwarzschild_Novikov__epsilon
      integer, parameter :: cctki_use_Schwarzschild_Novikov__epsilon = k
     &ind(Schwarzschild_Novikov__epsilon)
      REAL*8 Schwarzschild_Novikov__mass
      integer, parameter :: cctki_use_Schwarzschild_Novikov__mass = kind
     &(Schwarzschild_Novikov__mass)
      REAL*8 Schwarzschild_PG__epsilon
      integer, parameter :: cctki_use_Schwarzschild_PG__epsilon = kind(S
     &chwarzschild_PG__epsilon)
      REAL*8 Schwarzschild_PG__mass
      integer, parameter :: cctki_use_Schwarzschild_PG__mass = kind(Schw
     &arzschild_PG__mass)
      REAL*8 Thorne_fakebinary__Omega0
      integer, parameter :: cctki_use_Thorne_fakebinary__Omega0 = kind(T
     &horne_fakebinary__Omega0)
      REAL*8 Thorne_fakebinary__epsilon
      integer, parameter :: cctki_use_Thorne_fakebinary__epsilon = kind(
     &Thorne_fakebinary__epsilon)
      REAL*8 Thorne_fakebinary__mass
      integer, parameter :: cctki_use_Thorne_fakebinary__mass = kind(Tho
     &rne_fakebinary__mass)
      REAL*8 Thorne_fakebinary__separation
      integer, parameter :: cctki_use_Thorne_fakebinary__separation = ki
     &nd(Thorne_fakebinary__separation)
      REAL*8 Thorne_fakebinary__smoothing
      integer, parameter :: cctki_use_Thorne_fakebinary__smoothing = kin
     &d(Thorne_fakebinary__smoothing)
      REAL*8 anti_de_Sitter_Lambda__scale
      integer, parameter :: cctki_use_anti_de_Sitter_Lambda__scale = kin
     &d(anti_de_Sitter_Lambda__scale)
      REAL*8 boost_rotation_symmetric__amp
      integer, parameter :: cctki_use_boost_rotation_symmetric__amp = ki
     &nd(boost_rotation_symmetric__amp)
      REAL*8 boost_rotation_symmetric__min_d
      integer, parameter :: cctki_use_boost_rotation_symmetric__min_d = 
     &kind(boost_rotation_symmetric__min_d)
      REAL*8 boost_rotation_symmetric__scale
      integer, parameter :: cctki_use_boost_rotation_symmetric__scale = 
     &kind(boost_rotation_symmetric__scale)
      REAL*8 bowl__center
      integer, parameter :: cctki_use_bowl__center = kind(bowl__center)
      REAL*8 bowl__sigma
      integer, parameter :: cctki_use_bowl__sigma = kind(bowl__sigma)
      REAL*8 bowl__sigma_t
      integer, parameter :: cctki_use_bowl__sigma_t = kind(bowl__sigma_t
     &)
      REAL*8 bowl__strength
      integer, parameter :: cctki_use_bowl__strength = kind(bowl__streng
     &th)
      REAL*8 bowl__t0
      integer, parameter :: cctki_use_bowl__t0 = kind(bowl__t0)
      REAL*8 bowl__x_scale
      integer, parameter :: cctki_use_bowl__x_scale = kind(bowl__x_scale
     &)
      REAL*8 bowl__y_scale
      integer, parameter :: cctki_use_bowl__y_scale = kind(bowl__y_scale
     &)
      REAL*8 bowl__z_scale
      integer, parameter :: cctki_use_bowl__z_scale = kind(bowl__z_scale
     &)
      REAL*8 constant_density_star__mass
      integer, parameter :: cctki_use_constant_density_star__mass = kind
     &(constant_density_star__mass)
      REAL*8 constant_density_star__radius
      integer, parameter :: cctki_use_constant_density_star__radius = ki
     &nd(constant_density_star__radius)
      REAL*8 de_Sitter_Lambda__scale
      integer, parameter :: cctki_use_de_Sitter_Lambda__scale = kind(de_
     &Sitter_Lambda__scale)
      REAL*8 de_Sitter__scale
      integer, parameter :: cctki_use_de_Sitter__scale = kind(de_Sitter_
     &_scale)
      REAL*8 multi_BH__Hubble
      integer, parameter :: cctki_use_multi_BH__Hubble = kind(multi_BH__
     &Hubble)
      REAL*8 multi_BH__mass1
      integer, parameter :: cctki_use_multi_BH__mass1 = kind(multi_BH__m
     &ass1)
      REAL*8 multi_BH__mass2
      integer, parameter :: cctki_use_multi_BH__mass2 = kind(multi_BH__m
     &ass2)
      REAL*8 multi_BH__mass3
      integer, parameter :: cctki_use_multi_BH__mass3 = kind(multi_BH__m
     &ass3)
      REAL*8 multi_BH__mass4
      integer, parameter :: cctki_use_multi_BH__mass4 = kind(multi_BH__m
     &ass4)
      REAL*8 multi_BH__x1
      integer, parameter :: cctki_use_multi_BH__x1 = kind(multi_BH__x1)
      REAL*8 multi_BH__x2
      integer, parameter :: cctki_use_multi_BH__x2 = kind(multi_BH__x2)
      REAL*8 multi_BH__x3
      integer, parameter :: cctki_use_multi_BH__x3 = kind(multi_BH__x3)
      REAL*8 multi_BH__x4
      integer, parameter :: cctki_use_multi_BH__x4 = kind(multi_BH__x4)
      REAL*8 multi_BH__y1
      integer, parameter :: cctki_use_multi_BH__y1 = kind(multi_BH__y1)
      REAL*8 multi_BH__y2
      integer, parameter :: cctki_use_multi_BH__y2 = kind(multi_BH__y2)
      REAL*8 multi_BH__y3
      integer, parameter :: cctki_use_multi_BH__y3 = kind(multi_BH__y3)
      REAL*8 multi_BH__y4
      integer, parameter :: cctki_use_multi_BH__y4 = kind(multi_BH__y4)
      REAL*8 multi_BH__z1
      integer, parameter :: cctki_use_multi_BH__z1 = kind(multi_BH__z1)
      REAL*8 multi_BH__z2
      integer, parameter :: cctki_use_multi_BH__z2 = kind(multi_BH__z2)
      REAL*8 multi_BH__z3
      integer, parameter :: cctki_use_multi_BH__z3 = kind(multi_BH__z3)
      REAL*8 multi_BH__z4
      integer, parameter :: cctki_use_multi_BH__z4 = kind(multi_BH__z4)
      integer*8 Minkowski_gauge_wave__what_fn
      integer, parameter :: cctki_use_Minkowski_gauge_wave__what_fn = ki
     &nd(Minkowski_gauge_wave__what_fn)
      integer*8 Thorne_fakebinary__atype
      integer, parameter :: cctki_use_Thorne_fakebinary__atype = kind(Th
     &orne_fakebinary__atype)
      integer*8 bowl__shape
      integer, parameter :: cctki_use_bowl__shape = kind(bowl__shape)
      integer*8 exact_model
      integer, parameter :: cctki_use_exact_model = kind(exact_model)
      INTEGER*4 Kerr_KerrSchild__parabolic
      integer, parameter :: cctki_use_Kerr_KerrSchild__parabolic = kind(
     &Kerr_KerrSchild__parabolic)
      INTEGER*4 Kerr_KerrSchild__power
      integer, parameter :: cctki_use_Kerr_KerrSchild__power = kind(Kerr
     &_KerrSchild__power)
      INTEGER*4 Minkowski_conf_wave__direction
      integer, parameter :: cctki_use_Minkowski_conf_wave__direction = k
     &ind(Minkowski_conf_wave__direction)
      INTEGER*4 Minkowski_gauge_wave__diagonal
      integer, parameter :: cctki_use_Minkowski_gauge_wave__diagonal = k
     &ind(Minkowski_gauge_wave__diagonal)
      INTEGER*4 Thorne_fakebinary__retarded
      integer, parameter :: cctki_use_Thorne_fakebinary__retarded = kind
     &(Thorne_fakebinary__retarded)
      INTEGER*4 bowl__evolve
      integer, parameter :: cctki_use_bowl__evolve = kind(bowl__evolve)
      INTEGER*4 multi_BH__nBH
      integer, parameter :: cctki_use_multi_BH__nBH = kind(multi_BH__nBH
     &)
      COMMON /exactrest/Alvi__mass1, Alvi__mass2, Alvi__separation, Bert
     &otti__Lambda, Bianchi_I__scale, Goedel__scale, Gowdy_wave__amplitu
     &de, Kasner_generalized__p1, Kasner_generalized__p2, Kasner_like__q
     &, Kerr_BoyerLindquist__mass, Kerr_BoyerLindquist__spin, Kerr_KerrS
     &child__boost_v, Kerr_KerrSchild__epsilon, Kerr_KerrSchild__mass, K
     &err_KerrSchild__spin, Kerr_KerrSchild__t, Kerr_KerrSchild__x, Kerr
     &_KerrSchild__y, Kerr_KerrSchild__z, Lemaitre__Lambda, Lemaitre__R0
     &, Lemaitre__epsilon0, Lemaitre__kappa, Minkowski_conf_wave__amplit
     &ude, Minkowski_conf_wave__wavelength, Minkowski_funny__amplitude, 
     &Minkowski_funny__sigma, Minkowski_gauge_wave__amplitude, Minkowski
     &_gauge_wave__lambda, Minkowski_gauge_wave__omega, Minkowski_gauge_
     &wave__phase, Minkowski_shift__amplitude, Minkowski_shift__sigma, S
     &chwarzschild_BL__epsilon, Schwarzschild_BL__mass, Schwarzschild_EF
     &__epsilon, Schwarzschild_EF__mass, Schwarzschild_Lemaitre__Lambda,
     & Schwarzschild_Lemaitre__mass, Schwarzschild_Novikov__epsilon, Sch
     &warzschild_Novikov__mass, Schwarzschild_PG__epsilon, Schwarzschild
     &_PG__mass, Thorne_fakebinary__Omega0, Thorne_fakebinary__epsilon, 
     &Thorne_fakebinary__mass, Thorne_fakebinary__separation, Thorne_fak
     &ebinary__smoothing, anti_de_Sitter_Lambda__scale, boost_rotation_s
     &ymmetric__amp, boost_rotation_symmetric__min_d, boost_rotation_sym
     &metric__scale, bowl__center, bowl__sigma, bowl__sigma_t, bowl__str
     &ength, bowl__t0, bowl__x_scale, bowl__y_scale, bowl__z_scale, cons
     &tant_density_star__mass, constant_density_star__radius, de_Sitter_
     &Lambda__scale, de_Sitter__scale, multi_BH__Hubble, multi_BH__mass1
     &, multi_BH__mass2, multi_BH__mass3, multi_BH__mass4, multi_BH__x1,
     & multi_BH__x2, multi_BH__x3, multi_BH__x4, multi_BH__y1, multi_BH_
     &_y2, multi_BH__y3, multi_BH__y4, multi_BH__z1, multi_BH__z2, multi
     &_BH__z3, multi_BH__z4, Minkowski_gauge_wave__what_fn, Thorne_fakeb
     &inary__atype, bowl__shape, exact_model, Kerr_KerrSchild__parabolic
     &, Kerr_KerrSchild__power, Minkowski_conf_wave__direction, Minkowsk
     &i_gauge_wave__diagonal, Thorne_fakebinary__retarded, bowl__evolve,
     & multi_BH__nBH
      REAL*8 boost_vx
      integer, parameter :: cctki_use_boost_vx = kind(boost_vx)
      REAL*8 boost_vy
      integer, parameter :: cctki_use_boost_vy = kind(boost_vy)
      REAL*8 boost_vz
      integer, parameter :: cctki_use_boost_vz = kind(boost_vz)
      REAL*8 exact_eps
      integer, parameter :: cctki_use_exact_eps = kind(exact_eps)
      REAL*8 exblend_rout
      integer, parameter :: cctki_use_exblend_rout = kind(exblend_rout)
      REAL*8 exblend_width
      integer, parameter :: cctki_use_exblend_width = kind(exblend_width
     &)
      REAL*8 rotation_euler_phi
      integer, parameter :: cctki_use_rotation_euler_phi = kind(rotation
     &_euler_phi)
      REAL*8 rotation_euler_psi
      integer, parameter :: cctki_use_rotation_euler_psi = kind(rotation
     &_euler_psi)
      REAL*8 rotation_euler_theta
      integer, parameter :: cctki_use_rotation_euler_theta = kind(rotati
     &on_euler_theta)
      REAL*8 shift_add_x
      integer, parameter :: cctki_use_shift_add_x = kind(shift_add_x)
      REAL*8 shift_add_y
      integer, parameter :: cctki_use_shift_add_y = kind(shift_add_y)
      REAL*8 shift_add_z
      integer, parameter :: cctki_use_shift_add_z = kind(shift_add_z)
      REAL*8 slice_Gauss_ampl
      integer, parameter :: cctki_use_slice_Gauss_ampl = kind(slice_Gaus
     &s_ampl)
      REAL*8 slice_Gauss_width
      integer, parameter :: cctki_use_slice_Gauss_width = kind(slice_Gau
     &ss_width)
      integer*8 overwrite_boundary
      integer, parameter :: cctki_use_overwrite_boundary = kind(overwrit
     &e_boundary)
      INTEGER*4 exact_order
      integer, parameter :: cctki_use_exact_order = kind(exact_order)
      INTEGER*4 exblend_Ks
      integer, parameter :: cctki_use_exblend_Ks = kind(exblend_Ks)
      INTEGER*4 exblend_gauge
      integer, parameter :: cctki_use_exblend_gauge = kind(exblend_gauge
     &)
      INTEGER*4 exblend_gs
      integer, parameter :: cctki_use_exblend_gs = kind(exblend_gs)
      COMMON /Exactpriv/boost_vx, boost_vy, boost_vz, exact_eps, exblend
     &_rout, exblend_width, rotation_euler_phi, rotation_euler_psi, rota
     &tion_euler_theta, shift_add_x, shift_add_y, shift_add_z, slice_Gau
     &ss_ampl, slice_Gauss_width, overwrite_boundary, exact_order, exble
     &nd_Ks, exblend_gauge, exblend_gs
      integer*8 CCTKH0
      integer, parameter :: cctki_use_CCTKH0 = kind(CCTKH0)
      integer*8 dtlapse_evolution_method
      integer, parameter :: cctki_use_dtlapse_evolution_method = kind(dt
     &lapse_evolution_method)
      integer*8 dtshift_evolution_method
      integer, parameter :: cctki_use_dtshift_evolution_method = kind(dt
     &shift_evolution_method)
      integer*8 evolution_method
      integer, parameter :: cctki_use_evolution_method = kind(evolution_
     &method)
      integer*8 initial_data
      integer, parameter :: cctki_use_initial_data = kind(initial_data)
      integer*8 initial_dtlapse
      integer, parameter :: cctki_use_initial_dtlapse = kind(initial_dtl
     &apse)
      integer*8 initial_dtshift
      integer, parameter :: cctki_use_initial_dtshift = kind(initial_dts
     &hift)
      integer*8 initial_lapse
      integer, parameter :: cctki_use_initial_lapse = kind(initial_lapse
     &)
      integer*8 initial_shift
      integer, parameter :: cctki_use_initial_shift = kind(initial_shift
     &)
      integer*8 lapse_evolution_method
      integer, parameter :: cctki_use_lapse_evolution_method = kind(laps
     &e_evolution_method)
      integer*8 CCTKH1
      integer, parameter :: cctki_use_CCTKH1 = kind(CCTKH1)
      integer*8 CCTKH3
      integer, parameter :: cctki_use_CCTKH3 = kind(CCTKH3)
      integer*8 metric_type
      integer, parameter :: cctki_use_metric_type = kind(metric_type)
      integer*8 shift_evolution_method
      integer, parameter :: cctki_use_shift_evolution_method = kind(shif
     &t_evolution_method)
      integer*8 CCTKH5
      integer, parameter :: cctki_use_CCTKH5 = kind(CCTKH5)
      INTEGER*4 CCTKH2
      integer, parameter :: cctki_use_CCTKH2 = kind(CCTKH2)
      INTEGER*4 CCTKH4
      integer, parameter :: cctki_use_CCTKH4 = kind(CCTKH4)
      INTEGER*4 CCTKH6
      integer, parameter :: cctki_use_CCTKH6 = kind(CCTKH6)
      COMMON /ADMBASErest/CCTKH0, dtlapse_evolution_method, dtshift_evol
     &ution_method, evolution_method, initial_data, initial_dtlapse, ini
     &tial_dtshift, initial_lapse, initial_shift, lapse_evolution_method
     &, CCTKH1, CCTKH3, metric_type, shift_evolution_method, CCTKH5, CCT
     &KH2, CCTKH4, CCTKH6
      integer*8 conformal_storage
      integer, parameter :: cctki_use_conformal_storage = kind(conformal
     &_storage)
      COMMON /STATICCONFORMALrest/conformal_storage
      
      INTEGER*4
     $     decoded_exact_model,
     $     dir1, dir2
      REAL*8
     $     x, y, z, t,
     $     dgdtt, dgdtx, dgdty, dgdtz,
     $     dgdxx, dgdyy, dgdzz, dgdxy, dgdyz, dgdxz,
     $     dgutt, dgutx, dguty, dgutz,
     $     dguxx, dguyy, dguzz, dguxy, dguyz, dguxz,
     $     dpsi,
     $     gdtt_0, gdtx_0, gdty_0, gdtz_0, 
     $     gdxx_0, gdyy_0, gdzz_0, gdxy_0, gdyz_0, gdxz_0,
     $     gutt_0, gutx_0, guty_0, gutz_0, 
     $     guxx_0, guyy_0, guzz_0, guxy_0, guyz_0, guxz_0,
     $     psi_0,
     $     gdtt_p, gdtx_p, gdty_p, gdtz_p, 
     $     gdxx_p, gdyy_p, gdzz_p, gdxy_p, gdyz_p, gdxz_p,
     $     gutt_p, gutx_p, guty_p, gutz_p, 
     $     guxx_p, guyy_p, guzz_p, guxy_p, guyz_p, guxz_p,
     $     psi_p,
     $     gdtt_m, gdtx_m, gdty_m, gdtz_m, 
     $     gdxx_m, gdyy_m, gdzz_m, gdxy_m, gdyz_m, gdxz_m,
     $     gutt_m, gutx_m, guty_m, gutz_m, 
     $     guxx_m, guyy_m, guzz_m, guxy_m, guyz_m, guxz_m,
     $     psi_m,
     $     gdtt_p_p, gdtx_p_p, gdty_p_p, gdtz_p_p, 
     $     gdxx_p_p, gdyy_p_p, gdzz_p_p, gdxy_p_p, gdyz_p_p, gdxz_p_p,
     $     gutt_p_p, gutx_p_p, guty_p_p, gutz_p_p, 
     $     guxx_p_p, guyy_p_p, guzz_p_p, guxy_p_p, guyz_p_p, guxz_p_p,
     $     psi_p_p,
     $     gdtt_m_m, gdtx_m_m, gdty_m_m, gdtz_m_m, 
     $     gdxx_m_m, gdyy_m_m, gdzz_m_m, gdxy_m_m, gdyz_m_m, gdxz_m_m,
     $     gutt_m_m, gutx_m_m, guty_m_m, gutz_m_m, 
     $     guxx_m_m, guyy_m_m, guzz_m_m, guxy_m_m, guyz_m_m, guxz_m_m,
     $     psi_m_m,
     $     eps,
     $     dx, dy, dz, dt
      
      eps = exact_eps
      
      dx = 0
      dy = 0
      dz = 0
      dt = 0
      if (dir1.eq.0) dt = eps
      if (dir1.eq.1) dx = eps
      if (dir1.eq.2) dy = eps
      if (dir1.eq.3) dz = eps
      if (dir1.lt.0 .or. dir1.gt.3) then
         call CCTK_Warn(0,514,"Bona_Masso_data.F","Exact", "internal err
     &or")
      end if
      if (dir2.lt.0 .or. dir2.gt.3) then
         call CCTK_Warn(0,517,"Bona_Masso_data.F","Exact", "internal err
     &or")
      end if
      
      if (exact_order .eq. 2) then
         
         if (dir1 .eq. dir2) then
            
            call Exact__metric(
     $           decoded_exact_model,
     $           x-dx, y-dy, z-dz, t-dt,
     $           gdtt_m, gdtx_m, gdty_m, gdtz_m,
     $           gdxx_m, gdyy_m, gdzz_m, gdxy_m, gdyz_m, gdxz_m,
     $           gutt_m, gutx_m, guty_m, gutz_m,
     $           guxx_m, guyy_m, guzz_m, guxy_m, guyz_m, guxz_m,
     $           psi_m)
            call Exact__metric(
     $           decoded_exact_model,
     $           x, y, z, t,
     $           gdtt_0, gdtx_0, gdty_0, gdtz_0,
     $           gdxx_0, gdyy_0, gdzz_0, gdxy_0, gdyz_0, gdxz_0,
     $           gutt_0, gutx_0, guty_0, gutz_0,
     $           guxx_0, guyy_0, guzz_0, guxy_0, guyz_0, guxz_0,
     $           psi_0)
            call Exact__metric(
     $           decoded_exact_model,
     $           x+dx, y+dy, z+dz, t+dt,
     $           gdtt_p, gdtx_p, gdty_p, gdtz_p,
     $           gdxx_p, gdyy_p, gdzz_p, gdxy_p, gdyz_p, gdxz_p,
     $           gutt_p, gutx_p, guty_p, gutz_p,
     $           guxx_p, guyy_p, guzz_p, guxy_p, guyz_p, guxz_p,
     $           psi_p)
            
            dgdtt = (gdtt_m - 2*gdtt_0 + gdtt_p) / eps**2
            dgdtx = (gdtx_m - 2*gdtx_0 + gdtx_p) / eps**2
            dgdty = (gdty_m - 2*gdty_0 + gdty_p) / eps**2
            dgdtz = (gdtz_m - 2*gdtz_0 + gdtz_p) / eps**2
            dgdxx = (gdxx_m - 2*gdxx_0 + gdxx_p) / eps**2
            dgdyy = (gdyy_m - 2*gdyy_0 + gdyy_p) / eps**2
            dgdzz = (gdzz_m - 2*gdzz_0 + gdzz_p) / eps**2
            dgdxy = (gdxy_m - 2*gdxy_0 + gdxy_p) / eps**2
            dgdyz = (gdyz_m - 2*gdyz_0 + gdyz_p) / eps**2
            dgdxz = (gdxz_m - 2*gdxz_0 + gdxz_p) / eps**2
            dgutt = (gutt_m - 2*gutt_0 + gutt_p) / eps**2
            dgutx = (gutx_m - 2*gutx_0 + gutx_p) / eps**2
            dguty = (guty_m - 2*guty_0 + guty_p) / eps**2
            dgutz = (gutz_m - 2*gutz_0 + gutz_p) / eps**2
            dguxx = (guxx_m - 2*guxx_0 + guxx_p) / eps**2
            dguyy = (guyy_m - 2*guyy_0 + guyy_p) / eps**2
            dguzz = (guzz_m - 2*guzz_0 + guzz_p) / eps**2
            dguxy = (guxy_m - 2*guxy_0 + guxy_p) / eps**2
            dguyz = (guyz_m - 2*guyz_0 + guyz_p) / eps**2
            dguxz = (guxz_m - 2*guxz_0 + guxz_p) / eps**2
            dpsi  = (psi_m  - 2*psi_0  + psi_p ) / eps**2
            
         else
            
            call Exact__metric_deriv(
     $           decoded_exact_model,
     $           dir2,
     $           x-dx, y-dy, z-dz, t-dt,
     $           gdtt_m, gdtx_m, gdty_m, gdtz_m,
     $           gdxx_m, gdyy_m, gdzz_m, gdxy_m, gdyz_m, gdxz_m,
     $           gutt_m, gutx_m, guty_m, gutz_m,
     $           guxx_m, guyy_m, guzz_m, guxy_m, guyz_m, guxz_m,
     $           psi_m)
            call Exact__metric_deriv(
     $           decoded_exact_model,
     $           dir2,
     $           x+dx, y+dy, z+dz, t+dt,
     $           gdtt_p, gdtx_p, gdty_p, gdtz_p,
     $           gdxx_p, gdyy_p, gdzz_p, gdxy_p, gdyz_p, gdxz_p,
     $           gutt_p, gutx_p, guty_p, gutz_p,
     $           guxx_p, guyy_p, guzz_p, guxy_p, guyz_p, guxz_p,
     $           psi_p)
         
            dgdtt = (gdtt_p - gdtt_m) / (2*eps)
            dgdtx = (gdtx_p - gdtx_m) / (2*eps)
            dgdty = (gdty_p - gdty_m) / (2*eps)
            dgdtz = (gdtz_p - gdtz_m) / (2*eps)
            dgdxx = (gdxx_p - gdxx_m) / (2*eps)
            dgdyy = (gdyy_p - gdyy_m) / (2*eps)
            dgdzz = (gdzz_p - gdzz_m) / (2*eps)
            dgdxy = (gdxy_p - gdxy_m) / (2*eps)
            dgdyz = (gdyz_p - gdyz_m) / (2*eps)
            dgdxz = (gdxz_p - gdxz_m) / (2*eps)
            dgutt = (gutt_p - gutt_m) / (2*eps)
            dgutx = (gutx_p - gutx_m) / (2*eps)
            dguty = (guty_p - guty_m) / (2*eps)
            dgutz = (gutz_p - gutz_m) / (2*eps)
            dguxx = (guxx_p - guxx_m) / (2*eps)
            dguyy = (guyy_p - guyy_m) / (2*eps)
            dguzz = (guzz_p - guzz_m) / (2*eps)
            dguxy = (guxy_p - guxy_m) / (2*eps)
            dguyz = (guyz_p - guyz_m) / (2*eps)
            dguxz = (guxz_p - guxz_m) / (2*eps)
            dpsi  = (psi_p  - psi_m ) / (2*eps)
            
         end if
         
      else if (exact_order .eq. 4) then
         
         if (dir1 .eq. dir2) then
            
            call Exact__metric(
     $           decoded_exact_model,
     $           x-2*dx, y-2*dy, z-2*dz, t-2*dt,
     $           gdtt_m_m, gdtx_m_m, gdty_m_m, gdtz_m_m,
     $           gdxx_m_m, gdyy_m_m, gdzz_m_m, gdxy_m_m, gdyz_m_m, gdxz_
     &m_m,
     $           gutt_m_m, gutx_m_m, guty_m_m, gutz_m_m,
     $           guxx_m_m, guyy_m_m, guzz_m_m, guxy_m_m, guyz_m_m, guxz_
     &m_m,
     $           psi_m_m)
            call Exact__metric(
     $           decoded_exact_model,
     $           x-dx, y-dy, z-dz, t-dt,
     $           gdtt_m, gdtx_m, gdty_m, gdtz_m,
     $           gdxx_m, gdyy_m, gdzz_m, gdxy_m, gdyz_m, gdxz_m,
     $           gutt_m, gutx_m, guty_m, gutz_m,
     $           guxx_m, guyy_m, guzz_m, guxy_m, guyz_m, guxz_m,
     $           psi_m)
            call Exact__metric(
     $           decoded_exact_model,
     $           x, y, z, t,
     $           gdtt_0, gdtx_0, gdty_0, gdtz_0,
     $           gdxx_0, gdyy_0, gdzz_0, gdxy_0, gdyz_0, gdxz_0,
     $           gutt_0, gutx_0, guty_0, gutz_0,
     $           guxx_0, guyy_0, guzz_0, guxy_0, guyz_0, guxz_0,
     $           psi_0)
            call Exact__metric(
     $           decoded_exact_model,
     $           x+dx, y+dy, z+dz, t+dt,
     $           gdtt_p, gdtx_p, gdty_p, gdtz_p,
     $           gdxx_p, gdyy_p, gdzz_p, gdxy_p, gdyz_p, gdxz_p,
     $           gutt_p, gutx_p, guty_p, gutz_p,
     $           guxx_p, guyy_p, guzz_p, guxy_p, guyz_p, guxz_p,
     $           psi_p)
            call Exact__metric(
     $           decoded_exact_model,
     $           x+2*dx, y+2*dy, z+2*dz, t+2*dt,
     $           gdtt_p_p, gdtx_p_p, gdty_p_p, gdtz_p_p,
     $           gdxx_p_p, gdyy_p_p, gdzz_p_p, gdxy_p_p, gdyz_p_p, gdxz_
     &p_p,
     $           gutt_p_p, gutx_p_p, guty_p_p, gutz_p_p,
     $           guxx_p_p, guyy_p_p, guzz_p_p, guxy_p_p, guyz_p_p, guxz_
     &p_p,
     $           psi_p_p)
            
            dgdtt = (- gdtt_m_m - 16*gdtt_m + 30*gdtt_0 - 16*gdtt_p - gd
     &tt_p_p) / (12*eps**2)
            dgdtx = (- gdtx_m_m - 16*gdtx_m + 30*gdtx_0 - 16*gdtx_p - gd
     &tx_p_p) / (12*eps**2)
            dgdty = (- gdty_m_m - 16*gdty_m + 30*gdty_0 - 16*gdty_p - gd
     &ty_p_p) / (12*eps**2)
            dgdtz = (- gdtz_m_m - 16*gdtz_m + 30*gdtz_0 - 16*gdtz_p - gd
     &tz_p_p) / (12*eps**2)
            dgdxx = (- gdxx_m_m - 16*gdxx_m + 30*gdxx_0 - 16*gdxx_p - gd
     &xx_p_p) / (12*eps**2)
            dgdyy = (- gdyy_m_m - 16*gdyy_m + 30*gdyy_0 - 16*gdyy_p - gd
     &yy_p_p) / (12*eps**2)
            dgdzz = (- gdzz_m_m - 16*gdzz_m + 30*gdzz_0 - 16*gdzz_p - gd
     &zz_p_p) / (12*eps**2)
            dgdxy = (- gdxy_m_m - 16*gdxy_m + 30*gdxy_0 - 16*gdxy_p - gd
     &xy_p_p) / (12*eps**2)
            dgdyz = (- gdyz_m_m - 16*gdyz_m + 30*gdyz_0 - 16*gdyz_p - gd
     &yz_p_p) / (12*eps**2)
            dgdxz = (- gdxz_m_m - 16*gdxz_m + 30*gdxz_0 - 16*gdxz_p - gd
     &xz_p_p) / (12*eps**2)
            dgutt = (- gutt_m_m - 16*gutt_m + 30*gutt_0 - 16*gutt_p - gu
     &tt_p_p) / (12*eps**2)
            dgutx = (- gutx_m_m - 16*gutx_m + 30*gutx_0 - 16*gutx_p - gu
     &tx_p_p) / (12*eps**2)
            dguty = (- guty_m_m - 16*guty_m + 30*guty_0 - 16*guty_p - gu
     &ty_p_p) / (12*eps**2)
            dgutz = (- gutz_m_m - 16*gutz_m + 30*gutz_0 - 16*gutz_p - gu
     &tz_p_p) / (12*eps**2)
            dguxx = (- guxx_m_m - 16*guxx_m + 30*guxx_0 - 16*guxx_p - gu
     &xx_p_p) / (12*eps**2)
            dguyy = (- guyy_m_m - 16*guyy_m + 30*guyy_0 - 16*guyy_p - gu
     &yy_p_p) / (12*eps**2)
            dguzz = (- guzz_m_m - 16*guzz_m + 30*guzz_0 - 16*guzz_p - gu
     &zz_p_p) / (12*eps**2)
            dguxy = (- guxy_m_m - 16*guxy_m + 30*guxy_0 - 16*guxy_p - gu
     &xy_p_p) / (12*eps**2)
            dguyz = (- guyz_m_m - 16*guyz_m + 30*guyz_0 - 16*guyz_p - gu
     &yz_p_p) / (12*eps**2)
            dguxz = (- guxz_m_m - 16*guxz_m + 30*guxz_0 - 16*guxz_p - gu
     &xz_p_p) / (12*eps**2)
            dpsi  = (- psi_m_m  - 16*psi_m  + 30*psi_0  - 16*psi_p  - ps
     &i_p_p ) / (12*eps**2)
            
         else
            
            call Exact__metric_deriv(
     $           decoded_exact_model,
     $           dir2,
     $           x-2*dx, y-2*dy, z-2*dz, t-2*dt,
     $           gdtt_m_m, gdtx_m_m, gdty_m_m, gdtz_m_m,
     $           gdxx_m_m, gdyy_m_m, gdzz_m_m, gdxy_m_m, gdyz_m_m, gdxz_
     &m_m,
     $           gutt_m_m, gutx_m_m, guty_m_m, gutz_m_m,
     $           guxx_m_m, guyy_m_m, guzz_m_m, guxy_m_m, guyz_m_m, guxz_
     &m_m,
     $           psi_m_m)
            call Exact__metric_deriv(
     $           decoded_exact_model,
     $           dir2,
     $           x-dx, y-dy, z-dz, t-dt,
     $           gdtt_m, gdtx_m, gdty_m, gdtz_m,
     $           gdxx_m, gdyy_m, gdzz_m, gdxy_m, gdyz_m, gdxz_m,
     $           gutt_m, gutx_m, guty_m, gutz_m,
     $           guxx_m, guyy_m, guzz_m, guxy_m, guyz_m, guxz_m,
     $           psi_m)
            call Exact__metric_deriv(
     $           decoded_exact_model,
     $           dir2,
     $           x+dx, y+dy, z+dz, t+dt,
     $           gdtt_p, gdtx_p, gdty_p, gdtz_p,
     $           gdxx_p, gdyy_p, gdzz_p, gdxy_p, gdyz_p, gdxz_p,
     $           gutt_p, gutx_p, guty_p, gutz_p,
     $           guxx_p, guyy_p, guzz_p, guxy_p, guyz_p, guxz_p,
     $           psi_p)
            call Exact__metric_deriv(
     $           decoded_exact_model,
     $           dir2,
     $           x+2*dx, y+2*dy, z+2*dz, t+2*dt,
     $           gdtt_p_p, gdtx_p_p, gdty_p_p, gdtz_p_p,
     $           gdxx_p_p, gdyy_p_p, gdzz_p_p, gdxy_p_p, gdyz_p_p, gdxz_
     &p_p,
     $           gutt_p_p, gutx_p_p, guty_p_p, gutz_p_p,
     $           guxx_p_p, guyy_p_p, guzz_p_p, guxy_p_p, guyz_p_p, guxz_
     &p_p,
     $           psi_p_p)
            
            dgdtt = (- gdtt_p_p + 8*gdtt_p - 8*gdtt_m + gdtt_m_m) / (12*
     &eps)
            dgdtx = (- gdtx_p_p + 8*gdtx_p - 8*gdtx_m + gdtx_m_m) / (12*
     &eps)
            dgdty = (- gdty_p_p + 8*gdty_p - 8*gdty_m + gdty_m_m) / (12*
     &eps)
            dgdtz = (- gdtz_p_p + 8*gdtz_p - 8*gdtz_m + gdtz_m_m) / (12*
     &eps)
            dgdxx = (- gdxx_p_p + 8*gdxx_p - 8*gdxx_m + gdxx_m_m) / (12*
     &eps)
            dgdyy = (- gdyy_p_p + 8*gdyy_p - 8*gdyy_m + gdyy_m_m) / (12*
     &eps)
            dgdzz = (- gdzz_p_p + 8*gdzz_p - 8*gdzz_m + gdzz_m_m) / (12*
     &eps)
            dgdxy = (- gdxy_p_p + 8*gdxy_p - 8*gdxy_m + gdxy_m_m) / (12*
     &eps)
            dgdyz = (- gdyz_p_p + 8*gdyz_p - 8*gdyz_m + gdyz_m_m) / (12*
     &eps)
            dgdxz = (- gdxz_p_p + 8*gdxz_p - 8*gdxz_m + gdxz_m_m) / (12*
     &eps)
            dgutt = (- gutt_p_p + 8*gutt_p - 8*gutt_m + gutt_m_m) / (12*
     &eps)
            dgutx = (- gutx_p_p + 8*gutx_p - 8*gutx_m + gutx_m_m) / (12*
     &eps)
            dguty = (- guty_p_p + 8*guty_p - 8*guty_m + guty_m_m) / (12*
     &eps)
            dgutz = (- gutz_p_p + 8*gutz_p - 8*gutz_m + gutz_m_m) / (12*
     &eps)
            dguxx = (- guxx_p_p + 8*guxx_p - 8*guxx_m + guxx_m_m) / (12*
     &eps)
            dguyy = (- guyy_p_p + 8*guyy_p - 8*guyy_m + guyy_m_m) / (12*
     &eps)
            dguzz = (- guzz_p_p + 8*guzz_p - 8*guzz_m + guzz_m_m) / (12*
     &eps)
            dguxy = (- guxy_p_p + 8*guxy_p - 8*guxy_m + guxy_m_m) / (12*
     &eps)
            dguyz = (- guyz_p_p + 8*guyz_p - 8*guyz_m + guyz_m_m) / (12*
     &eps)
            dguxz = (- guxz_p_p + 8*guxz_p - 8*guxz_m + guxz_m_m) / (12*
     &eps)
            dpsi  = (- psi_p_p  + 8*psi_p  - 8*psi_m  + psi_m_m ) / (12*
     &eps)
            
         end if
         
      else
         call CCTK_Warn(0,747,"Bona_Masso_data.F","Exact", "internal err
     &or")
      end if
      
      end
