! --------------------------------------------------------------
! leapfrogA: use Leap Frog scheme to calculate A(n+1) from A(n) and A(n-1)
!
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

! Locals

  integer sh(3), nx,ny,nz
  real*8 dt,h

  sh = shape(an)

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

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

! Locals

  real*8 dt, h
  integer nx,ny,nz, sh(3)
  
  sh = shape(an)
  nx = sh(1)
  ny = sh(2)
  nz = sh(3)
  dt = params(P_TIMESTEP)
  h = params(P_GRIDSPACING)
  
  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(params, an, phip, phim)
  use input_params
  implicit none
  
  real*8, dimension(100) :: params
  real*8, dimension(:,:,:) :: an, phip, phim

! Locals	

  integer nx, ny, nz, sh(3)
  real*8 dt, h

  integer jerr, java_ser

  dt = params(P_TIMESTEP)
  h = params(P_GRIDSPACING)

  sh = shape(an)
  nx = sh(1)
  ny = sh(2)
  nz = sh(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(params, an, anp, bn, bnp, cn, cnp, dn, dnp, phip)
  use input_params
  implicit none

  real*8, dimension(100)   :: params
  real*8, dimension(:,:,:) :: an, anp, bn, bnp, &
     cn, cnp, dn, dnp, phip

  integer nx, ny, nz

  nx = int(params(P_NX))
  ny = int(params(P_NY))
  nz = int(params(P_NZ))

!   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(a)
  implicit none

  real*8, dimension(:,:,:) :: a

  integer, dimension(3) :: sh

  ! Locals

  integer nx,ny,nz

  sh = shape(a)

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

  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


subroutine swaptime(anm,an,anp, bnm,bn,bnp,&
     cnm,cn,cnp, dnm,dn,dnp, phim,phi,phip)
  implicit none
  
  real*8, dimension(:,:,:) :: anm,an,anp, bnm,bn,bnp, cnm,cn,cnp, &
       dnm,dn,dnp, phim,phi,phip
  
  ! Locals

  integer i,j,k, sh(3)

  real*8 r

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

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

  sh = shape(an)

!   Phim and phip have the same memory!
!   This swaps phip <-> phi and phim <-> phi

  do i=1, sh(1)
     do j=1, sh(2)
        do k=1, sh(3)
        
           r = phim(i,j,k)
           phim(i,j,k) = phi(i,j,k)
           phi(i,j,k) = r
           
        end do
     end do
  end do

  return  
end subroutine swaptime
