      program burgers
c---- 67-------------------------------------------------------flowtec-|
c     solves burgers equation u,t+u*u,x=0 numerically in
c     one dimension by a finite-volume  using:
c     ischeme           fluxes
c     
c     3                 roe      (lin. riemann solver, flux-diff-split)
c---- 67-------------------------------------------------------flowtec-|
c     30.05.01
c     BUGS: -?-
c
c     TODO:
c     (i)   second order scheme (muscl)
c     (ii)  add some source term 
c---- 67---------------------------------------------------------------|
      implicit double precision (a-h,o-z)
      parameter(nmax=400)
      dimension q(0:nmax+1),rhs(0:nmax+1) !counting 'ghost-cells'
      parameter(cfl=.4d0)
c---- 67---------------------------------------------------------------|
      parameter(ull=9.d0,urr=10.d0)
c     /* flx:   length of domain; tfin: final integration time */
c     /* nitmx: max. number of time steps                      */
      parameter(flx=1.d0,tfin=.1d0,nitmx=50000)
      parameter(ischeme=3,icase=2)
c---- 67---------------------------------------------------------------|
      time=0.d0
      nit=0
      dx=flx/dfloat(nmax)
      do i=0,nmax+1
         rhs(i)=0.d0
      enddo

      call set_initial_state(q,ull,urr,flx,nmax,icase)

c     /* main time loop */
 11   continue
      nit=nit+1
      dt=1.d35

      call set_boundary_conditions(q,ull,urr,flx,nmax,icase) 

c     /* flux balance for each cell-interface */
      do 20 icell=0,nmax

c     /* solve interfacial riemann problem analytically */
         if(ischeme.eq.1)then
            write(*,*)'not implemented'
            
         elseif(ischeme.eq.2)then
            write(*,*)'not implemented'               

         elseif(ischeme.eq.3)then

            call roeflux(q(icell),q(icell+1),flux,flamx)

            if(icell.gt.1) call centered_source(q(icell-1),q(icell),
     $           q(icell+1),source)

         else
            write(*,*)'small problem here...',ischeme
         endif

c     /* calculate maximum GLOBAL time step */
         if(flamx.gt.0.d0)then
            dt=dmin1(dt,dx*cfl/flamx)
         endif
         
c     /* accumulate flux difference */
         rhs(icell)=rhs(icell)+flux
         rhs(icell+1)=rhs(icell+1)-flux

c     /* add-on the source term locally */
         rhs(icell)=rhs(icell)+source

 20   continue

      time=time+dt
      write(*,*)'time = ',time, ' dt= ',dt, ' tfin = ',tfin
      if(time.gt.tfin.or.nit.gt.nitmx)then
         call output(q,time,flx,nmax)
         goto 10
      endif

c     /* update the variables */
      do ival=1,nmax
         q(ival)=q(ival)-dt/dx*rhs(ival)
         rhs(ival)=0.d0
      enddo

      goto 11
 10   continue
      stop 
      end
c----67-------------------------------------------------------flowtec-|
      subroutine output(q,time,flx,ni)
      implicit double precision (a-h,o-z)
      dimension q(ni)
      open(12,file='sol.dat')
      write(12,100)time,ni
      do i=1,ni
         write(12,101)dfloat(i-1)/dfloat(ni-1)*flx,
     $        q(i)
      enddo
      close(12)
 100  format('# 1:x 2:u ',x,2(E12.5,x))
 101  format(1x,2(E12.5,x))
      return
      end
c----67---------------------------------------------------------------|      
      subroutine roeflux(ul,ur,rflux,flamx)
c----67-------------------------------------------------------flowtec-|
c     solves the linearized riemann problem for st.venant equations 
c     (flat bottom) by using Roe's method (J.Comp.Physics, vol.43,1981)
c     * fluxes are passed upon exit in vector rflux(1:2)
c----67---------------------------------------------------------------|
c     21.05.99
c     BUGS: -?-
c----67---------------------------------------------------------------|
      implicit double precision (a-h,o-z)
      dimension a(1),cflux(1),uflux(1),diffv(1),rflux(1)
      logical lvacuum
      parameter (f23=2.d0/3.d0,f13=1.d0/3.d0,f19=1.d0/9.d0)
      parameter (small=1.d-6)
c----67---------------------------------------------------------------|      
      cl=ul
      cr=ur

c     /* roe's average for burgers */
      uc=(ul+ur)/.5d0
      cc=uc

c     /* diagonalization |A|=R*|lambda|*R1 */
      call aroe(uc,cc,a)

c     /* variable difference vector */
      diffv(1)=ur-ul

c     /* centered flux */
      cflux(1)=dflux(ur)+dflux(ul)

c     /* evaluate roe's flux: f=(fl+fr-|a|*(qr-ql))/2 */
      rflux(1)=(cflux(1)-a(1)*diffv(1))/2.d0

c     /* estimate the maximum signal velocity */
      flamx=dabs(cc)
c
      return
      end
c----67---------------------------------------------------------------|      
      subroutine aroe(uc,cc,a)
c----67-------------------------------------------------------flowtec-|
c     calculates roe's matrix
c----67---------------------------------------------------------------|
      implicit double precision (a-h,o-z)
      dimension a(1),flambda(1)
c----67---------------------------------------------------------------|      
c     /* eigenvalues */
      flambda(1)=uc

c     /* matrix*diag*matrix multiply */
      a=dabs(flambda(1))
c
      return
      end
c----67---------------------------------------------------------------|      
      double precision function dflux(u)
      implicit double precision (a-h,o-z)
c
      dflux=u*u/2.d0
c      
      end
c----67---------------------------------------------------------------|      
      subroutine centered_source(uim1,ui,uip1,source)
c----67-------------------------------------------------------flowtec-|
c     computes a centered source term for a scalar problem
c     u,t+c(u)*u,x=s(u)=d(u)*u,x
c     -->  \hat{s}=d(u_i)*[u_{i+1}-u_{i-1}]/(2dx)
c----67---------------------------------------------------------------|
c     21.05.99
c     BUGS: -?-
c----67---------------------------------------------------------------|
      implicit double precision (a-h,o-z)
c----67---------------------------------------------------------------|      
c     /* quasi-linear coefficient d(u_i) */
      call quasi_linear_source(ui,d) ! --> d

c     /* the numerical source term */
      source=d*(uip1-uim1)/2.d0 ! the 'dx' is actually factored later
c
      return
      end
c----67---------------------------------------------------------------|      
      subroutine quasi_linear_source(ui,d)
c----67-------------------------------------------------------flowtec-|
c     calculates the coefficient of a quasi-linear formulation of 
c     a scalar source term
c----67---------------------------------------------------------------|
      implicit double precision (a-h,o-z)
c----67---------------------------------------------------------------|      
      d=-ui
c
      return
      end
c----67---------------------------------------------------------------|      
      subroutine set_initial_state(u,ul,ur,flx,n,iflag)
c----67-------------------------------------------------------flowtec-|
c     set the initial state for the variable vector/scalar
c----67---------------------------------------------------------------|
      implicit double precision (a-h,o-z)
      dimension u(0:n+1)
c----67---------------------------------------------------------------|      
      if(iflag.eq.1)then        !riemann problem, i.e. discontinuity
         do i=1,n/2
            u(i)=ul
         enddo
         do i=n/2+1,n
            u(i)=ur
         enddo
      elseif(iflag.eq.2)then    !continuos initial data
         do i=1,n
            x=float(i-1)/float(n-1)*flx
            u(i)=ur
         enddo
      elseif(iflag.eq.3)then 
         do i=1,n/4
            u(i)=1
         enddo
         do i=n/4+1,n
            x=float(i-1)/float(n-1)*flx
            u(i)=x+1
         enddo
      elseif(iflag.eq.4)then 
         xc=flx/2.d0
         sigma=.05
         do i=1,n
            x=float(i-1)/float(n-1)*flx
            u(i)=exp(-((x-xc)/sigma)**2)
         enddo
      endif

c
      return
      end
c----67---------------------------------------------------------------|      
      subroutine set_boundary_conditions(u,ul,ur,flx,n,iflag)
c----67-------------------------------------------------------flowtec-|
c     set the boundary conditions for the variable vector/scalar
c     according to the condition in set_initial_state
c----67---------------------------------------------------------------|
      implicit double precision (a-h,o-z)
      dimension u(0:n+1)
c----67---------------------------------------------------------------|      
      if(iflag.eq.1)then        !riemann problem, i.e. discontinuity
         if(u(1).ge.0.0)then
            u(0)=ul
         else
            u(0)=u(1)
         endif
         if(u(n).le.0.0)then
            u(n+1)=ur
         else
            u(n+1)=u(n)
         endif
      elseif(iflag.eq.2)then    !step from boundary cond
         if(u(1).ge.0.0)then
            u(0)=ul
         else
            u(0)=u(1)
         endif
         if(u(n).le.0.0)then
            u(n+1)=ur
         else
            u(n+1)=u(n)
         endif
      elseif(iflag.eq.3)then 
         do i=1,n/4
            u(i)=1
         enddo
         do i=n/4+1,n
            x=float(i-1)/float(n-1)*flx
            u(i)=x+1
         enddo
      elseif(iflag.eq.4)then 
         xc=flx/2.d0
         sigma=.05
         do i=1,n
            x=float(i-1)/float(n-1)*flx
            u(i)=exp(-((x-xc)/sigma)**2)
         enddo
      endif

c
      return
      end
c----67---------------------------------------------------------------|      
