PROGRAM DIFFUSION

	implicit double precision (A-H,O-Z)

	dimension C(100,50)

	rl=20.

	NMAX=20

	JSAVE=5

      pi=3.141592654

	dx=1.

	LMAX=40

	DC=1.

	D=1.

	t=0.

	DT=.25

	do 40 J=1, 20

	x=0

	C(1,j)=1.

	do 10 I=2,nmax

	 X=X+DX

             TE1=DC*(1.-x/rl)

	         sum=0.
         
	        do 20 N=1,LMAX

                  dsum=DSIN(N*PI*X/RL)*Dexp(-N*N*pi*pi*D*T/RL**2.)/N


                    SUM=SUM+DSUM
20            continue

              c(i,J)=te1-2.*DC*sum/pi

	     
10        continue


        IF (J.EQ.JSAVE) THEN

	  open(10,FILE='c:\numerics\diffusio.res',status='unknown')
	  
	    write (10,62) T
		
 62        format (1x, 'T=',f10.2)
 
            x=0.
		  
		  do 27 i=1,nmax
		    write (10,63) X,C(i,J)
63          format(1x, f10.2,f10.5)

             

 27         continue
          close(10)
            stop

           endif



	t=T+DT

40     continue

	end