C=========================================================================
C     INTEGRATE 1ST-ORDER LINEAR ORDINARY DIFFERENTIAL EQUATION
C                                            BY IMPLICIT-FTCS METHOD     
C     dc/dt = (ce - c)/T
C          INITIAL-CONDITION  : c(0) = c0 
C=========================================================================
C------------------ MAIN PROGRAM -----------------------------------------
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (MXNT=1000) 
      DOUBLE PRECISION C(0:MXNT)
      CHARACTER CHARTRIAL*2
 
      CALL INPUT(C0,CE,T,TT,DT,THETA,CHARTRIAL)
      CALL SOLVER(MXNT,C0,CE,T,TT,DT,THETA,NT,C)
      CALL OUTPUT(MXNT,C0,CE,T,TT,DT,THETA,NT,C,CHARTRIAL)

      STOP 'NORMAL TERMINATION'
      END
C-------------------------------------------------------------------------
C------------------ SUB ROUTINE ------------------------------------------
      SUBROUTINE INPUT(C0,CE,T,TT,DT,THETA,CHARTRIAL)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      CHARACTER CHARTRIAL*2
      WRITE(*,*) 'Input initial concentration c0'
      READ(*,*) C0
      WRITE(*,*) 'Input equilibrium concentration ce'
      READ(*,*) CE
      WRITE(*,*) 'Input internal time constant T'
      READ(*,*) T
      WRITE(*,*) 'Input total time of calculation TT'
      READ(*,*) TT
      WRITE(*,*) 'Input time step dt'
      READ(*,*) DT
      WRITE(*,*) 'Input weighting parameter theta'
      READ(*,*) THETA
      IF((THETA.LT.0).OR.(THETA.GT.1)) STOP
     &'The range of theta is between 0 and 1'
      WRITE(*,*) 'Input 2 characters to append to filename'
      READ(*,'(A2)') CHARTRIAL
      RETURN
      END



      SUBROUTINE SOLVER(MXNT,C0,CE,T,TT,DT,THETA,NT,C)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION C(0:MXNT)
      NT = TT/DT
      IF(NT.GT.MXNT) STOP 'NT > MXNT'
C---- RESET
      DO N = 0,NT
         C(N) = 0.
      END DO
C---- COMPUTATION OF NUMERICAL SOLUTION BY IMPLICIT METHOD
      C(0) = C0
      DO N = 0,NT-1
         C(N+1) = (1-(1-THETA)*DT/T)/(1+THETA*DT/T)*(C(N)-CE)+CE
      END DO
      RETURN
      END



      SUBROUTINE OUTPUT(MXNT,C0,CE,T,TT,DT,THETA,NT,C,CHARTRIAL)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION C(0:MXNT)
      CHARACTER CHARTRIAL*2,INPDATA*13,OUTDATA*14
      INPDATA = './input'//CHARTRIAL//'.dat'
      OUTDATA = './output'//CHARTRIAL//'.dat'
C---- OUTPUT OF INPUT-PARAMETER
      OPEN(1,FILE=INPDATA,STATUS='REPLACE')
      WRITE(1,10) 'INITIAL CONCENTRATION     =   ', C0
      WRITE(1,10) 'EQUILIBRIUM CONCENTRATION =   ', CE
      WRITE(1,10) 'INTERNAL TIME CONSTANT    =   ', T
      WRITE(1,10) 'TOTAL TIME                =   ', TT
      WRITE(1,10) 'TIME STEP                 =   ', DT
      WRITE(1,10) 'WEIGHTING PARAMETER       =   ', THETA
10    FORMAT(A30, F15.8)
      CLOSE(1)
C---- OUTPUT OF NUMERICAL SOLUTION
      OPEN(2,FILE=OUTDATA,STATUS='REPLACE')
      DO N = 0,NT
         WRITE(2,20) N*DT, C(N)
20       FORMAT(2F20.8)
      END DO
      CLOSE(2)
      RETURN
      END