      SUBROUTINE INITAL
C...
C...  SIMULATION OF A LINEAR SYSTEM WITH A TIME DELAY
C...
C...  THE DIFFERENTIAL EQUATIONS FOR TWO CSTR*S WITH A LONG CONNECTING
C...  LINE ARE
C...
C...  DCA1/DT = (F0/V1)*CA0 - (F0/V1+K)*CA1, CA1(0) = 0        (1)(2)
C...
C...  DCA2/DT = (F0/V2)*CA1(T-D) - (F0/V2+K)*CA2, CA2(0) = 0   (3)(4)
C...
C...  WHERE
C...
C...     CA1     CONCENTRATION IN THE FIRST AND SECOND CSTR*S
C...     CA2     RESPECTIVELY
C...
C...     CA0     CONCENTRATION ENTERING THE FIRST CSTR
C...
C...     F0      FLOW RATE THROUGH THE SYSTEM
C...
C...     V1      VOLUMES OF THE FIRST AND SECOND CSTR*S, RESPECTIVELY
C...     V2      (V1 = V2 = V)
C...
C...     K       RATE CONSTANT FOR A FIRST-ORDER REACTION IN THE TWO
C...             CSTR*S
C...
C...     D       PURE TIME DELAY PRODUCED BY THE FLOW THROUGH THE LONG
C...             CONNECTING PIPE
C...
C...     T       TIME
C...
C...  THE NUMERICAL SOLUTION COMPUTED BY THE FOLLOWING CODE IS TO BE
C...  COMPARED WITH THE FOLLOWING ANALYTICAL SOLUTION TO EQUATIONS (1)
C...  TO (4) (WHICH CAN BE DERIVED BY THE LAPLACE TRANSFORM)
C...
C...  CA2(T) = C1*U(T-D)*(1 - EXP(-C2*(T-D)) - C3*EXP(-C2*(T-D)))
C...
C...     C1 = (F0/(F0+V*K))**2, C2 = (F0/V+K), C3 = (F0+V*K)
C...
C...            0, T LE 0
C...     U(T) =            (UNIT STEP FUNCTION)
C...            1, T GT 0
C...
      COMMON/T/     T, NSTOP, NORUN
     1      /Y/   CA1,   CA2
     2      /F/DCA1DT,DCA2DT
     3      /C/    C1,    C2,    C3,   F0V,   CA0,     D,     K,    IS
      REAL K
C...
C...  SET THE MODEL CONSTANTS
      F0=0.1
      V=1.
      D=5.
      K=0.1
      CA0=1.
      C1=(F0/(F0+V*K))**2
      C2= F0/V+K
      C3= F0+V*K
      F0V=F0/V
      IS=0
C...
C...  SET THE INITIAL CONDITIONS
      CA1=0.
      CA2=0.
C...
C...  INITIALIZE THE DERIVATIVES
      CALL DERV
      RETURN
      END
      SUBROUTINE DERV
      COMMON/T/     T, NSTOP, NORUN
     1      /Y/   CA1,   CA2
     2      /F/DCA1DT,DCA2DT
     3      /C/    C1,    C2,    C3,   F0V,   CA0,     D,     K,    IS
      REAL K
C...
C...  DIFFERENTIAL EQUATIONS
      DCA1DT=F0V*CA0-(F0V+K)*CA1
      DCA2DT=F0V*CA1D(T-D)-(F0V+K)*CA2
      RETURN
      END
      FUNCTION CA1D(DT)
C...
C...  FUNCTION CA1D SIMULATES THE PURE TIME DELAY
C...
      COMMON/DELAY/CA1S(2,161)
      COMMON/C/    C1,    C2,    C3,   F0V,   CA0,     D,     K,    IS
C...
C...  IF THE DELAYED TIME IS NEGATIVE, THE OUTPUT OF THE TIME DELAY
C...  IS THE INITIAL CONDITION (ZERO IN THE CASE OF EQUATIONS (2) AND
C...  (4))
      IF(DT.GT.0.)GO TO 1
      CA1D=0.
      RETURN
C...
C...  IF THE DELAYED TIME IS NONNEGATIVE, THE DELAYED VARIABLE IS
C...  IN ARRAY CA1S WHICH TO THIS POINT IN TIME HAS BEEN FILLED
C...  IN SUBROUTINE PRINT WITH THE SECOND SUBSCRIPT EQUAL TO IS, I.E.,
C...  THE CURRENT VALUE OF T CORRESPONDS TO IS FOR THE SECOND SUBSCRIPT
1     DO 2 I=1,IS
C...
C...  IF THE DELAYED TIME EXCEEDS THE STORED TIME CA1S(1,I), THE SEARCH
C...  IN THE ARRAY CA1S CONTINUES UNTIL THE DELAYED TIME IS LESS THAN
C...  THE STORED TIME.  FOR THIS LATTER CONDITION, THE DELAYED VARIABLE
C...  (CA1D) EQUALS THE STORED VARIABLE (CA1S(2,I))
      IF(DT.GE.CA1S(1,I))CA1D=CA1S(2,I)
C...
C...  SINCE THE DELAYED VARIABLE HAS BEEN DEFINED, THE SEARCH NEED NOT
C...  CONTINUE
      IF(DT.LT.CA1S(1,I))RETURN
2     CONTINUE
      RETURN
      END
      SUBROUTINE PRINT(NI,NO)
      COMMON/T/     T, NSTOP, NORUN
     1      /Y/   CA1,   CA2
     2      /F/DCA1DT,DCA2DT
     3      /C/    C1,    C2,    C3,   F0V,   CA0,     D,     K,    IS
      REAL K
      COMMON/DELAY/CA1S(2,161)
C...
C...  DIMENSION THE ARRAYS FOR PLOTTING
      DIMENSION TP(41), CP(2,41)
      DATA IP/0/
C...
C...  STORE THE CURRENT VALUES OF T AND CA1 IN ARRAY CA1S FOR USE IN
C...  FUNCTION CA1D
      IS=IS+1
      CA1S(1,IS)=T
      CA1S(2,IS)=CA1
C...
C...  PRINT THE NUMERICAL AND ANALYTICAL SOLUTIONS EVERY FOURTH CALL
C...  TO SUBROUTINE PRINT
      IF(((IS-1)/4*4).NE.(IS-1))RETURN
C...
C...  EVALUATE CA1(T-D) FOR PRINTING
      DELAY=CA1D(T-D)
C...
C...  EVALUATE THE ANALYTICAL SOLUTION TO EQUATIONS (1) TO (4) FOR
C...  PRINTING AND PLOTTING
      IF((T-D).LE.0.)CA2A=0.
      IF((T-D).GT.0.)CA2A=C1*(1.-EXP(-C2*(T-D))-C3*(T-D)*EXP(-C2*(T-D)))
C...
C...  PRINT THE NUMERICAL AND ANALYTICAL SOLUTIONS
      IF(IP.EQ.0)WRITE(NO,2)
2     FORMAT(9X,1HT,9X,3HCA1,5X,8HCA1(T-D),8X,3HCA2,4X,10HCA2 (ANAL))
      WRITE(NO,1)T,CA1,DELAY,CA2,CA2A
1     FORMAT(F10.2,4F12.4)
C...
C...  STORE THE NUMERICAL AND ANALYTICAL SOLUTIONS FOR PLOTTING
      IP=IP+1
      TP(IP)=T
      CP(1,IP)=CA2
      CP(2,IP)=CA2A
C...
C...  PLOT THE NUMERICAL AND ANALYTICAL SOLUTIONS FOR COMPARISON
      IF(IP.LT.41)RETURN
C...
C...  ******************************************************************
C...
C...  THE FOLLOWING CALLS TO A CONTINUOUS PLOTTER ARE MACHINE DEPENDENT
C...  AND WILL HAVE TO BE CHANGED FOR OTHER PLOTTERS OR CONVERTED TO
C...  COMMENTS
C...
C...  CALLS TO LEHIGH CALCOMP PLOTTER
C...  CALL QIKSAX(3,3)
C...  CALL QIKPLT(TP,CP,IP,3H*T*,16H*CA2, CA2(ANAL)*,2H**,2)
C...
C...  CALLS TO SSC LABORATORY PLOTTER
C...
C...  OPEN FILE FOR TOP DRAWER PLOTTING      
      OPEN(4,FILE='T.TOP',STATUS='NEW')
C...
C...  WRITE TOP DRAW FILE FOR PLOTTING ERRORS
      WRITE(4,3)
3     FORMAT(' SET LIMITS X FROM 0 TO 40 Y FROM 0 TO 0.3',/,
     1       ' SET FONT DUPLEX',/,
     2       ' TITLE LEFT "CA2(t), CA2(t) (anal)"',/,
     3       ' TITLE BOTTOM "   t"')
      WRITE(4,4)(TP(I),  CP(1,I),I=1,IP)
4     FORMAT(2F10.3)
      WRITE(4,5)
5     FORMAT(' JOIN ')
      WRITE(4,4)(TP(I),  CP(2,I),I=1,IP)
      WRITE(4,5)
C...
C...  ******************************************************************
C...
      IP=0
      RETURN
      END
