	implicit double precision (a-h,o-z)
	parameter(M=200)			! number of grid points
	parameter(EL=10.,D=2E-9*1.E+4*3600.*24.) ! use cm, day
	parameter(dz=EL/M,dt=0.01,xi=D*dt/2./dz/dz,CA0=0.01)
	parameter(Tmax=100,NT=9)
	dimension a(M),b(M),c(M),r(M),CA(0:M),CAold(0:M)
	dimension Tmark(NT),result(0:M,NT)
	data Tmark/0.,0.5,1.,5.,10.,20.,40.,60.,80./
c
c	Initial condition. Also solution on old time step
c	
	do i=1,M
	  CAold(i)=0.
	enddo
	CAold(0)=CA0
c
c	Build tri-diagonal matrix from Crank-Nicolson scheme
c	
c	print*,'D,dt,dz,xi=',D,dt,dz,xi

	do i=1,M
	  a(i)=xi
	  b(i)=-(1.+2.*xi)
	  c(i)=xi
	enddo
	a(M)=2.*xi					! Right BC
	
	do i=0,M
	  result(i,1)=CAold(i)
	enddo
	T=0.
	iT=2
	idt=0
	
	do while (T.lt.Tmax+dt/2.)
	
	r(1)=-2.*xi*CAold(0)+(2.*xi-1.)*CAold(1)-xi*CAold(2) ! left BC
	do i=2,M-1
	  r(i)=-xi*CAold(i-1)+(2.*xi-1.)*CAold(i)-xi*CAold(i+1)
	enddo
	r(M)=-xi*CAold(M-1)+(2.*xi-1.)*CAold(M)-xi*CAold(M-1) ! right BC
	
	call tridag(a,b,c,r,CA(1),M)
	CA(0)=CA0
	do i=0,M
	  CAold(i)=CA(i)
	enddo
	T=T+dt
	idt=idt+1
c	print*,'Time step idt=',idt,' T=',T
	
	if(abs(T-Tmark(iT)).lt.1.E-5) then
	  do i=0,M
	    result(i,iT)=CA(i)
	  enddo
	  iT=iT+1
	endif
	
	enddo
	
	open(12,file='diffus.res')
	do i=0,M
	  write(12,100) i*dz,(result(i,j),j=1,NT)
	enddo
	close(12)
100	format(<NT+1>E14.6)

	end
	
      SUBROUTINE tridag(a,b,c,r,u,n)
      implicit double precision (a-h,o-z)
      INTEGER n,NMAX
      double precision a(n),b(n),c(n),r(n),u(n)
      PARAMETER (NMAX=500)
      INTEGER j
      double precision bet,gam(NMAX)
      if(b(1).eq.0.)pause 'tridag: rewrite equations'
      bet=b(1)
      u(1)=r(1)/bet
      do 11 j=2,n
        gam(j)=c(j-1)/bet
        bet=b(j)-a(j)*gam(j)
        if(bet.eq.0.)pause 'tridag failed'
        u(j)=(r(j)-a(j)*u(j-1))/bet
11    continue
      do 12 j=n-1,1,-1
        u(j)=u(j)-gam(j+1)*u(j+1)
12    continue
      return
      END
