      SUBROUTINE INITAL
C...
C...  DYNAMICS OF A NONISOTHERMAL HOLDING TANK
C...
C...  THE DYNAMIC MODEL FOR A NONISOTHERMAL HOLDING TANK HEATED BY A
C...  STEAM COIL IS
C...
C...    TOTAL MATERIAL BALANCE
C...
C...    DV/DT = F0 - F                                              (1)
C...
C...    ENERGY BALANCE
C...
C...    D(V*T1)/DT = F0*T0 - F*T1 + (U*AT)/(RHO*CP)*(TS - T1)       (2)
C...
C...    VALVE EQUATION
C...
C...    F = CV*X*(H**0.5)                                           (3)
C...
C...  WHERE
C...
C...     V      VOLUME OF LIQUID IN THE TANK
C...
C...     T1     TEMPERATURE OF LIQUID IN THE TANK
C...
C...     T      TIME
C...
C...     F0     ENTERING LIQUID FLOW RATE = 5000
C...
C...     F      EXITING LIQUID FLOW RATE
C...
C...     H      HEIGHT OF LIQUID IN THE TANK = V/A
C...
C...     A      TANK CROSS SECTIONAL AREA = 10000
C...
C...     X      VALVE STEM POSITION = 1
C...
C...     CV     VALVE CONSTANT = 10
C...
C...     U      STEAM COIL HEAT TRANSFER COEFFICIENT = 100
C...
C...     AT     STEAM COIL HEAT TRANSFER AREA = 2000
C...
C...     RHO    LIQUID DENSITY = 1
C...
C...     CP     LIQUID HEAT CAPACITY = 1
C...
C...     TS     STEAM TEMPERATURE = 120
C...
C...  THE INITIAL CONDITIONS FOR EQUATIONS (1) AND (2) ARE
C...
C...     V(0) = H(0)*A = 50*10000 = 500000                         (4)
C...
C...     V(0)*T1(0) = 500000*25 = 12500000                         (5)
C...
C...  A NUMERICAL SOLUTION BY AN EULER, FIXED-STEP INTEGRATION IS ALSO
C...  PROGRAMMED IN SUBROUTINE PRINT FOR COMPARISON WITH THE DSS/2
C...  SOLUTION PRINTED IN SUBROUTINE PRINT.
C...
      COMMON/T/   TIME,  NSTOP,  NORUN
     1      /Y/      V,     VT
     2      /F/   DVDT,  DVTDT
     3      /C/      X,     CV,     F0,     TS,    RHO,     CP,
     4              T0,      A,     AT,      U,     C1,      T,
     5               H,     IP
C...
C...  SET THE PROBLEM PARAMETERS
      X=1.
      CV=10.
      F0=5000.
      TS=120.
      RHO=1.
      CP=1.
      T0=25.
      A=10000.
      AT=2000.
      U=100.
C...
C...  PRECOMPUTE THE CONSTANTS REQUIRED IN THE CALCULATION OF THE
C...  DERIVATIVES IN SUBROUTINE DERV
      C1=U*AT/(RHO*CP)
C...
C...  INITIALIZE THE INDEX FOR THE PLOTTED SOLUTION USED IN SUBROUTINE
C...  PRINT
      IP=0
C...
C...  SET THE MODEL INITIAL CONDITIONS
      V=500000.
      VT=V*25.
C...
C...  INITIALIZE THE MODEL DERIVATIVES
      CALL DERV
      RETURN
      END
      SUBROUTINE DERV
      COMMON/T/   TIME,  NSTOP,  NORUN
     1      /Y/      V,     VT
     2      /F/   DVDT,  DVTDT
     3      /C/      X,     CV,     F0,     TS,    RHO,     CP,
     4              T0,      A,     AT,      U,     C1,      T,
     5               H,     IP
C...
C...  MODEL ALGEBRA
      H=V/A
      F=CV*X*SQRT(H)
      T=VT/V
C...
C...  MODEL DIFFERENTIAL EQUATIONS
      DVDT=F0-F
      DVTDT=F0*T0-F*T+C1*(TS-T)
      RETURN
      END
      SUBROUTINE PRINT(NI,NO)
      COMMON/T/   TIME,  NSTOP,  NORUN
     1      /Y/      V,     VT
     2      /F/   DVDT,  DVTDT
     3      /C/      X,     CV,     F0,     TS,    RHO,     CP,
     4              T0,      A,     AT,      U,     C1,      T,
     5               H,     IP
C...
C...  DIMENSION THE ARRAYS REQUIRED FOR PLOTTING
      DIMENSION  TIMEP(51),    HP(51),    TP(51), DVDTP(51),DVTDTP(51)
C...
C...  PRINT A HEADING FOR THE SOLUTION
      IF(IP.EQ.0)WRITE(NO,1)
1     FORMAT(9X,1HT,9X,1HH,9X,1HT,6X,4HDVDT,5X,5HDVTDT)
      IP=IP+1
C...
C...  PRINT THE NUMERICAL SOLUTION EVERY FIFTH CALL TO SUBROUTINE PRINT
      IF(((IP-1)/5*5).EQ.(IP-1))THEN
         WRITE(NO,2)TIME,H,T,DVDT,DVTDT
2        FORMAT(4F10.2,E10.2)
      END IF
C...
C...  STORE THE SOLUTION FOR SUBSEQUENT PLOTTING
      TIMEP(IP)=TIME
      HP(IP)=H
      TP(IP)=T
      DVDTP(IP)=DVDT
      DVTDTP(IP)=DVTDT
C...
C...  PLOT THE SOLUTION AT THE END OF THE RUN
      IF(IP.LT.51)RETURN
      CALL SPLOTS(1,IP,TIMEP,HP)
      WRITE(NO,10)
10    FORMAT(1H ,//,' H VS TIME')
      CALL SPLOTS(1,IP,TIMEP,TP)
      WRITE(NO,11)
11    FORMAT(1H ,//,' T VS TIME')
      CALL SPLOTS(1,IP,TIMEP,DVDTP)
      WRITE(NO,12)
12    FORMAT(1H ,//,' DVDT VS TIME')
      CALL SPLOTS(1,IP,TIMEP,DVTDTP)
      WRITE(NO,13)
13    FORMAT(1H ,//,' DVTDT VS TIME')
C...
C...  PERFORM AN EULER INTEGRATION AT THE END OF THE RUN
      CALL EULER(NI,NO)
      RETURN
      END
      SUBROUTINE EULER(NI,NO)
C...
C...  SUBROUTINE EULER IS REALLY A *MINI DSS/2* IN THE SENSE THAT
C...  IT PERFORMS THE BASIC OPERATIONS OF THE DSS/2 MAIN PROGRAM,
C...  DSS2S, AND THE DSS/S INTEGRATION ROUTINE FOR ODES, INT1.
C...  OF COURSE, IT DOES NOT HAVE THE SOPHISTICATION OF INT1, E.G.,
C...  AUTOMATIC ERROR MONITORING AND CONTROL, SINCE IT PERFORMS ONLY
C...  A FIXED-STEP EULER INTEGRATION.  THEREFORE, THE NUMERICAL SOLU-
C...  TIONS GENERATED BY SUBROUTINE EULER MAY DIFFER SIGNIFICANTLY
C...  FROM THE SOLUTIONS GENERATED BY DSS/2 VIA SUBROUTINE INT1.
C...  HOWEVER, THE CALLS TO TWO USER-SUPPLIED SUBROUTINES, INITAL
C...  AND DERV, IN EULER PARALLEL THE CALLS TO THESE SUBROUTINES IN
C...  THE DSS/2 MAIN PROGRAM.  UNFORTUNATELY, CALLS TO SUBROUTINE
C...  PRINT FROM EULER ARE NOT POSSIBLE SINCE EULER WAS CALLED BY
C...  PRINT, I.E., RECURSIVE CALLS ARE NOT PERMITTED IN FORTRAN.
C...  THEREFORE, THE NUMERICAL SOLUTION PRODUCED BY EULER IS PRINTED
C...  BY WRITE STATEMENTS IN EULER RATHER THAN CALLS TO PRINT, AND
C...  IN THIS SENSE (DIFFERENCES IN PRINTING THE SOLUTION), THE ANALOGY
C...  BETWEEN DSS/2 AND EULER IS NOT COMPLETE.
C...
C...  ALSO, EULER DOES NOT READ DATA LINES LIKE DSS/2.  RATHER, THE
C...  BASIC PARAMETERS CONTROLLING THE ODE NUMERICAL INTEGRATION ARE
C...  SET EXPLICITLY IN EULER FOR THE PARTICULAR ODE PROBLEM.  HOWEVER,
C...  DATA PARALLELING THOSE OF DSS/2 COULD BE READ FROM EULER (VIA
C...  THE UNIT NUMBER NI).
C...
C...  ADDITIONAL INFORMATION IS AVAILABLE IN:
C...
C...     APPENDIX 1 - LISTING OF INTEGRATOR INT1
C...
C...     APPENDIX 2 - LISTING OF MAIN PROGRAM DSS2S
C...
C...  ARGUMENT LIST
C...
C...     NI     INPUT DEVICE (UNIT) NUMBER (INPUT)
C...
C...     NO    OUTPUT DEVICE (UNIT) NUMBER (INPUT)
C...
      COMMON/T/   TIME,  NSTOP,  NORUN
     1      /Y/      V,     VT
     2      /F/   DVDT,  DVTDT
     3      /C/      X,     CV,     F0,     TS,    RHO,     CP,
     4              T0,      A,     AT,      U,     C1,      T,
     5               H,     IP
C...
C...  PRINT A HEADING FOR THE SOLUTION
      WRITE(NO,1)
1     FORMAT(1H1,18H EULER INTEGRATION,//,
     1        9X,1HT,9X,1HH,9X,1HT,6X,4HDVDT,5X,5HDVTDT)
C...
C...  SET THE INITIAL CONDITIONS FOR THE EULER INTEGRATION
      TIME=0.
      CALL INITAL
C...
C...  PRINT THE INITIAL CONDITIONS
      WRITE(NO,2)TIME,H,T,DVDT,DVTDT
2     FORMAT(4F10.2,E10.2)
C...
C...  PERFORM THE EULER INTEGRATION OVER THE INTERVAL 0 LE TIME LE 50
C...  USING N STEPS
      N=50
      DT=50./FLOAT(N)
      DO 3 I=1,N
C...
C...  THE FOLLOWING PROGRAMMING OF THE EULER INTEGRATION IS SPECIFIC
C...  TO THIS TWO-ODE PROBLEM.  HOWEVER, THIS PROGRAMMING COULD EASILY
C...  BE GENERALIZED BY PUTTING THE EULER STEPPING FORMULA IN A DO
C...  LOOP THAT CYCLES THROUGH THE SET OF ODES.  THIS IS THE BASIC
C...  APPROACH FOLLOWED IN THE DSS/2 INTEGRATION ROUTINE INT1
      V=V+DVDT*DT
      VT=VT+DVTDT*DT
      TIME=TIME+DT
C...
C...  UPDATE THE DERIVATIVES
      CALL DERV
C...
C...  PRINT THE NUMERICAL SOLUTION
      IF((I/5*5).EQ.I)THEN
      WRITE(NO,2)TIME,H,T,DVDT,DVTDT
      END IF
3     CONTINUE
      RETURN
      END

