*DECK FLOW01
      SUBROUTINE INITAL
C...
C...  DYNAMICS OF LIQUID TRANSFER FROM A TANK CAR
C...
C...  THE LIQUID IN A TANK CAR, OF LENGTH L AND DIAMETER D, IS TO BE
C...  TRANSFERRED TO TWO VERTICAL STORAGE TANKS.  THE DELIVERY LINES
C...  DISCHARGE AT LEVELS Z1 AND Z2 ABOVE THE PUMP CENTER LINE.  THE
C...  PIPING SYSTEM IS SHOWN BELOW
C...
C...                                         .  .
C...                                         .  .
C...                                     .....  ..... STORAGE
C...            . . .                    .          .  TANK 1
C...         .         .                 .          .
C...       .             .               .          .          .   .
C...      .  CYLINDRICAL  .    ..........+          .          .   .
C...      .   TANK CAR    .    .         ............      .....   .....
C...       .             .     . L1      .          .      .           .
C...         .         .       .         ............      .           .
C...            . . .          ............................+
C...              .            .             L2            .           .
C...              +            .                           .............
C...            .....          .                   STORAGE .           .
C...       PUMP .   ............                    TANK 2 .............
C...            .....    L0
C...              .
C...             ...
C...
C...  THE MATHEMATICAL MODEL FOR THIS SYSTEM, WHICH CAN BE USED TO
C...  CALCULATE
C...
C...     (1)  THE TIME REQUIRED TO EMPTY THE TANK CAR
C...
C...     (2)  THE AMOUNT OF LIQUID PUMPED TO EACH STORAGE TANK (NOTE -
C...          DISCHARGES INTO THE STORAGE TANKS ARE NEVER SUBMERGED)
C...
C...  IS DETAILED BELOW.
C...
C...  THE CONTINUITY BALANCE ON THE TANK CAR IS
C...
C...     DV/DT = D(L*A)/DT = L*(DA/DZ)*(DZ/DT) = -SP*V0
C...
C...  WHERE
C...
C...     V     VOLUME OF LIQUID IN THE TANK CAR
C...
C...     T     TIME
C...
C...     L     LENGTH OF THE TANK CAR
C...
C...     A     CROSS SECTIONAL AREA OF THE LIQUID IN THE TANK CAR
C...
C...     Z     HEIGHT OF LIQUID IN THE TANK CAR
C...
C...     SP    CROSS SECTIONAL AREA OF THE EXIT PIPE FROM THE TANK CAR
C...
C...     V0    VELOCITY IN THE EXIT LINE
C...
C...  V, A AND Z ARE ALL FUNCTIONS OF TIME WHICH IS THE REASON FOR THE
C...  VARIOUS TIME DERIVATIVES IN THE PRECEDING DIFFERENTIAL EQUATION.
C...  A IS RELATED TO Z BY
C...
C...     A = PI*(D/2)**2 -(D**2/8)*(2*ARCCOS((Z - D/2)/(D/2))
C...
C...       - (8/D**2)*(Z - D/2)*SQRT(D**2/4 - (Z - D/2)**2))
C...
C...  FOR WHICH THE DERIVATIVE DA/DZ IS
C...
C...     DA/DZ = (D**2/4)*(2/D)*(1/SQRT(1 - ((Z - D/2)/(D/2))**2))
C...
C...           + SQRT(D**2/4 - (Z - D/2)**2)
C...
C...           - ((Z - D/2)**2)*(1/SQRT(D**2/4 - (Z - D/2)**2))
C...
C...  THE VELOCITIES IN THE LINES CAN BE CALCULATED BY EITHER MOMENTUM
C...  OR ENERGY BALANCES.  IN THIS CASE, WE CHOOSE ENERGY BALANCES.  THE
C...  ENERGY BALANCE FOR THE LINE OF LENGTH L0 IS
C...
C...     D(L0*SP*RHO*(1/(2*GC)*V0**2))/DT
C...
C...        = SP*(PR + (G/GC)*RHO*(Z + ZB))*V0
C...
C...        - SP*P*V0 - F0*V0*ABS(V0)*L0*RHO/(2*DP*GC)*V0*SP
C...
C...        + SP*V0*RHO*(-(G/GC)*ZB)
C...
C...  THE ENERGY BALANCE FOR THE LINE OF LENGTH L1 IS
C...
C...     D(L1*SP*RHO*(1/(2*GC)*V1**2))/DT
C...
C...        = SP*P*V1 - SP*PA*V1
C...
C...        - F1*V1*ABS(V1)*L1*RHO/(2*DP*GC)*V1*SP
C...
C...        + SP*V1*RHO*(G/GC)*(ZB - Z1)
C...
C...  THE ENERGY BALANCE FOR THE LINE OF LENGTH L2 IS
C...
C...     D(L2*SP*RHO*(1/(2*GC)*V2**2))/DT
C...
C...        = SP*P*V2 - SP*PA*V2 -
C...
C...        - F2*V2*ABS(V2)*L2*RHO/(2*DP*GC)*V2*SP
C...
C...        + SP*V2*RHO*(G/GC)*(ZB - Z2)
C...
C...  THE PRESSURE RISE ACROSS THE PUMP, PR, IS RELATED TO THE FLOW
C...  THROUGH THE PUMP, Q, BY THE QUADRATIC EQUATION
C...
C...     PR = A0 + A1*Q + A2*Q**2, Q = V0*SP
C...
C...  FROM CONTINUITY AT THE POINT WHERE THE THREE LINES CONNECT,
C...
C...     V0 = V1 + V2
C...
C...  WHICH CAN BE DIFFERENTIATED WITH RESPECT TO T
C...
C...     DV0/DT = DV1/DT + DV2/DT
C...
C...  EACH OF THESE THREE DERIVATIVES ALSO APPEARS IN ONE OF THE
C...  ENERGY BALANCES, I.E.,
C...
C...     DV0/DT = (GC/(L0*RHO))*(PR + (G/GC)*RHO*(Z + ZB) - P)
C...
C...        - F0*ABS(V0)/(2*DP)*V0 - G*ZB/L0
C...
C...     DV1/DT = (GC/(L1*RHO))*(P - PA)
C...
C...        - F1*ABS(V1)/(2*DP)*V1 - (G/L1)*(Z1 - ZB)
C...
C...     DV2/DT = (GC/(L2*RHO))*(P - PA)
C...
C...        - F2*ABS(V2)/(2*DP)*V2 - (G/L2)*(Z2 - ZB)
C...
C...  IF THE DERIVATIVES DV0/DT, DV1/DT AND DV2/DT ARE ELIMINATED
C...  BETWEEN THESE VARIOUS EQUATIONS
C...
C...     (GC/RHO)*(P/L1 + P/L2 + P/L0)
C...
C...     = (GC/RHO)*(PA/L1 + PA/L2 + PR/L0 + (G/GC)*RHO*(Z + ZB)/L0)
C...
C...     - F0*ABS(V0)/(2*DP)*V0 + F1*ABS(V1)/(2*DP)*V1
C...
C...     + F2*ABS(V2)/(2*DP)*V2 - G*(ZB/L0 - (Z1 - ZB)/L1
C...
C...     - (Z2 - ZB)/L2)
C...
C...  THIS EQUATION CAN THEN BE SOLVED FOR THE PRESSURE P AT THE INTER-
C...  SECTION OF THE THREE LINES.  THE RESULTING EQUATION FOR P CAN THEN
C...  BE USED TO CALCULATE P IN THE RHS OF THE ORDINARY DIFFERENTIAL
C...  EQUATIONS (ODES)
C...
C...     P = ((PR/L0 + PA/L1 + PA/L2 + (G/GC)*RHO*(Z + ZB)/L0)
C...
C...       - F0*ABS(V0)*RHO/(2*DP*GC)*V0
C...
C...       + F1*ABS(V1)*RHO/(2*DP*GC)*V1
C...
C...       + F2*ABS(V2)*RHO/(2*DP*GC)*V2
C...
C...       - (RHO*G/GC)*(ZB/L0 - (Z1 - ZB)/L1 - (Z2 - ZB)/L2))
C...
C...       /(1/L0 + 1/L1 + 1/L2)
C...
C...  THE DEPENDENT VARIABLES OF THE MODEL ARE Z (OR V), V0, V1, AND V2.
C...  ALGEBRAIC VARIABLES INCLUDE V (OR Z), P AND PR.  THUS, FOUR INITIAL
C...  CONDITIONS ARE REQUIRED (TWO MORE DEPENDENT VARIABLES, Q1T AND
C...  Q2T, ARE ADDED BELOW, SO A TOTAL OF SIX DEPENDENT INITIAL CONDI-
C...  TIONS WILL ACTUALLY BE REQUIRED)
C...
C...     Z(0), V0(0), V1(0), V2(0), Q1T(0), Q2T(0)
C...
C...  THE NUMERICAL PARAMETERS FOR THE MODEL ARE
C...
C...     RHO = 1000 KG/M**3
C...
C...     VIS = 0.001 KG/(M-SEC)
C...
C...     D = 2.5 M
C...
C...     L = 10 M
C...
C...     L0 = 30 M
C...
C...     ZB = Z2 = 10 M
C...
C...     Z1 = 20 M
C...
C...     L1 = 40 M
C...
C...     L2 = 30 M
C...
C...     DP = 0.075 M
C...
C...     PA = 1.0E+05 PASCALS (NEWTON/M**2 = N/M**2)
C...
C...     PR = 1.27E+06 + 5.27E+06*Q0 - 2.32E+08*Q0**2
C...
C...  SIMULATE THIS SYSTEM AND PLOT V0(T), V1(T), V2(T), Z(T), Q1T(T)
C...  Q2T(T) VS T.  Q1T(T) AND Q2T(T) ARE CUMULATIVE VOLUMES OF LIQUID
C...  FLOWING INTO THE STORAGE TANKS, I.E.,
C...
C...            T                T
C...     Q1T = INT Q1*DT, Q2T = INT Q2*DT
C...            0                0
C...
C...  THE INITIAL CONDITIONS ARE
C...
C...     Z(0) = 2.0 M, V0(0) = V1(0) = V2(0) = Q1T(0) = Q2T(0) = 0
C...
C...  SET THE MODEL NUMERICAL PARAMETERS
C...
C...  COMMON AREA
      COMMON/T/     T, NSTOP, NORUN
     1      /Y/     Z,    V0,    V1,    V2,   Q1T,   Q2T
     2      /F/  DZDT, DV0DT, DV1DT, DV2DT,    Q1,    Q2
     3      /C/   RHO,   VIS,     D,    DP,     L,    L0,    L1,    L2,
     4             ZB,    Z1,    Z2,    PA,     G,    GC,    SP,   C01,
     5            C11,   C21,   C02,   C12,   C22,    C1,    C2,   CP1,
     6            CP2,   CP3,   CP4,    IP
     7      /V/    F0,    F1,    F2,    PR,     P
      REAL          L,    L0,    L1,    L2
C...
C...  SET THE MODEL PARAMETERS
      RHO=1000.
      VIS=0.001
      D=2.5
      DP=0.075
      L=10.
      L0=30.
      L1=40.
      L2=30.
      ZB=10.
      Z1=20.
      Z2=10.
      PA=1.0E+05
      G=9.8
      GC=1.
C...
C...  COMPUTE SOME CONSTANTS
      PI=3.1415927
      SP=PI*(DP**2)/4.
      C01=GC/(L0*RHO)
      C11=GC/(L1*RHO)
      C21=GC/(L2*RHO)
      C02=(G/L0)*(0.-ZB)
      C12=(G/L1)*(Z1-ZB)
      C22=(G/L2)*(Z2-ZB)
      C1=(G/GC)*RHO
      C2=1./(2.*DP)
      CP1=(PA/L1+PA/L2)
      CP2=RHO/(2.*DP*GC)
      CP3=(ZB/L0-(Z1-ZB)/L1-(Z2-ZB)/L2)
      CP4=1./(1./L0+1./L1+1./L2)
C...
C...  MODEL INITIAL CONDITIONS
      Z=2.0
      V0=0.
      V1=0.
      V2=0.
      Q1T=0.
      Q2T=0.
C...
C...  COMPUTE THE INITIAL DERIVATIVES
      CALL DERV
      IP=0
      END
      SUBROUTINE DERV
      COMMON/T/     T, NSTOP, NORUN
     1      /Y/     Z,    V0,    V1,    V2,   Q1T,   Q2T
     2      /F/  DZDT, DV0DT, DV1DT, DV2DT,    Q1,    Q2
     3      /C/   RHO,   VIS,     D,    DP,     L,    L0,    L1,    L2,
     4             ZB,    Z1,    Z2,    PA,     G,    GC,    SP,   C01,
     5            C11,   C21,   C02,   C12,   C22,    C1,    C2,   CP1,
     6            CP2,   CP3,   CP4,    IP
     7      /V/    F0,    F1,    F2,    PR,     P
      COMMON/IO/   NI,    NO
      REAL          L,    L0,    L1,    L2
C...
C...  TEST FOR AN EMPTY TANK CAR
      IF(Z.GT.0.)THEN
C...
C...  MODEL ALGEBRA
C...
C...  FRICTION FACTORS FOR LINES OF L0, L1 AND L2
      CALL FF(V0,F0)
      CALL FF(V1,F1)
      CALL FF(V2,F2)
C...
C...  PUMP PRESSURE
      CALL PP
C...
C...  PRESSURE AT THE INTERSECTION OF THE THREE LINES
      CALL PC
C...
C...  MODEL DIFFERENTIAL EQUATIONS
      CALL DERV1
      CALL DERV2
      CALL DERV3
      CALL DERV4
      CALL DERV5
      CALL DERV6
      RETURN
C...
C...  TANK CAR IS EMPTY
      ELSE
      WRITE(NO,1)Z,T
1     FORMAT(1H ,//,5H Z = ,F8.4,8H AT T = ,F6.2)
C...
C...  WHEN THE SECOND ELEMENT IN COMMON/T/ (E.G., NSTOP) IS SET TO A
C...  NONZERO VALUE, THE INTEGRATION OF THE DERIVATIVES IN COMMON/F/
C...  IS STOPPED AND SUBROUTINE PRINT IS CALLED ONE MORE TIME TO PRINT
C...  THE SOLUTION AT THIS POINT.  DSS/2 THEN GOES ON TO THE NEXT RUN
C...  IF ONE IS PROGRAMMED (I.E., IF THREE MORE LINES OF DATA ARE PRO-
C...  VIDED)
      NSTOP=1
      END IF
      RETURN
      END
      SUBROUTINE DERV1
      COMMON/T/     T, NSTOP, NORUN
     1      /Y/     Z,    V0,    V1,    V2,   Q1T,   Q2T
     2      /F/  DZDT, DV0DT, DV1DT, DV2DT,    Q1,    Q2
     3      /C/   RHO,   VIS,     D,    DP,     L,    L0,    L1,    L2,
     4             ZB,    Z1,    Z2,    PA,     G,    GC,    SP,   C01,
     5            C11,   C21,   C02,   C12,   C22,    C1,    C2,   CP1,
     6            CP2,   CP3,   CP4,    IP
     7      /V/    F0,    F1,    F2,    PR,     P
      REAL          L,    L0,    L1,    L2
C...
C...  FUNCTION DERV1 EVALUATES THE DERIVATIVE DZ/DT
C...
      DZDT=-SP*V0/(L*DADZ(D,Z))
      RETURN
      END
      SUBROUTINE DERV2
      COMMON/T/     T, NSTOP, NORUN
     1      /Y/     Z,    V0,    V1,    V2,   Q1T,   Q2T
     2      /F/  DZDT, DV0DT, DV1DT, DV2DT,    Q1,    Q2
     3      /C/   RHO,   VIS,     D,    DP,     L,    L0,    L1,    L2,
     4             ZB,    Z1,    Z2,    PA,     G,    GC,    SP,   C01,
     5            C11,   C21,   C02,   C12,   C22,    C1,    C2,   CP1,
     6            CP2,   CP3,   CP4,    IP
     7      /V/    F0,    F1,    F2,    PR,     P
      REAL          L,    L0,    L1,    L2
C...
C...  SUBROUTINE DERV2 EVALUATES THE DERIVATIVE DV0/DT
C...
      DV0DT=C01*(PR+C1*(Z+ZB)-P)-F0*ABS(V0)*V0*C2+C02
      RETURN
      END
      SUBROUTINE DERV3
      COMMON/T/     T, NSTOP, NORUN
     1      /Y/     Z,    V0,    V1,    V2,   Q1T,   Q2T
     2      /F/  DZDT, DV0DT, DV1DT, DV2DT,    Q1,    Q2
     3      /C/   RHO,   VIS,     D,    DP,     L,    L0,    L1,    L2,
     4             ZB,    Z1,    Z2,    PA,     G,    GC,    SP,   C01,
     5            C11,   C21,   C02,   C12,   C22,    C1,    C2,   CP1,
     6            CP2,   CP3,   CP4,    IP
     7      /V/    F0,    F1,    F2,    PR,     P
      REAL          L,    L0,    L1,    L2
C...
C...  SUBROUTINE DERV3 EVALUATES THE DERIVATIVE DV1/DT
C...
      DV1DT=C11*(P-PA)-F1*ABS(V1)*V1*C2-C12
      RETURN
      END
      SUBROUTINE DERV4
      COMMON/T/     T, NSTOP, NORUN
     1      /Y/     Z,    V0,    V1,    V2,   Q1T,   Q2T
     2      /F/  DZDT, DV0DT, DV1DT, DV2DT,    Q1,    Q2
     3      /C/   RHO,   VIS,     D,    DP,     L,    L0,    L1,    L2,
     4             ZB,    Z1,    Z2,    PA,     G,    GC,    SP,   C01,
     5            C11,   C21,   C02,   C12,   C22,    C1,    C2,   CP1,
     6            CP2,   CP3,   CP4,    IP
     7      /V/    F0,    F1,    F2,    PR,     P
      REAL          L,    L0,    L1,    L2
C...
C...  SUBROUTINE DERV4 EVALUATES THE DERIVATIVE DV2/DT
C...
      DV2DT=C21*(P-PA)-F2*ABS(V2)*V2*C2-C22
      RETURN
      END
      SUBROUTINE DERV5
      COMMON/T/     T, NSTOP, NORUN
     1      /Y/     Z,    V0,    V1,    V2,   Q1T,   Q2T
     2      /F/  DZDT, DV0DT, DV1DT, DV2DT,    Q1,    Q2
     3      /C/   RHO,   VIS,     D,    DP,     L,    L0,    L1,    L2,
     4             ZB,    Z1,    Z2,    PA,     G,    GC,    SP,   C01,
     5            C11,   C21,   C02,   C12,   C22,    C1,    C2,   CP1,
     6            CP2,   CP3,   CP4,    IP
     7      /V/    F0,    F1,    F2,    PR,     P
      REAL          L,    L0,    L1,    L2
C...
C...  SUBROUTINE DERV5 EVALUATES THE DERIVATIVE Q1
C...
      Q1=V1*SP
      RETURN
      END
      SUBROUTINE DERV6
      COMMON/T/     T, NSTOP, NORUN
     1      /Y/     Z,    V0,    V1,    V2,   Q1T,   Q2T
     2      /F/  DZDT, DV0DT, DV1DT, DV2DT,    Q1,    Q2
     3      /C/   RHO,   VIS,     D,    DP,     L,    L0,    L1,    L2,
     4             ZB,    Z1,    Z2,    PA,     G,    GC,    SP,   C01,
     5            C11,   C21,   C02,   C12,   C22,    C1,    C2,   CP1,
     6            CP2,   CP3,   CP4,    IP
     7      /V/    F0,    F1,    F2,    PR,     P
      REAL          L,    L0,    L1,    L2
C...
C...  SUBROUTINE DERV6 EVALUATES THE DERIVATIVE Q2
C...
      Q2=V2*SP
      RETURN
      END
      SUBROUTINE FF(V,F)
      COMMON/T/     T, NSTOP, NORUN
     1      /Y/     Z,    V0,    V1,    V2,   Q1T,   Q2T
     2      /F/  DZDT, DV0DT, DV1DT, DV2DT,    Q1,    Q2
     3      /C/   RHO,   VIS,     D,    DP,     L,    L0,    L1,    L2,
     4             ZB,    Z1,    Z2,    PA,     G,    GC,    SP,   C01,
     5            C11,   C21,   C02,   C12,   C22,    C1,    C2,   CP1,
     6            CP2,   CP3,   CP4,    IP
     7      /V/    F0,    F1,    F2,    PR,     P
      REAL          L,    L0,    L1,    L2
C...
C...  SUBROUTINE FF CALCULATES THE FRICTION FACTOR
C...
C...  CALCULATE THE REYNOLDS NUMBER
      RE=RHO*ABS(V)*DP/VIS
C...
C...  IF THE REYNOLDS NUMBER IS LESS THAN 0.001, THE FRICTION FACTOR
C...  IS ARBITRARILY SET TO ZERO
      IF(RE.LT.0.001)F=0.
C...
C...  THE FRICTION FACTOR FROM POISEUILLE*S LAW (LAMINAR FLOW)
      IF((RE.GE.0.001).AND.(RE.LT.2000))F=1./RE
C...
C...  THE FRICTION FACTOR FOR TURBULENT FLOW
      IF(RE.GE.2000.)F=0.183/RE**0.2
      RETURN
      END
      SUBROUTINE PP
      COMMON/T/     T, NSTOP, NORUN
     1      /Y/     Z,    V0,    V1,    V2,   Q1T,   Q2T
     2      /F/  DZDT, DV0DT, DV1DT, DV2DT,    Q1,    Q2
     3      /C/   RHO,   VIS,     D,    DP,     L,    L0,    L1,    L2,
     4             ZB,    Z1,    Z2,    PA,     G,    GC,    SP,   C01,
     5            C11,   C21,   C02,   C12,   C22,    C1,    C2,   CP1,
     6            CP2,   CP3,   CP4,    IP
     7      /V/    F0,    F1,    F2,    PR,     P
      REAL          L,    L0,    L1,    L2
C...
C...  SUBROUTINE PP CALCULATES THE PUMP PRESSURE
C...
C...  CALCULATE THE VOLUMETRIC FLOW RATE FOR USE IN THE PUMP EQUATION
      Q0=V0*SP
C...
C...  CALCULATE THE PRESSURE FROM THE PUMP CHARACTERISTIC EQUATION
      PR=1.27E+06+5.27E+06*Q0-2.32E+08*Q0**2
      RETURN
      END
      SUBROUTINE PC
      COMMON/T/     T, NSTOP, NORUN
     1      /Y/     Z,    V0,    V1,    V2,   Q1T,   Q2T
     2      /F/  DZDT, DV0DT, DV1DT, DV2DT,    Q1,    Q2
     3      /C/   RHO,   VIS,     D,    DP,     L,    L0,    L1,    L2,
     4             ZB,    Z1,    Z2,    PA,     G,    GC,    SP,   C01,
     5            C11,   C21,   C02,   C12,   C22,    C1,    C2,   CP1,
     6            CP2,   CP3,   CP4,    IP
     7      /V/    F0,    F1,    F2,    PR,     P
      REAL          L,    L0,    L1,    L2
C...
C...  SUBROUTINE PC COMPUTES THE PRESSURE WHERE THE THREE LINES ARE
C...  CONNECTED
C...
C...  THE FOUR CONSTANTS, CP1, CP2, CP3 AND CP4 WERE PRECOMPUTED IN
C...  SUBROUTINE INITAL FOR EFFICIENCY (SINCE INITAL IS CALLED ONLY
C...  ONCE IN EACH RUN)
      P=(PR/L0+CP1+C1*(Z+ZB)/L0-F0*ABS(V0)*V0*CP2
     1                         +F1*ABS(V1)*V1*CP2
     2                         +F2*ABS(V2)*V2*CP2
     3                         -C1*CP3)*CP4
      RETURN
      END
      FUNCTION DADZ(D,Z)
C...
C...  FUNCTION DADZ COMPUTES THE DERIVATIVE OF THE WETTED CROSS
C...  SECTIONAL AREA OF THE TANK, A, WITH RESPECT TO THE LIQUID
C...  HEIGHT, Z
C...
C...  ONE HALF THE TANK DIAMETER IS USED IN SEVERAL PLACES IN THE EQUA-
C...  TION FOR DA/DZ, AND THEREFORE IT IS EVALUATED ONCE FOR REPEATED
C...  USE
      D2=D/2.
C...
C...  THE FOLLOWING THREE TERMS IN THE DERIVATIVE WERE DERIVED BY
C...  ANALYTICAL DIFFERENTIATION OF THE EQUATION FOR A(Z).  THE PROGRAM-
C...  MING CONTAINS NEEDLESSLY REPETITIVE ARITHMETIC, E.G., THE SQUARE
C...  ROOT IS COMPUTED THREE TIMES INSTEAD OF COMPUTING IT ONCE AND
C...  USING IT THREE TIMES.  THIS INEFFICIENT PROGRAMMING WAS DONE TO
C...  MORE CLEARLY DEMONSTRATE THE CORRESPONDENCE BETWEEN THE EQUATION
C...  FOR DA/DZ AND THE PROGRAMMING.  IN FACT, PROFESSOR C. A. SILEBI
C...  POINTED OUT THAT THE DERIVATIVE DA/DZ CAN BE WRITTEN DOWN FROM
C...  GEOMETRY IN ONE LINE AS DA = DZ*2*SQRT(R**2 - (Z-R)**2)) OR WITH
C...  MINOR REARRANGEMENT, DA/DZ = SQRT(D**2 - (2*Z - D)**2) WHICH IS
C...  IDENTICAL TO THE FOLLOWING EQUATIONS AFTER SIMPLIFICATION (R = D/2)
      TERM1=(D2**2)*(1./D2)*(1./SQRT(1.-((Z-D2)/D2)**2))
      TERM2=SQRT(D2**2-(Z-D2)**2)
      TERM3=-((Z-D2)**2)/SQRT(D2**2-(Z-D2)**2)
      DADZ=TERM1+TERM2+TERM3
      RETURN
      END
      SUBROUTINE PRINT(NI,NO)
      COMMON/T/     T, NSTOP, NORUN
     1      /Y/     Z,    V0,    V1,    V2,   Q1T,   Q2T
     2      /F/  DZDT, DV0DT, DV1DT, DV2DT,    Q1,    Q2
     3      /C/   RHO,   VIS,     D,    DP,     L,    L0,    L1,    L2,
     4             ZB,    Z1,    Z2,    PA,     G,    GC,    SP,   C01,
     5            C11,   C21,   C02,   C12,   C22,    C1,    C2,   CP1,
     6            CP2,   CP3,   CP4,    IP
     7      /V/    F0,    F1,    F2,    PR,     P
C...
C...  DIMENSION THE ARRAYS FOR PLOTTING
      DIMENSION ZP(51), VP(3,51), QP(2,51), TP(51)
      REAL          L,    L0,    L1,    L2
C...
C...  PRINT THE SOLUTION EVERY FIFTH CALL TO PRINT
      IP=IP+1
      IF(((IP-1)/5*5).NE.(IP-1))GO TO 3
      WRITE(NO,2)
2     FORMAT(9X,1HT,9X,1HZ,8X,2HV0,8X,2HV1,8X,2HV2,/,
     1      15X,5HDZ/DT,4X,6HDV0/DT,4X,6HDV1/DT,4X,6HDV2/DT,/)
C...
C...  PRINT THE NUMERICAL SOLUTION AND ITS DERIVATIVES
      WRITE(NO,1)T,     Z,    V0,    V1,    V2,
     1               DZDT, DV0DT, DV1DT, DV2DT
1     FORMAT(F10.2,4F10.3,/,10X,4F10.3,/)
C...
C...  STORE THE SOLUTION FOR PLOTTING
3     ZP(IP)=Z
      VP(1,IP)=V0
      VP(2,IP)=V1
      VP(3,IP)=V2
      QP(1,IP)=Q1T
      QP(2,IP)=Q2T
      TP(IP)=T
C...
C...  PLOT THE SOLUTION
      IF((NSTOP.EQ.0).AND.(IP.LT.51))RETURN
      CALL TPLOTS(1,IP,TP,ZP)
      WRITE(NO,10)
10    FORMAT(1H ,//,45H Z VS T                                       )
      CALL TPLOTS(3,IP,TP,VP)
      WRITE(NO,11)
11    FORMAT(1H ,//,45H 1 - V0, 2 - V1, 3 - V2 VS T                  )
      CALL TPLOTS(2,IP,TP,QP)
      WRITE(NO,12)
12    FORMAT(1H ,//,45H 1 - Q1T, 2 - Q2T VS T                        )
      RETURN
      END

