

subroutine initialize(params, b_box, shape, anm, an, anp, bnm, bn, bnp, &
     cnm, cn, cnp, dnm, dn, dnp, phi, temp)
  use input_params
  implicit none
  
  real*8, dimension(100):: params
  integer, dimension(3) :: shape
  real*8, dimension(6)  :: b_box
  real*8, dimension(:,:,:), target :: anm, an, anp, &
       bnm, bn, bnp, cnm, cn, cnp, dnm, dn, dnp, phi, temp

  ! Locals

  integer i, j, k, maxstp
  real*8 x, y, z, t, rhop, rhom
  real*8, dimension(:,:,:), pointer :: anpnew, bnpnew, cnpnew,&
       dnpnew, antemp
  real*8 anorm, bnorm, cnorm, dnorm
  real*8 h, dt, epsilon, vx, vy, vz, width, ampl

  real*8 DiffNorm, LeapFrogA, LeapFrogBCD

  h = params(P_GRIDSPACING)
  dt = params(P_TIMESTEP)
  epsilon = params(P_INITEPSILON)
  vx = params(P_VX) 
  vy = params(P_VY)
  vz = params(P_VZ)
  ampl = params(P_AMPL)
  width = params(P_WIDTH)

  t = 0.0d0

  do i=1, shape(1)
     do j=1, shape(2)
        do k=1, shape(3)
           
           x=h*(i-1) + b_box(1)
           y=h*(j-1) + b_box(3)
           z=h*(k-1) + b_box(5)

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

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

  do i=1, shape(1)
     do j=1, shape(2)
        do k=1, shape(3)
           
           x=h*(i-1) + b_box(1)
           y=h*(j-1) + b_box(3)
           z=h*(k-1) + b_box(5)

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

           an(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))

           bn(i,j,k) = -2.0 * (x-vx*t) * phi(i,j,k) / width**2
           
           cn(i,j,k) = -2.0 * (y-vy*t) * phi(i,j,k) / width**2

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

  maxstp = 100

  anpnew => an
  bnpnew => bn
  cnpnew => cn
  dnpnew => dn
  antemp => temp

  anm = an
  bnm = bn
  cnm = cn
  dnm = dn

  anp = an
  bnp = bn
  cnp = cn
  dnp = dn 

  do i = 1, maxstp

     an =  0.5 * (anm + anp)
     bn =  0.5 * (bnm + bnp)
     cn =  0.5 * (cnm + cnp)     
     dn =  0.5 * (dnm + dnp)

     antemp = an

     call leapfrogA(shape, dt, h, anm,an,anpnew, bn,cn,dn)

     call leapfrogBCD(shape, dt, h, antemp, bnm,bn,bnpnew, &
	cnm,cn,cnpnew, dnm,dn,dnpnew)

     anorm = DiffNorm (shape, anpnew, anp)
     bnorm = DiffNorm (shape, bnpnew, bnp)
     cnorm = DiffNorm (shape, cnpnew, cnp)
     dnorm = DiffNorm (shape, dnpnew, dnp)

     anp = anpnew
     bnp = bnpnew
     cnp = cnpnew
     dnp = dnpnew

     write (6,*) 'norms: ', anorm, bnorm, cnorm, dnorm

     if (sqrt(anorm ** 2 + bnorm ** 2 + cnorm ** 2 + dnorm ** 2) &
          .le. epsilon) exit

  end do

  an = anp
  bn = bnp
  cn = cnp
  dn = dnp

  return  
end subroutine initialize


! --------------------------------------------------------------
! DiffNorm: Take the 2-norm of the difference of two arrays.
!
real*8 function DiffNorm(shape, a, b)
  implicit none
  integer shape(3)
  real*8, dimension(shape(1),shape(2),shape(3)) :: a, b


  ! Locals

  integer i, j, k

  DiffNorm = 0.0d0

  do i=1, shape(1)
     do j=1, shape(2)
        do k=1, shape(3)

           DiffNorm = DiffNorm + (a(i,j,k) - b(i,j,k))**2

        end do
     end do
  end do

  DiffNorm = sqrt(DiffNorm / float(shape(1)*shape(2)*shape(3)))

  return
end

