program waveimplicit
C
C example of implicit solutions to  
C PDE  du/dx = - du/dt
C 
C
C  CWJ SDSU April 2005
C
      implicit none

C......... parameters of box..............

      real L  !  box is from 0 to L
      real dx !  discretiztion of box
      integer Nx !  number of points in box

C......... parameters of time.............

      real dt !  time step
      integer Nt !  number of time steps to take
      real t
      integer n

C.....parameters of wave..................

      integer size
      parameter(size=1000)
      real u1(size),u2(size)

C.......... parameters of initial gaussian

      real xc,xw    ! center and width
      integer j

C............READ IN PARAMETERS OF BOX...........

      print*,' Enter size of box L '
      read*,L
      print*,' Enter # of points '
      read*,Nx
      dx = L/float(nx)      

C...........read in initial parameters of gaussian wave 

      print*,' Enter center, width of initial gaussian '
      read*,xc,xw

C...........fill in gaussian.........................

      do j = 2,Nx-1
         u1(j) = exp(- (dx*j - xc)**2/(2*xw**2))
      enddo
      u1(1) = 0.
      u1(Nx) = 0.0

      open(unit=1,file='initial.dat',status='unknown')
      do j = 1,Nx
        write(1,*)j*dx,u1(j)
      enddo
      close(unit=1)

C...........read in parameters of time step.......

      print*,' Enter time step '
      print*,' (suggest less than ',dx,') '
      read*,dt       

      print*,' Enter time to travel  '
      read*,t
      Nt = int(t/dt)

      do n = 1,nt/2
         call implicit_step(size,Nx,dx,dt,u1,u2)
         call implicit_step(size,Nx,dx,dt,u2,u1)
      enddo
      
      open(unit=2,file='final.dat',status = 'unknown')
      do j = 1,Nx
        write(2,*)j*dx,u1(j)
      enddo
      close(unit=2)      
      open(unit=3,file='compare.dat',status = 'unknown')
      xc = xc-(nt/2)*2*dt
      do j = 1,Nx
        write(3,*)j*dx,exp(- (dx*j - xc)**2/(2*xw**2))
      enddo
      end
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      subroutine implicit_step(size,Nx,dx,dt,uold,unew)
C
C  does a single implicit step
C
      implicit none

      integer size
      real uold(size),unew(size)
   
      integer Nx
      real dx
      real dt

      integer size2
      parameter(size2=1000)
      real alo(size2),a0(size2),ahi(size2)
      real work(size2)

      integer j

      if(size2.ne.size)then
         print*,' mismatch in dimensions ',size,size2
         stop
      endif
      do j = 1,Nx
        a0(j) = 1.
        alo(j) = + dt/(2.*dx)
        ahi(j) = - dt/(2.*dx)
      enddo
      call tridiag(size,nx,alo,a0,ahi,uold,work,unew)

      unew(1) = 0. ! 
      unew(Nx) =0. ! 

      return
      end
  

      subroutine tridiag(size,n,aa,adiag,ac,b,work,x)
C
C  solves A x = b where A is tridiagonal
C  adiag = diagonal
C  aa  = lower off-diagonal
C  ac  = upper off diagonal
C
C  work(i) is temporary array
C
C  based on numerical recipes sect 2.4
C
C  

      implicit none
      integer size
      integer N
      real aa(size),adiag(size),ac(size)
      real b(size)
      real x(size)

      integer i
 
      real beta,work(size)

      if(adiag(1).eq.0.)then
         print*,1,adiag(1),' tridiag fails '
         stop
      endif

      beta = adiag(1)
      x(1) = b(1)/beta

      do i = 2,N
          work(i) =ac(i-1)/beta
          beta  = adiag(i)-aa(i)*work(i)
          if(beta.eq.0.)then
            print*,i,beta,' tridiag fails '
            stop
          endif
          x(i) =(b(i)-aa(i)*x(i-1))/beta
       enddo
C.............. backsubstitution..............

       do i = n-1,1,-1
           x(i) = x(i) - work(i+1)*x(i+1)
       enddo
       return
       end