program energy
c     Leap-frog advection with right boundary condition according
c     to Richtmyer and Morton, p. 141
c     Instability of linear extrapol. b.c. is only seen for very
c     long times, nt=3*ix/c=9/2*nt_outflow
      implicit NONE
      integer ix,nt,i,n,iop
      parameter (ix=50,nt=350)
      real um(ix),u(ix),up(ix),c
      open (unit=5,file='data',status='unknown')
c     iop=1 for Richtmyer boundary condition; =2 for linear extrapol.
      iop=2
c     Courant number c<1
      c=0.1
c     Initial conditions for square pulse
      do i=1,ix
      um(i)=0.
      u(i)=0.
      enddo
      do i=ix/3+1,2*ix/3
      um(i)=1.
      u(i)=1.
      enddo
c     Left boundary conditions
      u(1)=0.
      up(1)=0.
c     Time integration
      do n=0,nt
      if (n.eq.0) goto 1
c     Interior mesh
      do i=2,ix-1
      up(i)=um(i)-c*(u(i+1)-u(i-1))!+0.01*(u(i+1)-2.*u(i)+u(i-1))
      enddo
c     Right boundary condition
      if (iop.eq.1) then ! Richtmyer eq (6.37), shifted from J-1 to J
         up(ix)=((1.-c/2.)*um(ix)+c*u(ix-1))/(1.+c/2.)
      else               ! Linear extrapol
         up(ix)=2.*up(ix-1)-up(ix-2)
      endif
c     Updating
      do i=1,ix
      um(i)=u(i)
      u(i)=up(i)
      enddo
c     Output to file
 1    if (mod(n,nt/50).eq.0) then
      write (5,*) n, ' time steps'
      do i=1,ix
      write (5,*) i,u(i)
      enddo
      endif
      enddo ! Next time step
      close (5)
      stop
      end