/*@@
   @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_Evolution
   @date       
   @author     Tom Goodale
   @desc 
               Evolution for the wave equation
   @enddesc 
   @history 
 
   @endhistory 

@@*/

      subroutine WaveToyf90_Evolution(CCTK_ARGUMENTS)

      implicit none

      DECLARE_CCTK_ARGUMENTS
      DECLARE_CCTK_PARAMETERS

!     Declare local variables
      INTEGER  i,j,k
      INTEGER  istart, jstart, kstart, iend, jend, kend

      CCTK_REAL  dx,dy,dz,dt
      CCTK_REAL  dx2,dy2,dz2,dt2
      CCTK_REAL  dx2i,dy2i,dz2i

      CCTK_REAL  factor

!     Set up shorthands
!     -----------------
      dx = CCTK_DELTA_SPACE(1)
      dy = CCTK_DELTA_SPACE(2)
      dz = CCTK_DELTA_SPACE(3)
      dt = CCTK_DELTA_TIME

      dx2 = dx*dx
      dy2 = dy*dy
      dz2 = dz*dz
      dt2 = dt*dt

      dx2i = 1.0/dx2
      dy2i = 1.0/dy2
      dz2i = 1.0/dz2

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

      factor = 2*(1 - (dt2)*(dx2i + dy2i + dz2i))

!     Do the evolution
!     ----------------
      do k = kstart, kend
         do j = jstart, jend
            do i = istart, iend
               
               phi_n(i,j,k) = factor*phi(i,j,k) -
     1              phi_p(i,j,k) + (dt2) * 
     2              ((phi(i+1,j,k)+phi(i-1,j,k))*dx2i
     3              +(phi(i,j+1,k)+phi(i,j-1,k))*dy2i
     4              +(phi(i,j,k+1)+phi(i,j,k-1))*dz2i)

            end do
         end do
      end do
         
      end subroutine WaveToyF90_Evolution


 /*@@
   @routine    WaveToyF90_Boundaries
   @date       
   @author     Tom Goodale
   @desc 
               Boundary conditions for the wave equation
   @enddesc 
   @history 
 
   @endhistory 

@@*/

      subroutine WaveToyF90_Boundaries(CCTK_ARGUMENTS)

      implicit none

      DECLARE_CCTK_ARGUMENTS
      DECLARE_CCTK_PARAMETERS
      DECLARE_CCTK_FUNCTIONS

      integer i
      integer :: ierr=-1
      integer,dimension(3):: sw=1
      CCTK_REAL,parameter :: sval = 0.0
      CCTK_REAL,parameter :: rzero = 0.0
      CCTK_REAL,parameter :: rone = 1.0
      integer,parameter :: swdir = 1
      integer,parameter :: ione = 1,mione=-1
      integer,parameter :: itwo = 2,mitwo=-2
      integer,parameter :: ithree = 3,mithree=-3
      integer,parameter :: ifour = 4,izero = 0
      call CartSymGN(ierr,cctkGH,"wavetoy::scalarevolve")

      if (CCTK_EQUALS(bound,"flat")) then
         call BndFlatVN(ierr,cctkGH,sw,"wavetoy::phi")
      else if (CCTK_EQUALS(bound,"zero")) then
c         call BndScalarVN(ierr,cctkGH,sw,sval,"wavetoy::phi")
         call BndScalarDirVN(ierr,cctkGH,swdir,mione  ,sval,"wavetoy::phi")
         call BndScalarDirVN(ierr,cctkGH,swdir,ione   ,sval,"wavetoy::phi")
         call BndScalarDirVN(ierr,cctkGH,swdir,mitwo  ,sval,"wavetoy::phi")
         call BndScalarDirVN(ierr,cctkGH,swdir,itwo   ,sval,"wavetoy::phi")
         call BndScalarDirVN(ierr,cctkGH,swdir,mithree,sval,"wavetoy::phi")
         call BndScalarDirVN(ierr,cctkGH,swdir,ithree ,sval,"wavetoy::phi")
      else if (CCTK_EQUALS(bound,"radiation")) then
         call BndRadiativeVN(ierr,cctkGH,sw,rzero,rone,
     &        "wavetoy::phi","wavetoy::phi")
      end if

      if (ierr < 0)  then
        call CCTK_WARN(0,"WaveToyF90_Boundaries: Error in boundary routines - giving up!")
      end if

      end subroutine WaveToyF90_Boundaries