program leapfrog
      parameter(ni=21,ktime=20)
      dimension psiold(ni), psi(ni), psinew(ni), x(ni)
!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
!
!   Written by Dr. Pat Fitzpatrick
!   Assistant Professor of Meteorology
!   Department of Physics, Atmospheric Sciences, & General Science
!   Jackson State University
!   10/15/97
!
!   This program timesteps the 1D advection equation using the
!   leapfrog scheme for a given timestep, grid spacing, and
!   constant zonal wind
!
!   Variables:
!
!   ni is the number of grid points
!   u --- constant zonal wind speed
!   delx --- grid spacing
!   delt --- time step
!   ktime is the number of timesteps
!   psiold(i) is the old psi value
!   psi(i) is the current psi value
!   psinew(i) is the future psi value
!   x(i) --- x value at each grid point, defiend by (i-1)*delx
!&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
!
!      open(unit=20,file='fcst.dat',status='unknown')
!************************************************************
!   Define constant variables
!************************************************************
      pi=3.141592654
      delt=100
      delx=1000
      u=10
!************************************************************
!   Define x(i)
!************************************************************
      do i= 1,ni
         x(i)=delx*(i-1)
      enddo
!************************************************************
!  Initialize psi. Since we are using leapfrog, must know psiold
!  and psi to start time stepping
!************************************************************
      do i=1,ni
            psiold(i)=10.0*cos((2.0*pi*x(i))/(10.0*delx))
            psi(i)=10.0*cos((2.0*pi*(x(i)-u*delt))/(10.0*delx))
      enddo
!************************************************************
!  Output initial field (in integer format)
!************************************************************
      write(*,1000) (int(psiold(i)),i=1,ni-1)
      write(*,1000) (int(psi(i)),i=1,ni-1)
 1000 format(20(i3,1x))
!************************************************************
!  Begin timestepping
!************************************************************
      do 100 k=1,ktime
      do i=2,ni-1
         psinew(i)=psiold(i)-(u*delt/delx)*(psi(i+1)-psi(i-1))
      enddo
!************************************************************
!  Apply periodic boundary conditions
!************************************************************
         psinew(ni)=psiold(ni)-(u*delt/delx)*(psi(2)-psi(ni-1))
         psinew(1)=psinew(ni)
!************************************************************
!  Output field
!************************************************************
      write(*,1000) (int(psinew(i)),i=1,ni-1)
!************************************************************
!   Reset psi for next time step
!************************************************************
      do i=1,ni
         psiold(i)=psi(i)
         psi(i)=psinew(i)
      enddo
 100  continue
      end