C     ***************************************************
C     **  1-D DIFFUSION EQUATION (DU/DT+D/DX(DU/DX)=0) **
C     **                                               **
C     **  IMPLICIT METHOD                              **
C     **  EULER BACKWARD METHOD                        **
C     **                                               **
C     **                                               **
C     **  2000/01/12  VERSION 1.0                      **
C     **                                               **
C     **  CFD                                          **
C     ***************************************************
C
C     ******************
C     ** MAIN ROUTINE **
C     ******************
C
C
         CALL CONST
C
         CALL INIT
C
         CALL SOLVE
C
C
         STOP
         END
C
C
C     *************************************************
C     ** SOLVE DIFFUSION EQ. BY TIME MARCHING METHOD **
C     *************************************************
C
      SUBROUTINE SOLVE
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER (MX=256)
      COMMON /DIVIDE/ NX,NXM1
      COMMON /LENGTH/ DX,DXI
      COMMON /TMCND/ DT,NSTEP
      COMMON /RESULT/ X(MX),U(MX)
      DIMENSION A(MX),B(MX),C(MX),D(MX)
C
C     ***** IMPLICIT TIME MARCHING LOOP *****
C
      DO 500 ISTEP=1,NSTEP
C
      DO 100 I=2,NXM1
C
        A(I)=-DXI**2
        B(I)=2.D0*DXI**2+1.D0/DT
        C(I)=-DXI**2
  100   D(I)=U(I)/DT
C
C     *** BOUNDARY CONDITION ***
C
        A(1)=0.D0
        B(1)=1.D0
        C(1)=0.D0
        D(1)=0.D0
        A(NX)=0.D0
        B(NX)=1.D0
        C(NX)=0.D0
        D(NX)=0.D0
C
      WRITE(*,*)'           '
      WRITE(*,600) NSTEP,ISTEP
  600 FORMAT('  STEP=',2I5)
C
C
         CALL TDMA(NX,A,B,C,D)
C
      DO 110 I=1,NX
C
  110    U(I)=D(I)
C
  500 CONTINUE
C
          CALL OUTPUT
C
C
      RETURN
      END
C
C     *****************************
C     ** SUBROUTINE CONSTANT SET **
C     *****************************
C
      SUBROUTINE CONST
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER (MX=256)
      COMMON /DIVIDE/ NX,NXM1
      COMMON /TMCND/ DT,NSTEP
      COMMON /LENGTH/ DX,DXI
C
C     ***** NUMBER OF DIVIDE *****
C
         NX=11
         DX=1.D0/DBLE(NX-1)
         DXI=1.D0/DX
C
         NXM1=NX-1
C
         DT=0.96D0*DX**2
C
C     ***** TIME ITERATION *****
C
         NSTEP=20
C
C
      RETURN
      END
C
C     *****************************
C     ** SUBROUTINE INITIAL DATA **
C     *****************************
C
      SUBROUTINE INIT
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER (MX=256)
      COMMON /DIVIDE/ NX,NXM1
      COMMON /TMCND/ DT,NSTEP
      COMMON /LENGTH/ DX,DXI
      COMMON /RESULT/ X(MX),U(MX)
C
C     *****  INITIAL DISTRIBUTION  *****
C
         X(1)=0.D0
C
      DO 100 I=2,NX
C
        X(I)=X(I-1)+DX
C
      IF(I.LE.5) THEN
        U(I)=2.D0*X(I)
      ELSE
        U(I)=2.D0*(1.D0-X(I))
      END IF
C
  100 CONTINUE
C
C
      RETURN
      END
C
C     ***********************
C     ** SUBROUTINE OUTPUT **
C     ***********************
C
      SUBROUTINE OUTPUT
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      PARAMETER (MX=256)
      COMMON /DIVIDE/ NX,NXM1
      COMMON /LENGTH/ DX,DXI
      COMMON /RESULT/ X(MX),U(MX)
C
C
      WRITE(2,600)(X(I),U(I),I=1,NX)
  600 FORMAT(2F12.7)
C
C
      RETURN
      END
C
C     ***********************************
C     ** TRI-DIAGONAL MATRIX INVERSION **
C     ***********************************
C
      SUBROUTINE TDMA(N,A,B,C,D)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION A(N),B(N),C(N),D(N)
C
C
      DO 10 I=1,N-1
C
        C(I)=C(I)/B(I)
        D(I)=D(I)/B(I)
C
        B(I+1)=B(I+1)-C(I)*A(I+1)
        D(I+1)=D(I+1)-D(I)*A(I+1)
C
   10 CONTINUE
C
        C(N)=C(N)/B(N)
        D(N)=D(N)/B(N)
C
      DO 20 I=N-1,1,-1
C
   20    D(I)=D(I)-C(I)*D(I+1)
C
C
      RETURN
      END