path: root/src/Bona_Masso_data.F

c This routine calculates Bona-Masso initial data, making use of the
c subroutine Exact__metric() to calculate the spacetime metric and its
c inverse.  Note that this use of the Bona-Masso variables is independent
c of how (or even if) we are evolving the Einstein equations -- here
c the Bona-Masso variables are "just" used as intermediate variables.

c $Header$

#include "cctk.h"
#include "cctk_Parameters.h"

      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
      CCTK_INT decoded_exact_model
      logical psi_on
      CCTK_REAL 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!)

      CCTK_REAL
     $       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,
     $     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,
     $     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,
     $     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,
     $     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, 1,
     $        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, 2,
     $        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, 3,
     $        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, 2,
     $        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, 3,
     $        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, 3,
     $        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
      DECLARE_CCTK_PARAMETERS
      
      CCTK_INT
     $     decoded_exact_model,
     $     dir
      CCTK_REAL
     $     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 (CCTK_WARN_ABORT, "internal error")
      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
      DECLARE_CCTK_PARAMETERS
      
      CCTK_INT
     $     decoded_exact_model,
     $     dir1, dir2
      CCTK_REAL
     $     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 (CCTK_WARN_ABORT, "internal error")
      end if
      if (dir2.lt.0 .or. dir2.gt.3) then
         call CCTK_WARN (CCTK_WARN_ABORT, "internal error")
      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 - gdtt_p_p) / (12*eps**2)
            dgdtx = (- gdtx_m_m - 16*gdtx_m + 30*gdtx_0 - 16*gdtx_p - gdtx_p_p) / (12*eps**2)
            dgdty = (- gdty_m_m - 16*gdty_m + 30*gdty_0 - 16*gdty_p - gdty_p_p) / (12*eps**2)
            dgdtz = (- gdtz_m_m - 16*gdtz_m + 30*gdtz_0 - 16*gdtz_p - gdtz_p_p) / (12*eps**2)
            dgdxx = (- gdxx_m_m - 16*gdxx_m + 30*gdxx_0 - 16*gdxx_p - gdxx_p_p) / (12*eps**2)
            dgdyy = (- gdyy_m_m - 16*gdyy_m + 30*gdyy_0 - 16*gdyy_p - gdyy_p_p) / (12*eps**2)
            dgdzz = (- gdzz_m_m - 16*gdzz_m + 30*gdzz_0 - 16*gdzz_p - gdzz_p_p) / (12*eps**2)
            dgdxy = (- gdxy_m_m - 16*gdxy_m + 30*gdxy_0 - 16*gdxy_p - gdxy_p_p) / (12*eps**2)
            dgdyz = (- gdyz_m_m - 16*gdyz_m + 30*gdyz_0 - 16*gdyz_p - gdyz_p_p) / (12*eps**2)
            dgdxz = (- gdxz_m_m - 16*gdxz_m + 30*gdxz_0 - 16*gdxz_p - gdxz_p_p) / (12*eps**2)
            dgutt = (- gutt_m_m - 16*gutt_m + 30*gutt_0 - 16*gutt_p - gutt_p_p) / (12*eps**2)
            dgutx = (- gutx_m_m - 16*gutx_m + 30*gutx_0 - 16*gutx_p - gutx_p_p) / (12*eps**2)
            dguty = (- guty_m_m - 16*guty_m + 30*guty_0 - 16*guty_p - guty_p_p) / (12*eps**2)
            dgutz = (- gutz_m_m - 16*gutz_m + 30*gutz_0 - 16*gutz_p - gutz_p_p) / (12*eps**2)
            dguxx = (- guxx_m_m - 16*guxx_m + 30*guxx_0 - 16*guxx_p - guxx_p_p) / (12*eps**2)
            dguyy = (- guyy_m_m - 16*guyy_m + 30*guyy_0 - 16*guyy_p - guyy_p_p) / (12*eps**2)
            dguzz = (- guzz_m_m - 16*guzz_m + 30*guzz_0 - 16*guzz_p - guzz_p_p) / (12*eps**2)
            dguxy = (- guxy_m_m - 16*guxy_m + 30*guxy_0 - 16*guxy_p - guxy_p_p) / (12*eps**2)
            dguyz = (- guyz_m_m - 16*guyz_m + 30*guyz_0 - 16*guyz_p - guyz_p_p) / (12*eps**2)
            dguxz = (- guxz_m_m - 16*guxz_m + 30*guxz_0 - 16*guxz_p - guxz_p_p) / (12*eps**2)
            dpsi  = (- psi_m_m  - 16*psi_m  + 30*psi_0  - 16*psi_p  - psi_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 (CCTK_WARN_ABORT, "internal error")
      end if
      
      end