C   prog2b.f -- Solve viscous Burgers Equation
C   Written Apr 6, 2003 by Eric Olson for Mathematics 484
C
C   For information about the central finite difference scheme used
C   here see chapter 7 in Numerical Analysis for Applied Mathematics
C   Science and Engineering by Greenspan and Casulli.  Information
C   about the heat and Burgers equations may be found in 13 and 14 in
C   A First Course in the Numerical Analysis of Differential Equations
C   by Iserles.

      program main
      implicit none

      include "param2b.h"
C   param2b.h contains the following lines uncommented
C     integer N
C     parameter(N=16)

      real*8 u,t,dx,m_pi
      integer i
      common /soln/u(N)
      common /c_dx/dx
      data m_pi/3.14159265358979323D0/
      dx=2*m_pi/N
      write(*,10)
10    format(1X,'#prog2b.f -- Compute viscous Burger''s Equation')
      call doinit(u)
      call vectwr(u,0.0D0)
      t=0.0
      do i=1,8
        call euler(u,0.125D0,4096)
        call vectwr(u,i*0.125D0)
      enddo
      stop
      end

      subroutine doinit(u)
      implicit none
      include "param2b.h"
      real*8 u,dx,x,xx,m_pi
      integer i
      dimension u(N)
      data m_pi/3.14159265358979323D0/
      common /c_dx/dx
      do i=1,N
        x=dx*(i-1)
        xx=M_PI-x 
        xx=xx*xx
        u(i)=8*exp(-xx)
      enddo
      end

      subroutine force(f,u)
      implicit none
      include "param2b.h"
      real*8 f,u,dx
      common /c_dx/dx
      dimension f(N),u(N)
      integer i,ip1,im1
      do i=1,N
         ip1=mod(i,N)+1
         im1=mod(i+N-2,N)+1
         f(i)=(-u(i)*(u(ip1)-u(im1))+(u(im1)-2*u(i)+u(ip1))/dx)/dx
      enddo
      end

      subroutine euler(u,tn,m)
      implicit none
      include "param2b.h"
      real*8 u,f,tn,dt
      integer m,i,j
      dimension u(N)
      common /fvect/f(N)
      dt=tn/m
      do j=1,m
        call force(f,u)
        do i=1,N
            u(i)=u(i)+dt*f(i)
        enddo
      enddo
      end

      subroutine vectwr(u,t)
      implicit none
      include "param2b.h"
      real*8 u,t,dx
      integer i,k
      dimension u(N)
      common /c_dx/dx
      write(*,20) t
20    format(1X,//,1X,'#At time ',F8.5,' ...')
      do i=1,N
        write(*,30) dx*(i-1),u(i)
30      format(1X,2E14.5)
      enddo
      end
