
! --------------------------------------------------------------
! initial: Put down a Gaussian as initial data
!
subroutine initial(params, a0, b0, c0, d0, phi0)
  use input_params
  implicit none

  real*8, dimension(100)     :: params
  real*8, dimension (:,:,:) :: a0, b0, c0, d0, phi0
  INTENT(OUT) :: a0, b0, c0, d0, phi0

  ! Locals	
     
  integer nx, ny, nz
  real*8 h, vx, vy, vz, width, ampl, lb(3)

  integer i, j, k, sh(3)
  real*8 x, y, z, t, rhop, rhom

  sh = shape(a0)

  nx = sh(1)
  ny = sh(2)
  nz = sh(3)
  h = params(P_GRIDSPACING)
  vx = params(P_VX)
  vy = params(P_VY)
  vz = params(P_VZ)
  width = params(P_WIDTH)
  ampl = params(P_AMPL)
  lb(1) = params(P_LOCAL_LOWER_X)
  lb(2) = params(P_LOCAL_LOWER_Y)
  lb(3) = params(P_LOCAL_LOWER_Z)

  t = 0.0d0

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

           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) + params(P_XMIN)
           y=h*(j+lb(2)-1) + params(P_YMIN)
           z=h*(k+lb(3)-1) + params(P_ZMIN)

           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


! --------------------------------------------------------------
! initstep: average different time levels
!
subroutine initstep(params, anm, an, anp, bnm, bn, bnp, cnm, cn, cnp, &
     dnm, dn, dnp, temp)
  use input_params
  implicit none

  INTERFACE
     
     subroutine leapfrogA(params, anm,an,anp, bn,cn,dn)
       use input_params
       implicit none
       real*8, dimension(100) :: params
       real*8, dimension(:,:,:) :: anm,an,anp, bn,cn,dn
     end subroutine leapfrogA
     
     subroutine leapfrogBCD(params, an, bnm, bn, bnp, cnm, cn, cnp,&
          dnm, dn, dnp)
       use input_params
       implicit none
       real*8, dimension(100):: params
       real*8, dimension(:,:,:) :: an, bnm, bn, bnp, cnm, cn, cnp, dnm, dn, dnp
     end subroutine leapfrogBCD

     real*8 function DiffNorm(a, b)
       implicit none
       real*8, dimension(:,:,:) :: a, b
     end function DiffNorm

  END INTERFACE

  real*8, dimension (100) :: params
  real*8, dimension(:,:,:), target :: anm, an, anp, bnm, bn, bnp, &
       cnm, cn, cnp, dnm, dn, dnp, temp

  ! Locals

  integer maxstp, i
  real*8 anorm, bnorm, cnorm, dnorm
  integer nx, ny, nz, sh(3)
  real*8 vx, vy, vz, width, ampl, dt, h, epsilon
  real*8, dimension(:,:,:), pointer :: anpnew, bnpnew, cnpnew,&
       dnpnew, antemp

  sh = shape(an)
  nx = sh(1)
  ny = sh(2)
  nz = sh(3)
  h = params(P_GRIDSPACING)
  dt = params(P_TIMESTEP)
  vx = params(P_VX)
  vy = params(P_VY)
  vz = params(P_VZ)
  width = params(P_WIDTH)
  ampl = params(P_AMPL)
  epsilon = params(P_INITEPSILON)

  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(params, anm, an, anpnew, bn,cn,dn)

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

     anorm = DiffNorm (anpnew, anp)
     bnorm = DiffNorm (bnpnew, bnp)
     cnorm = DiffNorm (cnpnew, cnp)
     dnorm = DiffNorm (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

  ! Leave temp clean for the next subprogram to use it

  temp(:,:,:) = 0.0d0

  return
end subroutine initstep

	

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

  ! Locals

  integer i,j,k,sh(3)

  sh=shape(a)
  DiffNorm = 0.0d0

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

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

        end do
     end do
  end do

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

  return
end

