/*@@
   @file      WaveToyF90.F
   @date      
   @author    Tom Goodale
   @desc 
              Evolution routines for the wave equation solver
   @enddesc 
 @@*/

#include "cctk.h" 
#include "cctk_parameters.h"
#include "cctk_arguments.h"


 /*@@
   @routine    WaveToyF90_Boundaries
   @date       
   @author     Tom Goodale
   @desc 
               Boundary conditions for the wave equation
   @enddesc 
   @calls      ApplyFlatBC,ApplyRadiativeBC 
   @calledby   
   @history 
 
   @endhistory 

@@*/

      subroutine WaveToyF90_Boundaries(CCTK_FARGUMENTS)

      implicit none

      DECLARE_CCTK_FARGUMENTS
      DECLARE_CCTK_PARAMETERS

      INTEGER,DIMENSION(3):: sw=1
      CCTK_REAL,PARAMETER :: zero = 0.0
      INTEGER CCTK_Equals

      call ApplySymmetry(cctkGH,"wavetoy::scalarevolve")

      if (CCTK_Equals(bound,"flat")==1) then
         call ApplyFlatBC(cctkGH,sw,"wavetoy::phi")
      else if (CCTK_Equals(bound,"radiation")==1) then
         call ApplyRadiativeBC(cctkGH,zero,sw,"wavetoy::phi","wavetoyf90::phi_old")
      end if

      end subroutine wavetoyf90_boundaries

c --------------------------------------------------------------

 /*@@
   @routine    WaveToyF90_Evolution
   @date       
   @author     Tom Goodale
   @desc 
               Evolution for the wave equation
   @enddesc 
   @calls      CCTK_SyncGroup, wavetoy_boundaries
   @calledby   
   @history 
 
   @endhistory 

@@*/

      subroutine wavetoyf90_evolution(CCTK_FARGUMENTS)

      implicit none

      DECLARE_CCTK_FARGUMENTS
      DECLARE_CCTK_PARAMETERS

c     Declare local variables
      INTEGER :: i,j,k
      INTEGER :: istart, jstart, kstart, iend, jend, kend
      CCTK_REAL :: dx,dy,dz,dt

c     Set up shorthands
c     -----------------
      dx = cctk_delta_space(1)
      dy = cctk_delta_space(2)
      dz = cctk_delta_space(3)
      dt = cctk_delta_time

      istart = 2
      jstart = 2
      kstart = 2
      
      iend = cctk_lsh(1)-1
      jend = cctk_lsh(2)-1
      kend = cctk_lsh(3)-1

c     Do the evolution
c     ----------------
      do k = kstart, kend
         do j = jstart, jend
            do i = istart, iend
               
               phi_tmp(i,j,k) = 
     1		 2.0*(1.0 - (dt**2)*(1.0/dx**2 +   
     2              1.0/dy**2 +1.0/dz**2))*phi(i,j,k) -
     3              phi_old(i,j,k) + (dt**2) * 
     5              ((phi(i+1,j,k)+phi(i-1,j,k))/dx**2
     6              +(phi(i,j+1,k)+phi(i,j-1,k))/dy**2
     7              +(phi(i,j,k+1)+phi(i,j,k-1))/dz**2)

            end do
         end do
      end do
         
c     Update timeslices
c     -----------------
      phi_old = phi
      phi     = phi_tmp

c     Apply boundary conditions
c     -------------------------
      call wavetoyf90_boundaries(CCTK_FARGUMENTS)

c     Synchronize
c     -----------
      call CCTK_SyncGroup(cctkGH,"wavetoy::scalarevolve")

      end subroutine wavetoyf90_evolution


c --------------------------------------------------------------