! driver.f90
!
! Simulation of the wave equation in three dimensions 
! with a Gaussian as the initial data.

program wavedrive
  implicit none
  
  integer nx, ny, nz, maxitr, dosmooth, n, interval
  
  real*8, dimension(:,:,:),allocatable,target :: &
       anm, an, anp, bnm, bn, bnp, cnm, cn, cnp, &
       dnm, dn, dnp, phi, temp

!  real*8, dimension(:,:,:),allocatable,target :: temp
  
  real*8, dimension(:,:,:),pointer :: phim, phip
  
  real*8 width, ampl, vx, vy, vz, epsilon
  
  real*8 bbox(6), lb(3), cfl, dt, h, time
  
!  Read values from wave3d.in
 	
  call setup(nx, ny, nz, maxitr, dosmooth, width, ampl, &
       vx, vy, vz, epsilon, bbox, lb, cfl, dt, h, time, &
       interval) 
  
  allocate(anm(nx,ny,nz), an(nx,ny,nz), anp(nx,ny,nz), &
       bnm(nx,ny,nz), bn(nx,ny,nz), bnp(nx,ny,nz), &
       cnm(nx,ny,nz), cn(nx,ny,nz), cnp(nx,ny,nz), &
       dnm(nx,ny,nz), dn(nx,ny,nz), dnp(nx,ny,nz), &
       temp(nx,ny,nz), phi(nx,ny,nz))

! Only two levels at a time of phi are needed

  phim => temp
  phip => temp

! Put down a Gaussian curve

  call initial(nx,ny,nz, h, vx,vy,vz, width, ampl, bbox, lb, &
       an, bn, cn, dn, temp)
  
  call initstep(nx,ny,nz, dt, h, epsilon, bbox, anm, an, anp, &
       bnm, bn, bnp, cnm, cn, cnp, dnm, dn, dnp, temp)
  
  open(unit=2,access="sequential",&
       form="unformatted",file="wave3d.out")
  
  call startdump(nx,ny,nz, bbox, h, width, ampl, cfl, dosmooth)   
  
  do n = 1, maxitr
     
     if (mod(n, interval) .eq. 0) then
        call filedump(nx,ny,nz, time, an)
        call filedump(nx,ny,nz, time, bn)
        call filedump(nx,ny,nz, time, cn)
        call filedump(nx,ny,nz, time, dn)
        call filedump(nx,ny,nz, time, phi)       
!        call dump_2d_slice("phi", time, phi, nx,ny,nz, (ny+1)/2, &
!	bbox(1), bbox(3), h)
     end if
     
     time = dt*float(n)

     call boundupdt(nx,ny,nz, an,anp, bn,bnp, cn,cnp, dn,dnp, phip)

     call leapfrogA(nx,ny,nz, dt, h, anm,an,anp, bn,cn,dn)

     call leapfrogBCD(nx,ny,nz, dt, h, an, bnm,bn,bnp, cnm,cn,cnp, &
          dnm,dn,dnp)

     call updatephi(nx,ny,nz, dt, h, an, phip,phim)

     if (dosmooth .ne. 0) then
        call smooth(nx, ny, nz, anp)
        call smooth(nx, ny, nz, bnp)
        call smooth(nx, ny, nz, cnp)
        call smooth(nx, ny, nz, dnp)
        call smooth(nx, ny, nz, phip)
     end if
     
     anm = an
     bnm = bn
     cnm = cn
     dnm = dn
     
     an = anp
     bn = bnp
     cn = cnp
     dn = dnp

!    swaps phim with phi and phip with phi.  (phim == phip)

     call carefulswap(nx,ny,nz, phim, phi)
     
  end do
  
  close (unit=2)
  
  stop
end
