C-----FINITE DIFFERENCE SOLUTION OF ONE-PHASE DIFFUSIVITY EQUATION
C-----FOR CONSTANT LEFT SIDE AND RIGHT SIDE PRESSURES (PL AND PR)
C
C     P AND PNEW = PRESSURES
C     POR = POROSITY
C     PERM = PERMEABILITY
C     VISC = VISCOSITY
C     COMPR = COMPRESSIBILITY
C     L = LENGTH
C     PINIT= INITIAL PRESSURE
C     PL = LEFT SIDE PRESSURE
C     PR = RIGHT SIDE PRESSURE
C     DT = TIME STEP SIZE
C     TMAX = MAX TIME OF SIMULATION
C     N = NUMBER OF GRID BLOCKS
C     IPRINT = NUMBER OF TIME STEPS BETWEEN PRINTOUTS
C
C-----DECLARATIONS
C      
      REAL P(25),PNEW(25),X(25),A(25),B(25),C(25),D(25)
      REAL POR,PERM,VISC,COMPR,L,PL,PR,PINIT,DT,TMAX,PI,DX
      INTEGER I,J,K,N
C
C-----INITIALISATION OF PARAMETERS
C
      DATA POR/0.2/,PERM/1.0/,VISC/1.0/,COMPR/0.0001/,L/100./,PL/2.0/
     *,PR/1.0/PINIT/1.0/,DT/0.0005/,N/10/,TMAX/0.01/,IPRINT/1/
      DATA PI/3.14159/
C
   31 DX=L/N
      T=0.
      CONST=DT/DX/DX*PERM/POR/VISC/COMPR
      ALPHA=1./CONST
      ISW=0
      DO 5 I=1,N
      P(I)=PINIT
      PNEW(I)=PINIT
    5 X(I)=L*I/N-L/N/2.
C
C-----OPEN OUTPUT FILE
C
      OPEN(10,FILE=OUT,STATUS='UNKNOWN')
C
C-----TIME LOOP
C
      DO 200 J=1,1000
      ISW=ISW+1
      T=T+DT
C
C-----COEFFICIENTS FOR BLOCK 1
C
      A(1)=0.
      B(1)=-3.-ALPHA
      C(1)=1.
      D(1)=-ALPHA*P(1)-2.*PL
C
C-----COEFFICIENTS FOR BLOCKS 2 TO N-1
C
      DO 199 I=2,N-1
      A(I)=1.
      B(I)=-2.-ALPHA
      C(I)=1.
      D(I)=-ALPHA*P(I)
  199 CONTINUE
C
C-----COEFFICIENTS FOR BLOCK N
C
      A(N)=1.
      B(N)=-3.-ALPHA
      C(N)=0.
      D(N)=-ALPHA*P(N)-2.*PR
C
C-----GET NEW PRESSURES BY CALLING GAUSSIAN ELIMINATION ROUTINE
C
      CALL TRIDIA(N,A,B,C,D,PNEW)
C
C-----PRINT (?)
C
      IF(ISW.NE.IPRINT)GO TO 99
      ISW=0
      WRITE(10,100)T,PL,(PNEW(I),I=1,N),PR
  100 FORMAT(50F10.4)
   99 CONTINUE      
C
C-----END (?)
C
      IF(T.GE.TMAX)STOP
C
C-----UPDATING OF PRESSURES
C
   13 DO 14 I=1,N
   14 P(I)=PNEW(I)
  200 CONTINUE
      END
C
C----SUBROUTINE FOR GAUSSIAN ELIMINATION
       SUBROUTINE TRIDIA(N,A,B,C,D,P)
C
C------------------------------------------------------------------
C     THE SUBROUTINE USES GAUSSIAN ELIMINATION FOR SOLUTION OF THE
C     SET OF EQUATIONS
C   
C     A(I)*P(I-1) + B(I)*P(I) + C(I)*P(I+1) = D(I)
C
C     A(I),B(I),C(I),D(I)=MATRIX COEFFICIENTS
C     P=UNKNOWN PRESSURE
C     N=NUMBER OF EQUATIONS
C------------------------------------------------------------------	 
C
      REAL X
      REAL A(25),B(25),C(25),D(25),P(25),BB(25),DD(25)
      INTEGER I,K,N
      BB(1)=B(1)
      DD(1)=D(1)
      DO 60 I=2,N
      X=A(I)/BB(I-1)
      BB(I)=B(I)-X*C(I-1)
   60 DD(I)=D(I)-X*DD(I-1)
      P(N)=DD(N)/BB(N)
      DO 70 K=2,N
      I=N-K+1
   70 P(I)=(DD(I)-C(I)*P(I+1))/BB(I)
      RETURN
      END