
subroutine tstep(params, b_box, shape, t, anm, an, anp, bnm, bn, bnp, &
     cnm, cn, cnp, dnm, dn, dnp, phim, phi, phip, temp)
  use input_params
  implicit none

  real*8, dimension(100)     :: params
  integer, dimension(3) :: shape
  integer t
  real*8, dimension(6)  :: b_box
  real*8, dimension(:,:,:)  :: anm, an, anp, &
       bnm, bn, bnp, cnm, cn, cnp, dnm, dn, dnp, phim, phi, phip, temp


  integer nx,ny,nz
  real*8 dt,h

  nx = shape(1)
  ny = shape(2)
  nz = shape(3)
  dt = params(P_TIMESTEP)
  h = params(P_GRIDSPACING)

  write(6,*) '. '

  call dump_2d_slice("B", float(t), bn, nx,ny,nz, (ny+1)/2, &
       b_box(1), b_box(3), h)

  call boundupdt(shape, an,anp, bn,bnp, cn,cnp, dn,dnp, phip)

  call leapfrogA(shape, dt, h, anm,an,anp, bn,cn,dn)
  
  call leapfrogBCD(shape, dt, h, an, bnm,bn,bnp, cnm,cn,cnp, &
       dnm,dn,dnp)
  
  call updatephi(shape, dt, h, an, phip,phim)

  return  
end subroutine tstep

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

  integer shape(3)

  real*8 dt, h

  real*8, dimension(shape(1), shape(2), shape(3)) :: anm,an,anp, bn,cn,dn

  integer nx,ny,nz

  nx = shape(1)
  ny = shape(2)
  nz = shape(3)		

  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(shape, dt, h, an, bnm, bn, bnp, &
     cnm, cn, cnp, dnm, dn, dnp)			
  implicit none

  integer shape(3)

  real*8 dt, h

  real*8, dimension(shape(1), shape(2), shape(3)) :: an, bnm, bn, bnp, &
       cnm, cn, cnp, dnm, dn, dnp

  integer nx,ny,nz

  nx = shape(1)
  ny = shape(2)
  nz = shape(3)
  
  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(shape, dt, h, an, phip, phim)
  implicit none
  
  integer shape(3)

  real*8 dt, h

  real*8, dimension(shape(1), shape(2), shape(3)) :: an, phip, phim

  integer nx, ny, nz

  nx = shape(1)
  ny = shape(2)
  nz = shape(3)

  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(shape, an, anp, bn, bnp, cn, cnp, dn, dnp, phip)
  implicit none

  integer shape(3)

  real*8, dimension(shape(1), shape(2), shape(3)) :: an, anp, bn, bnp, &
     cn, cnp, dn, dnp, phip

  integer nx, ny, nz

  nx = shape(1)
  ny = shape(2)
  nz = shape(3)

!   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(:,:,:) :: 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
