program upwindtest
c        SOLVE du/dt + du/dx = 0
c        ON A PERIODIC DOMAIN

         implicit double precision(a-h,o-z)
         dimension u(101),uexact(101),x(101)
         common/grid/dt,dx,dtdx
         common/index/nx

         amplit=0.25
         pi=4.0d0*datan(1.0d0)
         tmax=1.d0
         nx=101
         mid=((nx-1)/2) + 1
         dx=1.0d0/(nx-1)
         do 01 i=1,nx
           x(i)=(i-1)*dx
01       continue

c        INITIAL DATA
         freq=10.0d0
         t=0.0d0
         do 10 i=1,nx
c          uexact(i)=amplit*dsin(2.0d0*pi*freq*x(i))
           uexact(i)=0.0d0
10       continue
         do 12 i=(3*(nx-1)/10)+1,(7*(nx-1)/10)+1
           uexact(i)=amplit
12       continue
         do 15 i=1,nx
           u(i)=uexact(i)
15       continue

c        CFL ESTIMATE
         dt=0.80d0*dx

           do 999 niter=1,5
           do 998 ninner=1,200

         tnew=t+dt
         if (tnew .gt. tmax) then
           dt=tmax-t
           tnew=t+dt
         endif
         t=tnew

         dtdx=dt/dx

         call upstep(u)

998        continue
999        continue

         write(*,*)"tfinal = ",t

c        WRITE DATA
         open(unit=1,file='upwind')
         do 1110 i=1,nx
           write(1,1066)x(i),u(i)
1110     continue
         close(1)
         open(unit=1,file='uexact')
         do 1120 i=1,nx
           write(1,1066) x(i),uexact(i)
1120     continue
         close(1)
1066     format('',2(f10.5,1x))

         stop
         end
c********************************************************************
         subroutine upstep(u)
         implicit double precision(a-h,o-z)
         dimension u(101),unew(101)
         common/grid/dt,dx,dtdx
         common/index/nx

c        UPDATE u
         do i=2,nx-1
           unew(i)=u(i)-dtdx*(u(i)-u(i-1))
         enddo
c        BOUNDARY  i=1 <--> i=nx
         unew(nx)=u(nx)-dtdx*(u(nx)-u(nx-1))
         unew(1)=unew(nx)

c        RELABEL
         do i=1,nx
           u(i)=unew(i)
         enddo

         return
         end