
! --------------------------------------------------------------
! initial: Put down a Gaussian as initial data
!
subroutine initial(nx,ny,nz, h, vx,vy,vz, width, ampl, bbox, lb, &
     a0, b0, c0, d0, phi0)
  implicit none
    
  integer nx, ny, nz

  real*8 h, vx, vy, vz, width, ampl, bbox(6), lb(3)

  real*8, dimension (nx,ny,nz) :: a0, b0, c0, d0, phi0

  ! Locals

  integer i, j, k

  real*8 x, y, z, t, rhop, rhom

  t = 0.0d0

  do i=1, nx
     do j=1, ny
        do k=1, nz
           
           x=h*(i+lb(1)-1) + bbox(1)
           y=h*(j+lb(2)-1) + bbox(2)
           z=h*(k+lb(3)-1) + bbox(3)

           rhop = sqrt( (x-vx*t)**2 + (y-vy*t)**2 + (z-vz*t)**2)

           phi0(i,j,k) = ampl * exp (-rhop**2 / width**2)
  
      end do
     end do
  end do

  do i=1, nx
     do j=1, ny
        do k=1, nz
           
           x=h*(i+lb(1)-1) + bbox(1)
           y=h*(j+lb(2)-1) + bbox(2)
           z=h*(k+lb(3)-1) + bbox(3)

           rhom = sqrt((x-vx*t)**2 +  (y-vy*t)**2 + (z-vz*t)**2)

           a0(i,j,k) = 2.0 * ampl / width**2 * ( &
                exp(-rhom**2/width**2) * &
                ((x-vx*t)*vx + (y-vy*t)*vy + (z-vz*y)*vz))

           b0(i,j,k) = -2.0 * (x-vx*t) * phi0(i,j,k) / width**2
           
           c0(i,j,k) = -2.0 * (y-vy*t) * phi0(i,j,k) / width**2

           d0(i,j,k) = -2.0 * (z-vz*t) * phi0(i,j,k) / width**2
       
        end do
     end do
  end do

  return
end subroutine initial

! --------------------------------------------------------------
! leapfrogA: use Leap Frog scheme to calculate A(n+1) from A(n) and A(n-1)
!
subroutine leapfrogA(nx,ny,nz, dt, h, anm,an,anp, bn,cn,dn)
  implicit none

  integer nx,ny,nz

  real*8 dt, h

  real*8, dimension(nx,ny,nz) :: anm,an,anp, bn,cn,dn

  anp(2:nx-1,2:ny-1,2:nz-1) = anm(2:nx-1,2:ny-1,2:nz-1) + &
       (dt/h) * (bn(3:nx,2:ny-1,2:nz-1) - bn(1:nx-2,2:ny-1,2:nz-1)) + &
       (dt/h) * (cn(2:nx,3:ny,  2:nz-1) - cn(2:nx-2,1:ny-2,2:nz-1)) + &
       (dt/h) * (dn(2:nx,2:ny-1,3:nz) -   dn(2:nx-2,2:ny-1,1:nz-2))

  return
end subroutine leapfrogA

! --------------------------------------------------------------
! leapfrogBCD: use Leap Frog scheme to calculate U(n+1) 
!              from U(n) and U(n-1)
!
subroutine leapfrogBCD(nx, ny, nz, dt, h, an, bnm, bn, bnp, &
     cnm, cn, cnp, dnm, dn, dnp)			
  implicit none

  integer nx, ny, nz

  real*8 dt, h

  real*8, dimension(nx,ny,nz) :: an, bnm, bn, bnp, cnm, cn, cnp, &
	dnm, dn, dnp
  
  bnp(2:nx-1,2:ny-1,2:nz-1) = bnm(2:nx-1,2:ny-1,2:nz-1) + &
       (dt/h) * (an(3:nx,2:ny-1,2:nz-1) - an(1:nx-2,2:ny-1,2:nz-1))

  cnp(2:nx-1,2:ny-1,2:nz-1) = cnm(2:nx-1,2:ny-1,2:nz-1) + &
       (dt/h) * (an(2:nx-1,3:ny,2:nz-1) - an(2:nx-1,1:ny-2,2:nz-1))
  
  dnp(2:nx-1,2:ny-1,2:nz-1) = dnm(2:nx-1,2:ny-1,2:nz-1) + &
       (dt/h) * (an(2:nx-1,2:ny-1,3:nz) - an(2:nx-1,2:ny-1,1:nz-2))
        

  return
end subroutine leapfrogBCD

! --------------------------------------------------------------
! updatephi: update phi (NB: phim and phip have the same storage)
!
subroutine updatephi(nx, ny, nz, dt, h, an, phip, phim)
  implicit none
  
  integer nx, ny, nz

  real*8 dt, h

  real*8, dimension(nx,ny,nz) :: an, phip, phim

  phip(1:nx,1:ny,1:nz) = phim(1:nx,1:ny,1:nz) + 2.0 * &
       dt * an(1:nx,1:ny,1:nz)

  return
end subroutine updatephi

! --------------------------------------------------------------
! boundupdt: update the boundaries using simple boundary condition 
!
subroutine boundupdt(nx, ny, nz, an, anp, bn, bnp, &
     cn, cnp, dn, dnp, phip)
  implicit none

  integer nx, ny, nz

  real*8, dimension(nx,ny,nz) :: an, anp, bn, bnp, &
     cn, cnp, dn, dnp, phip

!   phip(1,:,:) = 0.0d0
   anp(1,:,:)  = 0.0d0
   bnp(1,:,:)  = bn(1,:,:)
   cnp(1,:,:)  = cn(1,:,:)
   dnp(1,:,:)  = dn(1,:,:)
   
!   phip(nx,:,:) = 0.0d0
   anp(nx,:,:)  = 0.0d0
   bnp(nx,:,:)  = bn(nx,:,:)
   cnp(nx,:,:)  = cn(nx,:,:)
   dnp(nx,:,:)  = dn(nx,:,:)
   
!   phip(:,1,:) = 0.0d0
   anp(:,1,:)  = 0.0d0
   bnp(:,1,:)  = bn(:,1,:)
   cnp(:,1,:)  = cn(:,1,:)
   dnp(:,1,:)  = dn(:,1,:)
   
!   phip(:,ny,:) = 0.0d0
   anp(:,ny,:)  = 0.0d0
   bnp(:,ny,:)  = bn(:,ny,:)
   cnp(:,ny,:)  = cn(:,ny,:)
   dnp(:,ny,:)  = dn(:,ny,:)
   
!   phip(:,:,1) = 0.0d0
   anp(:,:,1)  = 0.0d0
   bnp(:,:,1)  = bn(:,:,1)
   cnp(:,:,1)  = cn(:,:,1)
   dnp(:,:,1)  = dn(:,:,1)
   
!   phip(:,:,nz) = 0.0d0
   anp(:,:,nz)  = 0.0d0
   bnp(:,:,nz)  = bn(:,:,nz)
   cnp(:,:,nz)  = cn(:,:,nz)
   dnp(:,:,nz)  = dn(:,:,nz)
      
   return
end subroutine boundupdt

! --------------------------------------------------------------
! smooth: use spatial averaging to smooth an array
!
subroutine smooth(nx,ny,nz, a)
  implicit none

  integer nx, ny, nz

  real*8, dimension(nx,ny,nz) :: a

  a(2:nx-1,2:ny-1,2:nz-1) = 0.5 * a(2:nx-1,2:ny-1,2:nz-1) + ( &
       a(3:nx  ,2:ny-1,2:nz-1) + a(1:nx-2,2:ny-1,2:nz-1) + &
       a(2:nx-1,3:ny  ,2:nz-1) + a(2:nx-1,1:ny-2,2:nz-1) + &
       a(2:nx-1,2:ny-1,3:nz  ) + a(2:nx-1,2:ny-1,1:nz-2) ) / 12.0

  return
end subroutine smooth
