C****************************************************************
C****                                           JK 19.01.01  ****
C****    SOLUTION TO ONE PHASE,                              ****         
C****    ONE DIMENSIONAL, HORIZONTAL FLOW EQUATION           ****
C****    EXPLICIT, IMPLICIT AND ANALYTICAL SOLUTIONS         ****
C****                                                        ****
C****    MAXIMUM GRID SIZE IS CURRENTLY SET TO 30            ****
C****************************************************************
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,N,IPRINT
C
C-----INITIALISATION OF PARAMETERS
C
      CHARACTER SOLUTION*4,OUTFILE*12,PARAM*6
      CHARACTER*10 PHEAD(28)
      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/
      DATA PHEAD/ "TIME      ","PL        ",
     *"P( 1)     ","P( 2)    ","P( 3)     ",
     *"P( 4)     ","P( 5)     ","P( 6)     ","P( 7)     ","P( 8)     ",
     *"P( 9)     ","P(10)     ","P(11)     ","P(12)     ","P(13)     ",
     *"P(14)     ","P(15)     ","P(16)     ","P(17)     ","P(18)     ",
     *"P(19)     ","P(20)     ","P(21)     ","P(22)     ","P(23)     ",
     *"P(24)     ","P(25)     ","PR        "/
C
C-----CORRECTIONS TO SYSTEM DATA MAY BE INPUTTED
C
      WRITE(*,*)'ONE-DIMENSIONAL, HORIZONTAL, ONE-PHASE FLOW IN'
      WRITE(*,*)'A POROUS ROD AT CONSTANT PRESSURE BOUNDARIES'
   20 WRITE(*,*)
      WRITE(*,*)'CURRENT SYSTEM PARAMETERS ARE:'
      WRITE(*,101)L
  101 FORMAT(' LENGTH, L=',F7.2)
      WRITE(*,102)POR
  102 FORMAT(' POROSITY, POR=',F6.2)
      WRITE(*,103)PERM
  103 FORMAT(' PERMEABILITY, PERM=',F6.2)
      WRITE(*,104)VISC
  104 FORMAT(' VISCOSITY, VISC=',F6.2)
      WRITE(*,105)COMPR
  105 FORMAT(' COMPRESSIBILITY, COMPR=',F9.5)
      WRITE(*,106)PINIT
  106 FORMAT(' INITIAL PRESSURE, PINIT=',F6.2)
      WRITE(*,107)PL
  107 FORMAT(' LEFT SIDE PRESSURE, PL=',F6.2)
      WRITE(*,108)PR
  108 FORMAT(' RIGHT SIDE PRESSURE, PR=',F6.2)
      WRITE(*,109)N
  109 FORMAT(' NUMBER OF GRID BLOCKS, N=',I3)
      WRITE(*,110)TMAX
  110 FORMAT(' MAX TIME OF SIMULATION, TMAX=',F8.4)
      WRITE(*,111)DT
  111 FORMAT(' TIME STEP SIZE, DT=',F9.5)
      WRITE(*,112)IPRINT
  112 FORMAT(' NO. OF TIME STEPS BETWEEN PRINTOUTS, IPRINT=',I3)
      WRITE(*,*)
      WRITE(*,*)'INPUT THE NAME OF A PARAMETER TO BE CHANGED',
     * ' (OR TYPE OK): '
      READ(*,*)PARAM
      IF(PARAM.EQ.'L     ')THEN
      WRITE(*,*)'INPUT NEW VALUE FOR LENGTH: '
      READ(*,*)L
      GO TO 20
      ELSEIF(PARAM.EQ.'POR   ')THEN
      WRITE(*,*)'INPUT NEW VALUE FOR POROSITY: '
      READ(*,*)POR
      GO TO 20
      ELSEIF(PARAM.EQ.'PERM  ')THEN
      WRITE(*,*)'INPUT NEW VALUE FOR PERMEABILITY: '
      READ(*,*)PERM
      GO TO 20
      ELSEIF(PARAM.EQ.'VISC  ')THEN
      WRITE(*,*)'INPUT NEW VALUE FOR VISCOSITY: '
      READ(*,*)VISC
      GO TO 20
      ELSEIF(PARAM.EQ.'COMPR ')THEN
      WRITE(*,*)'INPUT NEW VALUE FOR COMPRESSIBILITY: '
      READ(*,*)COMPR
      GO TO 20
      ELSEIF(PARAM.EQ.'PINIT ')THEN
      WRITE(*,*)'INPUT NEW VALUE FOR INITIAL PRESSURE: '
      READ(*,*)PINIT
      GO TO 20
      ELSEIF(PARAM.EQ.'PL    ')THEN
      WRITE(*,*)'INPUT NEW VALUE FOR LEFT SIDE PRESSURE: '
      READ(*,*)PL
      GO TO 20
      ELSEIF(PARAM.EQ.'PR    ')THEN
      WRITE(*,*)'INPUT NEW VALUE FOR RIGHT SIDE PRESSURE: '
      READ(*,*)PR
      GO TO 20
      ELSEIF(PARAM.EQ.'N    ')THEN
      WRITE(*,*)'INPUT NEW VALUE FOR NUMBER OF GRID BLOCKS: '
      READ(*,*)N
      GO TO 20
      ELSEIF(PARAM.EQ.'TMAX  ')THEN
      WRITE(*,*)'INPUT NEW VALUE FOR MAXIMUM TIME: '
      READ(*,*)TMAX
      GO TO 20
      ELSEIF(PARAM.EQ.'DT    ')THEN
      WRITE(*,*)'INPUT NEW VALUE FOR TIME STEP SIZE: '
      READ(*,*)DT
      GO TO 20
      ELSEIF(PARAM.EQ.'IPRINT ')THEN
      WRITE(*,*)'INPUT NEW VALUE FOR PRINT INTERVAL: '
      READ(*,*)IPRINT
      GO TO 20
      ELSEIF(PARAM.EQ.'OK    ')THEN
      ENDIF
C
   31 DX=L/N
      PHEAD(N+3)="PR        "
      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-----INPUT IF IMPLICIT, EXPLICIT OR ANALYTICAL SOLUTION IS WANTED
C
      WRITE(*,*)' SPECITY EXPL, IMPL OR ANAL SOLUTION: '
      READ(*,*)SOLUTION
C
C-----INPUT NAME OF OUTPUT FILE
C
      WRITE(*,*)' SPECIFY NAME OF OUTPUT FILE '
      READ(*,*)OUTFILE
C
C-----OPEN OUTPUT FILE
C
      OPEN(10,FILE=OUTFILE,STATUS='UNKNOWN')
C
C-----WRITE SYSTEM PARAMETERS, HEADINGS AND INITIAL PRESSUERS
C-----TO THE OUTPUT FILE
C
      WRITE(10,*)'ONE-DIMENSIONAL, HORIZONTAL, ONE-PHASE FLOW IN'
      WRITE(10,*)'A POROUS ROD AT CONSTANT PRESSURE BOUNDARIES'
      WRITE(10,*)
      WRITE(10,*)'CURRENT SYSTEM PARAMETERS ARE:'
      WRITE(10,101)L
      WRITE(10,102)POR
      WRITE(10,103)PERM
      WRITE(10,104)VISC
      WRITE(10,105)COMPR
      WRITE(10,106)PINIT
      WRITE(10,107)PL
      WRITE(10,108)PR
      WRITE(10,109)N
      WRITE(10,110)TMAX
      WRITE(10,111)DT
      WRITE(10,112)IPRINT
      WRITE(10,*)
      WRITE(10,*)' SOLUTION TYPE IS: ',SOLUTION
      WRITE(10,*)
      WRITE(10,*)' COMPUTED PRESSURES ARE: '
      WRITE(10,*)
      WRITE(10,98)(PHEAD(I),I=1,N+3)
   98 FORMAT( 4X,50A10)
      WRITE(10,94)T,PL,(P(I),I=1,N),PR
   94 FORMAT(50F10.4)
C
C-----TIME LOOP
C
      DO 200 J=1,1000
      ISW=ISW+1
      T=T+DT
C
    4 IF(SOLUTION.EQ.'EXPL')THEN
C
C-----EXPLICIT SOLUTION
C
      PNEW(1)=P(1)+CONST*(P(2)-3.0*P(1)+2.0*PL)
      PNEW(N)=P(N)+CONST*(2.0*PR-3.0*P(N)+P(N-1))
      DO 7 I=2,N-1
    7 PNEW(I)=P(I)+CONST*(P(I+1)-2.0*P(I)+P(I-1))
C
C-----IMPLICIT SOLUTION
C
      ELSEIF(SOLUTION.EQ.'IMPL')THEN
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-----ANALYTICAL SOLUTION
C
      ELSEIF(SOLUTION.EQ.'ANAL')THEN
C
    3 DO 8 I=1,N
      PNEW(I)=PL+(PR-PL)*X(I)/L
      DO 11 K=1,1000
      DP=(PR-PL)*2./PI/K*EXP(-PI*PI*K*K/L/L*PERM*T/POR/VISC/COMPR)
     * *SIN(X(I)*K*PI/L)
      PNEW(I)=PNEW(I)+DP
      IF(K.GT.10.AND.ABS(DP).LT.0.0000001)GOTO 8
   11 CONTINUE
    8 CONTINUE      
      ENDIF
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,N,K
      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

