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

! Using Cactus infrastructure
#include "cctk.h" 

! Using Cactus parameters
#include "cctk_parameters.h"

! Using Cactus arguments lists
#include "cctk_arguments.h"


 /*@@
   @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

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

!     Set up shorthands
!     -----------------
      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

!     Do the evolution
!     ----------------
      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
         
!     Synchronize
!     -----------
      call CCTK_SyncGroup(cctkGH,"wavetoyf90::temps")

!     Apply boundary conditions
!     -------------------------
      call WaveToyF90_Boundaries(CCTK_PASS_FTOF)

!     Update timeslices
!     -----------------
      phi_old = phi
      phi     = phi_tmp


      end subroutine WaveToyF90_Evolution


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

@@*/

      subroutine WaveToyF90_Boundaries(CCTK_FARGUMENTS)

      implicit none

      DECLARE_CCTK_FARGUMENTS
      DECLARE_CCTK_PARAMETERS
      DECLARE_CCTK_FUNCTIONS

      integer :: ierr
      integer,dimension(3):: sw=1
      CCTK_REAL,parameter :: zero = 0.0
      CCTK_REAL,parameter :: one = 1.0
      
      call CartSymBCGroup(ierr,cctkGH,"wavetoyf90::temps")

      if (CCTK_EQUALS(bound,"flat")) then
         call FlatBCVar(ierr,cctkGH,sw,"wavetoyf90::phi_tmp")
      else if (CCTK_EQUALS(bound,"radiation")) then
         call RadiativeBCVar(ierr,cctkGH,zero,one,sw,
     &        "wavetoyf90::phi_tmp","wavetoy::phi")
      end if

      if (ierr < 0)  then
        call CCTK_WARN(0,"Boundary conditions not applied - giving up!")
      end if

      end subroutine WaveToyF90_Boundaries