      SUBROUTINE INITAL
C...
C...  DYNAMICS OF TWO TANKS CONNECTED BY A LONG LINE WITH FEED AND
C...  DISCHARGE LINES
C...
C...  THE FOLLOWING SYSTEM OF EQUATIONS MODELS TWO HOLDING TANKS
C...  CONNECTED BY A LONG LINE WITH A FEED LINE TO THE FIRST TANK
C...  AND A DISCHARGE LINE WITH A VALVE FROM THE SECOND TANK
C...
C...  MATERIAL BALANCE FOR TANK 1
C...
C...     DH1/DT = (1/A1)*(F0 - F1), F1 = V1*AP
C...
C...  MATERIAL BALANCE FOR TANK 2
C...
C...     DH2/DT = (1/A2)*(F1 - F2), F2 = CV*X*(H2**0.5)
C...
C...  MOMENTUM BALANCE FOR THE FLUID IN THE CONNECTING LINE
C...
C...     (1/GC)*D(AP*L*RHO*V1)/DT =
C...
C...         RHO*(G/GC)*H1*AP - RHO*(G/GC)*H2*AP - FK
C...
C...     FK = F*RHO*L*AP*V1*ABS(V1)*V1/(2*GC*DP)
C...
C...     F = 64/RE, RE LE 2000, F = 0.18/RE**0.2, RE GT 2000
C...
C...     RE = DP*ABS(V1)*RHO/VIS
C...
C...  THE INITIAL CONDITIONS FOR THIS THIRD-ORDER SYSTEM ARE
C...
C...     H1(0) = 2 M, H2(0) = 5 M, V1 = 0 M/S
C...
C...  THE MODEL VARIABLES AND PARAMETERS ARE
C...
C...     H1      HEIGHT OF LIQUID IN THE FIRST TANK (M)
C...
C...     H2      HEIGHT OF LIQUID IN THE SECOND TANK (M)
C...
C...     V1      AVERAGE FLUID VELOCITY IN THE CONNECTING LINE (M/S)
C...
C...     T       TIME (S)
C...
C...     F0      ENTERING FLOW RATE TO THE FIRST TANK (0.2 M**3/S)
C...
C...     F1      FLOW RATE FROM THE FIRST TANK (M**3/S)
C...
C...     F2      FLOW RATE FROM THE SECOND TANK (M**3/S)
C...
C...     AP      CONNECTING LINE CROSS SSCTIONAL AREA (M**2)
C...
C...     DP      CONNECTING LINE DIAMETER (0.3 M)
C...
C...     D1      DIAMETER OF THE FIRST TANK (2.5 M)
C...
C...     D2      DIAMETER OF THE SECOND TANK (3.0 M)
C...
C...     A1      CROSS SECTIONAL AREA OF THE FIRST TANK (M**2)
C...
C...     A2      CROSS SECTIONAL AREA OF THE SECOND TANK (M**2)
C...
C...     RHO     FLUID DENSITY (1000 KG/M**3)
C...
C...     VIS     FLUID VISOCITY (0.001 KG/M-S)
C...
C...     CV      VALVE CONSTANT (0.1 M**3/S-M**0.5)
C...
C...     X       VALVE STEM POSITION (1 DIMENSIONLESS)
C...
C...     L       CONNECTING LINE LENGTH (100 M)
C...
C...     G       ACCELERATION OF GRAVITY (9.8 M/S**2)
C...
C...     GC      CONVERSION FACTOR (1 KG-M/N-S**2)
C...
C...     F       FRICTION FACTOR (DIMENSIONLESS)
C...
C...  THE NUMERICAL SOLUTION TO THIS SYSTEM OF ALGEBRAIC AND DIFFER-
C...  ENTIAL EQUATIONS IS PRINTED AND PLOTTED IN SUBROUTINE PRINT
C...
      COMMON/T/     T
     1      /Y/    H1,    H2,    V1
     2      /F/ DH1DT, DH2DT, DV1DT
     3      /C/    D1,    D2,    A1,    A2,    F0,    F1,    F2,
     4             CV,     X,    AP,     L,    DP,   RHO,   VIS,
     5              G,    GC,     F,    RE,    C1,    C2,    IP
      REAL L
C...
C...  DEFINE THE MODEL PARAMETERS
      D1=2.5
      D2=3.0
      F0=0.2
      CV=0.1
      X=1.0
      L=100.
      DP=0.30
      RHO=1000.
      VIS=0.001
      G=9.8
      GC=1.0
C...
C...  PRECOMPUTE SOME CONSTANTS
      PI=3.1415927
      A1=PI*(D1**2)/4.
      A2=PI*(D2**2)/4.
      AP=PI*(DP**2)/4.
      C1=G/L
      C2=1./(2.*DP)
C...
C...  SET THE MODEL INITIAL CONDITIONS
      H1=2.
      H2=5.
      V1=0.
C...
C...  INITIALIZE A COUNTER FOR THE PLOTTED SOLUTION
      IP=0
C...
C...  INITIALIZE THE MODEL DERIVATIVES
      CALL DERV
      RETURN
      END
      SUBROUTINE DERV
      COMMON/T/     T
     1      /Y/    H1,    H2,    V1
     2      /F/ DH1DT, DH2DT, DV1DT
     3      /C/    D1,    D2,    A1,    A2,    F0,    F1,    F2,
     4             CV,     X,    AP,     L,    DP,   RHO,   VIS,
     5              G,    GC,     F,    RE,    C1,    C2,    IP
      REAL L
C...
C...  MODEL ALGEBRA
      RE=DP*ABS(V1)*RHO/VIS
      IF( RE.LT.0.01)F=6400.
      IF((RE.GE.0.01).AND.(RE.LE.2000.))F=64./RE
      IF( RE.GT.2000.)F=0.18/RE**0.2
      F1=AP*V1
      F2=CV*X*SQRT(H2)
C...
C...  MODEL DIFFERENTIAL EQUATIONS
      DH1DT=(1./A1)*(F0-F1)
      DH2DT=(1./A2)*(F1-F2)
      DV1DT=C1*(H1-H2)-C2*F*V1*ABS(V1)
      RETURN
      END
      SUBROUTINE PRINT(NI,NO)
      COMMON/T/     T
     1      /Y/    H1,    H2,    V1
     2      /F/ DH1DT, DH2DT, DV1DT
     3      /C/    D1,    D2,    A1,    A2,    F0,    F1,    F2,
     4             CV,     X,    AP,     L,    DP,   RHO,   VIS,
     5              G,    GC,     F,    RE,    C1,    C2,    IP
      REAL L
C...
C...  DIMENSION THE ARRAYS FOR THE PLOTTED SOLUTION
      DIMENSION  TP(51), H1P(51), H2P(51), V1P(51), REP(51)
C...
C...  PRINT A HEADING FOR THE NUMERICAL SOLUTION
      IP=IP+1
      IF(IP.EQ.1)WRITE(NO,1)
1     FORMAT(9X,1HT,8X,2HH1,8X,2HH2,8X,2HV1,8X,2HRE,8X,2HF1,8X,2HF2)
C...
C...  PRINT THE NUMERICAL SOLUTON EVERY FIFTH CALL TO SUBROUTINE PRINT
      IF(((IP-1)/5*5).EQ.(IP-1))THEN
         WRITE(NO,2)T,H1,H2,V1,RE,F1,F2
2        FORMAT(F10.1,3F10.3,E10.3,2F10.4)
      END IF
C...
C...  STORE THE SOLUTION FOR SUBSEQUENT PLOTTING
      TP(IP)=T
      H1P(IP)=H1
      H2P(IP)=H2
      V1P(IP)=V1
      REP(IP)=RE
C...
C...  PLOT THE NUMERICAL SOLUTION
      IF(IP.LT.51)RETURN
      CALL SPLOTS(1,IP,TP,H1P)
      WRITE(NO,11)
11    FORMAT(1H ,//,8H H1 VS T)
      CALL SPLOTS(1,IP,TP,H2P)
      WRITE(NO,12)
12    FORMAT(1H ,//,8H H2 VS T)
      CALL SPLOTS(1,IP,TP,V1P)
      WRITE(NO,13)
13    FORMAT(1H ,//,8H V1 VS T)
      CALL SPLOTS(1,IP,TP,REP)
      WRITE(NO,14)
14    FORMAT(1H ,//,8H RE VS T)
C...
C...  ****************************************************************
C...
C...  THE FOLLOWING CALLS TO A CONTINUOUS PLOTTER ARE MACHINE DEPENDENT,
C...  AND MUST CONVERTED FOR OTHER COMPUTERS OR CHANGED TO COMMENTS
C...
C...  CALL QIKSAX(3,3)
C...  CALL QIKPLT(TP,H1P,IP,3H*T*,7H*H1(T)*,12H*H1(T) VS T*,1)
C...  CALL QIKSAX(3,3)
C...  CALL QIKPLT(TP,H2P,IP,3H*T*,7H*H2(T)*,12H*H2(T) VS T*,1)
C...  CALL QIKSAX(3,3)
C...  CALL QIKPLT(TP,V1P,IP,3H*T*,7H*V1(T)*,12H*V1(T) VS T*,1)
C...  CALL QIKSAX(3,3)
C...  CALL QIKPLT(TP,REP,IP,3H*T*,4H*RE*,9H*RE VS T*,1)
C...
C...  ****************************************************************
C...
      RETURN
      END
