program sample2

c  Kreider 5/12/08

c  add in file IO commands

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,0:1),x(1001),t(100001),dx,dt,c,f,h
      double precision Tmax,r,ux
      integer i,n,m,k

      character*1 cflag
      character*13 flname
      character*3 chan3
      integer np

      write(6,*) 'Enter id flag (1 character)'
      read(5,77) cflag
77    format(a1)

      np = 10000

      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
      k = 0
      call ooutfn('smp2',chan3(k),cflag,1)
      do i=1,n
         u(i,0) = h(x(i))
         write(1,*) x(i),u(i,0)
      end do
      close(1)

c   main loop
c   u(i,0) is old value u_i^k
c   u(i,1) is new value u_i^{k+1}
      do k=1,m-1
         u(1,1) = f(t(k+1))
         do i=2,n-1
            ux = u(i,0) - u(i-1,0)
            u(i,1) = u(i,0) - r*ux
         end do
         u(n,1) = 2.d0*u(n-1,1) - u(n-2,1)
         if (np*((k+1)/np) .eq. k+1) then
            call ooutfn('smp2',chan3((k+1)/np),cflag,1)
            write(6,*) k+1,'   ',chan3((k+1)/np)
            do i=1,n
               write(1,*) x(i),u(i,1)
            end do
            close(1)
         end if
         do i=1,n
            u(i,0) = u(i,1)
         end do
      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

      subroutine ooutfn(fpr,fid,cflag,lun)
c  usage
c      type *,' Enter file id'
c      accept (2a), id
c      call ooutfn(id,'i',1) creates output file i'id'.dat
c      taken from Nazanin Imani

      character*3 fid
      character*1 cflag
      character*4 fpr
      integer lun
      character*20 flname

      flname = fpr // cflag // fid
      open(unit=lun,file=flname,err=150)
      return
150     write(6,*),'*** ERROR IN OPENING FILE ***'
      return
      end

      function chan3(M)
c     Stephen Cardarelli, Oct 2003
c     provides time step files 000 to 999
      implicit none
      character*3 chan3
      integer i,j,k,M,countout

      countout = M
      i = countout/100
      countout = countout - (i*100)
      j = countout/10
      k = countout - (j*10)

      chan3 = char(i+48)//char(j+48)//char(k+48)
      return
      end