program sample1

c  Kreider 5/12/08
c  sample code to show typical structures
c  problem:  u_t + c u_x = 0
c            u(x,0) = h(x)
c            u(0,t) = f(t)

      implicit none
      double precision u(1001,100001),x(1001),t(100001),dx,dt,c,f,h
      double precision Tmax,r,ux
      integer i,n,m,k

      c = .5d0
      dx = 1.d-2
      dt = 1.d-5
      r = c*dt/dx
      Tmax = 1.d0

c        BELOW
c        add 1.1 to account for occasional roundoff error
c        the right side is computed in double precision, then
c        converted to an integer
      n = (1.d0-0.d0)/dx + 1.1
      m = Tmax/dt + 1.1

      write(6,*) n,m,r

c   grids
      do i=1,n
         x(i) = 0.d0 + dfloat(i-1)*dx
      end do
      do k=1,m
         t(k) = 0.d0 + dfloat(k-1)*dt
      end do

c   initial condition
      do i=1,n
         u(i,1) = h(x(i))
      end do

c   main loop
      do k=1,m-1
         u(1,k+1) = f(t(k+1))
         do i=2,n-1
            ux = u(i,k) - u(i-1,k)
            u(i,k+1) = u(i,k) - r*ux
         end do
         u(n,k+1) = 2.d0*u(n-1,k+1) - u(n-2,k+1)
      end do

c   write out first 10 for assignment in 9 and all for graphing in 10
      open(9,file='sample1.out')
      open(10,file='sample1.all')
      do i=1,10
         write(9,*) x(i),u(i,m)
      end do
      do i=1,n
         write(10,*) x(i),u(i,m)
      end do

      stop
      end

      function f(t)
      implicit none
      double precision f,t
      f = 1.d0 + sin(t)
      return
      end

      function h(x)
      implicit none
      double precision h,x
      h = exp(-4.d0*(x-.2d0)**2)
      return
      end