c     $Header$

#include "cctk.h"
#include "cctk_Arguments.h"
#include "cctk_Functions.h"
#include "cctk_Parameters.h"

      subroutine WaveToyStaggered_Evolve (CCTK_ARGUMENTS)
      implicit none
      DECLARE_CCTK_ARGUMENTS
      DECLARE_CCTK_FUNCTIONS
      DECLARE_CCTK_PARAMETERS
      
      CCTK_REAL dt, dx(3)
      integer i, j, k
      integer d
      integer ierr
      
      dt = CCTK_DELTA_TIME
      do d = 1, 3
         dx(d) = CCTK_DELTA_SPACE(d)
      end do
      
      do k = 2, cctk_lsh(3) - 1
         do j = 2, cctk_lsh(2) - 1
            do i = 2, cctk_lsh(1) - 1
               phi(i,j,k) =
     $              ( (psix_p(i,j,k) - psix_p(i-1,j,k)) / dx(1)
     $              + (psiy_p(i,j,k) - psiy_p(i,j-1,k)) / dx(2)
     $              + (psiz_p(i,j,k) - psiz_p(i,j,k-1)) / dx(3))
            end do
         end do
      end do
      call CCTK_SyncGroup (ierr, cctkGH, "WaveToyStaggered::scalarevolve")
      if (ierr.ne.0) call CCTK_WARN (0, "internal error")
      
      do k = 2, cctk_lsh(3) - 1
         do j = 2, cctk_lsh(2) - 1
            do i = 2, cctk_lsh(1) - 1
               psix(i,j,k) = (phi_p(i+1,j,k) - phi_p(i,j,k)) / dx(1)
               psiy(i,j,k) = (phi_p(i,j+1,k) - phi_p(i,j,k)) / dx(2)
               psiz(i,j,k) = (phi_p(i,j,k+1) - phi_p(i,j,k)) / dx(3)
            end do
         end do
      end do
      call CCTK_SyncGroup (ierr, cctkGH, "WaveToyStaggered::scalarevolvex")
      if (ierr.ne.0) call CCTK_WARN (0, "internal error")
      call CCTK_SyncGroup (ierr, cctkGH, "WaveToyStaggered::scalarevolvey")
      if (ierr.ne.0) call CCTK_WARN (0, "internal error")
      call CCTK_SyncGroup (ierr, cctkGH, "WaveToyStaggered::scalarevolvez")
      if (ierr.ne.0) call CCTK_WARN (0, "internal error")
      
      end