program burgnc

c     Solve the Burger's equation problem:
c     u_t + u u_x = 0     u(x,0) = 1 + sin(x)    u(x,t) = u(x + 2 pi,t)

c     This program uses a nonconservative upwind scheme

c     Compile using: f77 -o burgnc burgnc.f burgex.f

c     Run using: burgnc

c     Plot results using:
c         ~dommelen/bin/xyplot burgnc.plt 2 ' ' ' ' ' ' 2 '0,100,0,100,.4'
c     which creates a printable postscript file named dump.ps

c     To use your own plotting software, change the format of file burgnc.plt
c     below as needed.

      implicit none

c Declarations

c     average initial velocity
      real u0

c     time step index n
      integer n
c     index at which the velocities u^n are actually stored
      integer n_s
c     index at which the velocities u^{n+1} are actually stored
      integer np1_s
c     next value of n at which to print out diagnostics
      integer nprint

c     time
      real t

c     declared number of mesh cells
      integer j_dim
      parameter (j_dim=1000)
c     actual number of mesh cells
      integer j_max
c     mesh point index
      integer j

c     x-positions of the mesh points
      real x(0:j_dim)

c     velocity values.  only two time levels are stored at any point
      real u(-1:j_dim,0:1)

c     mesh cell size Delta x and time step Delta t
      real dx,dt

c     maximum Courant number C = u Delta t / Delta x
      real C

c     output time interval
      real t_out
c     output time index and number of output times
      integer out,outmx

c     exact solution, for comparing to
      real burgex
      external burgex

c     pointwise and rms errors in the solution
      real err,errrms

c     minimum and maximum velocity
      real u_min,u_max

c     constants
      real pi

      pi=4.*atan(1.)

c executable statements

c open files

c     input file
      open(1,file='burg.in',status='old')

c     output file
      open(2,file='burgnc.out',status='unknown')

c     plot file
      open(3,file='burgnc.plt',status='unknown')

c read input file and write to output

      read(1,*)u0
      write(2,*)'average initial velocity: ',u0

      read(1,*)j_max
      write(2,*)'number of mesh cells j_max: ',j_max
      if(j_max.gt.j_dim)then
         print*,'*** declared dimension exceeded'
         stop
      endif
      if(j_max.lt.2)then
         print*,'*** number of mesh cells should be at least 2'
         stop
      endif

      read(1,*)C
      write(2,*)'maximum Courant number u Delta t/Delta x: ',C
      if(C.le.0.)then
         print*,'*** Courant number must be positive'
         stop
      endif

      read(1,*)t_out
      write(2,*)'output time interval: ',t_out
      if(t_out.le.0.)then
         print*,'*** output time interval must be positive'
         stop
      endif

      read(1,*)outmx
      write(2,*)'number of output time intervals to compute: ',outmx
      if(outmx.le.0)then
         print*,'*** number of output times must be positive'
         stop
      endif

      close(1)

c initialize the parameters

c     mesh size
      dx=2.*pi/j_max

c     output count
      out=1

c     print value of n
      nprint=1

c initial velocities

c     set all initial velocities and x-positions
      do 100 j=0,j_max
         x(j)=j*dx
         u(j,0)=u0+sin(x(j))
 100  continue

c     initialize time
      t=0.

c main loop over all time steps

      do 2000 n=0,123456789

c        where to store the velocities at n and n+1
         n_s=mod(n,2)
         np1_s=mod(n+1,2)

c        perform output

 1010    if(t.lt.out*t_out) goto 1090

c        write computed velocities to the plot file
         write(3,*)-j_max-1
         write(3,1050)(x(j),u(j,n_s),j=0,j_max)
 1050    format(f12.7,',',f12.7)

c        write the exact solution to the plot file as a 201 point curve
         write(3,*)201,'    ',j_max,t
         write(3,1050)(0.005*2.*pi*j,
     &    burgex(0.005*2.*pi*j,t,u0),
     &    j=0,200)

         out=out+1
         if(t.ge.out*t_out)goto 1010
 1090    if(out.gt.outmx)goto 2010

c        time step
         dt=1.e30
         do 1200 j=0,j_max-1
            if(u(j,n_s).ne.0.)dt=min(dt,C*dx/abs(u(j,n_s)))
 1200    continue

c        increment time
         t=t+dt

c        set the j=-1 value of u using periodicity
         u(-1,n_s)=u(j_max-1,n_s)

c        compute the velocities
         do 1400 j=0,j_max-1
            if(u(j,n_s).ge.0.)then
               u(j,np1_s)=u(j,n_s)-dt/dx*u(j,n_s)*
     &          (u(j,n_s)-u(j-1,n_s))
            else
               u(j,np1_s)=u(j,n_s)-dt/dx*u(j,n_s)*
     &          (u(j+1,n_s)-u(j,n_s))
            endif
 1400    continue

c        set the boundary velocity at j_max using periodicity
         u(j_max,np1_s)=u(0,np1_s)

c        print diagnostics only occasionally
         if(n+1.ne.nprint)goto 2000
         nprint=nprint*2

c        compute the rms error and minimum and maximum velocity
         errrms=0.
         u_min=u(0,np1_s)
         u_max=u(0,np1_s)
         do 1500 j=0,j_max-1
            err=abs(u(j,np1_s)-burgex(x(j),t,u0))
            errrms=errrms+err**2
            u_min=min(u_min,u(j,np1_s))
            u_max=max(u_max,u(j,np1_s))
 1500    continue
         errrms=sqrt(errrms/j_max)
         print1510,n+1,t,u_min,u_max,errrms
         write(2,1510)n+1,t,u_min,u_max,errrms
 1510    format('n =',i3,'   t =',f7.4,
     &    '   u-range:',f9.5,' to',f8.5,'   error:',e10.3)

 2000 continue
 2010 continue

      stop
      end