/*@@
   @file      WaveToyFortran.F90
   @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    WaveToyFortran_Evolution
   @date       
   @author     Tom Goodale
   @desc 
               Evolution for the wave equation
   @enddesc 
   @calls      CCTK_SyncGroup, wavetoy_boundaries
   @calledby   
   @history 
 
   @endhistory 

@@*/

subroutine WaveToyFreeF90_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) =                         &
                2.0*(1.0 - (dt**2)*(1.0/dx**2 +     &
                1.0/dy**2 +1.0/dz**2))*phi(i,j,k) - &
                phi_old(i,j,k) + (dt**2) *          &
                ((phi(i+1,j,k)+phi(i-1,j,k))/dx**2  &
                +(phi(i,j+1,k)+phi(i,j-1,k))/dy**2  &
                +(phi(i,j,k+1)+phi(i,j,k-1))/dz**2)
           
        end do
     end do
  end do
     
  ! Synchronize
  ! -----------
  call CCTK_SyncGroup(cctkGH,"wavetoyFreeF90::temps")
  
  ! Apply boundary conditions
  ! -------------------------
  call WaveToyFreeF90_Boundaries(CCTK_PASS_FTOF)
  
  ! Update timeslices
  ! -----------------
  phi_old = phi
  phi     = phi_tmp
  

end subroutine WaveToyFreeF90_Evolution


 /*@@
   @routine    WaveToyFortran_Boundaries
   @date       
   @author     Tom Goodale
   @desc 
               Boundary conditions for the wave equation
   @enddesc 
   @calls      FlatBC,RadiativeBC 
   @calledby   
   @history 
 
   @endhistory 

@@*/

subroutine WaveToyFreeF90_Boundaries(CCTK_FARGUMENTS)
  
  implicit none
  
  DECLARE_CCTK_FARGUMENTS
  DECLARE_CCTK_PARAMETERS
  DECLARE_CCTK_FUNCTIONS
  
  integer :: ierr
  integer, dimension(3):: sw(3)=1
  CCTK_REAL, parameter :: zero = 0.0
  CCTK_REAL, parameter :: one = 1.0
  
  call CartSymBCGroup(ierr,cctkGH,"wavetoyfreef90::temps")
  
  if (CCTK_EQUALS(bound,"flat")) then
     call ApplyFlatBC(ierr,cctkGH,sw,"wavetoyfreef90::phi_tmp")
  else if (CCTK_EQUALS(bound,"radiation")) then
     call RadiativeBCVar(ierr,cctkGH,zero,one,sw,   &
          "wavetoyfreef90::phi_tmp","wavetoy::phi")
  end if
  
  if (ierr < 0)  then
     call CCTK_WARN(0,"Boundary conditions not applied - giving up!")
  end if
  
end subroutine WaveToyFreeF90_Boundaries