C     SUBROUTINE  **READ**                                                      
      SUBROUTINE READ                                                           
      INTEGER    A                                                              
      DIMENSION  ZZ(16000),I(3),E(2),J(2),K(2),A(7)                             
C                                                                               
      EQUIVALENCE (ZZ(11053),I(1)),                (ZZ(11056),E(1))             
     X,          (ZZ(11058),J(1)),                (ZZ(11060),K(1))              
     X,          (ZZ(11062),A(1) ),               (ZZ(11069),LFLAG)             
     X,          (ZZ(11090),IXX),                 (ZZ(11099),IZ)                
     X,          (I(1),I1),(I(2),I2),(I(3),I3)                                  
     X,          (J(1),J1),(J(2),J2)                                            
      COMMON  ZZ                                                                
C                                                                               
      IF(LFLAG.EQ.1) GO TO 10                                                   
      READ(5,1)AI1,I1,AI2,I2,I3,E,AJ1,J1,AJ2,J2,K,A                             
      GO TO 14                                                                  
   10 READ(5,12)AI1,I1,AI2,I2,I3,IZ,E,AJ1,J1,AJ2,J2,K,A                         
      IF(E(2).EQ.0.)E(2)=1.                                                     
   12 FORMAT(2(A1,I3),I4,I3,F9.0,F12.0,2(A1,I3),2I4,7A4)                        
   14 CONTINUE                                                                  
      I2S= I2                                                                   
      IF(LFLAG.EQ.-1) I3= MAX0(I3,1)                                            
      IF(LFLAG.EQ.1)IZ= MAX0(IZ,1)                                              
      IX= IXX+1                                                                 
      LP=0                                                                      
    1 FORMAT(2(A1,I3),I4,2F12.0,2(A1,I3),2I4,7A4)                               
    2 FORMAT (6X,2(A1,I3,1X),I4,2E15.6,2(1X,A1,I3),2I5,2X,5A4,2X,2A4)           
    9 FORMAT(1X,I4,1X,2(A1,I3,1X),2I4,1PE11.4,E15.6,2(1X,A1,I3),                
     12I5,2X,5A4,2X,2A4)                                                        
      IP=I1                                                                     
      I1 = LINK(AI1,I1)                                                         
      IF(I1.GE.0.AND.LFLAG.EQ.7) GO TO 6                                        
      IF(LFLAG.EQ.1.AND.I1.GE.0) LP=1                                           
      IF(LP.GT.0)WRITE(6,9)IX,AI1,IP,AI2,I2,I3,IZ,E,AJ1,J1,AJ2,J2,K,A           
      IF(LP.LE.0)WRITE(6,2)   AI1,IP,AI2,I2,I3,E,AJ1,J1,AJ2,J2,K,A              
      LP= 0                                                                     
      IF (I1) 6,3,3                                                             
    3 I2 = LINK(AI1,I2)                                                         
      IF(I2.LT.2999)I2= LINK(AI2,I2S)                                           
C     LFLAG = + FOR FUNIN                                                       
      IF (LFLAG) 6,6,4                                                          
    4 J1 = LINK(AJ1,J1)                                                         
      J2 = LINK(AJ2,J2)                                                         
    6 RETURN                                                                    
      END                                                                       
      SUBROUTINE SPILL                                                          
      DIMENSION ZZ(16000),DYTCQ(11030),CL(4,12)                                 
      EQUIVALENCE (ZZ(0002),DYTCQ(1))                                           
     X,          (ZZ(11062),N2)                                                 
      COMMON  ZZ                                                                
      DATA CL(1,1),CL(2,1),CL(3,1),CL(4,1)/4HADMI,4HTTAN,4HCES ,4H    /
     1,    CL(1,2),CL(2,2),CL(3,2),CL(4,2)/4HTEMP,4HERAT,4HURES,4H    /
     2,    CL(1,3),CL(2,3),CL(3,3),CL(4,3)/4HCAPA,4HCITI,4HES  ,4H    /
     3,    CL(1,4),CL(2,4),CL(3,4),CL(4,4)/4HGENE,4HRATI,4HON R,4HATES/
     4,    CL(1,5),CL(2,5),CL(3,5),CL(4,5)/4HCONS,4HTANT,4HS  (,4HD)  /
     5,    CL(1,6),CL(2,6),CL(3,6),CL(4,6)/4HCONS,4HTANT,4HS  (,4HA)  /
     6,    CL(1,7),CL(2,7),CL(3,7),CL(4,7)/4HCONS,4HTANT,4HS  (,4HB)  /
     7,    CL(1,8),CL(2,8),CL(3,8),CL(4,8)/4HINER,4HTL. ,4HVALU,4HES-L/
     8,    CL(1,9),CL(2,9),CL(3,9),CL(4,9)/4H*030,4HGAIN,4HS  (,4HG)  /
     9,CL(1,10),CL(2,10),CL(3,10),CL(4,10)/4HSPEC,4H.CON,4HSTAN,4HTS  /
      MAXSP= 1999                                                               
      ID=1                                                                      
      ISPIL= 0                                                                  
      WRITE (6,36)                                                              
   36 FORMAT (1H1,42X,11H DYTCQ DUMP/) 
    3 ISPIL = ISPIL+1                                                           
      IF(ISPIL.GT.10) GO TO 33                                                  
      MAXSP= MAXSP+ 1000                                                        
      IF(MAXSP.GT.3000)ID= MAXSP-998                                            
      IF(ISPIL.EQ.10)MAXSP=11020                                                
      IF(ISPIL.EQ.10)ID= 11001                                                  
      IF(N2.EQ.4.AND.ISPIL.EQ.8) GO TO 114                                      
   14 WRITE(6,10)   (CL(I,ISPIL),I=1,4)                                         
      GO TO 115                                                                 
  114 WRITE(6,110)                                                              
  110 FORMAT(9H0   LOC. /5X,'NUMBER',10X,'STEADY STATE RESIDUALS')              
  115 CONTINUE                                                                  
   10 FORMAT(9H0   LOC. / 10H    NUMBER,10X, 4A4/)                              
      IZFLG= 0                                                                  
      DTEST= DYTCQ(ID)                                                          
   15 IXL = ID + 9                                                              
      IF (MAXSP-ID) 28,16,16                                                    
   16 IF(DYTCQ(ID).EQ.DTEST)GO TO 20                                            
   17 IF (IZFLG) 18,18,31                                                       
   18 IMOD= ID                                                                  
      IF(ISPIL.GT.1)IMOD= MOD(ID,1000)                                          
      WRITE (6,19) IMOD,(DYTCQ(IX),IX=ID,IXL)                                   
   19 FORMAT(1H ,2X,I4,2X,1P10E10.3)
      ID = IXL+1                                                                
      DTEST= DYTCQ(ID)                                                          
      GO TO 15                                                                  
   20 IZ = ID+1                                                                 
   21 IF(DYTCQ(IZ).NE.DTEST) GO TO 17                                           
   22 IF (IZ-IXL) 23,24,24                                                      
   23 IZ = IZ+1                                                                 
      GO TO 21                                                                  
   24 IF (IZFLG) 25,25,26                                                       
   25 IZFLG = ID                                                                
   26 ID = IXL+1                                                                
      GO TO 15                                                                  
   28 IF (IZFLG) 3,3,29                                                         
   29 JMOD = MOD(IZFLG,1000)                                                    
      LMOD = MOD(MAXSP,1000)                                                    
      IF(ISPIL.EQ.1) JMOD=IZFLG                                                 
      IF(ISPIL.EQ.1) LMOD= MAXSP                                                
      WRITE(6,32)JMOD,LMOD,DTEST                                                
      GO TO 3                                                                   
   31 IZE= IXL-10                                                               
      JMOD = MOD(IZFLG,1000)                                                    
      KMOD = MOD(IZE,1000)                                                      
      IF(ISPIL.EQ.1) JMOD= IZFLG                                                
      IF(ISPIL.EQ.1) KMOD= IZE                                                  
      WRITE(6,32) JMOD,KMOD,DTEST                                               
   32 FORMAT(1H0,34X,'LOCATIONS ',I4,' THROUGH ',I4,' EQUAL ',1PE12.5 /)        
      IZFLG = 0                                                                 
      GO TO 18                                                                  
   33 RETURN                                                                    
      END                                                                       
C     SUBROUTINE **LATENT**                                                     
      SUBROUTINE LATENT(NTEST)                                                  
      INTEGER NY1(3000),NY2(1000),NZA(1000),INXP(1000),INWP(1000)               
     X, NXP(1000),KEY(1000),IPR(500),IPH(100),NTHEAD(100),NTTAIL(100)           
     X,NFLAG(1000),IZ(1000)                                                     
      DIMENSION  ZZ(16000),T(1000),TOLD(1000),C(1000)                           
     X,PHT1(100),PHH1(100),PHT2(100),PHH2(100),QCH1(100),QCH2(100)              
     X,Q(999)       ,SUMY(999),SUMYT(999)                                       
        COMMON ZZ                                                               
       COMMON/I2BYT/NY1,NY2,NZA,NXP,INXP,INWP,KEY,IPR,IPH                       
     X,NTHEAD,NTTAIL,NFLAG,IZ                                                   
      EQUIVALENCE(ZZ(3002),T(1)),(ZZ(01),TIME),   (ZZ(4002),C(1))               
     X,          (ZZ(11102),PHT1(1)),             (ZZ(11072),DTHETA)            
     X,          (ZZ(11081),IPHMAX),              (ZZ(11082),IPRMAX)            
     X,          (ZZ(11202),PHH1(1) ),            (ZZ(11302),PHT2(1) )          
     X,          (ZZ(11402),PHH2(1) ),            (ZZ(11502),QCH1(1) )          
     X,          (ZZ(11602),QCH2(1) )                                           
     X,          (ZZ(5002),Q(1)),                 (ZZ(11100),ODTIME)            
     X,          (ZZ(12002),SUMY(1)),             (ZZ(13002),TOLD(1))           
     X,          (ZZ(14002),SUMYT(1)),            (ZZ(11083),IPFLAG)            
     X,          (ZZ(11062),NT),                  (ZZ(11005),CON4)              
      DATA COOL,HEAT/4HCOOL,4HHEAT/                                             
    1 DO 90 I =1,IPHMAX                                                         
        K = IPH(I)                                                              
      IF(NFLAG(K).NE.1) GO TO 90                                                
      IF(NTEST.EQ.1) GO TO 200                                                  
      IF(PHH1(I).EQ.0.) GO TO 90                                                
    2   N2 = 1                                                                  
      DT= T(K)-TOLD(K)                                                          
        IF (DT) 6,90,7                                                          
    6 A = COOL                                                                  
         TXS = T(K) -PHT2(I)                                                    
      TXSP = TOLD(K)- PHT2(I)                                                   
         N3 = 2                                                                 
         GO TO 10                                                               
    7 A = HEAT                                                                  
         TXS = T(K) -PHT1(I)                                                    
      TXSP = TOLD(K)- PHT1(I)                                                   
         N3 =1                                                                  
   10   TEST = TXS*TXSP                                                         
        IF (TEST) 11,11,80                                                      
   11   GO TO (12,13),N3                                                        
   12    QCH = QCH1(I)                                                          
         PHT = PHT1(I)                                                          
         PHH = PHH1(I)                                                          
         GO TO 14                                                               
   13    QCH = QCH2(I)                                                          
         PHT = PHT2(I)                                                          
         PHH = PHH2(I)                                                          
   14   IF(TEST) 15,16,80                                                       
   15 TIMEP = TIME - DTHETA*(TXS/DT)                                            
      WRITE (6,100) K,N3,A,TIMEP,PHT                                            
  100 FORMAT (7H0  NODE,I4, 14H, PHASE CHANGE,I2,2H, ,A4,20HING BEGINS A        
     1T TIME =,F8.2,8H, TEMP =,F8.2)                                            
   16   QCH = QCH + C(K)*TXS                                                    
        IF (DT) 18,90,17                                                        
   17    IF (QCH -PHH) 19,20,20                                                 
   18    IF (QCH) 21,21,19                                                      
   19     T(K)= PHT                                                             
          N2  = 2                                                               
          GO TO 79                                                              
   20     DTP = (QCH -PHH)/C(K)                                                 
          QCH = PHH                                                             
          GO TO 78                                                              
   21     DTP =(QCH)/C(K)                                                       
          QCH =0.0                                                              
   78 TIMEP = TIME - DTHETA*(DTP/DT)                                            
         T(K)  = PHT + DTP                                                      
      WRITE (6,101) K,N3,A,TIMEP,PHT,QCH                                        
  101 FORMAT (7H0  NODE,I4,14H, PHASE CHANGE,I2,2H, ,A4,20HING ENDS   AT        
     1 TIME =,F8.2,8H, TEMP =,F8.2, 22H, INTEGRATED HEATING =,E11.4)            
   79 IF (N3 .EQ. 2) GO TO 62                                                   
   61     QCH1(I) =QCH                                                          
          GO TO 80                                                              
   62     QCH2(I) =QCH                                                          
   80 IF (N2 .EQ. 2) GO TO 90                                                   
   81 IF (N3 .EQ. 2) GO TO 83                                                   
   82     TXS = T(K) -PHT2(I)                                                   
      TXSP = TOLD(K)- PHT2(I)                                                   
          N3  = 2                                                               
          GO TO 84                                                              
   83     TXS = T(K) -PHT1(I)                                                   
      TXSP = TOLD(K)- PHT1(I)                                                   
           N3 = 1                                                               
   84      N2 = 2                                                               
           GO TO 10                                                             
  200 IF(PHH1(I).NE.0.) GO TO 90                                                
      TNEW= PHT1(I)                                                             
      TEST3= ABS(T(K)-TNEW)                                                     
      DTT= DTHETA+ODTIME                                                        
      IF(TEST3.LT.CON4) GO TO 90                                                
      TR= 1.+ ODTIME/DTHETA                                                     
      DENOM= TR*C(K)+ODTIME*SUMY(K)                                             
      TNEW2=(TR*C(K)*T(K)+DTT*SUMYT(K)-DTHETA*TOLD(K)*SUMY(K))/DENOM            
      TEST3 =(TNEW2- TNEW)*(TNEW-T(K))                                          
                                                                                
      IF(TEST3.LE.0.) GO TO 90                                                  
      DEN=TOLD(K)* SUMY(K)-SUMYT(K)                                             
      B1=C(K)*(TNEW-T(K))-ODTIME*(SUMYT(K)-TNEW*SUMY(K))                        
      B1= B1/DEN                                                                
      C1= ODTIME*C(K)*(TNEW-T(K))/DEN                                           
      TEST2= B1**2-4.*C1                                                        
      IF(TEST2.LT.0.) GO TO 90                                                  
      TSQRT= SQRT(TEST2)                                                        
      DTIME = 0.                                                                
      TG1= .5*(TSQRT-B1)                                                        
      TG2=-.5*(B1+TSQRT)                                                        
      IF(TG1.GT.0..AND.TG1.LT.DTHETA)DTIME = TG1                                
      IF(TG2.GT.0..AND.TG2.LT.DTHETA)DTIME = AMAX1(TG2,DTIME)                   
      IF(NT.EQ.2) DTIME =C(K)*(TNEW-T(K))/(SUMYT(K)-T(K)*SUMY(K))               
      IF(DTIME.LE.0..OR.DTIME.GT.DTHETA) GO TO 90                               
      TIME= TIME -DTHETA + DTIME                                                
      DTHETA = DTIME                                                            
      WRITE(6,210) K,TNEW,TIME                                                  
  210 FORMAT(7H0  NODE,I4,' REACHED A VALUE OF ',G12.5,' AT TIME=',             
     1G12.5)                                                                    
      IPFLAG= 1                                                                 
   90 CONTINUE                                                                  
      RETURN                                                                    
      END                                                                       
C     FUNCTION  **LINK**                                                        
C                                                                               
C THIS FUNCTION IDENTIFIES THE TYPE 0F INFORMATION READ IN - D,Y,T,C,Q          
C                                                                               
      FUNCTION LINK(LA,LB)                                                      
      DIMENSION ITEST(13)                                                       
C THE INPUT DATA PREFIXES D,Y,T,C,Q ARECONVERTED TO NUMERICAL PREFIXES          
C 0,1,2,3,4 FOR THE PROGRAM                                                     
C                                                                               
      DATA ITEST/1H*,1HY,1H0,1H ,1H1,1H2,1HT,1HC,1HQ,1HD,1HA,1HB,1HL/           
C                                                                               
      LAB=0                                                                     
      DO 1 I=1,13                                                               
      KTEST = ITEST(I)                                                          
      IF(LA.EQ.KTEST)GO TO(4,5,5,5,5,5,6,8,9,15,13,14,16),I                     
    1 CONTINUE                                                                  
      CALL ERROR                                                                
      GO TO 10                                                                  
C     PREFIX IS AN * (CONTROL CARD)                                             
    4 LAB = -LB                                                                 
      GO TO 10                                                                  
C     PREFIX IS Y,0, ,1,OR2 (ADMITTANCE VALUES)                                 
    5 LAB = LB                                                                  
      IF(I.GT.4)LAB=LAB+1000*(I-4)                                              
      GO TO 10                                                                  
C     PREFIX IS  T (TEMPERATURES) OR P (PRESSURES)                              
    6 LAB= 3000+ LB                                                             
      GO TO 10                                                                  
C PREFIX IS C (CAPACITANCE)                                                     
    8 LAB= 4000+ LB                                                             
      GO TO 10                                                                  
C PREFIX IS Q (HEAT GENERATION RATE)                                            
    9 LAB= 5000+ LB                                                             
      GO TO 10                                                                  
C     PREFIX IS A  CONSTANTS FOR CALCULATION                                    
   13 LAB=7000 + LB                                                             
      GO TO 10                                                                  
C     PREFIX IS B  CONSTANTS FOR CALCULATIONS                                   
   14 LAB=8000 + LB                                                             
      GO TO 10                                                                  
C     PREFIX IS D  DUMMY CONSTANTS FOR CALCULATIONS                             
   15 LAB=6000 + LB                                                             
      GO TO 10                                                                  
C     PREFIX IS U (MAXIMUM TEMPERATURES)                                        
   16 LAB= 9000 + LB                                                            
   10 LINK = LAB                                                                
      RETURN                                                                    
      END                                                                       
C     SUBROUTINE **TABXX**                                                      
C     THIS SUBROUTINE USES LINEAR INTERPOLATION IN EVALUATING TABLE             
C     FUNCTIONS                                                                 
      SUBROUTINE TABXX(Z,X,NTAB,NPRE,KX)                                        
      INTEGER NY1(3000),NY2(1000),NZA(1000),INXP(1000),INWP(1000)               
     X, NXP(1000),KEY(1000),IPR(500),IPH(100),NTHEAD(100),NTTAIL(100)           
     X,NFLAG(1000),IZ(1000),NPREX(100)                                          
      DIMENSION  ZZ(16000),TABL(1000)                                           
       COMMON/I2BYT/NY1,NY2,NZA,NXP,INXP,INWP,KEY,IPR,IPH                       
     X,NTHEAD,NTTAIL,NFLAG,IZ                                                   
      EQUIVALENCE(ZZ(11069), LFLAG),              (ZZ(11070), IERFL)            
     X,          (ZZ(001),TIME),                  (ZZ(15002),TABL(1))           
     X,          (ZZ(11087),NARITH),              (ZZ(11090),IXX)               
      COMMON ZZ                                                                 
      NTH=NTHEAD(NTAB)                                                          
      NTT=NTTAIL(NTAB)                                                          
      IF(LFLAG.GE.10) GO TO 20                                                  
      DO 25  J=1,100                                                            
   25 NPREX(J)= NTHEAD(J)                                                       
      LFLAG = 10                                                                
   20 NT = NPREX(NTAB)                                                          
      IF (X-TABL(NT)) 7,1,2                                                     
    1 Z = TABL(NT+1)                                                            
      GO TO 16                                                                  
    2 IF (NT-NTT) 3,12,12                                                       
    3 IF (X-TABL(NT+2))6,4,5                                                    
    4 NT=NT+2                                                                   
      GO TO 1                                                                   
    5 NT=NT+2                                                                   
      GO TO 2                                                                   
    6 SLOPE=(TABL(NT+3)-TABL(NT+1))/(TABL(NT+2)-TABL(NT))                       
      Z = TABL(NT+1)+SLOPE*(X-TABL(NT))                                         
      GO TO 16                                                                  
    7 IF (NT-NTH) 12,12,8                                                       
    8 IF (X-TABL(NT-2)) 10,9,11                                                 
    9 NT = NT-2                                                                 
      GO TO 1                                                                   
   10 NT = NT-2                                                                 
      GO TO 7                                                                   
   11 SLOPE=(TABL(NT+1)-TABL(NT-1))/(TABL(NT)-TABL(NT-2))                       
      Z = TABL(NT+1)-SLOPE*(TABL(NT)-X)                                         
      GO TO 16                                                                  
   12 IERFL=0                                                                   
      IF (KX) 13,17,13                                                          
   17 WRITE(6,18)IXX,X,NTAB                                                     
   18 FORMAT(1H0,12HIN OPERATION,I4,1H,,E12.5,'LIES OUTSIDE LIMITS OF ',
     X'TABLE ',I4)
       GO TO 16
   13 WRITE(6,14)IXX,KX,X,NTAB                                                  
   14 FORMAT(1H0,12HIN OPERATION,I4,1H,,I6,1H=,E12.5, 
     129H LIES OUTSIDE LIMITS OF TABLE,I4)  
   16 NPREX(NTAB) = NT                                                          
      RETURN                                                                    
      END                                                                       
C     SUBROUTINE TO INTEGRATE BETWEEN VALUES IN A TABLE TO OBTAIN A             
C       MEAN VALUE                                                              
      SUBROUTINE INTGRT(XMIN,XMAX,NDX,NTAB,Z)                                   
C                                                                               
      DIMENSION  A(100),X(100)                                                  
      INTEGER NY1(3000),NY2(1000),NZA(1000),INXP(1000),INWP(1000)               
     X, NXP(1000),KEY(1000),IPR(500),IPH(100),NTHEAD(100),NTTAIL(100)           
     X,NFLAG(1000),IZ(1000)                                                     
       COMMON/I2BYT/NY1,NY2,NZA,NXP,INXP,INWP,KEY,IPR,IPH                       
     X,NTHEAD,NTTAIL,NFLAG,IZ                                                   
C                                                                               
   10 ANDX = NDX                                                                
      K = NDX + 1                                                               
      NPREV = NTHEAD(NTAB)                                                      
C                                                                               
   20 DELTAX = (XMAX-XMIN)/ANDX                                                 
      X(1) = XMIN                                                               
      DO 30 I=2,K                                                               
      X(I) = X(I-1) + DELTAX                                                    
   30 CONTINUE                                                                  
C                                                                               
      DO 50 I1=1,K                                                              
      KX = 30000 + I1                                                           
      XX = X(I1)                                                                
      CALL TABXX(V,XX,NTAB,NPREV,KX)                                            
      A(I1) = V                                                                 
   50 CONTINUE                                                                  
C                                                                               
      SUMA = 0.0                                                                
      DO 100 I2=2,K                                                             
      DELA =((A(I2)+A(I2-1))/2.0) * DELTAX                                      
      SUMA = SUMA + DELA                                                        
  100 CONTINUE                                                                  
C                                                                               
      Z = SUMA/(XMAX-XMIN)                                                      
      RETURN                                                                    
      END                                                                       
C     SUBROUTINE  **HEADER**                                                    
      SUBROUTINE HEADER                                                         
      COMMON /LAB4/SPROG(2)                                                     
      COMMON /LAB2/DATE(3)                                                      
   10 WRITE (6,20)                                                              
   20 FORMAT(1H1/1H-/1H-,36X,35HROCKETDYNE THERMAL ANALYZER PROGRAM) 
      WRITE(6,30)SPROG                                                          
   30 FORMAT(1H0,41X,2A4,' AUGUST,1969' )                                       
      WRITE(6,40)                                                               
40     FORMAT(1H0,38X,'MODIFIED FOR U.S.C.SAP USERS GROUP 5/30/79')             
      WRITE(6,60)                                                               
   60 FORMAT(1H0,//19X,'THIS PROGRAM USES A THREE TIME LEVEL INTEGRATION        
     1 METHOD. THE THREE'/19X,'LEVEL METHOD CAN BE USED BY SETTIN               
     2G CON(12)=1.')                                                            
      WRITE(6,70) DATE                                                          
   70 FORMAT(1H0,//41X,'THIS RUN WAS MADE ON ',3A4 )                            
      RETURN                                                                    
      END                                                                       
C     SUBROUTINE  **BIVAR**                                                     
C           BIVARIATE TABLES                                                    
      SUBROUTINE BIVAR                                                          
      INTEGER NY1(3000),NY2(1000),NZA(1000),INXP(1000),INWP(1000)               
     X, NXP(1000),KEY(1000),IPR(500),IPH(100),NTHEAD(100),NTTAIL(100)           
     X,NFLAG(1000),IZ(1000)                                                     
       COMMON/I2BYT/NY1,NY2,NZA,NXP,INXP,INWP,KEY,IPR,IPH                       
     X,NTHEAD,NTTAIL,NFLAG,IZ                                                   
      DIMENSION  ZZ(16000),TABL(1000)                                           
      COMMON ZZ                                                                 
      EQUIVALENCE(ZZ(001),TIME),                  (ZZ(15002),TABL(1))           
      EQUIVALENCE(ZZ(11069), LFLAG),              (ZZ(11070), IERFL)            
     X,          (ZZ(11091),NTLAST),              (ZZ(11092),NWW)               
     X,          (ZZ(11093),NXX),                 (ZZ(11094),NZZ)               
     X,          (ZZ(11095), W ),                 (ZZ(11096), X )               
     X,          (ZZ(11097), Z ) ,                (ZZ(11098), NTAB)             
      DATA  ITDTAG/4HTIME/                                                      
      XX= X                                                                     
    1 NTABB = NTAB                                                              
      NTH = NTHEAD(NTABB)                                                       
      NTT = (NTTAIL(NTABB))-1                                                   
      IERFL = 1                                                                 
    2 IF (W-TABL(NTH)) 16,13,3                                                  
    3 IF (W-TABL(NTT))4,14,16                                                   
    4 NT = NTH+2                                                                
      DO 6 I=NT,NTT,2                                                           
    5 IF(W-TABL(I))7,50,6                                                       
    6 CONTINUE                                                                  
   50 NT = I                                                                    
      GO TO 15                                                                  
    7 NT=I                                                                      
      ASSIGN 10 TO IBV                                                          
      NTAB2 = TABL(NT-1)+0.1                                                    
      W2 = TABL(NT-2)                                                           
      W1 = TABL(NT)                                                             
    8 NTAB = TABL(NT+1)+0.1                                                     
      IF (NTHEAD(NTAB)) 31,31,29                                                
   29 IX=IXX                                                                    
      NN = NTHEAD(NTAB)                                                         
      CALL TABXX(ZX,XX,NTAB,NN,NWW)                                             
      Z = ZX                                                                    
      IF(IERFL.EQ.0)  IERFL = 2                                                 
      GO TO (9,27), IERFL                                                       
    9 GO TO IBV,(28,10)                                                         
   10 Z1 = Z                                                                    
      NTAB = NTAB2                                                              
      IF (NTHEAD(NTAB)) 31,31,30                                                
   30 IX=IXX                                                                    
      NN = NTHEAD(NTAB)                                                         
      CALL TABXX(ZX,XX,NTAB,NN,NWW)                                             
      Z = ZX                                                                    
      IF(IERFL.EQ.0)  IERFL = 2                                                 
      IF (IERFL-1)12,12,27                                                      
   12 Z = ((Z1-Z)*(W-W2)/(W1-W2))+Z                                             
      GO TO 28                                                                  
   13 NT=NTH                                                                    
   15 ASSIGN 28 TO IBV                                                          
      GO TO 8                                                                   
   14 NT=NTT                                                                    
      GO TO 15                                                                  
   16 IF (NWW) 19,25,19                                                         
  19   ITDTAG=NWW                                                               
  25   ITDVTG=NZZ                                                               
      WRITE (6,18) ITDVTG, ITDTAG, ITDTAG, W, NTAB                              
   18 FORMAT ('0',2X,'***IN THE EVALUATION OF ',A4,', AS A FUNCTION OF '        
     1,A4,', ',A4,' =',1E12.5,', WHICH EXCEEDS LIMITS OF TABLE',I3,'**')        
      GO TO 27                                                                  
   31 WRITE (6,32) NTAB                                                         
   32 FORMAT (29H0 A NONEXISTENT TABLE, NUMBER,I3, 17H, HAS BEEN CALLED)
   27 IERF = 2                                                                  
      ISFLAG = 1                                                                
   28 RETURN                                                                    
      END                                                                       
C     SUBROUTINE  **SETUP**                                                     
C           **** MODIFIED FEB 1968****                                          
      SUBROUTINE SETUP(JSTART)                                                  
      DIMENSION ZZ(16000),Q(1000),Y(3000),T(1000),C(1000),CON(20)               
     X,PHT1(100),PHT2(100),QCH1(100),PHH1(100),QCH2(100),PHH2(100)              
     X,PRTIM(10),IP(3000),DYTCQ(1)                                              
        INTEGER NY1(3000),NY2(1000),NZA(1000),INXP(1000),INWP(1000)             
      INTEGER NXP(1000),KEY(1000),NFLAG(1000),IZ(1000),IPR(500),IPH(100)        
       INTEGER NTHEAD(100),NTTAIL(100)                                          
       COMMON /I2BYT/NY1,NY2,NZA,NXP,INXP,INWP,KEY,IPR,IPH                      
     X,NTHEAD,NTTAIL,NFLAG,IZ                                                   
      COMMON /LAB1/ PLOT(100,10),IPLOT(10),TMAX(1000)                           
      COMMON /LAB2/DATE(2)                                                      
      EQUIVALENCE (ZZ(001),TIME),                 (ZZ(002),Y(1),DYTCQ(1)        
     X) ,        (ZZ(3002),T(1)) ,                (ZZ(4002),C(1))               
     X,          (ZZ(5002),Q(1))                                                
     X,      (ZZ(10001),LYMIN),     (ZZ(10501),LYMAX)                           
     X,          (ZZ(11002),CON(1)),              (ZZ(11010),CON9,SIGMA)        
     X,          (ZZ(11011),CON10,TZERO),         (ZZ(11032),PRTIM(1))          
      EQUIVALENCE(ZZ(11069), LFLAG),              (ZZ(11070), IERFL)            
     X,          (ZZ(11071),PRTMX),               (ZZ(11072),DTHETA)            
     X,          (ZZ(11073),RCMIN),               (ZZ(11074), JSAVE)            
     X,          (ZZ(11075),DELTAT),              (ZZ(11076),DTMAX )            
     X,          (ZZ(11077), IMIN),               (ZZ(11078), IMAX )            
     X,          (ZZ(11079), IXMIN),              (ZZ(11080), IXMAX)            
     X,          (ZZ(11081),IPHMAX),              (ZZ(11082),IPRMAX)            
     X,          (ZZ(11083),IPFLAG),              (ZZ(11084),ISFLAG)            
     X,          (ZZ(11085),MAXNOT),              (ZZ(11086),MAXNOY)            
     X,          (ZZ(11087),NARITH),              (ZZ(11088),NCRIT )            
      EQUIVALENCE (ZZ(11089),NTIME),              (ZZ(11090), NT )              
     X,          (ZZ(11091),NTLAST),              (ZZ(11092),NWW)               
     X,          (ZZ(11093),NXX),                 (ZZ(11094),NZZ)               
     X,          (ZZ(11095), W ),                 (ZZ(11096), X )               
     X,          (ZZ(11097), Z ) ,                (ZZ(11098), NTAB)             
     X,          (ZZ(11099),CFLAG),               (ZZ(11100),ODTIME)            
     X,          (ZZ(11101),DTIME ),              (ZZ(11102),PHT1(1) )          
     X,          (ZZ(11202),PHH1(1) ),            (ZZ(11302),PHT2(1) )          
     X,          (ZZ(11402),PHH2(1) ),            (ZZ(11502),QCH1(1) )          
     X,          (ZZ(11602),QCH2(1) )                                           
     X,          (ZZ(12002),IP(1)),     (ZZ(11056),NYMIN)                       
      COMMON ZZ                                                                 
C     SET LATENT HEAT DATA                                                      
  200 DO 206 I = 1,IPHMAX                                                       
        K = IPH(I)                                                              
        IF(T(K) -PHT1(I)) 201,203,202                                           
  201   QCH1(I) = 0.0                                                           
        GO TO 203                                                               
  202   QCH1(I) = PHH1(I)                                                       
  203   IF (T(K) -PHT2(I)) 204,206,205                                          
  204     QCH2(I) =0.0                                                          
          GO TO 206                                                             
  205     QCH2(I) = PHH2(I)                                                     
  206     CONTINUE                                                              
  310 IF((JSTART.EQ.0).AND.(CON(11).NE.0.)) GO TO 750                           
      LIMY=1                                                                    
        MEGA=1000000                                                            
      JLIMY=1                                                                   
      DO 350 J=1,MAXNOY                                                         
  350 IP(J)=0                                                                   
  410 WRITE (6,420)                                                             
  420 FORMAT(1H1,36X,28HSUMMARY OF NODE CONNECTIONS //
     1,'  NODE  INITIAL   NODE  CAPACITY    GENERATION              NODE
     2  NODE   CONNECTOR        INITIAL '/
     3,'    NO    VALUE    FLAG                 RATE                  NO'
     4 FLAG    NUMBER         ADMITTANCE '//)
  500 DO 690 J1=1,MAXNOT                                                        
      IF(NFLAG(J1).EQ.0) GO TO 690                                              
      DO 570 J2=LIMY,MAXNOY                                                     
      JT=NY1(J2)/MEGA                                                           
      IF(IP(J2).GE.2)GO TO 570                                                  
      IF(JT.GT.0)IP(J2)=2                                                       
      JL= MOD(NY1(J2),MEGA)/1000                                                
      JM= MOD(NY1(J2),1000)                                                     
      IF(JL.EQ.J1.OR.JM.EQ.J1)GO TO 580                                         
  570 CONTINUE                                                                  
      GO TO 690                                                                 
  580 LMIN=J2+1                                                                 
      JT=MOD(NY1(J2),MEGA)/1000                                                 
      IF(JT.EQ.J1)JT=MOD(NY1(J2),1000)                                          
      WRITE(6,610)J1,T(J1),NFLAG(J1),C(J1),Q(J1),JT,NFLAG(JT),J2,Y(J2)          
      IP(J2)=IP(J2)+1                                                           
  610 FORMAT(1H ,I5,1X,E12.4,I4,2E12.4,12X,2I5,5X,I6,7X,1PE12.4)    
      DO 630 L2=LMIN,MAXNOY                                                     
      JT= NY1(L2)/MEGA                                                          
      IF(IP(L2).GE.2)GO TO 630                                                  
      IF(JT.GT.0)IP(L2)=0                                                       
      JT1=MOD(NY1(L2),MEGA)/1000                                                
      JT2= MOD(NY1(L2),1000)                                                    
      IF(J1.EQ.JT1.OR.J1.EQ.JT2)IP(L2)=IP(L2)+1                                 
      IF(JT1.EQ.J1) WRITE(6,620)JT2,NFLAG(JT2),L2,Y(L2)                         
      IF(JT2.EQ.J1) WRITE(6,620)JT1,NFLAG(JT1),L2,Y(L2)                         
  620 FORMAT(60X,I4,I5,5X,I6,7X,1PE12.4)                                        
  630 CONTINUE                                                                  
      DO 650 L2=LIMY,MAXNOY                                                     
      IF(IP(L2).LT.2)GO TO 660                                                  
  650 JLIMY= L2                                                                 
  660 LIMY=JLIMY                                                                
  690 CONTINUE                                                                  
  700 IF(LYMAX.LE.0)GO TO 750                                                   
      WRITE(6,800)                                                              
        IF(NYMIN.LE.0)WRITE(6,820)                                              
        IF(NYMIN.GT.0)WRITE(6,810)                                              
       WRITE(6,830)                                                             
  800 FORMAT(1H1,36X,'SUMMARY OF FLOW NODE CONNECTORS '/)                       
810     FORMAT(1H0/                                                             
     1'  CONNECTOR      INITIAL        DOWNSTREAM    NODE      INITIAL          
     2  UPSTREAM    NODE      INITIAL')                                         
820     FORMAT(1H0/                                                             
     1'  CONNECTOR       INITIAL       UPSTREAM     NODE        INITIAL         
     2  DOWNSTREAM   NODE  INITIAL')                                            
830     FORMAT(1H0/                                                             
     1'   NUMBER         VALUE           NODE        FLAG       VALUE           
     4  NODE        FLAG       VALUE '//)                                       
      DO 850 J1=1,LYMAX                                                         
      J2=NY2(J1)                                                                
      IF(J2.EQ.0) GO TO 850                                                     
      JT1= MOD(NY1(J2),MEGA)/1000                                               
      JT2= MOD(NY1(J2),1000)                                                    
      WRITE(6,860)J2,Y(J2),JT1,NFLAG(JT1),T(JT1),JT2,NFLAG(JT2),T(JT2)          
  850 CONTINUE                                                                  
  860 FORMAT(1H ,I8,4X,1PE12.4,2X,2I11,6X,G12.4,I4,I11,6X,G12.4)                
  750 JSTART = 0                                                                
      RETURN                                                                    
      END                                                                       
C     SUBROUTINE  **TIMER**                                                     
      SUBROUTINE TIMER                                                          
      INTEGER NY1(3000),NY2(1000),NZA(1000),INXP(1000),INWP(1000)               
      INTEGER NXP(1000),KEY(1000),NFLAG(1000),IZ(1000),IPR(500),IPH(100)        
     X,NTHEAD(100),NTTAIL(100)                                                  
      DIMENSION  ZZ(16000),DYTCQ(11030),Y(3000),T(1000)                         
     X,C(1000)        ,Q(1000)        ,GAIN(1000) ,A(1000),B(1000)              
     X,D(1000),SUMY(1000),SUMYT(1000), TOLD(1000),L(1000)                       
     X,PHT1(100),PHH1(100),PHT2(100),PHH2(100),QCH1(100)                        
     X,QCH2(100),CON(20),PRTIM(10),DELTIM(10),TITLE(17,10),TABL(1000)           
     X,TPROG(2),AA(7)                                                           
       COMMON/I2BYT/NY1,NY2,NZA,NXP,INXP,INWP,KEY,IPR,IPH                       
     X,NTHEAD,NTTAIL,NFLAG,IZ                                                   
      EQUIVALENCE (ZZ(001),TIME),         (ZZ(002),Y(1),DYTCQ(1))               
     X,          (ZZ(3002),T(1)),                 (ZZ(4002),C(1))               
     X,          (ZZ(5002),Q(1)),                 (ZZ(6002),D(1) )              
     X,          (ZZ(7002),A(1) ),              (ZZ(8002), B(1) )               
     X,          (ZZ(9002),L(1) ),                (ZZ(10002),GAIN(1) )          
     X,          (ZZ(11002),CON(1)),              (ZZ(11010),CON9,SIGMA)        
     X,          (ZZ(11011),CON10,TZERO), (ZZ(11032),PRTIM(1))                  
     X,          (ZZ(11042),DELTIM(1)),           (ZZ(11053),IREAD)             
      EQUIVALENCE(ZZ(11060),MIN1T),               (ZZ(11061),MAX1T)             
     X,          (ZZ(11062),N2)                                                 
      EQUIVALENCE(ZZ(11069), LFLAG),              (ZZ(11070), IERFL)            
     X,          (ZZ(11071),PRTMX),               (ZZ(11072),DTHETA)            
     X,          (ZZ(11073),RCMIN),               (ZZ(11074), JSAVE)            
     X,          (ZZ(11075),DELTAT),              (ZZ(11076),DTMAX )            
     X,          (ZZ(11077), IMIN),               (ZZ(11078), IMAX )            
     X,          (ZZ(11079), IXMIN),              (ZZ(11080), IXMAX)            
     X,          (ZZ(11081),IPHMAX),              (ZZ(11082),IPRMAX)            
     X,          (ZZ(11083),IPFLAG),              (ZZ(11084),ISFLAG)            
     X,          (ZZ(11085),MAXNOT),              (ZZ(11086),MAXNOY)            
     X,          (ZZ(11087),NARITH),              (ZZ(11088),NCRIT )            
      EQUIVALENCE (ZZ(11089),NTIME),              (ZZ(11090), NT )              
     X,          (ZZ(11091),NTLAST),              (ZZ(11092),NWW)               
     X,          (ZZ(11093),NXX),                 (ZZ(11094),NZZ)               
     X,          (ZZ(11095), W ),                 (ZZ(11096), X )               
     X,          (ZZ(11097), Z ) ,                (ZZ(11098), NTAB)             
     X,          (ZZ(11099),OTIME2),              (ZZ(11100),ODTIME)            
     X,          (ZZ(11101),DTIME),               (ZZ(11102),PHT1(1))           
     X,          (ZZ(11202),PHH1(1)),             (ZZ(11302),PHT2(1))           
     X,          (ZZ(11402),PHH2(1)),             (ZZ(11502),QCH1(1))           
     X,          (ZZ(11602),QCH2(1))                                            
      EQUIVALENCE(ZZ(12002),SUMY(1)),             (ZZ(13002),TOLD(1))           
     X,          (ZZ(14002),SUMYT(1)),            (ZZ(15002),TABL(1))           
      COMMON  ZZ                                                                
      IF(IERFL.EQ.2) GO TO 50                                                   
        RC=1.E36                                                                
      IF(TIME.EQ.PRTIM(1))POTIME=PRTIM(1)                                       
      IF(N2.EQ.4) GO TO 100                                                     
      IF(DELTAT.LE.0.)DELTAT=DELTIM(1)                                          
      MINDO= MIN1T                                                              
      MAXDO= MAX1T                                                              
      ITEST= 10                                                                 
      IF(CON(14).GT.0.)ITEST= CON(14)                                           
      RCMIN= 1.0E36                                                             
      TTEST=0.1                                                                 
      IF(CON(8).GT.0.0)TTEST=CON(8)                                             
      ICOUNT= 0                                                                 
      DO 10  J= MIN1T,MAX1T                                                     
      IF(SUMY(J).LE.0.) GO TO 10                                                
        IF(NFLAG(J).LT.1)GO TO 10                                               
      IF(C(J).GT.0.) RC= C(J)/SUMY(J)                                           
    7 IF(RC.GE.RCMIN) GO TO 10                                                  
      NCRIT = J                                                                 
      RCMIN= AMIN1(RC,RCMIN)                                                    
   10 CONTINUE                                                                  
      DELTAT= RCMIN*CON(3)                                                      
      GO TO 30                                                                  
   15 CONTINUE                                                                  
      IF(IPFLAG.LE.0.AND.ICOUNT.LE.ITEST) GO TO 30                              
      DELTAT = 2.*AMAX1(RCMIN,DELTAT)                                           
      ICOUNT= 0                                                                 
      DTT= ODTIME+DELTAT                                                        
      RCMIN= 1.0E36                                                             
      TR= 1.+ ODTIME/DELTAT                                                     
      DO 20 J = MINDO,MAXDO                                                     
        IF(NFLAG(J).LT.1)GO TO 20                                               
      IF(SUMY(J).LE.0.) GO TO 20                                                
      IF(C(J).LE.0.) GO TO 16                                                   
      RC = C(J)/SUMY(J)                                                         
      IF(RC.LT.RCMIN)NCRIT= J                                                   
   16 CONTINUE                                                                  
      RCMIN= AMIN1(RC,RCMIN)                                                    
   17 CONTINUE                                                                  
      DENOM=TR*C(J)+ODTIME*SUMY(J)                                              
      TNEW=TR*C(J)*T(J)+DTT*SUMYT(J)-DELTAT*SUMY(J)*TOLD(J)                     
   18 TNEW = TNEW/DENOM                                                         
      DENOM=C(J)+ODTIME*SUMY(J)                                                 
      TDUM=(DTT*SUMYT(J)-TOLD(J)*(DELTAT*SUMY(J)-C(J)))/DENOM                   
      ERTMP=ABS(TDUM-TNEW)*(1.+DELTAT/RC)                                       
      ERTMP=ERTMP*(1.+DELTAT/RC)/TR                                             
      IF(ERTMP.LE.TTEST) GO TO 20                                               
      EMAX=0.5*RCMIN*CON(3)                                                     
      IF(DELTAT.LT.EMAX)GO TO 20                                                
      DELTAT= DELTAT/2.                                                         
      DTT= ODTIME + DELTAT                                                      
      TR= 1.+ ODTIME/DELTAT                                                     
      GO TO 17                                                                  
   20 CONTINUE                                                                  
   25 CONTINUE                                                                  
   30 ICOUNT = ICOUNT + 1                                                       
      COMP = POTIME -TIME                                                       
      IF(COMP.LE.0.)COMP= DELTIM(NTIME)                                         
      IC= COMP/DELTAT                                                           
      IF(IC.LT.1) IC=1                                                          
      DIC= IC                                                                   
      DELTAT = COMP/DIC                                                         
      IF(N2.EQ.2.AND.DELTAT.GT.RCMIN)DELTAT= COMP/(DIC+1)                       
      DELTAT= AMIN1(CON(1),DELTAT)                                              
   40 IF(DELTAT.LT.CON(2)) IERFL=2                                              
      IF(IERFL.LT.2) GO TO 100                                                  
      WRITE(6,41)DELTAT,CON(2)                                                  
   41 FORMAT(1H ,10X,'DELTAT=',E12.4,' AND IS LESS THAN THE MINIMUM VALUE       
     1 ALLOWED (',E12.4,'),RUN TERMINATED')                                     
   50 RETURN                                                                    
  100 IF(IPFLAG.NE.1) GO TO 140                                                 
      IPFLAG= 0                                                                 
      IF(TIME.LT.POTIME) GO TO 140                                              
      POTIME= POTIME + DELTIM(NTIME)                                            
      CHECK= PRTIM(NTIME+1)-0.001*DELTIM(NTIME)                                 
      IF(POTIME.LT.CHECK) GO TO 140                                             
      NTIME=NTIME+1                                                             
                                                                                
      POTIME= PRTIM(NTIME)                                                      
      IF(DELTIM(NTIME).LE.0.) DELTIM(NTIME)= DELTIM(NTIME-1)                    
  140 TIMEX = TIME + DELTAT                                                     
      TEST= POTIME - 0.10*DELTAT                                                
      IF(TIMEX.LT.TEST) GO TO 150                                               
      DELTAT=  POTIME - TIME                                                    
      TIMEX= POTIME                                                             
      IPFLAG=1                                                                  
  150 DTHETA = DELTAT                                                           
      TIME= TIMEX                                                               
      IF(TIME.GE.PRTMX)ISFLAG=1                                                 
      RETURN                                                                    
      END                                                                       
        SUBROUTINE DUMMY                                                        
       COMMON/LAB2/PDATE(3),CPU                                                 
       ENTRY ERF                                                                
       ENTRY ERFC                                                               
       ENTRY TIMEV(TIME2)                                                       
       RETURN                                                                   
       ENTRY COUNTV                                                             
       RETURN                                                                   
       ENTRY SUB98                                                              
       ENTRY USER1                                                              
       ENTRY USER2                                                              
       ENTRY CDATEV(DATE)
       ENTRY PLOTPT                                                             
       RETURN                                                                   
       END                                                                      
