***************************************************************

      PROGRAM  PROG4

********************************************************************************

*           PROGRAM FOR PUMP RUNDOWN FROM POWER FAILURE

*           AND SUBSEQUENT RESTART. PROVISION IS MADE FOR

*           COLUMN SEPARATION, AIR-VACUUM VALVES, AND

*           ONE-WAY VACUUM VALVES. A BYPASS LINE IS IN PLACE

*           AROUND THE PUMPS.

*

********************************************************************************

*

      DIMENSION X[ALLOCATABLE](:,:),VL[ALLOCATABLE](:,:),

     $VR[ALLOCATABLE](:,:),H[ALLOCATABLE](:,:),HEAD[ALLOCATABLE](:,:),

     $HHIGH[ALLOCATABLE](:,:),HLOW[ALLOCATABLE](:,:),

     $VNEWL[ALLOCATABLE](:,:),VNEWR[ALLOCATABLE](:,:),

     $HNEW[ALLOCATABLE](:,:),HSTEAD[ALLOCATABLE](:,:),

     $THIGH[ALLOCATABLE](:,:),TLOW[ALLOCATABLE](:,:),

     $CAV[ALLOCATABLE](:,:),HCAV[ALLOCATABLE](:,:),

     $CAVMAX[ALLOCATABLE](:,:),TCAV[ALLOCATABLE](:,:)

      DIMENSION L[ALLOCATABLE](:),D[ALLOCATABLE](:),

     $A[ALLOCATABLE](:),PIPEZ[ALLOCATABLE](:),

     $F[ALLOCATABLE](:),VZERO[ALLOCATABLE](:),

     $NPAR[ALLOCATABLE](:),C[ALLOCATABLE](:),AK[ALLOCATABLE](:),

     $SINE[ALLOCATABLE](:),AREA[ALLOCATABLE](:),

     $DELTT[ALLOCATABLE](:),RATIO[ALLOCATABLE](:),

     $AIRVAC[ALLOCATABLE](:),CTZERO[ALLOCATABLE](:),

     $BUBBLE[ALLOCATABLE](:)

      DIMENSION TSTART(10),TRAMP(10),RPMSTRT(10),RPMEND(10),RPMS(10),

     $QP(10)

      DIMENSION QN(6),HNSQ(6),TNSQ(6),IOU(6),TIOU(6)

     $,QNN(6),HNN(6),C7(6),C8(6),IZ(6)

      LOGICAL FAIL,BYPASS,SAVE,PPLOT,HVPRNT,RERUN,RSTART

      LOGICAL BUBBLE,GRAPH,PFILE,NOGOOD

      INTEGER PIPE,AIRVAC

      REAL LA,L,NEXP,KLB

      CHARACTER TITLE*80,TITLE2*80,VACYES*3,VACNO*3,AV*3

      CHARACTER CH,CH12,CH78,CH79,CHN,FNAME*12

      DATA SAVE,PPLOT,HVPRNT,RERUN,PFILE,GRAPH/6*.FALSE./

*

      NAMELIST /SPECS/ NPIPES,NPARTS,HRES,HSUMP,ZEND,HATM,QTRY,QACC,

     $                 TMAX,DTNEW,HVPRNT,PPLOT,PFILE,GRAPH,RERUN,

     $                 DB,KLB

      NAMELIST /PUMPS/ NPUMPS,NSTAGE,RPM,RPMZ,WRSQ,NSTART,QN,HNSQ,TNSQ

      NAMELIST /RESTART/ TSTART,TRAMP,RPMSTRT,RPMEND

      NAMELIST /OUTCTRL/ IOU,TIOU

*

*------------------------------------------------------------------------------

      WRITE(*,791)

      READ(*,790) FNAME

      OPEN(5,FILE=FNAME)

      READ(5,100) TITLE

      READ(5,100) TITLE2

      READ(5,SPECS)

      NP=NPIPES

      ALLOCATE (L(NP),D(NP),A(NP),PIPEZ(NP),F(NP),VZERO(NP),NPAR(NP),

     $C(NP),AK(NP),SINE(NP),AREA(NP),DELTT(NP),RATIO(NP),AIRVAC(NP),

     $BUBBLE(NP),CTZERO(NP))

      READ(5,*) (PIPE,D(I),L(I),F(I),A(I),PIPEZ(I),AIRVAC(I),I=1,NPIPES)

      READ(5,PUMPS)

      IF(NSTART.GT.0) READ(5,RESTART)

      READ(5,OUTCTRL)

      WRITE(*,792)

      READ(*,790) FNAME

      OPEN(6,FILE=FNAME,STATUS='NEW')

*------------------------------------------------------------------------------

*

      III=2

      ICOUNT=0

      IOUT=IOU(1)

      CH12=CHAR(12)

      CH=CHAR(27)

      CH78=CHAR(78)

      CH79=CHAR(79)

      CHN=CHAR(3)

      DELTN=100000.

      RPMZ=RPM

      VACYES='YES'

      VACNO=' NO'

      RSTART=.FALSE.

      FAIL=.FALSE.

      BYPASS=.FALSE.

      IF(PPLOT.OR.HVPRNT.OR.GRAPH) SAVE=.TRUE.

      TTTT=10000.

      DO 303 I=1,NSTART

      RPMS(I)=0.

      QP(I)=0.

      IF(TSTART(I).LT.TTTT) TTTT=TSTART(I)

303 CONTINUE

      PINERT=WRSQ/32.2

      PI=3.141592

      COEF=RPMZ/RPM

      DO 302 I=1,6

      QNN(I)=QN(I)*COEF

302 HNN(I)=HNSQ(I)*COEF*COEF

      NEXP=1.0

      IF(F(1).GT.10.) NEXP=0.85

*

* COMPUTE STEADY STATE DISCHARGE

      IF(QTRY.LT.0.001) QTRY=QN(4)

      COEF=12.*144.*144.*16.*NPUMPS*NPUMPS/(64.4*449.*449.*PI*PI)

      IF(F(1).GT.10.) COEF=3.03*(12.**4.87)*(NPUMPS/(449.*.7854))**1.85

      EX=2.0

      IF(F(1).GT.10.) EX=1.85

      DO 320 NZ=1,20

      DO 300 I=1,5

      IF(QTRY.GT.QN(I).AND.QTRY.LE.QN(I+1)) GO TO 310

300 CONTINUE

      WRITE(6,250)

      STOP

310 HN=HNSQ(I)+(QTRY-QN(I))*(HNSQ(I+1)-HNSQ(I))/(QN(I+1)-QN(I))

      SUM=0.

      DO 301 J=1,NPIPES

      IF(F(1).GT.10.) SUM=SUM+L(J)*QTRY**1.85/((F(J)**1.85)*

     $(D(J)**4.87))

      IF(F(1).LT.10.) SUM=SUM+F(J)*L(J)*QTRY*QTRY/D(J)**5

301 CONTINUE

      FUNCT=HSUMP-HRES+NSTAGE*HN-COEF*SUM

      FPRIME=NSTAGE*(HNSQ(I+1)-HNSQ(I))/(QN(I+1)-QN(I))-2.*COEF*SUM/QTRY

      QNEXT=QTRY-FUNCT/FPRIME

      IF(ABS(QNEXT-QTRY).LT.QACC) GO TO 350

      QTRY=QNEXT

320 CONTINUE

      WRITE(6,251)

      STOP

350 QLINE=QNEXT*NPUMPS

      Q=QNEXT

      HPUMP=NSTAGE*HN+HSUMP

      COEF=4.*144./(449.*PI)

      DO 2 I=1,NPIPES

      VZERO(I)=QLINE*COEF/D(I)**2

      AREA(I)=0.25*PI*D(I)*D(I)/144.

      C(I)=32.2/A(I)

      DELTT(I)=L(I)/(NPARTS*(2.*VZERO(I)+A(I)))

      IF(I.EQ.NPIPES) SINE(I)=(ZEND-PIPEZ(I))/L(I)

      IF(I.EQ.NPIPES) GO TO 2

      SINE(I)=(PIPEZ(I+1)-PIPEZ(I))/L(I)

   2 CONTINUE

*

* ** COMPUTE PUMP BYPASS PARAMETERS **

*

      AB=.25*PI*DB**2/144.

      CCA=(KLB*AREA(1)**2)/(64.4*AB**2)

      CCB=1.0/C(1)

*

* ** COMPUTE MINIMUM DELTA T **

*

      DELT=DELTT(1)

      KMIN=1

      DO 12 I=2,NPIPES

      IF(DELTT(I).GT.DELT) GO TO 12

      DELT=DELTT(I)

      KMIN=I

 12 CONTINUE

      IF(RERUN) DELT=DTNEW

      DO 13 I=1,NPIPES

      IF(I.EQ.KMIN) GO TO 14

      ANPAR=L(I)/(DELT*(VZERO(I)+A(I)))

      NPAR(I)=ANPAR

      TEST=NPAR(I)

      IF(ABS(TEST+1-ANPAR).LT..0001) NPAR(I)=NPAR(I)+1

      GO TO 166

 14 NPAR(I)=NPARTS

166 IF(NPAR(I).GT.NPARR) NPARR=NPAR(I)

 13 CONTINUE

      IX=NPARR+1

      IZZ=IX+1

      ALLOCATE(X(NP,IX),H(NP,IX),HLOW(NP,IX),HHIGH(NP,IX),HEAD(NP,IZZ),

     $HNEW(NP,IX),HSTEAD(NP,IX),VL(NP,IX),VR(NP,IX),VNEWL(NP,IX),

     $VNEWR(NP,IX),CAV(NP,IX),HCAV(NP,IX),CAVMAX(NP,IX),THIGH(NP,IX),

     $TLOW(NP,IX),TCAV(NP,IX))

      DO 177 I=1,NPIPES

      K=NPAR(I)+1

177 HEAD(I,K+1)=0.

*

      WRITE(6,1011) CH,CH78,CHN

      WRITE(6,201) TITLE

      WRITE(6,2011) TITLE2

      WRITE(6,202)

      WRITE(6,101) NPIPES,NPARTS,HRES,HSUMP,HATM,ZEND,TMAX,QACC,DB,KLB

      WRITE(6,1010) (IOU(I),I=1,6),(TIOU(I),I=1,6)

      WRITE(6,104)

      WRITE(6,105) NPUMPS,NSTAGE,RPMZ,WRSQ,RPM

      WRITE(6,106)

      WRITE(6,107) (QN(I),HNSQ(I),TNSQ(I),I=1,6)

      WRITE(6,108) (I,TSTART(I),RPMSTRT(I),RPMEND(I),TRAMP(I),

     $I=1,NSTART)

      HPUMPP=HPUMP-HSUMP

      WRITE(6,1088) QLINE,HPUMPP

      WRITE(6,1011) CH12

      IF(F(1).LT.10.) WRITE(6,102)

      IF(F(1).GT.10.) WRITE(6,1022)

      DO 3 I=1,NPIPES

      IF(F(1).LT.10.) WRITE(6,103) I,D(I),L(I),A(I),PIPEZ(I),F(I),

     $VZERO(I)

      IF(F(1).GT.10.) WRITE(6,1033) I,D(I),L(I),A(I),PIPEZ(I),F(I),

     $VZERO(I)

  3 CONTINUE

      WRITE(6,1020)

      DO 4 I=1,NPIPES

      AV=VACNO

      IF(AIRVAC(I).GE.1) AV=VACYES

      LA=L(I)/A(I)

      AINT=1.0-DELT*NPAR(I)/LA

  4 WRITE(6,1030) I,DELTT(I),NPAR(I),SINE(I),LA,AINT,AV

      WRITE(6,1040) DELT

*

*  CONVERT TO Q/N, H/NSQ, T/NSQ

      COEF=60.*550./(2.*PI)

      DO 1 I=1,6

      QN(I)=QN(I)/(449.*RPM)

      TNSQ(I)=NSTAGE*COEF*TNSQ(I)/RPM**3

   1 HNSQ(I)=HNSQ(I)/RPM**2

*

*

* ** SET UP CONDITIONS FOR STEADY STATE (T=0) **

*

      H(1,1)=HPUMP

      VL(1,1)=100000.

      VNEWL(1,1)=VL(1,1)

      VR(1,1)=VZERO(1)

      DO 20 I=1,NPIPES

      DELL=L(I)/NPAR(I)

      IF(F(I).LT.10.) AK(I)=12.*F(I)*DELT/(2.*D(I))

      IF(F(I).GT.10.) AK(I)=12.*DELT*195./(2.*D(I)*(F(I)**1.85)*

     $(D(I)/12.)**.17)

      IF(F(I).GT.10.) F(I)=195./((F(I)**1.85)*(VZERO(I)**.15)*

     $((D(I)/12.)**.17))

      DELHF=12.*F(I)*DELL*VZERO(I)**2/(64.4*D(I))

      K=NPAR(I)+1

      X(I,1)=0.

      HEAD(I,1)=H(I,1)-PIPEZ(I)

      HCAV(I,1)=PIPEZ(I)-HATM

      RATIO(I)=DELT/DELL

      CAV(I,1)=0.

      CAVMAX(I,1)=0.

      TLOW(I,1)=0.

      THIGH(I,1)=0.

      TCAV(I,1)=0.

      BUBBLE(I)=.FALSE.

      DO 21 J=2,K

      X(I,J)=(J-1)*DELL/L(I)

      H(I,J)=H(I,1)-DELHF*(J-1)

      HEAD(I,J)=H(I,J)-(PIPEZ(I)+X(I,J)*L(I)*SINE(I))

      HCAV(I,J)=PIPEZ(I)+X(I,J)*L(I)*SINE(I)-HATM

      CAV(I,J)=0.

      CAVMAX(I,J)=0.

      TLOW(I,J)=0.

      THIGH(I,J)=0.

      TCAV(I,J)=0.

      VL(I,J)=VZERO(I)

      VR(I,J)=VZERO(I)

      IF(J.EQ.K) VR(I,J)=100000.

      IF(J.EQ.K) VNEWR(I,J)=VR(I,J)

 21 CONTINUE

      IF(I.EQ.NPIPES) GO TO 20

      H(I+1,1)=H(I,K)

      VR(I+1,1)=VZERO(I+1)

      VL(I+1,1)=100000.

      VNEWL(I+1,1)=VL(I+1,1)

 20 CONTINUE

      DO 23 I=1,NPIPES

      IF(AIRVAC(I).GE.1) HCAV(I,1)=PIPEZ(I)

      IF(I.GT.1) KK=NPAR(I-1)+1

      IF(I.GT.1.AND.AIRVAC(I).GE.1) HCAV(I-1,KK)=PIPEZ(I)

      K=NPAR(I)+1

      DO 23 J=1,K

      HSTEAD(I,J)=H(I,J)

      HLOW(I,J)=H(I,J)

 23 HHIGH(I,J)=H(I,J)

      PHMAX=-100.

      PHMIN=100000.

      DO 11 I=1,NPIPES

      K=NPAR(I)+1

      DO 11 J=1,K

      IF(HEAD(I,J).LT.PHMAX) GO TO 16

      IPMAX=I

      XMAX=X(I,J)

      PHMAX=HEAD(I,J)

      TTMAX=T

      GO TO 17

 16 IF(HEAD(I,J).GT.PHMIN) GO TO 17

      IPMIN=I

      XMIN=X(I,J)

      PHMIN=HEAD(I,J)

      TTMIN=T

 17 CONTINUE

 11 CONTINUE

*

* ** WRITE OUT STEADY STATE CONDITIONS **

*

      T=0.0

      IENTRY=0

      INDEX=TMAX/DELT+1

      WRITE(6,1011) CH12

      WRITE(6,204)

      WRITE(6,205) T

      DO 22 I=1,NPIPES

      K=NPAR(I)+1

      WRITE(6,206) I,(X(I,J),HEAD(I,J),H(I,J),VL(I,J),

     $VR(I,J),CAV(I,J),J=1,K)

 22 CONTINUE

      II=1

      IF(SAVE) CALL GRAFCAV(II,NP,IX,INDEX,DELT,VL,VR,HEAD,CAV,

     $NPAR,IENTRY)

      IENTRY=1

      HPUMP=HPUMP-HSUMP

      WRITE(6,212) RPM,Q,HPUMP

      Q=Q/449.

*

*-------------------------------------------------------------------------------

*                  BEGIN TRANSIENT ANALYSIS

*-------------------------------------------------------------------------------

*

      DO 99 II=1,INDEX

      T=T+DELT

*

* ** COMPUTE H AND V AT INTERIOR NODES **

*

      DO 30 I=1,NPIPES

      K=NPAR(I)

      DO 30 J=2,K

      RA=RATIO(I)*A(I)

      VLEFT=VL(I,J)-RA*(VL(I,J)-VR(I,J-1))

      VRITE=VR(I,J)-RA*(VR(I,J)-VL(I,J+1))

      HLEFT=H(I,J)-RA*(H(I,J)-H(I,J-1))

      HRITE=H(I,J)-RA*(H(I,J)-H(I,J+1))

      IF (CAV(I,J).GT.0.) GO TO 31

      VNEWL(I,J)=0.5*(VLEFT+VRITE+C(I)*(HLEFT-HRITE)+C(I)*DELT*SINE(I)

     $*(VLEFT-VRITE)-AK(I)*(VLEFT*ABS(VLEFT)**NEXP+VRITE*

     $ABS(VRITE)**NEXP))

      HNEW(I,J)=0.5*(HLEFT+HRITE+(VLEFT-VRITE)/C(I)+DELT*SINE(I)*

     $(VLEFT+VRITE)-(AK(I)/C(I))*(VLEFT*ABS(VLEFT)**NEXP-VRITE*

     $ABS(VRITE)**NEXP))

      VNEWR(I,J)=VNEWL(I,J)

      IF (HNEW(I,J).LT.HCAV(I,J)) GO TO 31

      GO TO 30

 31 C1=VLEFT+C(I)*HLEFT-AK(I)*VLEFT*ABS(VLEFT)**NEXP+C(I)*DELT*VLEFT*

     $SINE(I)

      C3=VRITE-C(I)*HRITE-AK(I)*VRITE*ABS(VRITE)**NEXP-C(I)*DELT*VRITE*

     $SINE(I)

      VNEWL(I,J)=C1-C(I)*HCAV(I,J)

      VNEWR(I,J)=C3+C(I)*HCAV(I,J)

      HNEW(I,J)=HCAV(I,J)

      CAV(I,J)=(VNEWR(I,J)-VNEWL(I,J))*DELT*AREA(I)+CAV(I,J)

      IF (CAV(I,J).GT. 0.) GO TO 30

      DELTAH=ABS(0.5*A(I)*(VNEWR(I,J)-VNEWL(I,J))/32.2)

      HNEW(I,J)=HCAV(I,J)+DELTAH

      VNEWL(I,J)=0.5*(VNEWL(I,J)+VNEWR(I,J))

      VNEWR(I,J)=VNEWL(I,J)

      CAV(I,J)=0.

 30 CONTINUE

*

* ** COMPUTE H AND V AT INTERIOR JUNCTIONS **

*

      KK=NPIPES-1

      DO 40 I=1,KK

      K=NPAR(I)+1

      VLEFT=VL(I,K)-RATIO(I)*A(I)*(VL(I,K)-VR(I,K-1))

      HLEFT=H(I,K)-RATIO(I)*A(I)*(H(I,K)-H(I,K-1))

      CCC=VLEFT+C(I)*HLEFT+C(I)*VLEFT*DELT*SINE(I)-AK(I)*VLEFT*

     $ABS(VLEFT)**NEXP

      VRITE=VR(I+1,1)-RATIO(I+1)*A(I+1)*(VR(I+1,1)-VL(I+1,2))

      HRITE=H(I+1,1)-RATIO(I+1)*A(I+1)*(H(I+1,1)-H(I+1,2))

      CC=VRITE-C(I+1)*HRITE-C(I+1)*DELT*VRITE*SINE(I+1)-AK(I+1)*VRITE*

     $ABS(VRITE)**NEXP

      IF (CAV(I,K).GT.0.) GO TO 41

      HNEW(I,K)=(-CC*AREA(I+1)+CCC*AREA(I))/(C(I+1)*AREA(I+1)+C(I)*

     $AREA(I))

      IF(HNEW(I,K).LT.HCAV(I,K)) GO TO 41

      HNEW(I+1,1)=HNEW(I,K)

      VNEWL(I,K)=CCC-C(I)*HNEW(I,K)

      VNEWR(I+1,1)=CC+C(I+1)*HNEW(I+1,1)

      GO TO 40

 41 VNEWL(I,K)=CCC-C(I)*HCAV(I,K)

      VNEWR(I+1,1)=CC+C(I+1)*HCAV(I+1,1)

      HNEW(I,K)=HCAV(I,K)

      HNEW(I+1,1)=HCAV(I+1,1)

      DELCAV=DELT*(VNEWR(I+1,1)*AREA(I+1)-VNEWL(I,K)*AREA(I))

      IF(DELCAV.GT.0.) GO TO 43

      IF(AIRVAC(I+1).EQ.2.AND.HCAV(I,K).LE.PIPEZ(I+1).AND.

     $DELCAV.LT.0.) CTZERO(I+1)=CAV(I,K)

      IF(CTZERO(I+1).LT.(5.*AREA(I+1))) GO TO 43

      IF(AIRVAC(I+1).EQ.2.AND.HCAV(I,K).LE.PIPEZ(I+1).AND.

     $DELCAV.LT.0.) BUBBLE(I+1)=.TRUE.

 43 CAV(I,K)=CAV(I,K)+DELCAV

      CAV(I+1,1)=CAV(I,K)

      IF(.NOT. BUBBLE(I+1)) GO TO 42

      HCAV(I,K)=PIPEZ(I+1)-HATM+HATM*(CTZERO(I+1)/CAV(I,K))**1.2

      IF(HCAV(I,K).LT.PIPEZ(I+1)) BUBBLE(I+1)=.FALSE.

      IF(HCAV(I,K).LT.PIPEZ(I+1)) HCAV(I,K)=PIPEZ(I+1)

      HCAV(I+1,1)=HCAV(I,K)

      IF(BUBBLE(I+1)) GO TO 40

 42 IF(CAV(I,K).GT.0.) GO TO 40

      DELH2=ABS((0.5*A(I+1)/32.2)*((AREA(I)/AREA(I+1))*VNEWL(I,K)-

     $VNEWR(I+1,1)))

      HNEW(I,K)=HCAV(I,K)+DELH2*2.0*AREA(I+1)*A(I)/(A(I+1)*AREA(I)

     $+AREA(I+1)*A(I))

      HNEW(I+1,1)=HNEW(I,K)

      VNEWL(I,K)=VNEWL(I,K)-(32.2/A(I))*(HNEW(I,K)-HCAV(I,K))

      VNEWR(I+1,1)=VNEWL(I,K)*AREA(I)/AREA(I+1)

      CAV(I,K)=0.

      CAV(I+1,1)=0.

 40 CONTINUE

*

* ** COMPUTE H AND V AT DOWNSTREAM END **

*     --THIS BOUNDARY CONDITION IS FOR A CONSTANT HEAD RESERVOIR--

*

      NV=NPAR(NPIPES)+1

      VLEFT=VL(NPIPES,NV)-RATIO(NPIPES)*A(NPIPES)*(VL(NPIPES,NV)-

     $VR(NPIPES,NV-1))

      HLEFT=H(NPIPES,NV)-RATIO(NPIPES)*A(NPIPES)*(H(NPIPES,NV)-

     $H(NPIPES,NV-1))

      CC=VLEFT+C(NPIPES)*HLEFT+C(NPIPES)*VLEFT*DELT*SINE(NPIPES)-

     $AK(NPIPES)*VLEFT*ABS(VLEFT)**NEXP

      HNEW(NPIPES,NV)=HRES

      VNEWL(NPIPES,NV)=CC-C(NPIPES)*HNEW(NPIPES,NV)

*

*  ** COMPUTE H AND V AT UPSTREAM END **

*

*           BOUNDARY CONDITIONS AT THE PUMP

*

      IF(T.GT.TTTT) RSTART=.TRUE.

      IF(RSTART) BYPASS=.FALSE.

      IF(RSTART) GO TO 554

      IF(BYPASS) GO TO 554

556 QQ=Q/RPM

      IF(QQ.LE.0.) GO TO 552

      DO 550 I=1,5

      IF(QQ.GT.QN(I).AND.QQ.LE.QN(I+1)) GO TO 551

550 CONTINUE

      ICOUNT=ICOUNT+1

      IF(ICOUNT.LT.5) WRITE(6,260)

      I=5

      TN2=TNSQ(6)

      GO TO 553

551 IF(RSTART) GO TO 554

      TN2=TNSQ(I)+(QQ-QN(I))*(TNSQ(I+1)-TNSQ(I))/(QN(I+1)-QN(I))

      GO TO 553

552 TN2=TNSQ(1)

      I=1

553 RPM=RPM-30.*TN2*DELT/(PI*PINERT)*RPM*RPM

      IF(RPM.LT.0.) BYPASS=.TRUE.

554 VRITE=VR(1,1)-RATIO(1)*A(1)*(VR(1,1)-VL(1,2))

      HRITE=H(1,1)-RATIO(1)*A(1)*(H(1,1)-H(1,2))

      C1=VRITE-C(1)*HRITE-C(1)*DELT*VRITE*SINE(1)-AK(1)*VRITE*

     $ABS(VRITE)**NEXP

*

*  IF PUMP RESTART HAS BEGUN, EXIT THE POWER FAILURE MODE OF

*  OPERATION

*

      IF(RSTART) GO TO 650

      IF(.NOT.BYPASS) GO TO 555

      CCC=-HSUMP-C1/C(1)

      VPP=(-CCB+SQRT(CCB**2-4.*CCA*CCC))/(2.*CCA)

      IF(VPP.GT.0.) GO TO 557

      VPP=0.

      HPP=-C1/C(1)

      GO TO 530

557 HPP=(VPP-C1)/C(1)

      GO TO 530

555 DO 500 NZ=1,20

      HA=HNSQ(I)

      HB=HNSQ(I+1)

      QA=QN(I)

      QB=QN(I+1)

      C55=NSTAGE*RPM*RPM

      C77=(HA-HB)/(QA-QB)

      C88=HB-QB*C77

      C66=HSUMP+C55*C88+C55*C77*C1*AREA(1)/(NPUMPS*RPM)

      C99=C55*C77*C(1)*AREA(1)/(NPUMPS*RPM)

      HPP=C66/(1.-C99)

      IF(HPP.GT.HSUMP) GO TO 558

      CCC=-HSUMP-C1/C(1)

      VPP=(-CCB+SQRT(CCB**2-4.*CCA*CCC))/(2.*CCA)

      IF(VPP.GT.0.) GO TO 559

      VPP=0.

      HPP=-C1/C(1)

      GO TO 530

559 HPP=(VPP-C1)/C(1)

      BYPASS=.TRUE.

      GO TO 530

558 VPP=C1+HPP*C(1)

      IF(VPP.LE.0.) VPP=0.

      IF(VPP.LE.0.) HPP=-C1/C(1)

      QPP=VPP*AREA(1)/(NPUMPS*RPM)

      DO 510 J=1,5

      IF(QPP.GT.QN(J).AND.QPP.LE.QN(J+1)) GO TO 520

510 CONTINUE

      IF(QPP.GT.0..AND.I.EQ.5) GO TO 530

      IF(QPP.GT.0..AND.I.LT.5) I=I+1

      IF(QPP.LE.0..AND.I.EQ.1) GO TO 530

      IF(QPP.LE.0..AND.I.GT.1) I=I-1

      GO TO 500

520 IF(J.EQ.I) GO TO 530

      I=J

500 CONTINUE

      WRITE(6,270)

      STOP

530 HNEW(1,1)=HPP

      VNEWR(1,1)=VPP

      IF(BYPASS) GO TO 531

      IF(VPP.GT.0.) HPUMP=HPP-HSUMP

      IF(HPP.LT.HSUMP) HPUMP=0.

      IF(HPP.LT.HSUMP) BYPASS=.TRUE.

      IF(VPP.LE.0.) HPUMP=HNSQ(1)*RPM*RPM*NSTAGE

      Q=QPP*RPM

      GO TO 531

*......................................................................

*   PUMPS RESTARTING

*......................................................................

650 CONTINUE

      ICHEK=0

*

* COMPUTE SPEED OF EACH STARTUP PUMP

*

      DO 651 I=1,NSTART

      IF(T.LT.TSTART(I)) GO TO 651

      NACT=I

      IF(T.GT.(TSTART(I)+TRAMP(I))) GO TO 652

      RPMS(I)=RPMSTRT(I)+(T-TSTART(I))*(RPMEND(I)-RPMSTRT(I))/TRAMP(I)

      GO TO 651

652 RPMS(I)=RPMEND(I)

651 CONTINUE

*

* COMPUTE LOCATION ON PUMP CURVES FOR FIRST TRY

*

      DO 655 J=1,NSTART

      IF(RPMS(J).LE.0.) GO TO 655

      QP(J)=QP(J)/RPMS(J)

      DO 656 I=1,5

      IF(QP(J).GE.QN(I).AND.QP(J).LE.QN(I+1)) GO TO 657

656 CONTINUE

      I=5

      IF(QP(J).LE.0.) I=1

657 IZ(J)=I

      HA=HNSQ(I)

      HB=HNSQ(I+1)

      QA=QN(I)

      QB=QN(I+1)

      C7(J)=(HA-HB)/(QA-QB)

      C8(J)=HB-QB*C7(J)

655 CONTINUE

*

* COMPUTE TERMS REQUIRED TO GENERATE SOLUTIONS FOR PUMP STARTUP

*

658 SUM7=0.

      SUM8=0.

      DO 660 I=1,NACT

      IF(RPMS(I).LE.0.) GO TO 660

      SUM7=SUM7+1.0/(C7(I)*RPMS(I))

      SUM8=SUM8+RPMS(I)*C8(I)/C7(I)

660 CONTINUE

*

* FIND SOLUTION FOR PUMP STARTUP

*

      C15=C(1)*AREA(1)-SUM7/NSTAGE

      C16=-C1*AREA(1)-HSUMP*SUM7/NSTAGE-SUM8

      HPP=C16/C15

      VPP=C1+C(1)*HPP

      NACT2=NACT

      DO 665 I=1,NACT

      IF(RPMS(I).LE.0.) GO TO 665

      QP(I)=(RPMS(I)/C7(I))*((HPP-HSUMP)/(NSTAGE*RPMS(I)**2)

     $-C8(I))

      IF(QP(I).LE.0.) NACT2=NACT2-1

      IF(QP(I).LE.0.) QP(I)=0.

665 CONTINUE

      IF(NACT2.EQ.0) GO TO 654

      IF(NACT2.EQ.NACT) GO TO 654

      NACT=NACT2

      GO TO 658

654 IF(VPP.LE.0.) VPP=0.

      IF(VPP.LE.0.) HPP=-C1/C(1)

*

* CHECK TO SEE IF Q/N ON PROPER SEGMENT

*

689 CONTINUE

      NOGOOD=.FALSE.

      DO 685 IZZ=1,NSTART

      IF(RPMS(IZZ).LE.0.) GO TO 685

      QPP=QP(IZZ)/RPMS(IZZ)

      I=IZ(IZZ)

      DO 686 J=1,5

      IF(QPP.GT.QN(J).AND.QPP.LE.QN(J+1)) GO TO 688

686 CONTINUE

      IF(QPP.GT.0..AND.I.EQ.5) GO TO 685

      IF(QPP.GT.0..AND.I.LT.5) IZ(IZZ)=I+1

      IF(QPP.LE.0..AND.I.EQ.1) GO TO 685

      IF(QPP.LE.0..AND.I.GT.1) IZ(IZZ)=I-1

      GO TO 684

688 IF(J.EQ.I) GO TO 685

      IZ(IZZ)=J

684 NOGOOD=.TRUE.

685 CONTINUE

      ICHEK=ICHEK+1

      IF(ICHEK.GT.20) WRITE(6,270)

      IF(ICHEK.GT.20) STOP

      IF(NOGOOD) GO TO 691

      IF(HPP.GT.HSUMP) GO TO 699

      CCC=-HSUMP-C1/C(1)

      VPP=(-CCB+SQRT(CCB**2-4.*CCA*CCC))/(2.*CCA)

      HPP=(VPP-C1)/C(1)

      BYPASS=.TRUE.

      GO TO 699

*

* RECOMPUTE SOLUTION FOR REVISED POINTS ON CURVE

*

691 DO 690 J=1,NSTART

      IF(RPMS(J).LE.0.) GO TO 690

      QPP=QP(J)/RPMS(J)

      I=IZ(J)

      HA=HNSQ(I)

      HB=HNSQ(I+1)

      QA=QN(I)

      QB=QN(I+1)

      C7(J)=(HA-HB)/(QA-QB)

      C8(J)=HB-QB*C7(J)

690 CONTINUE

      GO TO 658

*

* SOLUTION IS SUCCESSFUL. COMPUTE FINAL VALUES.

*

699 HNEW(1,1)=HPP

      VNEWR(1,1)=VPP

      IF(BYPASS) GO TO 696

      IF(VPP.GT.0.) HPUMP=HPP-HSUMP

*

* ** CHECK CONTINUITY AT PUMPS **

*

      QPUMPS=0.

      DO 695 I=1,NSTART

      IF(RPMS(I).LE.0.) GO TO 695

      QPUMPS=QPUMPS+QP(I)

695 CONTINUE

      QLINE=VPP*AREA(1)

      IF(ABS(QLINE).LE..0001.AND.ABS(QPUMPS).LE..0001) GO TO 696

      IF(ABS((QPUMPS-QLINE)/QLINE).LT.0.01) GO TO 696

      WRITE(*,271) T

      STOP

696 CONTINUE

*

* ** LOCATE MAXIMUM AND MINIMUM HEADS AND H-VALUES **

*

531 DO 61 I=1,NPIPES

      K=NPAR(I)+1

      DO 61 J=1,K

      IF(J.EQ.1) GO TO 532

      DELTX=L(I)/(NPAR(I)*(ABS(VNEWL(I,J))+A(I)))

      IF(DELTX.LT.DELTN) DELTN=DELTX

532 CONTINUE

      IF(J.EQ.K) GO TO 533

      DELTX=L(I)/(NPAR(I)*(ABS(VNEWR(I,J))+A(I)))

      IF(DELTX.LT.DELTN) DELTN=DELTX

533 CONTINUE

      IF(HNEW(I,J).LT.HLOW(I,J)) TLOW(I,J)=T

      IF(HNEW(I,J).LT.HLOW(I,J)) HLOW(I,J)=HNEW(I,J)

      IF(HNEW(I,J).GT.HHIGH(I,J)) THIGH(I,J)=T

      IF(HNEW(I,J).GT.HHIGH(I,J)) HHIGH(I,J)=HNEW(I,J)

      IF(CAV(I,J).GT.CAVMAX(I,J)) TCAV(I,J)=T

      IF(CAV(I,J).GT.CAVMAX(I,J)) CAVMAX(I,J)=CAV(I,J)

      HEAD(I,J)=HNEW(I,J)-(PIPEZ(I)+X(I,J)*L(I)*SINE(I))

      IF(HEAD(I,J).LT.PHMAX) GO TO 66

      IPMAX=I

      XMAX=X(I,J)

      PHMAX=HEAD(I,J)

      TTMAX=T

      GO TO 67

 66 IF(HEAD(I,J).GT.PHMIN) GO TO 67

      IPMIN=I

      XMIN=X(I,J)

      PHMIN=HEAD(I,J)

      TTMIN=T

 67 CONTINUE

 61 CONTINUE

*

* ** WRITE OUT H,V,AND HEAD VALUES AND TEST FOR TMAX **

*

      IF(FAIL) GO TO 71

      IF(T.LT.TIOU(III)) GO TO 70

      IOUT=IOU(III)

      III=III+1

 70 CONTINUE

      IF(MOD(II,IOUT).NE.0) GO TO 72

 71 WRITE(6,205) T

      DO 73 I=1,NPIPES

      K=NPAR(I)+1

      WRITE(6,206) I,(X(I,J),HEAD(I,J),HNEW(I,J),VNEWL(I,J),

     $VNEWR(I,J),CAV(I,J),J=1,K)

 73 CONTINUE

      IF(BYPASS.AND.VNEWR(1,1).GT.0.) WRITE(6,215)

      IF(BYPASS) GO TO 72

      IF(RSTART) GO TO 75

      QQQ=449.*Q

      WRITE(6,212) RPM,QQQ,HPUMP

      GO TO 72

 75 DO 76 I=1,NSTART

      IF(RPMS(I).LE.0.) GO TO 76

      QQQ=449.*QP(I)

      WRITE(6,212) RPMS(I),QQQ,HPUMP

 76 CONTINUE

 72 CONTINUE

      IF(SAVE) CALL GRAFCAV(II,NP,IX,INDEX,DELT,VNEWL,VNEWR,HEAD,CAV,

     $NPAR,IENTRY)

      IF(FAIL) GO TO 400

*

* ** PREPARE FOR NEXT TIME STEP COMPUTATION **

*

      DO 80 I=1,NPIPES

      K=NPAR(I)+1

      DO 80 J=1,K

      VL(I,J)=VNEWL(I,J)

      VR(I,J)=VNEWR(I,J)

 80 H(I,J)=HNEW(I,J)

      IF(T.GT.TMAX) GO TO 400

 99 CONTINUE

*

*-------------------------------------------------------------------------------

*             END OF TRANSIENT ANALYSIS

*-------------------------------------------------------------------------------

*

* ** CHECK STABILITY CRITERIA **

*

400 IF(DELTN.GE.(0.9999*DELT)) GO TO 402

      DELT=0.999*DELTN

      WRITE(*,216) DELT

      WRITE(6,216) DELT

      GO TO 98

*

* ** WRITE OUT MAXIMUM AND MINIMUM HEADS AND H-VALUES **

*

402 CONTINUE

      WRITE(6,1011) CH12

      WRITE(6,207)

      DO 401 I=1,NPIPES

      WRITE(6,208) I

      K=NPAR(I)+1

      DO 401 J=1,K

      HEADMX=HHIGH(I,J)-(PIPEZ(I)+X(I,J)*L(I)*SINE(I))

      HEADMN= HLOW(I,J)-(PIPEZ(I)+X(I,J)*L(I)*SINE(I))

401 WRITE(6,209) X(I,J),HEADMX,THIGH(I,J),HEADMN,TLOW(I,J),

     $HHIGH(I,J),HLOW(I,J),CAVMAX(I,J),TCAV(I,J)

      PMAX=PHMAX/2.31

      PMIN=PHMIN/2.31

      WRITE(6,210) PHMAX,PMAX,IPMAX,XMAX,TTMAX

      WRITE(6,211) PHMIN,PMIN,IPMIN,XMIN,TTMIN

      IF(PFILE) WRITE(6,1011) CH12

      IF(PFILE) CALL PROFILE(NPIPES,L,NPAR,PIPEZ,SINE,X,

     $HLOW,HHIGH,HSTEAD)

      IENTRY=2

      IF(HVPRNT) WRITE(6,1011) CH12

      IF(HVPRNT) CALL GRAFCAV(II,NP,IX,INDEX,DELT,VNEWL,VNEWR,HEAD,CAV,

     $NPAR,IENTRY)

      IENTRY=3

      WRITE(6,1011) CH,CH79

      IF(PPLOT) WRITE(6,1011) CH12

      IF(PPLOT) CALL GRAFCAV(II,NP,IX,INDEX,DELT,VNEWL,VNEWR,HEAD,CAV,

     $NPAR,IENTRY)

      IENTRY=4

      IF(GRAPH) WRITE(6,1011) CH12

      IF(GRAPH) CALL GRAFCAV(II,NP,IX,INDEX,DELT,VNEWL,VNEWR,HEAD,CAV,

     $NPAR,IENTRY)

 98 CONTINUE

*

*               ******* FORMAT STATEMENTS *******

*

100 FORMAT(A)

1011 FORMAT(3A)

101 FORMAT(5X,'NUMBER OF PIPES =',I3/5X,'MINIMUM NO. OF PARTS PIPE IS

     $DIVIDED INTO =',I3//5X,'RESERVOIR ELEVATION AT DOWNSTREAM END =',

     $F7.1,' FT'/5X,'ELEVATION OF PUMP SUMP =',F7.1,' FT'//

     $5X,'ATMOSPHERIC PRESSURE HEAD =',F5.1,' FT OF WATER'

     $//5X,'ELEVATION OF DOWNSTREAM END OF LAST PIPE =',F7.1,' FT'//5X,

     $'MAXIMUM REAL TIME OF ANALYSIS =',F7.1,' SEC'//

     $5X,'ACCURACY OF ITERATION FOR STEADY STATE DISCHARGE =',F5.1,

     $' GPM'//5X,'PUMP BYPASS DIAMETER = ',F6.2,' IN'/5X,

     $'PUMP BYPASS LOSS COEFFICIENT =',F6.2)

1010 FORMAT(/5X,'PRINT OUTPUT SCHEDULE IS IOUT =',I8,I5,4I8/

     $5X,'CORRESPONDING INTERVALS (SEC) =',F4.1,5F8.2)

102 FORMAT(//22X,'PIPE INPUT DATA'/22X,15('-')//5X,'PIPE DIAM-IN LENGT

     $H-FT WAVE SPD-FPS PIPEZ-FT F-VALUE VEL-FPS'/

     $5X,4('-'),1X,7('-'),1X,9('-'),1X,12('-'),1X,8('-'),1X,7('-'),

     $1X,7('-'))

1022 FORMAT(//22X,'PIPE INPUT DATA'/22X,15('-')//5X,'PIPE DIAM-IN LENGT

     $H-FT WAVE SPD-FPS PIPEZ-FT C-VALUE VEL-FPS'/

     $5X,4('-'),1X,7('-'),1X,9('-'),1X,12('-'),1X,8('-'),1X,7('-'),

     $1X,7('-'))

1020 FORMAT(/5X,'PIPE DELT-SEC PARTS   SINE  L/A-SEC INTERPOLATION',1X,

     $'A-V VALVE'/5X,4('-'),1X,8('-'),1X,5('-'),2X,6('-'),1X,7('-'),

     $1X,13('-'),1X,9('-'))

103 FORMAT(5X,I3,2X,F6.2,1X,F9.1,6X,F6.0,4X,F6.0,3X,F5.4,3X,

     $F5.2)

1033 FORMAT(5X,I3,2X,F6.2,1X,F9.1,6X,F6.0,4X,F6.0,3X,F5.0,3X,F5.2)

1030 FORMAT(5X,I3,2X,F6.3,3X,I3,2X,F7.5,1X,F6.2,5X,F5.3,9X,A)

1040 FORMAT(//5X,'DELTA T =',F9.5,' SEC'//)

104 FORMAT(///15X,'* PUMP DATA *'/)

 105 FORMAT(5X,'NUMBER OF PUMPS IN PARALLEL =',I3/5X,'EACH PUMP HAS',I2

     $,' STAGES'/5X,'STEADY STATE PUMP SPEED =',F7.1,' RPM'/5X,'EACH PUM

     $P AND MOTOR UNIT HAS A MOMENT OF INERTIA OF',F6.1,' LB-FT SQ'//5X

     $'PUMP SPEED FOR DATA SHOWN BELOW IS',F7.1,' RPM')

106 FORMAT(/8X,'Q-GPM  H/STAGE-FT  BHP/STAGE'/8X,5('-'),2X,10('-'),2X,

     $9('-'))

107 FORMAT(6X,F7.1,2X,F8.1,4X,F8.2)

108 FORMAT(/3X,'PUMP NO.',I2,' STARTS AT',F8.2,' SEC AT',F6.0,

     $' RPM AND RAMPS TO',F6.0,' RPM IN',F7.2,' SEC')

 1088 FORMAT(//5X,'* STEADY STATE LINE FLOW RATE =',F9.1,

     $' GPM AT A PUMP HEAD =',F7.1,' FT *'/)

 201 FORMAT(//10X,A)

2011 FORMAT(10X,A//)

  202 FORMAT(21X,14('*')/21X,'* INPUT DATA *'/21X,14('*')/)

  204 FORMAT(/'  PRESSURE HEADS, H-VALUES AND VELOCITIES AS FUNCTIONS

     $ OF TIME'/2X,60('-'))

205 FORMAT(//20X,4X,'  X   HEAD-FT  H-FT   VL-FPS VR-FPS CAV-FT3'/

     $' TIME =',F9.3,' SEC',4X,'----- ------- -----   ------ ------',

     $1X,'-------')

206 FORMAT(/14X,'PIPE',I2,4X,F5.3,2F7.0,2X,2F7.2,F8.2/

     $(24X,F5.3,2F7.0,2X,2F7.2,F8.2))

207 FORMAT(//29X,27('*')/29X,'* TABLE OF EXTREME VALUES *'/29X,

     $27('*')//11X,'X    MAX HEAD  TIME   MIN HEAD  TIME    '

     $'MAX H    MIN H   CAVMAX  TIME'/9X,5('-'),2X,8('-'),2X,

     $4('-'),3X,8('-'),2X,4('-'),3X,6('-'),3X,6('-'),3X,6('-'),

     $2X,4('-'))

208 FORMAT(1X,'PIPE',I3)

209 FORMAT(9X,F5.3,2X,F7.1,1X,F6.1,3X,F7.1,1X,F6.1,2X,F7.1,2X,

     $F7.1,2X,F6.1,1X,F6.1)

210 FORMAT(//' MAXIMUM HEAD =',F6.1,' FT (',F5.1,' PSI) IN PIPE',

     $I3,' AT X =',F5.3,' AT TIME =',F7.2,' SEC')

211 FORMAT(/' MINIMUM HEAD =',F6.1,' FT (',F5.1,' PSI) IN PIPE',

     $I3,' AT X =',F5.3,' AT TIME =',F7.2,' SEC')

212 FORMAT(/5X,'PUMP SPEED =',F7.1,' RPM',5X,'PUMP DISCHARGE =',F7.1

     $,' GPM EACH',5X,'PUMP HEAD =',F6.1,' FT')

213 FORMAT(//2X,'THE DIMENSIONED SIZE ',I4,' OF THE ARRAYS X(,), V(,),

     $... IS TOO SMALL FOR THE NPARTS YOU SELECTED.'/2X,'YOU MUST DECREA

     $SE NPARTS OR INCREASE THE SIZE OF THE ARRAYS TO AT LEAST ',I4)

214 FORMAT(////2X,'COLUMN SEPARATION HAS OCCURRED AT',F6.2,' SEC IN PI

     $PE',I3,' AT LOCATION',F5.3)

215 FORMAT(/15X,'PUMP HEAD HAS DROPPED TO ZERO AND BYPASS IS OPEN')

216 FORMAT(/' THE CHECK FOR STABILITY CRITERIA HAS FOUND A VIOLATION.'

     $/' THE ANALYSIS MUST BE RERUN WITH A SMALLER DELTA T.'/

     $' CHANGE DATA FILE TO SET DTNEW =',F7.4,' SEC & RERUN = "T"')

250 FORMAT(/////10X,'YOU HAVE EXCEEDED MAXIMUM DISCHARGE VALUE INPUT.

     $EXECUTION IS TERMINATED.')

251 FORMAT(/////10X,'ITERATION TO COMPUTE STEADY STATE DISCHARGE WAS U

     $NSUCCESSFUL AFTER 20 TRIES')

260 FORMAT(10X,'YOU HAVE EXCEEDED MAXIMUM BHP VALUE INPUT.')

270 FORMAT(/////10X,'THE SEARCH FOR A VALUE OF Q HAS CYCLED 20 TIMES W

     $ITHOUT SUCCESS. EXECUTION IS TERMINATED.')

271 FORMAT(//' CONTINUITY CHECK AT THE PUMPS HAS FAILED AT T ='

     $,F10.4,' SEC')

790 FORMAT(A)

791 FORMAT(/' ENTER THE NAME OF YOUR INPUT DATA FILE: '\)

792 FORMAT(/' ENTER THE NAME OF THE FILE ON WHICH YOU WANT YOUR OUTPUT

     $ WRITTEN: '\)

      END