      SUBROUTINE KIP1 (NSP,NCEQ,NHCP,HCP,SK,NAE,AEA,                     
     &                 KINEV,TEMP,RK,GAS1,GAS2,KFAIL)                    
C*********************************************************************** 
C                                                                      * 
C DATE OF LATEST CHANGE: JANUARY 3, '86                                * 
C                                                                      * 
C*********************************************************************** 
C                                                                      * 
C  SUBROUTINE KIP1 CALCULATES THE REACTION CONSTANTS                   * 
C                                                                      * 
C*********************************************************************** 
C                                                                      * 
C  THE REACTION CONSTANT IS COMPUTED ACCORDING TO:                     * 
C                                                                      * 
C  A. (IF 1 PARAMETER FOR ONE CHEM.EQUATION IS GIVEN)                  * 
C     REACTION CONSTANT IS SET TO THE VALUE OF THE GIVEN PARAMETER     * 
C                                                                      * 
C  B. (IF 2 PARAMETERS FOR ONE CHEM.EQUATION ARE GIVEN)                * 
C     ARRHENIUS LAW:                                                   * 
C     RK := A * EXP(-E/(GAS2*TEMP))                                    * 
C     WITH: PARAMETERS ASSUMED TO BE:  LN(A) ,  E                      * 
C                                                                      * 
C  C. (IF 3 PARAMETERS FOR ONE CHEM.EQUATION ARE GIVEN)                * 
C     MODIFIED ARRHENIUS LAW:                                          * 
C     RK := A * EXP(-E/(GAS2*TEMP)) * TEMP**ALFA                       * 
C     WITH: PARAMETERS ASSUMED TO BE:  LN(A) ,  E , ALFA               * 
C                                                                      * 
C  D. (IF NO PARAMETERS FOR ONE CHEM.EQUATION ARE GIVEN)               * 
C     THIS IS ONLY POSSIBLE FOR REVERSE REACTIONS, IF KINETIC PARA-    * 
C     METERS FOR ASSOCIATED FORWARD REACTIONS ARE GIVEN. THE REACTION  * 
C     CONSTANTS ARE THEN COMPUTED BY MEANS OF MASS-ACTION-LAW,         * 
C     INVOLVING ENTHALPIES AND ENTROPIES OF SPECIES PARTICIPATING IN   * 
C     THESE REACTIONS.                                                 * 
C                                                                      * 
C*********************************************************************** 
C                                                                      * 
C  INPUT PARAMETERS:                                                   * 
C-------------------                                                   * 
C                                                                      * 
C  NSP             NUMBER OF CHEMICAL SPECIES                          * 
C  NCEQ            NUMBER OF CHEMICAL EQUATIONS                        * 
C  NAE             = 3*NCEQ                                            * 
C  TEMP            TEMPERATURE                                         * 
C  GAS1            UNIVERSAL GAS CONSTANT (FOR UNIVERSAL GAS LAW)      * 
C  GAS2            UNIVERSAL GAS CONSTANT (FOR ARRHENIUS LAW)          * 
C  AEA(1),..,AEA(NCEQ):                                                * 
C                  FIRST (KINETIC) PARAMETER OF EQUATION 1,...,NCEQ    * 
C  AEA(NCEQ+1),..,AEA(2*NCEQ):                                         * 
C                  SECOND (KINETIC) PARAMETER OF EQUATION 1,...,NCEQ   * 
C  AEA(2*NCEQ+1),..,AEA(3*NCEQ):                                       * 
C                  THIRD (KINETIC) PARAMETER OF EQUATION 1,...,NCEQ    * 
C  HCP(NHCP)       MOLAR ENTHALPIES AND HEAT CAPACITIES OF SPECIES     * 
C                  (OUTPUT OF SUBROUTINE THERMO)                       * 
C  SK(NSP)         MOLAR ENTROPIES OF SPECIES                          * 
C                  (OUTPUT OF SUBROUTINE ENTROP)                       * 
C  KINEV(NCEQ)     POINTER FOR COMPUTATION OF RATE CONSTANTS           * 
C                  = 1 : CONSTANT TYPE                                 * 
C                  = 2 : ARRHENIUS TYPE                                * 
C                  = 3 : MODIFIED ARRHENIUS TYPE                       * 
C                  = 4 : ZERO ACTIVATION ENERGY                        * 
C                  = 5 : REVERSE RATE CONSTANT DUE TO FREE ENTHALPY    * 
C                  = 6 : ARRHENIUS TYPE OF REVERSE REACTION, PARA-     * 
C                        METERS DETERMINED IN SIMULA (KREV=2)          * 
C                  = 7 : MODIFIED ARRHENIUS TYPE OF REVERSE REACTION,  * 
C                        PARAMETERS DETERMINED IN SIMULA (KREV=2)      * 
C                                                                      * 
C  COMMON BLOCKS:                                                      * 
C----------------                                                      * 
C                                                                      * 
C  COMMON / LARK1 / LR(NLR)                                            * 
C  COMMON / LARK2 / PLR(2  *NCEQ)                                      * 
C         POINTER TO THE RIGHT HAND SIDE OF DIFFERENTIAL EQUATIONS     * 
C         (SEE DESCRIPTION IN DRIVER)                                  * 
C                                                                      * 
C                                                                      * 
C  OUTPUT PARAMETERS:                                                  * 
C--------------------                                                  * 
C                                                                      * 
C  RK(1),...,RK(NCEQ):  COMPUTED RATE CONSTANTS                        * 
C  KFAIL                INTERNAL FLAG                                  * 
C                       = 0 SUCCESSFUL CALL OF KIP1                    * 
C                       = 1 FAIL: X TOO LARGE FOR EXP(X)               * 
C                       = 2 FAIL: NEGATIVE OR ZERO TEMPERATURE         * 
C                                                                      * 
C*********************************************************************** 
C                                                                        
CI4       INTEGER KINEV(NCEQ),LR,PLR                                     
      INTEGER*2 KINEV(NCEQ),LR,PLR                                       
      DOUBLE PRECISION AEA(NAE)                                          
CSP       REAL AEA(NAE)                                                  
      DOUBLE PRECISION RK(NCEQ),TEMP,ZERO,EXPO,GAS1,GAS2,RCORR,RLSE      
CSP       REAL RK(NCEQ),TEMP,ZERO,EXPO,GAS1,GAS2,RCORR,RLSE              
      DOUBLE PRECISION EXPTST,TLOG,RT,EXPOF,EXPOFH,EXPOR,EXPORH,RKFH     
CSP       REAL EXPTST,TLOG,RT,EXPOF,EXPOFH,EXPOR,EXPORH,RKFH             
      DOUBLE PRECISION DUM1,DUM2,DELTAS,DELTAH,SK(NSP),HCP(NHCP)         
CSP       REAL DUM1,DUM2,DELTAS,DELTAH,SK(NSP),HCP(NHCP)                 
      DOUBLE PRECISION EXPTT                                             
CSP       REAL EXPTT                                                     
      INTEGER I,I1,I2,II,IL1,IL2,IR2,J,KIN,LRI,NII                       
C                                                                        
      COMMON / MACHIN / DUM1,DUM2,EXPTT                                  
C                                                                        
      COMMON / LARK1  / LR(5000)                                         
      COMMON / LARK2  / PLR(2000)                                        
      COMMON / COUNT  / NSTEP,NFCN,NANFA,NTFAC,NFAC,NSOL,NTHERM,NKIP
      SAVE /MACHIN/, /LARK1/, /LARK2/, /COUNT/
C                                                                        
      DATA  ZERO/0.D0/ , RLSE/3.D1/                                      
CSP       DATA  ZERO/0.E0/ , RLSE/3.E1/                                  
C                                                                        
      KFAIL=0                                                            
      NKIP=NKIP+1                                                        
      EXPTST=EXPTT-RLSE                                                  
      IF(TEMP.LE.ZERO) GOTO 99                                           
      RT=GAS2*TEMP                                                       
      RCORR=GAS2/GAS1                                                    
      TLOG=DLOG(TEMP)                                                    
CSP       TLOG=ALOG(TEMP)                                                
C                                                                        
      DO 10 I=1,NCEQ                                                     
      I1=NCEQ+I                                                          
      I2=NCEQ+I1                                                         
      KIN=KINEV(I)                                                       
C                                                                        
      GOTO (1,2,3,4,5,6,7) , KIN                                         
6     GOTO 2                                                             
7     GOTO 3                                                             
C                                                                        
C  CONSTANT TYP                                                          
1     RK(I)=AEA(I)                                                       
      GOTO 10                                                            
C  ARRHENIUS TYP                                                         
2     EXPOF=AEA(I)-AEA(I1)/RT                                            
      IF(EXPOF.GT.EXPTST) GOTO 99                                        
22    RK(I)=DEXP(EXPOF)                                                  
CSP       RK(I)=EXP(EXPOF)                                               
      GOTO 10                                                            
C  MODIFIED ARRHENIUS TYP                                                
3     EXPOF=AEA(I) + AEA(I2)*TLOG - AEA(I1)/RT                           
      IF(EXPOF.GT.EXPTST) GOTO 99                                        
33    RK(I)=DEXP(EXPOF)                                                  
CSP       RK(I)=EXP(EXPOF)                                               
      GOTO 10                                                            
C  NO EXPONENTIAL                                                        
4     EXPOF=AEA(I) + AEA(I2)*TLOG                                        
      IF(EXPOF.GT.EXPTST) GOTO 99                                        
44    RK(I)=DEXP(EXPOF)                                                  
CSP       RK(I)=EXP(EXPOF)                                               
      GOTO 10                                                            
C  REACTION CONSTANT FOR REVERSE REACTIONS                               
5     II=I-1                                                             
      DELTAS=ZERO                                                        
      DELTAH=ZERO                                                        
      NII=NCEQ+II                                                        
      IL1=PLR(NII-1)+1                                                   
      IF(II.EQ.1) IL1=1                                                  
      IL2=PLR(II)                                                        
      IR2=PLR(NII)                                                       
      DO 70 J=IL1,IR2                                                    
      LRI=LR(J)                                                          
      IF(J.GT.IL2) GOTO 60                                               
      DELTAS=DELTAS - SK(LRI)                                            
      DELTAH=DELTAH - HCP(LRI)                                           
      GOTO 70                                                            
60    DELTAS=DELTAS + SK(LRI)                                            
      DELTAH=DELTAH + HCP(LRI)                                           
70    CONTINUE                                                           
C  THIS MULTIPLICATION BY 'RCORR' IS NECESSARY FOR GAS1.NE.GAS2          
      DELTAH=DELTAH*RCORR                                                
      DELTAS=DELTAS*RCORR                                                
      AEA(I1)=AEA(I1-1) - DELTAH                                         
      EXPO=AEA(I) + AEA(I2)*TLOG - (AEA(I1)+TEMP*DELTAS)/RT              
      IF(EXPO.GT.EXPTST) GOTO 99                                         
      RK(I)=DEXP(EXPO)                                                   
CSP       RK(I)=EXP(EXPO)                                                
10    CONTINUE                                                           
C                                                                        
C  NORMAL RETURN , NO ERROR                                              
C  ------------------------                                              
      RETURN                                                             
C                                                                        
C  OVERFLOW RETURN                                                       
C  ---------------                                                       
99    KFAIL=1                                                            
      IF(TEMP.LE.ZERO) KFAIL=2                                           
      RETURN                                                             
C                                                                        
C----------------------------------------------------------------------- 
C                  END OF SUBROUTINE KIP1                                
C----------------------------------------------------------------------- 
C                                                                        
      END 
      SUBROUTINE SIMULA (RW,NRW,RWC,NRWC,IW,NIW,IWP,NIWP,CW,NCW,KFLAG)   
C*********************************************************************** 
C                                                                      * 
C DATE OF LATEST CHANGE: JANUARY 3, '86                                * 
C                                                                      * 
C*********************************************************************** 
C                                                                      * 
C  INITIATION ROUTINE FOR SIMULATION                                   * 
C  READS ALL SIMULATION INPUT FROM DATABASE 'INVAL' , TRANSFERS PART   * 
C  OF IT TO DATABASE 'OUT' (SIMULATION OUTPUT) AND STARTS INTEGRATOR.  * 
C  SPECIAL DEVICES FOR COMPUTATION OF REVERSE REACTION CONSTANTS.      * 
C                                                                      * 
C*********************************************************************** 
C                                                                      * 
C  INPUT PARAMETERS                                                    * 
C------------------                                                    * 
C                                                                      * 
C   RW(NRW)    REAL WORK SPACE                                         * 
C   RWC(NRWC)  REAL WORK SPACE                                         * 
C   NRW,NRWC   DIMENSIONS OF REAL WORK SPACES                          * 
C   IW(NIW)    INTEGER WORK SPACE                                      * 
C   IWP(NIWP)  INTEGER WORK SPACE                                      * 
C   NIW,NIWP   DIMENSIONS OF INTEGER WORK SPACES                       * 
C   CW         CHARACTER WORK SPACE                                    * 
C   NCW        DIMENSION OF CHARACTER WORK SPACE                       * 
C                                                                      * 
C   KFLAG      PRINT PARAMETER FOR INTEGRATOR OUTPUT                   * 
C              .EQ.0: ERROR MESSAGES                                   * 
C              .EQ.1: INTEGRATION MONITOR                              * 
C              .EQ.2: ENHANCED INTEGRATION MONITOR                     * 
C                                                                      * 
C                                                                      * 
C  OUTPUT PARAMETERS                                                   * 
C-------------------                                                   * 
C    KFLAG     ERROR FLAG                                              * 
C              .EQ.0: SUCCESSFULL SIMULATION                           * 
C              .EQ.-1: SIMULATION FAILED                               * 
C                                                                      * 
C  COMMON BLOCKS                                                       * 
C---------------                                                       * 
C                                                                      * 
C     INTEGER          ITIN,ITOUT,CIN,HDS,POUT,SIN,DOUT,THDAT *          
C     COMMON / LUNIT / ITIN,ITOUT,CIN,HDS,POUT,SIN,DOUT,THDAT *          
C                                                                        
C  LOGICAL UNITS ( SEE DESCRIPTION OF LARKIN DRIVING ROUTINE )         * 
C                                                                      * 
C                                                                      * 
C   COMMON / LARK1 / LR(NLR)                                           * 
C   COMMON / LARK2 / PLR(2*NCEQ)                                       * 
C  INTEGER ARRAYS CONTAINING POINTERS TO THE ODE-SYSTEM AND SPARSE     * 
C                                            JACOBIAN                  * 
C   COMMON / KINPAR / RK(NCEQ)                                         * 
C  REAL ARRAY FOR COMPUTED KINETIC PARAMETERS                          * 
C                                                                      * 
C   COMMON / RECOND  / TEMP,DENSY,PRESS,COLUME,XMASS                   * 
C  REACTION CONDITIONS OF THE TEST GAS MIXTURE. (TRANSIENT PARAMETERS) * 
C                                                                      * 
C   COMMON / MACHIN / DUM1,SMALL,DUM2                                  * 
C  SEE DESCRIPTION IN DRIVING ROUTINE                                  * 
C                                                                      * 
C   COMMON / RATE1 / R(NCEQ)                                           * 
C  REAL ARRAY FOR COMPUTED RATES (USED IN SUBROUTINE META FOR OUTPUT   * 
C  PURPSES ONLY)                                                       * 
C                                                                      * 
C   COMMON / DELHR / DHR(NCEQ)                                         * 
C  REAL ARRAY FOR ENTHALPY PER REACTION (USED IN SUBROUTINE META FOR   * 
C  OUTPUT PURPOSES ONLY)                                               * 
C                                                                      * 
C  REQUESTED DIMENSIONS FOR SIMULATION                                 * 
C       RW  : REAL WORK SPACE OF DIMENSION NRW                         * 
C             NRW   = 18*NSP + LDIM + NLU                              * 
C       RWC : REAL WORK SPACE OF DIMENSION NRWC                        * 
C             NRWC  = 14*NSP + 3*NCEQ                                  * 
C                                                                      * 
C       IW : INTEGER WORK SPACE OF DIMENSION NIW                       * 
C            NIW  = 13*NSP + NLU                                       * 
C       IWP: INTEGER WORK SPACE OF DIMENSION NIWP                      * 
C            NIWP = 3*NSP + 2*NCEQ + LDIM                              * 
C                                                                      * 
C            NSP  =  NUMBER OF SPECIES                                 * 
C            LDIM =  NUMBER OF NONZEROS IN JACOBIAN                    * 
C            NLU  =  ESTIMATED MAXIMUM NUMBER OF NONZEROS IN           * 
C                    LU-DECOMPOSITION (NLU = LDIM + #FILL-IN ENTRIES)  * 
C            NCEQ =  NUMBER CHEMICAL EQUATIONS                         * 
C            NLR  =  NUMBER OF ENTRIES IN ARRAY LR                     * 
C       CW : CHARACTER WORK SPACE OF DIMENSION NCW                     * 
C            NCW  = MAX (NAMEL*NSP , 360)                              * 
C            NAMEL= MAXIMUM LENGTH OF SPECIES NAMES (.LE.10 CHARS)     * 
C                                                                      * 
C       TP(50)      : ARRAY FOR PRESCRIBED OUTPUT POINTS (.LE.50)      * 
C       ATWS(10)    : ARRAY FOR ELEMENT ATOMIC WEIGHT (.LE.10)         * 
C       ELEM(10,5)  : ARRAY FOR ELEMENT COMPOSITION OF SPECIES         * 
C                    .LE.10 ELEMENTS WITH NAME LENGTH OF .LE.5 DIGITS  * 
C       LENTB(5)    : LENTB(K) = 'LENGTH' OF K-TH THIRD BODY           * 
C       NAMETB(25)  : NAMES OF THIRD BODIES                            * 
C                                                                      * 
C                                                                      * 
C  SUBROUTINES USED (PROVIDED BY LARKIN)                               * 
C------------------                                                    * 
C       META1                     INTEGRATION ROUTINE                  * 
C       THERMO , ENTROP , KIP1    (SEE DESCRIPTION IN META1)           * 
C                                                                      * 
C*********************************************************************** 
C                                                                        
      INTEGER LENTB(5)                                                   
CI4       INTEGER IW(NIW),IWP(NIWP),LR,PLR                               
      INTEGER*2 IW(NIW),IWP(NIWP),LR,PLR                                 
      CHARACTER*1 ELEM(10,5),NAMETB(25),CW(NCW)                          
C                                                                        
      INTEGER I,I1,I2,ICO,IDISC,IEL,IH,IH1,IH2,IHEAD,IKIN,INDEX,         
     & IPREV,IPRINT,ITB,ITPM,IUNIT,J,JFLAG,JH,JJ1,JJ2,K,K1,K2,K3,        
     & KD1,KD2,KFAIL,KK1,KK2,KREV,L,L1,LDIM,LH1,LH2,LJ,MCH1,MCH2,        
     & MCH3,MCH4,MH1,MH2,MH3,MODEL,MSCW,MSIW,MSIWP,MSRW,MSRWC,           
     & N,NAE,NCEQ,NCEQ2,NCEQ3,NCH,NCO,NCO12,NCO2,NCTB,NEN,NFREE,         
     & NH1,NH2,NHCP,NIWH,NIWPH,NLR,NP1,NRWCH,NRWH,NSK,NSP,NTB            
      REAL R                                                             
      DOUBLE PRECISION RWC(NRWC),RW(NRW),TP(50),ATWS(10)                 
CSP       REAL RWC(NRWC),RW(NRW),TP(50),ATWS(10)                         
      DOUBLE PRECISION DENSY,DIV,DIV1,DUM1,DUM2,EPS,GAS1,GAS2,HMAX       
CSP       REAL DENSY,DIV,DIV1,DUM1,DUM2,EPS,GAS1,GAS2,HMAX               
      DOUBLE PRECISION ONE,PRESS,RGAS,RK,RKTEST,SMALL,SMALLH             
CSP       REAL ONE,PRESS,RGAS,RK,RKTEST,SMALL,SMALLH                     
      DOUBLE PRECISION TEMP,THIGH                                        
CSP       REAL TEMP,THIGH                                                
      DOUBLE PRECISION TLOW,U,VOLUME,XLTLOW,XLTT,XMASS,ZERO              
CSP       REAL TLOW,U,VOLUME,XLTLOW,XLTT,XMASS,ZERO                      
C                                                                        
      INTEGER          ITIN,ITOUT,MOUT,EOUT,CIN,HDS,POUT,SIN,DOUT,THDAT  
      COMMON / LUNIT / ITIN,ITOUT,MOUT,EOUT,CIN,HDS,POUT,SIN,DOUT,THDAT  
C                                                                        
C                                                                        
      COMMON / MACHIN / DUM1,SMALL,DUM2                                  
C                                                                        
      COMMON / LARK1  / LR(5000)                                         
C                       LR(NLR)                                          
      COMMON / LARK2  / PLR(2000)                                        
C                       PLR(2*NCEQ)                                      
      COMMON / KINPAR / RK(1000)                                         
C                       RK(NCEQ)                                         
      COMMON / RATE1  / R(1000)                                          
C                       R(NCEQ)                                          
      COMMON / COUNT  / NSTEP,NFCN,NANFA,NTFAC,NFAC,NSOL,NTHERM,NKIP     
C                                                                        
      COMMON / RECOND / TEMP,DENSY,PRESS,VOLUME,XMASS 
      SAVE /LUNIT/, /MACHIN/, /LARK1/, /LARK2/, /KINPAR/, /RATE1/,
     & /COUNT/
C                                                                        
      DATA ZERO/0.D0/ ,ONE/1.D0/ , RKTEST/-1.D+31/                       
CSP        DATA ZERO/0.E0/ , ONE/1.E0/ , RKTEST/-1.E+31/                 
C                                                                        
      SMALLH=SMALL*SMALL                                                 
C                                                                        
C REWIND IN-/OUTPUT FILES                                                                                                       
      REWIND SIN 
      REWIND DOUT                                                   
C                                                                        
C  1. CHECK FOR MINIMUM STORAGE                                          
C  ----------------------------                                          
      MSCW=360-NCW                                                       
      IF(360.GT.NCW) GOTO 9809                                           
C                                                                        
C----------------------------------------------------------------------- 
C  SET STANDARD VALUE FOR COMPUTATION OF REVERSE REACTION CONSTANTS      
C  ( SEE DESCRIPTION IN DRIVING ROUTINE)                                 
C  KREV = 1 : RECOMPUTATION FOR EVERY CHANGE IN TEMPERATURE WITH         
C             RESPECT TO ENTHALPIES AND ENTROPIES                        
C       = 2 : COMPUTATION ONLY  O N C E  FOR A GIVEN TEMPERATURE RANGE   
C             DEFINE LOW TEMP. 'TLOW' AND HIGH TEMP. 'THIGH' BELOW       
C----------------------------------------------------------------------- 
      KREV=IW(1)                                                         
      TLOW=RW(1)                                                         
      THIGH=RW(2)                                                        
C                                                                        
C----------------------------------------------------------------------- 
C  READ SIMULATION INPUT                                                 
C----------------------------------------------------------------------- 
C                                                                        
C  READ HEAD                                                             
      READ(SIN,9000) IHEAD                                               
      IHEAD=IHEAD*72                                                     
      IF(IHEAD.LE.0) GOTO 5                                              
      READ(SIN,9005) (CW(L),L=1,IHEAD)                                   
C                                                                        
C  WRITE HEAD FOR INTERACTIVE SESSION                                    
5     IF(KFLAG.GE.1) WRITE(ITOUT,9003)                                   
      IF(KFLAG.GE.1.AND.IHEAD.GT.0) WRITE(ITOUT,8001) (CW(L),L=1,IHEAD)  
      IF(KFLAG.GE.1.AND.KREV.EQ.2) WRITE(ITOUT,9002) TLOW,THIGH          
      IF(KFLAG.GE.1.AND.KREV.EQ.1) WRITE(ITOUT,9004)                     
C                                                                        
C  WRITE HEAD IN OUTPUT DATABASE 'OUT'                                   
C----------------------------------------------------------------------- 
C                                                                        
      IH=IHEAD/72                                                        
      WRITE(DOUT,8000) IH                                                
      IF(IHEAD.GT.0) WRITE(DOUT,8001) (CW(L),L=1,IHEAD)                  
      WRITE(DOUT,8020)                                                   
C                                                                        
C                                                                        
C  READ MODEL PARAMETER                                                  
      READ(SIN,9000) MODEL                                               
      IF(KFLAG.LT.1) GOTO 9                                              
      IF(MODEL.EQ.1) WRITE(ITOUT,9100)                                   
      IF(MODEL.EQ.2) WRITE(ITOUT,9150)                                   
      IF(MODEL.EQ.3) WRITE(ITOUT,9200)                                   
      IF(MODEL.EQ.4) WRITE(ITOUT,9250)                                   
      IF(MODEL.EQ.5) WRITE(ITOUT,9300)                                   
      IF(MODEL.EQ.6) WRITE(ITOUT,9350)                                   
C                                                                        
C  READ ELEMENT LIST                                                     
9     READ(SIN,9000) IEL                                                 
      IF(IEL.LE.0) GOTO 11                                               
      DO 10 I=1,IEL                                                      
10    READ(SIN,9007) (ELEM(I,L),L=1,5),ATWS(I)                           
C                                                                        
C  READ DIMENSIONS                                                       
11    READ(SIN,9000) NSP                                                 
      N=NSP                                                              
      IF(MODEL.GE.3) N=NSP+1                                             
      IF(MODEL.GE.5) N=NSP+2                                             
C                                                                        
C  2. CHECK FOR MINIMUM STORAGE                                          
C  ----------------------------                                          
      IH=10*NSP                                                          
      MSIW=IH-NIW                                                        
      IF(IH.GE.NIW) GOTO 9801                                            
      MSCW=IH-NCW                                                        
      IF(IH.GE.NCW) GOTO 9809                                            
      IH=2*N                                                             
      MSRW=IH-NRW                                                        
      IF(IH.GE.NRW) GOTO 9803                                            
C                                                                        
      IH1=10*NSP                                                         
      MH1=NRW-NSP                                                        
      MH2=NIW-NSP*(5+IEL)                                                
      MH3=NIW-IEL*NSP                                                    
C                                                                        
C  READ ARRANGEMENT,NAMES,INITIAL VALUES AND MOLECULAR WEIGHTS           
C  AND (FOR IEL.GT.0) ELEMENT COMPOSITION OF SPECIES                     
      DO 20 I=1,NSP                                                      
      IF(IEL.GT.0) READ(SIN,9009) IH,(CW(10*(IH-1)+L),L=1,10),RW(IH),    
     &             RW(MH1+IH),(IW(MH3+IEL*(IH-1)+L),L=1,IEL)             
20    IF(IEL.LE.0) READ(SIN,9010) IH,(CW(10*(IH-1)+L),L=1,10),RW(IH)     
C                                                                        
C  READ TEMPERATURE ( = INITIAL VALUE FOR Y(NSP+1) FOR MODEL>1 )         
      READ(SIN,9015) TEMP                                                
      IF(MODEL.GT.1) RW(NSP+1)=TEMP                                      
C                                                                        
C  READ PRESSURE ( = INITIAL VALUE FOR Y(NSP+2) FOR MODEL=3 )            
      READ(SIN,9015) PRESS                                               
      IF(MODEL.EQ.3 .OR. MODEL.EQ.4) RW(NSP+2) = PRESS                   
C                                                                        
C  READ DENSITY ( = INITIAL VALUE FOR Y(NSP+2) FOR MODEL=4 )             
      READ(SIN,9015) DENSY                                               
      IF(MODEL.GE.5) RW(NSP+2) = DENSY                                   
C                                                                        
C  READ UNIVERSAL GAS CONSTANTS AND IDENTIFIER FOR UNIT SYSTEM           
      READ(SIN,9016) GAS1,GAS2,IUNIT                                     
C                                                                        
C  READ THIRD-BODY-COMPOSITIONS                                          
      READ(SIN,9000) NTB                                                 
      JH=0                                                               
      DO 68 I=1,5                                                        
68    LENTB(I)=0                                                         
      IF(NTB.EQ.0) GOTO 75                                               
      DO 70 I=1,NTB                                                      
      IH1=(I-1)*10 + 1                                                   
      IH2=I*10                                                           
      READ(SIN,9020) LENTB(I),(NAMETB(L),L=IH1,IH2)                      
      ITB=LENTB(I)                                                       
      IF(ITB.LT.0) GOTO 70                                               
      DO 71 J=1,ITB                                                      
      JH=JH+1                                                            
71    READ(SIN,9030) IW(JH),RW(N+JH)                                     
70    CONTINUE                                                           
75    CONTINUE                                                           
      NCTB=JH                                                            
      IF(NCTB.EQ.0) NCTB=1                                               
C                                                                        
C  READ KINETIC PARAMETERS                                               
      READ(SIN,9000) NCEQ                                                
      NCEQ2=NCEQ+NCEQ                                                    
      NCEQ3=NCEQ2+NCEQ                                                   
C  CHECK FOR MINIMUM STORAGE                                             
      MSRWC=NCEQ3-NRWC                                                   
      IF(NCEQ3.GT.NRWC) GOTO 9805                                        
      MSIWP=NCEQ2-NIWP                                                   
      IF(NCEQ2.GT.NIWP) GOTO 9807                                        
      DO 30 I=1,NCEQ                                                     
      I1=NCEQ+I                                                          
      I2=NCEQ+I1                                                         
      READ(SIN,9040) RWC(I),RWC(I1),RWC(I2)                              
      IF(RWC(I).LE.RKTEST) RWC(I)=ZERO                                   
      IF(RWC(I1).LE.RKTEST) RWC(I1)=ZERO                                 
      IF(RWC(I2).LE.RKTEST) RWC(I2)=ZERO                                 
30    CONTINUE                                                           
C                                                                        
C  READ IDENTIFIER FOR THIRD-BODIES (ICOLLI)                             
      READ(SIN,9070) (IWP(L),L=1,NCEQ)                                   
      NCH=NCTB+NCEQ                                                      
C                                                                        
C  READ IDENTIFIER FOR EVALUATION OF KINETIC PARAMETERS (KINEV)          
      READ(SIN,9070) (IWP(NCEQ+L),L=1,NCEQ)                              
C                                                                        
C READ COEFFICIENTS FOR ENTHALPIES (ENTH) AND ENTROPIES (SK0)            
      READ(SIN,9000) ICO                                                 
      NCO=ICO                                                            
      IF(NCO.EQ.0) NCO=1                                                 
      NCO12=12*NCO                                                       
      MCH1=NCEQ3+NCO12                                                   
      IF(ICO.EQ.0) GOTO 5544                                             
C  CHECK FOR MINIMUM STORAGE                                             
      IH=NCEQ3+14*NCO                                                    
      MSRWC=IH-NRWC                                                      
      IF(IH.GT.NRWC) GOTO 9805                                           
C  SET ENTHALPIES AND ENTROPIES TO ZERO                                  
      DO 33 I=1,NCO12                                                    
33    RWC(NCEQ3+I)=ZERO                                                  
      NCO2=NCO+NCO                                                       
      DO 34 I=1,NCO2                                                     
34    RWC(MCH1+I)=ZERO                                                   
C                                                                        
      DO 9521 I=1,ICO                                                    
      READ(SIN,9001) INDEX                                               
      IWP(NCEQ2+I)=INDEX                                                 
      KK2=NCEQ3 + I*12                                                   
      KK1=KK2-11                                                         
      LH2=MCH1 + I*2                                                     
      LH1=LH2-1                                                          
      READ(SIN,9060) (RWC(L),L=KK1,KK2),(RWC(LJ),LJ=LH1,LH2)             
9521  CONTINUE                                                           
C                                                                        
C  READ SYMBOLIC REPRESENTATION OF FCN AND JACOBIAN                      
5544  CONTINUE                                                           
      READ(SIN,9000) NLR                                                 
      READ(SIN,9080) (PLR(J),J=1,NCEQ)                                   
      JJ1=NCEQ+1                                                         
      JJ2=2*NCEQ                                                         
      READ(SIN,9080) (PLR(J),J=JJ1,JJ2)                                  
      READ(SIN,9080) (LR(J),J=1,NLR)                                     
C                                                                        
C  READ SPARSE PATTERN (JA,IA,IDA)                                       
      READ(SIN,9000) LDIM                                                
C                                                                        
C  3. CHECK FOR MINUMUM STORAGE                                          
C  ----------------------------                                          
      IH=13*N + 2*LDIM + NCTB                                            
      MSIW=IH-NIW                                                        
      IF(IH.GT.NIW) GOTO 9801                                            
      IH=IH+N                                                            
      IF(ICO.NE.0)   IH=IH+4*NSP                                         
      IF(MODEL.GE.5) IH=IH+2*NSP                                         
      MSRW=IH-NRW                                                        
      IF(IH.GT.NRW) GOTO 9803                                            
      IH=NCEQ2+NCO+LDIM+2*N+1                                            
      MSIWP=IH-NIWP                                                      
      IF(IH.GT.NIWP) GOTO 9807                                           
      MCH4=NCEQ2+NCO                                                     
      MSCW=360-NCW                                                       
      IF(360.GT.NCW) GOTO 9809                                           
C                                                                        
      READ(SIN,9080) (IWP(MCH4+J),J=1,LDIM)                              
      NP1=N+1                                                            
      MCH2=MCH4+LDIM                                                     
      MCH3=MCH2+N+1                                                      
      READ(SIN,9080) (IWP(MCH2+J),J=1,NP1)                               
      READ(SIN,9080) (IWP(MCH3+J),J=1,N)                                 
C                                                                        
C  READ IDENTIFIER FOR DISCRETISATION METHOD                             
      READ(SIN,9000) IDISC                                               
      IF(IDISC.GT.1) IDISC=1                                             
      IF(IDISC.LT.0) IDISC=0                                             
C                                                                        
C  READ RELATIVE PRECISION EPS                                           
      READ(SIN,9015) EPS                                                 
C                                                                        
C  READ OUTPUT POINTS                                                    
      READ(SIN,9000) ITPM                                                
      READ(SIN,9050) (TP(I),I=1,ITPM) 
C                                                                        
C  READ PRINT PARAMETER                                                  
      READ(SIN,9000) IPRINT                                              
C                                                                        
C  REACTION IS ASSUMED TO HAPPEN AT AN INITIAL VOLUME OF 1 (IN APPRO-    
C  PRIATE UNITS), SO XMASS(TSTART)=XMASS(TEND)=DENSY*VOLUME=DENSY        
      VOLUME=1.D0                                                        
CSP       VOLUME=1.                                                      
      XMASS=DENSY                                                        
C                                                                        
C----------------------------------------------------------------------- 
C  COMPUTE KINETIC PARAMETERS  O N C E  FOR MODEL=1                      
C----------------------------------------------------------------------- 
      NTHERM=0                                                           
      KD1=2*N                                                            
      KD2=KD1+NSP                                                        
      NHCP=3*NSP                                                         
      NEN=12*NCO                                                         
      NSK=2*NCO                                                          
      NAE=3*NCEQ                                                         
      IF((MODEL.GT.1.OR.ICO.LE.0).AND.IPRINT.GT.0) GOTO 69               
      DO 65 I=1,NCEQ                                                     
65    IF(IWP(NCEQ+I).EQ.5) GOTO 61                                       
      GOTO 69                                                            
61    CALL THERMO (NSP,NHCP,RW(KD2+1),NEN,RWC(NCEQ3+1),NCO,IWP(NCEQ2+1), 
     &             TEMP,GAS1,0)                                          
      CALL ENTROP (NSP,TEMP,GAS1,RW(KD1+1),NEN,NSK,RWC(NCEQ3+1),         
     &             RWC(MCH1+1),NCO,IWP(NCEQ2+1),JFLAG)                   
69    NKIP=0                                                             
      IF(MODEL.EQ.1.OR.IPRINT.EQ.0) CALL KIP1 (NSP,NCEQ,NHCP,RW(KD2+1)   
     &      ,RW(KD1+1),NAE,RWC(1),IWP(NCEQ+1),TEMP,RK,GAS1,GAS2,KFAIL)   
      IF(MODEL.EQ.1) GOTO 100                                            
C                                                                        
C----------------------------------------------------------------------- 
C  COMPUTATION OF REVERSE KINETIC PARAMETERS FOR KREV=2                  
C----------------------------------------------------------------------- 
      IF(KREV.NE.2) GOTO 100                                             
      DO 92 I=1,NCEQ                                                     
92    IF(IWP(NCEQ+I).EQ.5) GOTO 91                                       
      GOTO 100                                                           
91    IPREV=2*N+NCTB+4*NSP                                               
C  KINETIC PARAMETERS FOR TEMP=TLOW (RK(1),...,RK(NCEQ))                 
      CALL THERMO (NSP,NHCP,RW(KD2+1),NEN,RWC(NCEQ3+1),NCO,IWP(NCEQ2+1), 
     &             TLOW,GAS1,0)                                          
      CALL ENTROP (NSP,TLOW,GAS1,RW(KD1+1),NEN,NSK,RWC(NCEQ3+1),         
     &             RWC(MCH1+1),NCO,IWP(NCEQ2+1),JFLAG)                   
      CALL KIP1 (NSP,NCEQ,NHCP,RW(KD2+1),RW(KD1+1),NAE,                  
     &         RWC(1),IWP(NCEQ+1),TLOW,RK,GAS1,GAS2,KFAIL)               
C  KINETIC PARAMETERS FOR TEMP=THIGH (RW(IPREV+1),...,RW(IPREV(NCEQ))    
      CALL THERMO (NSP,NHCP,RW(KD2+1),NEN,RWC(NCEQ3+1),NCO,IWP(NCEQ2+1), 
     &             THIGH,GAS1,0)                                         
      CALL ENTROP (NSP,THIGH,GAS1,RW(KD1+1),NEN,NSK,RWC(NCEQ3+1),        
     &             RWC(MCH1+1),NCO,IWP(NCEQ2+1),JFLAG)                   
      CALL KIP1 (NSP,NCEQ,NHCP,RW(KD2+1),RW(KD1+1),NAE,                  
     &           RWC(1),IWP(NCEQ+1),THIGH,RW(IPREV+1),GAS1,GAS2,KFAIL)   
C  SET CONSTANTS                                                         
      XLTT=DLOG(TLOW/THIGH)                                              
CSP       XLTT=ALOG(TLOW/THIGH)                                          
      RGAS=GAS2                                                          
      DIV=(ONE/THIGH-ONE/TLOW)/RGAS                                      
      XLTLOW=DLOG(TLOW)                                                  
CSP       XLTLOW=ALOG(TLOW)                                              
      DIV1=ONE/(RGAS*TLOW)                                               
      DO 90 K=1,NCEQ                                                     
      K1=NCEQ+K                                                          
      IF(IWP(K1).NE.5) GOTO 90                                           
      K2=NCEQ+K1                                                         
      K3=IPREV+K                                                         
CC      RWC(K2)=RWC(K2-1)                                                
      IF(RK(K).LT.SMALLH) RK(K)=SMALLH                                   
      IF(RW(K3).LT.SMALLH) RW(K3)=SMALLH                                 
      U=DLOG(RK(K)/RW(K3))                                               
CSP       U=ALOG(RK(K)/RW(K3))                                           
      RWC(K1)=(U-RWC(K2)*XLTT)/DIV                                       
      U=DLOG(RK(K))                                                      
CSP       U=ALOG(RK(K))                                                  
      RWC(K)=U-RWC(K2)*XLTLOW+DIV1*RWC(K1)                               
      IKIN=7                                                             
      IF(RWC(K2).EQ.ZERO) IKIN=6                                         
      IWP(K1)=IKIN                                                       
C  WRITE KINETIC PARAMETERS ON DATASET 'HDS'                             
      U=DEXP(RWC(K))                                                     
CSP       U=EXP(RWC(K))                                                  
90    WRITE(HDS,7000) K,RWC(K),RWC(K1),RWC(K2)                           
7000  FORMAT(1X,'(',I4,')',4X,3D20.10)                                   
CSP 7000  FORMAT(1X,'(',I4,')',4X,3E20.10)                               
C                                                                        
100   CONTINUE                                                           
C                                                                        
C----------------------------------------------------------------------- 
C  WRITE IN OUTPUT DATABASE 'OUT'                                        
C----------------------------------------------------------------------- 
C                                                                        
C  DIMENSIONS                                                            
      WRITE(DOUT,8002)                                                   
      WRITE(DOUT,8003) NSP,MODEL,NCEQ,IEL,NTB,NLR,IUNIT,ITPM             
      WRITE(DOUT,8020)                                                   
C                                                                        
C  SYMBOLIC REPRESENTATION OF FCN                                        
      IF(IPRINT.EQ.4) GOTO 35                                            
      IF(IPRINT.LT.2) GOTO 35                                            
      WRITE(DOUT,8014)                                                   
      WRITE(DOUT,8016) (LR(L),L=1,NLR)                                   
      WRITE(DOUT,8020)                                                   
      WRITE(DOUT,8015)                                                   
      WRITE(DOUT,8016) (PLR(L),L=1,NCEQ2)                                
      WRITE(DOUT,8020)                                                   
C                                                                        
C  OUTPUT POINTS                                                         
35    WRITE(DOUT,8004)                                                   
      WRITE(DOUT,8005) (TP(L),L=1,ITPM)                                  
      WRITE(DOUT,8020)                                                   
C                                                                        
C  ELEMENT LIST                                                          
      IF(IEL.LE.0) GOTO 40                                               
      WRITE(DOUT,8006)                                                   
      WRITE(DOUT,8007) ((ELEM(I,L),L=1,5),I=1,IEL)                       
      WRITE(DOUT,8020)                                                   
40    CONTINUE                                                           
C                                                                        
C  ELEMENT ATOMIC WEIGHT                                                 
      IF(IEL.LE.0) GOTO 50                                               
      WRITE(DOUT,8008)                                                   
      WRITE(DOUT,8013) (ATWS(L),L=1,IEL)                                 
      WRITE(DOUT,8020)                                                   
50    CONTINUE                                                           
C                                                                        
C  SPECIES LIST                                                          
      WRITE(DOUT,8009)                                                   
      IH=NSP*10                                                          
      IH1=NTB*10                                                         
      IF(IH1.LE.0) WRITE(DOUT,8010) (CW(L),L=1,IH)                       
      IF(IH1.GT.0) WRITE(DOUT,8010) (CW(L),L=1,IH),                      
     &                              (NAMETB(L1),L1=1,IH1)                
      WRITE(DOUT,8020)                                                   
C                                                                        
C  SPECIES MOLECULAR WEIGHT                                              
      IF(IEL.LE.0) GOTO 60                                               
      WRITE(DOUT,8011)                                                   
      WRITE(DOUT,8013) (RW(MH1+I),I=1,NSP)                               
      WRITE(DOUT,8020)                                                   
60    CONTINUE                                                           
C                                                                        
C  THIRD BODIES                                                          
      IF(NTB.EQ.0) GOTO 80                                               
      WRITE(DOUT,8017) (LENTB(L),L=1,5)                                  
      WRITE(DOUT,8016) (IW(J),J=1,NCTB)                                  
      WRITE(DOUT,8020)                                                   
      WRITE(DOUT,8019)                                                   
      WRITE(DOUT,8013) (RW(N+L),L=1,NCTB)                                
      WRITE(DOUT,8020)                                                   
80    CONTINUE                                                           
C                                                                        
C  KINETIC PARAMETERS FOR MODEL=1                                        
      IF(MODEL.GT.1.AND.IPRINT.GT.0) GOTO 81                             
      WRITE(DOUT,8012)                                                   
      WRITE(DOUT,8018) TP(1)                                             
      WRITE(DOUT,8013) (RK(I),I=1,NCEQ)                                  
      WRITE(DOUT,8020)                                                   
81    CONTINUE                                                           
C                                                                        
C----------------------------------------------------------------------- 
C  PERFORM SIMULATION (INTEGRATION) OF THE KINETIC REACTION MODEL        
C----------------------------------------------------------------------- 
C                                                                        
C  SET INPUT PARAMETERS                                                  
      NRWH=NRW-N-NCTB                                                    
      NIWH=NIW-NCTB                                                      
      NH1=N+NCTB+1                                                       
      NH2=NCTB+1
      HMAX=TP(ITPM)-TP(1)                                                
      NRWCH=NCEQ3+2                                                      
      IF(ICO.GT.0) NRWCH=NCEQ3+14*NCO                                    
      NIWPH=NCEQ2+NCO+LDIM+N+N+1                                         
C                                                                        
C  CALL INTEGRATION SUBROUTINE                                           
      CALL META1 (N,NSP,NCEQ,LDIM,NTB,NCTB,NRWH,NIWH,RW(1),RW(N+1),      
     & IW(1),IWP,NIWPH,TP,ITPM,HMAX,EPS,RW(NH1),RWC,NRWCH,IW(NH2),GAS1,  
     & GAS2,MODEL,IDISC,LENTB,ICO,KFLAG,IPRINT)                          
C                                                                        
C  CHECK FOR POSSIBLE REDUCTION OF WORKSPACE (RWC AND IWP)               
      NFREE=NRWC-NRWCH                                                   
      IF(NFREE.GT.0) WRITE(ITOUT,9500) NFREE                             
9500  FORMAT(/,' REAL WORK SPACE RWC MAY BE REDUCED BY',I7,              
     &         ' STORAGE LOCATIONS')                                     
      IF(360.GT.NIWPH) NIWPH=360                                         
      NFREE=NIWP-NIWPH                                                   
      IF(NFREE.GT.0) WRITE(ITOUT,9501) NFREE                             
9501  FORMAT(/,' INTEGER WORK SPACE IWP MAY BE REDUCED BY',I7,           
     &         ' STORAGE LOCATIONS')                                     
C----------------------------------------------------------------------- 
C                 NORMAL RETURN, NO ERROR                                
C----------------------------------------------------------------------- 
C                                                                        
      RETURN                                                             
C----------------------------------------------------------------------- 
C                   ERROR RETURN                                         
C----------------------------------------------------------------------- 
9801  WRITE(ITOUT,9802) MSIW                                             
9802  FORMAT(/,' INTEGER WORK SPACE EXHAUSTED; ENLARGE IW BY AT',        
     &  ' LEAST',/,I6,' STORAGE LOCATIONS')                              
      GOTO 9999                                                          
9803  WRITE(ITOUT,9804) MSRW                                             
9804  FORMAT(/,' REAL WORK SPACE EXHAUSTED; ENLARGE RW BY AT',           
     &  ' LEAST',/,I6,' STORAGE LOCATIONS')                              
      GOTO 9999                                                          
9805  WRITE(ITOUT,9806) MSRWC                                            
9806  FORMAT(/,' REAL WORK SPACE EXHAUSTED; ENLARGE RWC BY AT',          
     &  ' LEAST',/,I6,' STORAGE LOCATIONS')                              
      GOTO 9999                                                          
9807  WRITE(ITOUT,9808) MSIWP                                            
9808  FORMAT(/,' INTEGER WORK SPACE EXHAUSTED; ENLARGE IWP BY AT',       
     &  ' LEAST',/,I6,' STORAGE LOCATIONS')                              
      GOTO 9999                                                          
9809  WRITE(ITOUT,9810) MSCW                                             
9810  FORMAT(/,' CHARACTER WORK SPACE EXHAUSTED; ENLARGE CW BY AT',      
     &  ' LAEST',/,I6,' STORAGE LOCATIONS')                              
9999  KFLAG=-1                                                           
      RETURN                                                             
C                                                                        
C----------------------------------------------------------------------- 
C               FORMATS (FOR FILE 'SIN' AND TERMINAL)                    
C----------------------------------------------------------------------- 
9000  FORMAT(19X,I5)                                                     
9001  FORMAT(I5)                                                         
9002  FORMAT(/,' REVERSE REACTION CONSTANTS ONLY VALID FOR A USER',      
     & ' PRESCRIBED',/,' TEMPERATURE RANGE FROM',F8.2,' TO',F8.2,        
     & ' IF INTERNAL COMPUTATION DEMANDED')                              
9003  FORMAT(//,19X,'ACTUAL MODEL:',/,19X,13('-'),/)                     
9004  FORMAT(/,' REVERSE REACTION CONSTANTS ARE RECOMPUTED BY',          
     &  ' MEANS OF',/,' MASS-ACTION-LAW AT EVERY TEMPERATURE CHANGE',    
     & ', IF NOT GIVEN BY THE USER')                                     
9005  FORMAT(72A1)                                                       
9007  FORMAT(5A1,3X,D21.12)                                              
CSP 9007  FORMAT(5A1,3X,E21.12)                                          
9009  FORMAT(I4,1X,10A1,1X,D21.12,1X,D21.12,1X,2(I8),/,8(I8))            
CSP 9009  FORMAT(I4,1X,10A1,1X,E21.12,1X,E21.12,1X,2(I8),/,8(I8))        
9010  FORMAT(I5,1X,10A1,2X,D21.12)                                       
CSP 9010  FORMAT(I5,1X,10A1,2X,E21.12)                                   
9015  FORMAT(19X,D21.12)                                                 
CSP 9015  FORMAT(19X,E21.12)                                             
9016  FORMAT(19X,2D21.12,I3)                                             
CSP 9016  FORMAT(19X,2E21.12,I3)                                         
9020  FORMAT(9X,I5,3X,10A1)                                              
9030  FORMAT(5X,I5,D21.12)                                               
CSP 9030  FORMAT(5X,I5,E21.12)                                           
9040  FORMAT(3D24.14)                                                    
CSP 9040  FORMAT(3D24.12)                                                
9050  FORMAT(D21.12)                                                     
CSP 9050  FORMAT(E21.12)                                                 
9060  FORMAT(5D15.8,/,5D15.8,/,4D15.8)                                   
CSP 9060  FORMAT(5E15.8,/,5E15.8,/,4E15.8)                               
9070  FORMAT(36I2)                                                       
9080  FORMAT(9I8)                                                        
9100  FORMAT(/,10X,' *** ISOTHERMIC MODEL ***')                          
9150  FORMAT(/,10X,' *** PRESCRIBED TEMPERATURE MODEL ***')              
9200  FORMAT(/,10X,' *** CONSTANT DENSITY MODEL ***')                    
9250  FORMAT(/,10X,' *** PRESCRIBED DENSITY MODEL ***')                  
9300  FORMAT(/,10X,' *** CONSTANT PRESSURE MODEL ***')                   
9350  FORMAT(/,10X,' *** PRESCRIBED PRESSURE MODEL ***')                 
C                                                                        
C----------------------------------------------------------------------- 
C                 FORMATS (FOR FILE 'DOUT')                              
C----------------------------------------------------------------------- 
C                                                                        
8000  FORMAT(' &HEAD',9X,I3)                                             
8001  FORMAT(1X,72A1)                                                    
8002  FORMAT(' &DIM')                                                    
8003  FORMAT(' NSP=',I6,', MODEL=',I1,', NCEQ=',I6,', IEL=',I6,',',/,    
     & ' NTB=',I1,', NLR=',I6,', IU=',I1,', ITPM=',I3,',')               
8004  FORMAT(' &TOUT',/,', TP=')                                         
8005  FORMAT(3(1X,D21.12,','))                                           
CSP 8005  FORMAT(3(1X,E21.12,','))                                       
8006  FORMAT(' &ELEMN')                                                  
8007  FORMAT(7(1X,5A1,'   ,'))                                           
8008  FORMAT(' &EATW',/,' RK=')                                          
8009  FORMAT(' &NAME')                                                   
8010  FORMAT(5(1X,10A1,'  ,'))                                           
8011  FORMAT(' &SATW',/,' RK=')                                          
8012  FORMAT(' &KIPA')                                                   
8018  FORMAT(' T=' ,D21.12,', RK=')                                      
CSP 8018  FORMAT(' T=' ,E21.12,', RK=')                                  
8013  FORMAT(5(1X,D13.6,','))                                            
CSP 8013  FORMAT(5(1X,E13.6,','))                                        
8014  FORMAT(' &FCN1',/,' LR=')                                          
8015  FORMAT(' &FCN2',/,' PLR=')                                         
8016  FORMAT(12(I5,','))                                                 
8017  FORMAT(' &TBO1',/,' LR=',5(I5,','))                                
8019  FORMAT(' &TBO2',/,' C=')                                           
8020  FORMAT(' &END')                                                    
C                                                                        
C----------------------------------------------------------------------- 
C                 END OF SUBROUTINE SIMULA                               
C----------------------------------------------------------------------- 
C                                                                        
      END                                                                
      SUBROUTINE SIMDAT (N,NSP,NCEQ,MODEL,CONC,GAS1,T,IPRINT,LPRINT)     
C                                                                        
C*********************************************************************** 
C                                                                      * 
C  DATE OF LATEST CHANGE: MARCH 1, '86                                 * 
C                                                                      * 
C*********************************************************************** 
C                                                                      * 
C  SUBROUTINE SIMDAT PERFORMS OUTPUT HANDLING OF THE COMPUTED          * 
C  SOLUTION (CONCENTRATIONS, TEMPERATURE, DENSITY, PRESSURE) AND       * 
C  ADDITIONAL CHEMICAL INFORMATION                                     * 
C                                                                      * 
C*********************************************************************** 
C                                                                      * 
C  INPUT PARAMETERS                                                    * 
C------------------                                                    * 
C                                                                      * 
C      NSP               NUMBER OF FIRST ORDER DIFFERENTIAL EQUATIONS  * 
C                        (NUMBER OF CHEMICAL SPECIES)                  * 
C      NCEQ              NUMBER OF CHEMICAL EQUATIONS                  * 
C      MODEL             MODEL PARAMETER (SEE DRIVER)                  * 
C      CONC(NSP)         CONCENTRATIONS AT T                           * 
C      TOUT(IOUT)        PRESCRIBED OUTPUT POINTS                      * 
C      IOUT              NUMBER OF OUTPUT POINTS                       * 
C      T                 ACTUAL INTEGRATION POINT                      * 
C      TOLD              LAST INTEGRATION POINT                        * 
C      IPRINT            PRINT PARAMETER FOR ADDITIONAL CHEMICAL       * 
C                        OUTPUT                                        * 
C                        SEE DESCRIPTION BELOW                         * 
C      LPRINT            LOGICAL PRINT PARAMETER                       * 
C                        TRUE : IF TP(I).EQ.T                          * 
C                        FALSE : ELSE                                  * 
C                                                                      * 
C  COMMON BLOCKS                                                       * 
C---------------                                                       * 
C                                                                      *  
C                                                                      * 
C     INTEGER          ITIN,ITOUT,MOUT,EOUT,CIN,HDS,POUT,SIN,DOUT,THDAT* 
C     COMMON / LUNIT / ITIN,ITOUT,MOUT,EOUT,CIN,HDS,POUT,SIN,DOUT,THDAT* 
C                                                                      * 
C   LOGICAL UNITS:                                                     * 
C   ITIN,ITOUT,MOUT,EOUT,CIN,HDS,POUT,SIN,THDAT: UNITS NOT USED        * 
C   DOUT: LOGICAL UNIT FOR DATABASE 'OUT'                              * 
C                                                                      * 
C     COMMON / KINPAR / RK(NCEQ)                                       * 
C   KINETIC PARAMETERS                                                 * 
C                                                                      * 
C     COMMON / RATE1 / R(NCEQ)                                         * 
C   R(I): RATE OF EQUATION I AT TIME T                                 * 
C                                                                      * 
C     COMMON / DELHR / DHR(NCEQ)                                       * 
C   DHR(J): ENTHALPY PER J-TH REACTION                                 * 
C                                                                      * 
C----------------------------------------------------------------------* 
C  ACTUAL IMPLEMENTED OUTPUT HANDLING:                                 * 
C  (MAY BE ALTERED BY THE SKILLFULL USER)                              * 
C----------------------------------------------------------------------* 
C                                                                      * 
C  AT EACH INTEGRATION POINT (I.E. CALL OF SIMDAT) THE FOLLOWING DATA  * 
C  ARE WRITTEN ON OUTPUT DATASET 'OUT':                                * 
C  1.) VALUE OF T                                                      * 
C  2.) VALUES OF CONC(I),I=1,NSP                                       * 
C  3.) VALUE OF T                                                      * 
C  4.) VALUE OF P,  P = SUM OF CONC(I),I=1,NSP                         * 
C  5.) VALUE OF T                                                      * 
C  6.) REACTION CONDITIONS (TEMPERATURE,DENSITY,PRESSURE,VOLUME)       * 
C                                                                      * 
C  ADDITIONALLY, FOR                                                   * 
C  IPRINT = 1  AND  T=TOUT(I),I=1,IOUT  (I.E. LPRINT = TRUE):          * 
C    7.) T                                                             * 
C    8.) RK(I),I=1,NCEQ (KINETIC PARAMETERS)                           * 
C  IPRINT = 2  AND  T=TOUT(I),I=1,IOUT  (I.E. LPRINT = TRUE):          * 
C    7.) T                                                             * 
C    8.) R(I),I=1,NCEQ   (RATES)                                       * 
C    9.) T                                                             * 
C   10.) DHR(I),I=1,NCEQ  (ENTHALPY PER REACTION)                      * 
C  IPRINT = 3  AND  T=TOUT(I),I=1,IOUT  (I.E. LPRINT = TRUE):          * 
C    7.) T                                                             * 
C    8.) RK(I),I=1,NCEQ                                                * 
C    9.) T                                                             * 
C   10.) R(I),I=1,NCEQ                                                 * 
C   11.) T                                                             * 
C   12.) DHR(I),I=1,NCEQ  (ENTHALPY PER REACTION)                      * 
C                                                                      * 
C                                                                      * 
C  FOR:                                                                * 
C  IPRINT = 4:                                                         * 
C    7.) T                                                             * 
C    8.) RK(I),I=1,NCEQ (KINETIC PARAMETERS)                           * 
C  IPRINT = 5:                                                         * 
C    7.) T                                                             * 
C    8.) R(I),I=1,NCEQ   (RATES)                                       * 
C    9.) T                                                             * 
C   10.) DHR(I),I=1,NCEQ  (ENTHALPY PER REACTION)                      * 
C  IPRINT = 6:                                                         * 
C    7.) T                                                             * 
C    8.) RK(I),I=1,NCEQ                                                * 
C    9.) T                                                             * 
C   10.) R(I),I=1,NCEQ                                                 * 
C   11.) T                                                             * 
C   12.) DHR(I),I=1,NCEQ  (ENTHALPY PER REACTION)                      * 
C                                                                      * 
C*********************************************************************** 
C                                                                        
C                                                                        
      DOUBLE PRECISION CONC(N),P,RK,T,ZERO,TEMP,DENSY,PRESS,VOLUME       
CSP       REAL CONC(N),P,RK,T,ZERO,TEMP,DENSY,PRESS,VOLUME               
      DOUBLE PRECISION GAS1,XMASS                                        
CSP       REAL GAS1,XMASS                                                
C                                                                        
      REAL R,DHR                                                         
      INTEGER I,IPR                                                      
C                                                                        
      LOGICAL LPRINT                                                     
C                                                                        
      INTEGER          ITIN,ITOUT,MOUT,EOUT,CIN,HDS,POUT,SIN,DOUT,THDAT  
      COMMON / LUNIT / ITIN,ITOUT,MOUT,EOUT,CIN,HDS,POUT,SIN,DOUT,THDAT  
C                                                                        
C                                                                        
      COMMON / RECOND / TEMP,DENSY,PRESS,VOLUME,XMASS                    
C                                                                        
      COMMON / KINPAR / RK(1000)                                         
      COMMON / RATE1  / R(1000)                                          
      COMMON / DELHR  / DHR(1000) 
      SAVE /LUNIT/, /RECOND/, /KINPAR/, /RATE1/, /DELHR/
C                                                                        
      DATA ZERO/0.D0/                                                    
CSP       DATA ZERO/0.E0/                                                
C                                                                        
C  WRITE CONCENTRATIONS                                                  
      WRITE(DOUT,9001)                                                   
      WRITE(DOUT,9005) T                                                 
      WRITE(DOUT,9010) (CONC(I),I=1,NSP)                                 
      WRITE(DOUT,9000)                                                   
C                                                                        
C  COMPUTE SUM OF CONCENTRATIONS                                         
      P=ZERO                                                             
      DO 11 I=1,NSP                                                      
11    P=P+CONC(I)                                                        
C                                                                        
C  WRITE SUM                                                             
      WRITE(DOUT,9100)                                                   
      WRITE(DOUT,9105) T,P                                               
      WRITE(DOUT,9000)                                                   
C                                                                        
C  WRITE REACTION CONDITIONS                                             
      IF(MODEL.LE.2) PRESS=P*TEMP*GAS1                                   
      WRITE(DOUT,9500)                                                   
      WRITE(DOUT,9505) T,TEMP,DENSY,PRESS,VOLUME                         
      WRITE(DOUT,9000)                                                   
C                                                                        
      IF(IPRINT.EQ.0) RETURN                                             
      IF(IPRINT.LE.3.AND..NOT.LPRINT) RETURN                             
C                                                                        
C  WRITE ADDITIONAL CHEMICAL OUTPUT ACCORDING TO IPRINT                  
      IPR=IPRINT                                                         
      IF(IPRINT.GT.3) IPR=IPR-3                                          
      IF(IPR.NE.1.AND.IPR.NE.3) GOTO 20                                  
C                                                                        
C  WRITE KINETIC PARAMETERS                                              
      IF(MODEL.EQ.1) GOTO 20                                             
      WRITE(DOUT,9400)                                                   
      WRITE(DOUT,9305) T                                                 
      WRITE(DOUT,9010) (RK(I),I=1,NCEQ)                                  
      WRITE(DOUT,9000)                                                   
20    IF(IPR.EQ.1) RETURN                                                
C                                                                        
C  WRITE RATES                                                           
      WRITE(DOUT,9300)                                                   
      WRITE(DOUT,9305) T                                                 
      WRITE(DOUT,9010) (R(I),I=1,NCEQ)                                   
      WRITE(DOUT,9000)                                                   
C                                                                        
      IF(MODEL.LE.2) RETURN                                              
C                                                                        
C  WRITE ENTHALPIES OF REACTIONS                                         
      WRITE(DOUT,9600)                                                   
      WRITE(DOUT,9011) T                                                 
      WRITE(DOUT,9010) (DHR(I),I=1,NCEQ)                                 
      WRITE(DOUT,9000)                                                   
C                                                                        
      RETURN                                                             
C                                                                        
C  THE FOLLOWING FORMATS ARE CHOSEN TO ALLOW THE READ WITH NAMELISTS     
C IN SUBROUTINE LIST                                                     
C                                                                        
9000  FORMAT(' &END')                                                    
9001  FORMAT(' &CONC')                                                   
9005  FORMAT(' T=' ,D21.12,', C=')                                       
CSP 9005  FORMAT(' T=' ,E21.12,', C=')                                   
9010  FORMAT(5(1X,D13.6,','))                                            
CSP 9010  FORMAT(5(1X,E13.6,','))                                        
9100  FORMAT(' &SUMC')                                                   
9105  FORMAT(' T=' ,D21.12,', SUM=',D13.6)                               
CSP 9105  FORMAT(' T=' ,E21.12,', SUM=',E13.6)                           
9300  FORMAT(' &RATE')                                                   
9305  FORMAT(' T=' ,D21.12,', RK=')                                      
CSP 9305  FORMAT(' T=' ,E21.12,', RK=')                                  
9400  FORMAT(' &KIPA')                                                   
9500  FORMAT(' &RECO')                                                   
9505  FORMAT(' T=',D21.12,',',/,' TEMP= ',D14.7,', DENSY= ',             
CSP 9505  FORMAT(' T=',E21.12,',',/,' TEMP= ',E14.7,', DENSY= ',         
     & D14.7,',',/,' PRESS= ',D14.7,', VOLUME= ',D14.7,',')              
CSP      & E14.7,',',/,' PRESS= ',E14.7,', VOLUME= ',E14.7,',')          
9600  FORMAT(' &DHRE')                                                   
9011  FORMAT(' T=' ,D21.12,', DHR=')                                     
CSP 9011  FORMAT(' T=' ,E21.12,', DHR=')                                 
C                                                                        
C----------------------------------------------------------------------- 
C  END OF SUBROUTINE SIMDAT                                              
C----------------------------------------------------------------------- 
      END                                                                
      SUBROUTINE META1 (N,NSP,NCEQ,LDIM,NTB,NCTB,NW,NIW,Y,COE,ICOE,      
     &  JW,NJW,TP,ITPM,HMAX,EPS,W,WC,NWC,IW,GAS1,GAS2,MODEL,IDISC,LENTB,  
     &   ICO,KFLAG,IPRINT)                                                 
C*********************************************************************** 
C                                                                      * 
C DATE OF LATEST CHANGE: MARCH 1, '86                                  * 
C                                                                      * 
C*********************************************************************** 
C                                                                      * 
C  INTERFACE ROUTINE FOR INTEGRATION METHOD                            * 
C                                                                      * 
C*********************************************************************** 
C                                                                      * 
C                                                                      * 
C   INPUT PARAMETERS OF SUBROUTINE META1:                              * 
C----------------------------------------                              * 
C    (* MARKS TRANSIENT PARAMETERS )                                   * 
C                                                                      * 
C    N                  NUMBER OF FIRST-ORDER DIFFERENTIAL EQUATIONS   * 
C    NSP                NUMBER OF CHEMICAL SPECIES                     * 
C    NCEQ               NUMBER OF CHEMICAL REACTION EQUATIONS          * 
C    LDIM               NUMBER OF NONZERO ELEMENTS IN JACOBIAN         * 
C    NW                 DIMENSION OF REAL WORK SPACE W                 * 
C    NIW                DIMENSION OF INTEGER WORK SPACE IW             * 
C  * Y(N)               INITIAL VALUES (CONCENTRATIONS) Y(1),...,Y(N)  * 
C  * TP(ITPM)           REQUESTED OUTPUT POINTS                        * 
C    ITPM               NUMBER OF OUTPUT POINTS                        * 
C                       TP(1):    STARTING POINT OF INTEGRATION        * 
C                       TP(ITPM): END POINT OF INTEGRATION             * 
C    HMAX               USER PRESCRIBED MAXIMUM PERMITTED STEPSIZE     * 
C                       ( MAY CAUSE A NON-OPTIMAL INTEGRATION )        * 
C    EPS                PRESCRIBED RELATIVE PRECISION                  * 
C    MODEL              CHOSEN TYPE OF CHEMICAL SIMULATION             * 
C                       1 : ISOTHERMIC MODEL                           * 
C                       2 : PRESCRIBED TEMPERATURE PROFILE             * 
C                       3 : CONSTANT DENSITY MODEL                     * 
C                       4 : PRESCRIBED DENSITY PROFILE                 * 
C                       5 : CONSTANT PRESSURE MODEL                    * 
C                       6 : PRESCRIBED PRESSURE PROFILE                * 
C    IDISC              INTEGRATION METHOD FLAG                        * 
C                       0 : SEMI-IMPLICIT EULER METHOD                 * 
C                       1 : SEMI-IMPLICIT MID-POINT-RULE               * 
C    GAS1               GAS CONSTANT FOR IDEAL GAS LAW                 * 
C    GAS2               GAS CONSTANT FOR ARRHENIUS LAW                 * 
C  * KFLAG              PRINT PARAMETER FOR INTEGRATOR OUTPUT          * 
C                      -1 : NO OUTPUT                                  * 
C                       0 : ERROR MESSAGES                             * 
C                       1 : INTEGRATION MONITOR                        * 
C                       2 : EXTENDED INTEGRATION MONITOR               * 
C    IPRINT             PRINT PARAMETER FOR ADDITIONAL CHEMICAL DATA   * 
C    IW                 INTEGER WORK SPACE                             * 
C    W                  REAL WORK SPACE                                * 
C                                                                      * 
C                                                                      * 
C                                                                      * 
C   OUTPUT PARAMETERS                                                  * 
C--------------------                                                  * 
C                                                                      * 
C    TP(ITPM)           FINAL POINT OF ACTUALLY PERFORMED INTEGRATION  * 
C    Y(N)               FINAL VALUES AT T(ITPM)                        * 
C    H                  STEPSIZE PROPOSAL FOR NEXT INTEGRATION STEP    * 
C    KFLAG              ERROR FLAG                                     * 
C                       GE.0 : INTEGRATION SUCCESFULLY COMPLETED       * 
C                       LT.0 : INTEGRATION FAILED                      * 
C                                                                      * 
C   INTERNAL PARAMETERS                                                * 
C----------------------                                                * 
C                                                                      * 
C    ICO                NUMBER OF SPECIES, FOR WHICH THERMODYNAMICAL   * 
C                       COEFFICIENTS ARE SUPPLIED                      * 
C    NTB                NUMBER OF THIRD BODIES                         * 
C    NCTB               DIMENSION OF WORK SPACES COE AND ICOE          * 
C    COE(NCTB)          !  INTERNAL CODING OF THE                      * 
C    ICOE(NCTB)         !  THIRD BODIES USED IN REACTION SYSTEM        * 
C    NWC                DIMENSION OF REAL WORK SPACE WC                * 
C    NJW                DIMENSION OF INTEGER WORK SPACE JW             * 
C    JW                 INTEGER WORK SPACE                             * 
C    WC                 REAL WORK SPACE                                * 
C    THE WORK SPACES JW AND WC CONTAIN CHEMICAL INFORMATION,           * 
C    READ FROM DATASET 'INVAL' IN SUBROUTINE SIMULA.                   * 
C    LENTB(5)           INTEGER FIELD, LENTB(K) IS THE 'LENGTH' OF     * 
C                       THE K-TH THIRD BODY                            * 
C                                                                      * 
C                                                                      * 
C   EXTERNAL SUBROUTINES ( SUPPLIED BY 'LARKIN' )                      * 
C-----------------------                                               * 
C                                                                      * 
C    FCN , JACOBI , SIMDAT , MA30LA , MA30LB , MA30LC , MA30LM         * 
C    THERMO , ENTROP , KIP1                                            * 
C                                                                      * 
C   COMMON BLOCKS                                                      * 
C----------------                                                      * 
C                                                                      * 
CCKIR COMMON / LUNIT / ITIN,ITOUT,MOUT,EOUT,CIN,HDS,POUT,SIN,DOUT      * 
C                                                                        
C     INTEGER          ITIN,ITOUT,MOUT,EOUT,CIN,HDS,POUT,SIN,DOUT,THDAT  
C     COMMON / LUNIT / ITIN,ITOUT,MOUT,EOUT,CIN,HDS,POUT,SIN,DOUT,THDAT  
C                                                                        
C                                                                      * 
C     COMMON / COUNT / NSTEP,NFCN,NANFA,NTFAC,NFAC,NSOL                * 
C       COUNT PARAMETERS                                               * 
C         NSTEP :  NUMBER OF INTEGRATION STEPS                         * 
C         NFCN  :  NUMBER OF FCN-EVALUATIONS                           * 
C         NANFA :  NUMBER OF ANALYSE/FACTOR CALL                       * 
C         NTFAC :  NUMBER OF FACTOR CALL WITH PIVOT RATIO TEST         * 
C         NFAC  :  NUMBER OF SIMPLE FACTOR CALL                        * 
C         NSOL  :  NUMBER OF SOLVED LINER SYSTEMS                      * 
C                                                                      * 
C                                                                      * 
C                                                                      * 
C    FOR MORE DETAILS SEE DESCRIPTION OF SUBROUTINE METAS1             * 
C                                                                      * 
C*********************************************************************** 
C                                                                        
C                                                                        
CI4       INTEGER IW(NIW),JW(NJW),ICOE(NCTB)                             
      INTEGER*2 IW(NIW),JW(NJW),ICOE(NCTB)                               
C                                                                        
      CHARACTER*20 HEAD6                                                 
C                                                                        
      DOUBLE PRECISION W(NW),WC(NWC)                                     
CSP       REAL W(NW),WC(NWC)                                             
      DOUBLE PRECISION Y(N),COE(NCTB),TP(ITPM),ONE                       
CSP       REAL Y(N),COE(NCTB),TP(ITPM),ONE                               
      DOUBLE PRECISION CON1,CON2,EPS,EPSH,HI,HMAX,GAS1,GAS2              
CSP       REAL CON1,CON2,EPS,EPSH,HI,HMAX,GAS1,GAS2                      
      DOUBLE PRECISION TEMP,ZERO                                         
CSP       REAL TEMP,ZERO                                                 
C                                                                        
      INTEGER LENTB(5)                                                   
      INTEGER ITPM, JM, KF, LD, LDIM, M2, M4, M6, MIW, MIWH, MW, MWH, N1,
     & N10, N11, N2, N3, N4, N5, N6, N7, N8, N9, NAE, NB1, NC1, NC2, NCO,
     & NEN, NH, NHCP, NI1, NI2, NI3, NI4, NI5, NP13, NP5, NP7, NSK                
C                                                                        
      INTEGER          ITIN,ITOUT,MOUT,EOUT,CIN,HDS,POUT,SIN,DOUT,THDAT  
      COMMON / LUNIT / ITIN,ITOUT,MOUT,EOUT,CIN,HDS,POUT,SIN,DOUT,THDAT  
C                                                                        
C                                                                        
      COMMON / COUNT / NSTEP,NFCN,NANFA,NTFAC,NFAC,NSOL,NTHERM,NKIP  
      SAVE /LUNIT/, /COUNT/
C                                                                        
      DATA ZERO/0.D0/ , ONE/1.D0/                                        
CSP       DATA ZERO/0.E0/ , ONE/1.E0/                                    
      DATA HEAD6/'IBM-MESSUNG: METAS1'/                                  
C                                                                        
      KF=KFLAG                                                           
C                                                                        
C  PREPARATIONS TO CALL INTEGRATOR                                       
C                                                                        
C  SET COUNT PARAMETERS TO ZERO                                          
      NSTEP=0                                                            
      NFCN=0                                                             
      NANFA=0                                                            
      NTFAC=0                                                            
      NFAC=0                                                             
      NSOL=0                                                             
      NTHERM=0                                                           
C                                                                        
C  SET POINTERS (DIVIDE WORK SPACES)                                     
      JM=5                                                               
      IF(IDISC.LE.0) JM=7                                                
      NB1=1                                                              
      IF(MODEL.GE.5) NB1=NSP                                             
      NHCP=1                                                             
      IF(ICO.GT.0) NHCP=3*NSP                                            
      N1=N+1                                                             
      N2=N1+N                                                            
      N3=N2+N                                                            
      N4=N3+N                                                            
      N5=N4+NB1                                                          
      N6=N5+NB1                                                          
      N7=N6+N                                                            
      N8=N7+N                                                            
      N9=N8+JM*N                                                         
      N10=N9+NHCP                                                        
      N11=N10+NSP                                                        
      NH=N11+LDIM                                                        
      NP5=5*N+1                                                          
      NP7=7*N+1                                                          
      NP13=13*N+1                                                        
C                                                                        
      NCO=ICO                                                            
      IF(NCO.EQ.0) NCO=1                                                 
      NAE=3*NCEQ                                                         
      NEN=1                                                              
      IF(ICO.GT.0) NEN=12*NCO                                            
      NSK=1                                                              
      IF(ICO.GT.0) NSK=2*NCO                                             
C                                                                        
      NC1=1+NAE                                                          
      NC2=NC1+NEN                                                        
C                                                                        
      NI1=1+NCEQ                                                         
      NI2=NI1+NCEQ                                                       
      NI3=NI2+NCO                                                        
      NI4=NI3+LDIM                                                       
      NI5=NI4+N+1                                                        
C                                                                        
C  CALCULATE FREE SPACE IN WORK ARRAYS FOR STORAGE OF LU-DECOMPOSITION   
      MW=NW-NH+1                                                         
      MWH=MW                                                             
      LD=LDIM                                                            
      MIW=NIW-13*N                                                       
      MIWH=MIW                                                           
      M2=N+N                                                             
      M4=M2+M2                                                           
      M6=M4+M2                                                           
C                                                                        
C  SET INITIAL STEPSIZE GUESS HI                                         
C  ( FOR HI.EQ.0 , HI IS COMPUTETD INTERNALLY BY SUBROUTINE METAS1 )     
C                                                                        
      HI=ZERO                                                            
C                                                                        
C  SET PARAMETERS FOR SCALING (SEE DESCRIPTION OF SUBROUTINE METAS1)     
C                                                                        
C  FOR INCREASING COMPONENTS                                             
      CON1=1.D-12                                                        
CSP       CON1=1.E-12                                                    
C  FOR DECREASING COMPONENTS                                             
      CON2=CON1                                                          
      IF(EPS.GT.ZERO) CON2=ONE                                           
      IF(EPS.LT.ZERO) EPS=-EPS                                           
C  (FOR OSCILLATING SYSTEMS SET CON2:=CON1, INTERNALLY DONE FOR EPS<0)   
C                                                                        
C                                                                        
C     CALL MONINI (6,HEAD6,N)                                            
      CALL METAS1 (N,NSP,Y,TP,ITPM,EPS,CON1,CON2,HMAX,HI,KFLAG,IPRINT,   
     & W(1),W(N1),W(N2),W(N3),W(N4),W(N5),W(N6),W(N7),W(N8),W(N9),       
     & W(N10),W(N11),IW(1),IW(N1),IW(NP5),IW(NP7),IW(NP13),NHCP,M2,M4,   
     & M6,MWH,LDIM,MIWH,NCEQ,W(NH),GAS1,GAS2,MODEL,NB1,LENTB,IDISC,NTB,  
     & NCTB,N1,NSK,NAE,NEN,NCO,COE,ICOE,JW(1),JW(NI1),JW(NI2),JW(NI3),
     & JW(NI4),JW(NI5),WC(1),WC(NC1),WC(NC2),JM)
C     CALL MONEND                                                        
C                                                                        
C  PRINT STATISTICS                                                      
C      IF(KF.LT.2) GOTO 50                                               
      WRITE(ITOUT,8050) NSTEP,NFCN,NANFA,NTFAC,NFAC,NSOL,NTHERM,NKIP     
C                                                                        
C50    RETURN                                                            
      RETURN                                                             
8050  FORMAT(/,' ************',/,' STATISTICS ,'/,' ************',/,     
     &        ' STEPS    ',I10,/,                                         
     &        ' FCN-EVAL ',I10,/,                                         
     &        ' ANAL/FAC ',I10,/,                                         
     &        ' TEST-FAC ',I10,/,                                         
     &        ' FACTOR   ',I10,/,                                         
     &        ' SUBST    ',I10,/,                                         
     &        ' THERMO   ',I10,/,                                         
     &        ' KIN.PAR. ',I10)                                           
C                                                                        
C----------------------------------------------------------------------- 
C                END OF SUBROUTINE META1                                 
C----------------------------------------------------------------------- 
      END                                                                
      SUBROUTINE METAS1 (N,NSP,Y,TP,ITPM,EPS,CON1,CON2,HMAX,HI,KFLAG,    
     & IPRINT,YM,DELQ,DZ,DEL,BV0,BVK,SH,SM,DT,HCP,SK,A,IS,IKEEP,IWO,     
     & IW1,ILEQ,NHCP,M2,M4,M6,NWTMP,LDIM,NLEQ,NCEQ,WTMP,GAS1,GAS2,MODEL, 
     & NB1,LENTB,IDISC,NTB,NCTB,NIA,NSK,NAE,NEN,NCO,COE,ICOE,ICOLLI,     
     & KINEV,IPTD,JA,IA,IDA,AEA,ENTH,SK0,JM)                             
C                                                                        
C---------------------------------------------------------------------   
C  SUBROUTINE METAS1 IS A SPECIAL VERSION OF SUBROUTINE METAN1           
C  ADAPTED FOR LARGE SYSTEMS OF REACTION KINETICS                        
C---------------------------------------------------------------------   
C                                                                        
C                                                                        
C  METAS1 (IN CONTRAST TO METAN1) USES SPECIAL DEVICES FOR:              
C                                                                        
C        - LU-DECOMPOSITION OF SPARSE MATRIX                             
C        - SCALING OF THE ODE-SYSTEM COMING FROM REACTON KINETICS        
C        - OUTPUT PERFORMANCE                                            
C                                                                        
C                                                                        
C  REFERENCES:                                                           
C-------------                                                           
C /1/  G.BADER, P. DEUFLHARD:                                            
C      A SEMI-IMPLICIT MIDPOINT RULE FOR STIFF SYSTEMS OF ORDINARY       
C      DIFFERENTIAL EQUATIONS                                            
C      UNIV. HEIDELBERG, SFB 123: TECHN. REP. 114 (1981)                 
C                                                                        
C /2/  P. DEUFLHARD:                                                     
C      ORDER AND STEPSIZE CONTROL IN EXTRAPOLATION METHODS               
C      UNIV. HEIDELBERG, SFB 123: TECHN. REP. 93 (1980)                  
C                                                                        
C /3/  P. DEUFLHARD, G. BADER, U. NOWAK:                                 
C      LARKIN - A SOFTWARE PACKAGE FOR THE SIMULATION OF LARGE           
C      SYSTEMS ARISING IN CHEMICAL REACTION KINETICS                     
C      (UNIV. HEIDELBERG, SFB 123: TECHN. REP. 100 (1980))               
C      IN:                                                               
C      K.H. EBERT, P. DEUFLHARD, W. JAEGER (ED.):                        
C      MODELLING OF CHEMICAL REACTION SYSTEMS.                           
C      SPRINGER SERIES CHEM. PHYS. 18 (1981)                             
C                                                                        
C /4/  I. S. DUFF, U. NOWAK:                                             
C      ON SPARSE MATRIX TECHNIQUES IN A STIFF INTEGRATOR                 
C      OF EXTRAPOLATION TYP                                              
C      UNIV. HEIDELBERG, SFB 123: TECHN. REP. (1982)                     
C                                                                        
C                                                                        
C*********************************************************************** 
C                                                                      * 
C DATE OF LATEST CHANGE:  MARCH 1, '86                                 * 
C                                                                      * 
C*********************************************************************** 
C                                                                      * 
C                                                                      * 
C   INPUT PARAMETERS OF SUBROUTINE METAS1:                             * 
C-----------------------------------------                             * 
C    ( * MARKS TRANSIENT PARAMETERS )                                  * 
C                                                                      * 
C    N                  NUMBER OF FIRST-ORDER DIFFERENTIAL EQUATIONS   * 
C    NSP                NUMBER OF CHEMICAL SPECIES                     * 
C  * Y(N)               INITIAL VALUES Y(1),...,Y(N)                   * 
C  * TP(I)              REQUESTED OUTPUT POINTS                        * 
C    ITPM               NUMBER OF OUTPUT POINTS                        * 
C                       TP(1):    STARTING POINT OF INTEGRATION        * 
C                       TP(ITPM): FINAL POINT OF INTEGRATION           * 
C    EPS                PRESCRIBED RELATIVE PRECISION                  * 
C    CON1               PARAMETER FOR SCALING FOR RISING COMPONENTS:   * 
C                       Y(I) IS COMPUTED DUE TO A STRICT RELATIVE      * 
C                       ERROR CRITERION IF Y(I).GT.YMAX  *CON1         * 
C                       ( YMAX:= MAXIMUM VALUE OF INITIAL CONC. )      * 
C    CON2               PARAMETER FOR SCALING OF DECREASING COMP.      * 
C                       Y(I) IS COMPUTED DUE TO A STRICT RELATIVE      * 
C                       ERROR CRITERION IF Y(I).GT.YMAX(I)  *(CON2/EPS)* 
C    HMAX               MAXIMUM PERMITTED STEPSIZE                     * 
C                       ( MAY CAUSE A NON-OPTIMAL INTEGRATION )        * 
C  * HI                 USER PRESCRIBED INITIAL STEPSIZE GUESS         * 
C                       IF HI.LE.0, HI IS COMPUTED INTERNALLY          * 
C  * KFLAG              PRINT PARAMETER FOR INTEGRATOR OUTPUT          * 
C                       -1 : NO OUTPUT                                 * 
C                        0 : ERROR MESSAGES                            * 
C                        1 : INTEGRATION MONITOR                       * 
C                        2 : EXTENDED INTEGRATION MONITOR              * 
C    IPRINT             PRINT PARAMETER FOR ADDITIONAL CHEMICAL OUTP.  * 
C                       ( .GT.0: RATES AND/OR INTEGRATED RATES ARE     * 
C                         INTERNALLY COMPUTED, SEE DESCRIPTION OF      * 
C                         SUBROUTINE SIMDAT )                          * 
C                                                                      * 
C  * NWTMP               DIMENSION OF ARRAY WTMP (ESTIMATED NUMBER OF  * 
C                        STORAGE LOCATIONS FOR LU-DECOMPOSITION)       * 
C                                                                      * 
C   PARAMETERS FOR INTERNAL USE ONLY                                   * 
C-----------------------------------                                   * 
C   INTEGER ARRAYS: IS(N) , IKEEP(M4) , IWO(M2) , IW1(M6) ,            * 
C                   ILEQ(NLEQ) , IPTD(NCO) , ICOE(NCTB) ,              * 
C                   ICOLLI(NCEQ) , KINEV(NCEQ)                         * 
C   FOR SPARSE PATTERN OF JACOBIAN: JA(LDIM) , IA(NIA) , IDA(N)        * 
C                                                                      * 
C   REAL ARRAYS: YM(N) , DELQ(N) , DZ(N) , DEL(N) ,                    * 
C                SH(N) , SM(N) , DT(N,JM) , SK(NSP) , BV0(NB1) ,       * 
C                BVK(NB1) , HCP(NHCP) , COE(NCTB) , AEA(NAE) ,         * 
C                ENTH(NEN) , SK0(NSK) , COLLI(5)                       * 
C                A(LDIM) , WTMP(NWTMP)                                 * 
C                                                                      * 
C   ASSOCIATED DIMENSIONS:                                             * 
C   M2=2*N , M4=4*N , M6=6*N                                           * 
C   NAE=3*NCEQ , NEN=12*NCO , NHCP=3*NSP , NSK=2*NSP , NIA=N+1         * 
C   NCEQ = NUMBER OF CHEMICAL EQUATIONS                                * 
C   LDIM = DIMENSION OF ARRAY A WHICH CONTAINS THE NONZERO VALUES      * 
C          OF JACOBIAN                                                 * 
C   NWTMP = DIMENSION OF ARRAY WTMP WHICH IS USED FOR STORAGE OF       * 
C           LU DECOMPOSITION OF JACOBIAN AND ADDDITIONALLY PROVIDES    * 
C           A WORK ARRAY OF LENGTH N                                   * 
C   NLEQ = DIMENSION OF ARRAY ILEQ WHICH IS USED FOR STORAGE OF        * 
C          COLUMN INDICES OF DECOMPOSED JACOBIAN AND ADDITIONALLY      * 
C          PROVIDES A WORK ARRAY OF LENGTH LDIM+N                      * 
C   NB1  = NSP, IF MODEL.GE.5                                          * 
C        = 1    ELSE                                                   * 
C   NCTB = SUM(1,..K) LENTB(K), IF NTB.GT.0  (SEE META)                * 
C        = 1    ELSE                                                   * 
C   NCO  = NUMBER OF SPECIES, FOR WHICH THERMODYNAMICAL DATA IS        * 
C          SUPPLIED BY THE USER                                        * 
C                                                                      * 
C   FOR FURTHER DESCRIPTION SEE SUBROUTINE META1                       * 
C                                                                      * 
C   OUTPUT PARAMETERS                                                  * 
C--------------------                                                  * 
C                                                                      * 
C    T(ITPM)            FINAL POINT OF ACTUALLY PERFORMED INTEGRATION  * 
C    Y(N)               FINAL VALUES AT T                              * 
C    H                  STEPSIZE PROPOSAL FOR NEXT INTEGRATION STEP    * 
C    KFLAG              ERROR FLAG                                     * 
C                       GE.0 : INTEGRATION SUCCESFULLY COMPLETED       * 
C                       LE.0 : INTEGRATION FAILED                      * 
C                       -1 : STEPSIZE TOO SMALL                        * 
C                       -2 : MORE THAN ISMAX BASIC INTEGRATION STEPS   * 
C                       -3 : NOT ENOUGH SPACE FOR TO PERFORM           * 
C                            LU-DECOMPOSITION (ENLARGE WTMP,ILEQ)      * 
C                       -4 : NEGATIVE INITIAL CONCENTRATION INDICATED  * 
C                       -5 : MORE THAN JRMAX STEPSIZE REDUCTIONS IN    * 
C                            ONE BASIC INTEGRATION STEP                * 
C                       -6 : PRECRIBED OUTPUT POINTS INVALID           * 
C                                                                      * 
C    NWTMP               MAXIMUM NUMBER OF STORAGE LOCATIONS USED FOR  * 
C                        A LU-DECOMPOSITION                            * 
C                                                                      * 
C   EXTERNAL SUBROUTINES ( SUPPLIED BY 'LARKIN' )                      * 
C-----------------------                                               * 
C                                                                      * 
C    FCN ............... RIGHT-HAND SIDE OF FIRST-ORDER                * 
C                        DIFFERENTIAL EQUATIONS                        * 
C        FCN CALLS THE SUBROUTINES :                                   * 
C        THERMO , ENTROP , KIP1     (SUPPLIED BY LARKIN)               * 
C        HEAT , DPROF , PPROF       (SUPPLIED BY USER FOR MODEL=2,4,6) * 
C                                                                      * 
C    JACOBI..............JACOBIAN OF THE DIFFERENTIAL EQUATION         * 
C                                                                      * 
C    SIMDAT..............STORES OUTPUT-INFORMATION ON DATA-SET 'OUT'   * 
C                        AT EACH INTEGRATION POINT                     * 
C                                                                      * 
C    THERMO..............COMPUTATION OF MOLAR ENTHALPY, MOLAR          * 
C                        HEAT CAPACITY AND (OPTIONAL) TEMPERATURE-     * 
C                        DERIVATIVE OF HEAT CAPACITY OF ALL SPECIES    * 
C                                                                      * 
C    ENTROP..............COMPUTATION OF MOLAR ENTROPY OF ALL SPECIES   * 
C                                                                      * 
C    KIP1................COMPUTATION OF RATE CONSTANTS OF ALL          * 
C                        REACTIONS (IF MODEL.GT.2)                     * 
C                                                                      * 
C     MA30LA, MA30LB, MA30LC, MA30LD, MA30LM:                          * 
C                        SUBROUTINES FOR SPARSE LU-DECOMPOSITION AND   * 
C                        FOR/BACKWARD SUBSTITUTION                     * 
C                                                                      * 
C                                                                      * 
C--------------------------------------------------------------------  * 
C                                                                      * 
C                                                                      * 
C   COMMON BLOCKS                                                      * 
C----------------                                                      * 
C                                                                      * 
C                                                                      *  
C                                                                      * 
C     INTEGER          ITIN,ITOUT,MOUT,EOUT,CIN,HDS,POUT,SIN,DOUT,THDAT* 
C     COMMON / LUNIT / ITIN,ITOUT,MOUT,EOUT,CIN,HDS,POUT,SIN,DOUT,THDAT* 
C                                                                      * 
C  LOGICAL UNITS ( SEE DESCRIPTION OF LARKIN DRIVING ROUTINE )         * 
C                                                                      * 
C     COMMON / COUNT / NSTEP,NFCN,NANFA,NTFAC,NFAC,NSOL                * 
C  COUNT PARAMETERS ( SEE DESCRIPTION OF SUBROUTINE META )             * 
C                                                                      * 
C     COMMON / MACHIN / EPMA , SMALLH , EXPTST                         * 
C  MACHINE DEPENDENT CONSTANTS:                                        * 
C      EPMA: RELATIVE MACHINE PRECISION                                * 
C      SMALLH: SQUARE ROOT OF SMALLES POSITIVE MACHINE NUMBER          * 
C      EXPTST: BIGGEST REAL NUMBER X ALLOWED FOR USE IN EXP(X)         * 
C                                                                      * 
C   COMMON / RATE1 / R(NCEQ)                                           * 
C  REAL ARRAY FOR COMPUTATION OF RATES                                 * 
C                                                                      * 
C   COMMON / DELHR / DHR(NCEQ)                                         * 
C  REAL ARRAY FOR COMPUTED ENTHALPY PER REACTION                       * 
C                                                                      * 
C   COMMON / MA30LE / LP,BD1,BD2,BD3                                   * 
C   COMMON / MA30LF / IRNCP,ICNCP,IRANK,MINIRN,MINICN                  * 
C   COMMON / MA30LG / DEPS,RMIN                                        * 
C  FOR COMMUNICATION WITH SPARSE SOLVER                                * 
C                                                                      * 
C  THE FOLLOWING COMMON - BLOCKS ARE ONLY USED IN SUBROUTINES CALLED   * 
C  BY META1. (DEFINED HERE FOR SEAK OF COMPLETENESS)                   * 
C                                                                      * 
C   COMMON / KINPAR / RK(NCEQ)                                         * 
C  REAL ARRAY FOR COMPUTED REACTION CONSTANTS                          * 
C                                                                      * 
C   COMMON / LARK1 / LR(NLR)                                           * 
C   COMMON / LARK2 / PLR(2*NCEQ)                                       * 
C  INTEGER ARRAYS FOR POINTER STRUCTURE OF ODE'S.                      * 
C                                                                      * 
C   COMMON / KINPAR / RK(NCEQ)                                         * 
C  REAL ARRAY FOR COMPUTED REACTION CONSTANTS                          * 
C                                                                      * 
C*********************************************************************** 
C                                                                      * 
C  INTERNAL REAL WORK ARRAYS:                                          * 
C----------------------------                                          * 
C                                                                      * 
C    ENTH(12*(I-1)+1,...,12*I)    ENTHALPY COEFFICIENTS OF THE I-TH    * 
C                                 CHEMICAL SPECIES.                    * 
C    AEA(3*NCEQ)                  KINETIC PARAMETERS (SEE KIP1)        * 
C    SK0(2*NSP)                   STANDARD MOLAR ENTROPIES FOR SPECIES * 
C                                 FOR HIGH AND LOW TEMPERATURES.       * 
C    HCP(1,...,NSP)               MOLAR ENTHALPIES OF SPECIES          * 
C    HCP(NSP+1,...,2*NSP)         MOLAR HEAT CAPACITIES OF SPECIES     * 
C    HCP(2*NSP+1,...,3*NSP)       TEMP. DERIVATIVE OF HEAT CAPACITIES  * 
C    SK(NSP)                      MOLAR ENTROPIES OF SPECIES           * 
C    BV0(NSP)                     VARIABLE LEFT HAND MATRIX OF THE ODE * 
C                                 SYSTEM FOR CONSTANT PRESSURE MODEL   * 
C                                 AT TIME POINT T0.                    * 
C    BVK(NSP)                     LIKE BV0, BUT AT TIME POINT TK.      * 
C    SH(N)                        SCALING DIAGONAL MATRIX.             * 
C    SM(N)                        MAXIMUM VALUES OF COMPONENTS         * 
C    COLLI(5)                     VALUES OF THIRD BODIES               * 
C    DT(N,JM)                     EXTRAPOLATION TABLEAU                * 
C    A(LDIM)                      ROW ORIENTED SPARSE MATRIX           * 
C    COE(NCTB)                    EFFICIENCIES OF THIRD BODIES         * 
C                                                                      * 
C  INTERNAL INTEGER WORK ARRAYS:                                       * 
C-------------------------------                                       * 
C                                                                      * 
C    ICOE(NCTB)                   POINTER TO SPECIES, ASSOCIATED TO COE* 
C    IPTD(NCO)                    POINTER TO SPECIES, FOR WHICH THERMO-* 
C                                 DYNAMICAL DATA IS SUPPLIED BY USER.  * 
C    ICOLLI(NCEQ)                 ICOLLI(K) = NUMBER OF THIRD BODY     * 
C                                 USED IN REACTION K.                  * 
C    KINEV(NCEQ)                  SEE SUBROUTINE KIP1.                 * 
C    JA(LDIM)                     COLUMN POINTER NON ZERO ENTRIES OF A.* 
C    IA(N+1)                      IA(L) POINTS TO THE FIRST NON-ZERO   * 
C                                 ELEMENT IN ROW L OF SPARSE MATRIX A. * 
C    IDA(N)                       IDENTIFIES DIAGONAL ENTRIES OF A.    * 
C                                                                      * 
C*********************************************************************** 
C                                                                        
C----------------------------------------------------------------------- 
C  DECLARATIONS                                                          
C----------------------------------------------------------------------- 
C 
      DOUBLE PRECISION Y(N),YM(N),DELQ(N),DZ(N),DEL(N)                   
CSP       REAL Y(N),YM(N),DELQ(N),DZ(N),DEL(N)                           
      DOUBLE PRECISION ETA(360),SK(NSP)                                  
CSP       REAL ETA(360),SK(NSP)                                          
      DOUBLE PRECISION BV0(NB1),BVK(NB1),SH(N),SM(N)                     
CSP       REAL BV0(NB1),BVK(NB1),SH(N),SM(N)                             
      DOUBLE PRECISION DT(N,JM),HCP(NHCP)                                
CSP       REAL DT(N,JM),HCP(NHCP)                                        
      DOUBLE PRECISION COE(NCTB),AEA(NAE),ENTH(NEN),SK0(NSK)             
CSP       REAL COE(NCTB),AEA(NAE),ENTH(NEN),SK0(NSK)                     
      DOUBLE PRECISION A(LDIM),WTMP(NWTMP)                               
CSP       REAL A(LDIM),WTMP(NWTMP)                                       
C                                                                        
C FOR HANDLING OF MORE THAN 32000 NONZEROS WITH INTEGER*2 VERSION        
C CHANGE ARRAY IWO TO STANDARD INTEGER                                   
C                                                                        
CI4       INTEGER IS(N),IKEEP(M4),IW1(M6),ILEQ(NLEQ),IWO(M2),IPTD(NCO)   
      INTEGER*2 IS(N),IKEEP(M4),IW1(M6),ILEQ(NLEQ),IWO(M2),IPTD(NCO)     
CI4       INTEGER ICOE(NCTB),ICOLLI(NCEQ),KINEV(NCEQ),LR,PLR             
      INTEGER*2 ICOE(NCTB),ICOLLI(NCEQ),KINEV(NCEQ),LR,PLR               
CI4       INTEGER JA(LDIM),IA(NIA),IDA(N)                                
      INTEGER*2 JA(LDIM),IA(NIA),IDA(N)                                  
C                                                                        
      INTEGER IDISP(2),LENTB(5)                                          
      INTEGER I, IDI, IDISP2, IELBO2, IELBOW, IERRMX, IEST, IFCNJ,       
     &  IFLAG, IH, II, IL1, IL2, IMAX, IMIN, IR2, ISMAX, ISTEP, ITP, IX,  
     &  J, J1, JK, JO, JOH, JRED, JRMAX, K, KFAIL, KH1, KINRED, KM, KO,   
     &  KOH, KONV, L, LICN, LICNP1, LICNR, LIRN, LRI, M, MFILL, N2, N3,   
     &  N4, N5, NFREE, NFREEH, NP1, NP2, NP3, NP4, NP5, NSP1, NSTC,       
     &  NZAMAX, NZV                                                       
C                                                                        
      INTEGER NJ(7),INCR(7) , NRED(6)                                    
C             NJ(JM),INCR(JM),NRED(KM)                                   
      REAL R,DHR                                                         
C                                                                        
      DOUBLE PRECISION TP(ITPM),  D(7,7), AJ(7), AL(7,7) , COLLI(5)      
CSP       REAL TP(ITPM), D(7,7),  AJ(7),  AL(7,7) , COLLI(5)             
C                        D(JM,JM),AJ(JM),AL(JM,JM)                       
C                                                                        
      DOUBLE PRECISION WZ,WY,SS,EPMIN,EPDIFF,ETADIF,ETAD                 
CSP       REAL WZ,WY,SS,EPMIN,EPDIFF,ETADIF,ETAD                         
      DOUBLE PRECISION GAS1,GAS2                                         
CSP       REAL GAS1,GAS2                                                 
      DOUBLE PRECISION TIMESC                                            
CSP       REAL TIMESC                                                    
      DOUBLE PRECISION B1,BPR1,BPR2,C,CONTRA,CON1,CON1H,CON2,CON2H       
CSP       REAL B1,BPR1,BPR2,C,CONTRA,CON1,CON1H,CON2,CON2H               
      DOUBLE PRECISION COSTF,COSTJ,COSTLR,COSTS,DABS,DENSY,DDOT          
CSP       REAL COSTF,COSTJ,COSTLR,COSTS,ABS,DENSY,DDOT                   
      DOUBLE PRECISION DEPS,DBLE,DM,DMA,DMH,DELN,DELQN,DUM1              
CSP       REAL DEPS,FLOAT,DM,DMA,DMH,DELN,DELQN,DUM1                     
      DOUBLE PRECISION DSQRT,EPH,EPMA,EPMACH,EMIN,EPS,ERR,FC,FCM,FCO     
CSP       REAL SQRT,EPH,EPMA,EPMACH,EMIN,EPS,ERR,FC,FCM,FCO              
      DOUBLE PRECISION FJ,FJ1,FMIN,FN,FRATIO                             
CSP       REAL FJ,FJ1,FMIN,FN,FRATIO                                     
      DOUBLE PRECISION G,GD,GH,H,HALF,HMAX,HMAXU                         
CSP       REAL G,GD,GH,H,HALF,HMAX,HMAXU                                 
      DOUBLE PRECISION HMIN,HI,HR,H1,OMJ,OMJO                            
CSP       REAL HMIN,HI,HR,H1,OMJ,OMJO                                    
      DOUBLE PRECISION ONE,ONE1,ONE2,PRESS,QUART                         
CSP       REAL ONE,ONE1,ONE2,PRESS,QUART                                 
      DOUBLE PRECISION RED,RK,RMIN,RO,ROOT                               
CSP       REAL RED,RK,RMIN,RO,ROOT                                       
      DOUBLE PRECISION SAFE,SAFEDM,SCALET,SHH,SHI                        
CSP       REAL SAFE,SAFEDM,SCALET,SHH,SHI                                
      DOUBLE PRECISION SCCP,SUMC,SUMDC,SB,SBY,TEMP,TEMPH                 
CSP       REAL SCCP,SUMC,SUMDC,SB,SBY,TEMP,TEMPH                         
      DOUBLE PRECISION SMALL,SMALLH,T,TA,TEN,TENTH,TEND                  
CSP       REAL SMALL,SMALLH,T,TA,TEN,TENTH,TEND                          
      DOUBLE PRECISION TGROW,TH,RMAX                                     
CSP       REAL TEND,TGROW,TH,RMAX                                        
      DOUBLE PRECISION THRESH,TN,TOLD,TWO,U,V,VOLUME,W,XKAMIN,XKAACT     
CSP       REAL THRESH,TN,TOLD,TWO,U,V,VOLUME,W,XKAMIN,XKAACT             
      DOUBLE PRECISION XMASS,YMAX,YMAXH,ZERO                             
CSP       REAL XMASS,YMAX,YMAXH,ZERO                                     
C                                                                        
      LOGICAL PRERR,PRM,PRM1,PRM2,WINDOW,LPRINT,LEUL,LMPR                
      LOGICAL BD1,BD2,BD3,BANFA                                          
      LOGICAL LENTRO,LHCP,LDHPR                                          
C                                                                        
      INTEGER          ITIN,ITOUT,MOUT,EOUT,CIN,HDS,POUT,SIN,DOUT,THDAT  
      COMMON / LUNIT / ITIN,ITOUT,MOUT,EOUT,CIN,HDS,POUT,SIN,DOUT,THDAT  
C                                                                        
      COMMON / MACHIN / EPMA,SMALLH,DUM1                                 
      COMMON / COUNT  / NSTEP,NFCN,NANFA,NTFAC,NFAC,NSOL,NTHERM,NKIP     
C                                                                        
      COMMON / LARK1  / LR(5000)                                         
      COMMON / LARK2  / PLR(2000)                                        
C                                                                        
      COMMON / MA30LE / LP,BD1,BD2,BD3                                   
      COMMON / MA30LF / IRNCP,ICNCP,IRANK,MINIRN,MINICN                  
      COMMON / MA30LG / DEPS,RMIN                                        
C                                                                        
      COMMON / RATE1  / R(1000)                                          
      COMMON / DELHR  / DHR(1000)                                        
      COMMON / KINPAR / RK(1000)                                         
C                                                                        
      COMMON / RECOND / TEMP,DENSY,PRESS,VOLUME,XMASS  
      SAVE /LUNIT/, /MACHIN/, /COUNT/, /LARK1/, /LARK2/, 
     & /MA30LE/, /MA30LF/, /MA30LG/, /RATE1/, /DELHR/, /KINPAR/, 
     & /RECOND/
C                                                                        
      DATA  ZERO/0.D0/ , FMIN/1.D-2/ , QUART/0.25D0/ , HALF/0.5D0/       
CSP       DATA  ZERO/0.E0/ , FMIN/1.E-2/ , QUART/0.25E0/ , HALF/0.5E0/   
      DATA  SAFE/0.5D0/ , ONE/1.D0/ , ONE1/1.01D0/                       
CSP       DATA  SAFE/0.5E0/ , ONE/1.E0/ , ONE1/1.01E0/                   
      DATA ONE2 /1.1D0/ , TWO/2.D0/ , TEN/1.D1/ , TENTH/1.D-1/           
CSP       DATA ONE2/1.1E0/ ,TWO/2.E0/ , TEN/1.E1/ , TENTH/1.E-1/         
      DATA SAFEDM/0.8D0/ , RO/0.25D0/ , ETADIF/1.D-6/                    
CSP       DATA SAFEDM/0.8E0/ , RO/0.25E0/ , ETADIF/1.E-6/                
      DATA EPMIN/1.D-10/ , RMAX/0.75D0/ , TIMESC/3.86D4/                 
CSP       DATA EPMIN/1.E-10/ , RMAX/0.75E0/ , TIMESC/3.86E4/             
C                                                                        
C  DEFINE LOGICAL UNITS FOR OUTPUT                                       
C  ITUT: ERROR/WARNING MESSAGES                                          
      EOUT=ITOUT                                                        
C  MOUT: INTEGRATION MONITOR                                             
      MOUT=ITOUT                                                        
C  LP: SPARSE SOLVER MESSAGES                                            
      LP=EOUT                                                           
C                                                                        
C  MACHINE DEPENDENT NUMBERS PROVIDED BY COMMON BLOCK MACHIN             
      EPMACH=EPMA*TEN                                                   
      SMALL=SMALLH*TEN**4                                               
      EPDIFF=DSQRT(EPMA)                                                
CSP        EPDIFF=SQRT(EPMA)                                             
C                                                                        
C----------------------------------------------------------------------- 
C  SET INTERNAL PARAMETERS                                               
C----------------------------------------------------------------------- 
C  (STANDARD VALUES FIXED BELOW                                          
C   TO BE ALTERED, IF NECESSARY, BY THE SKILLFUL USER)                   
C                                                                        
C  MAXIMUM PERMITTED NUMBER OF INTEGRATION STEPS PER INTERVAL            
      ISMAX=1000                                                         
C                                                                        
C  MAXIMUM PERMITTED NUMBER OF STEPSIZE REDUCTIONS PER STEP              
      JRMAX=12                                                           
C                                                                        
C  PARAMETERS FOR SPARSE LINEAR EQUATION SOLVER                          
C  THRESHOLD VALUES FOR CONDITIONAL PIVOTING:                            
C  FOR ANALYSE/FACTOR                                                    
      THRESH=1.D-2                                                       
CSP       THRESH=1.E-2                                                   
C                                                                        
      IF(EPMACH.GT.1.D-10) THRESH=THRESH*5.D0                            
CSP       IF(EPMACH.GT.1.E-10) THRESH=THRESH*5.                          
      IF(EPMACH.GT.1.D-5) THRESH=0.5D0                                   
CSP       IF(EPMACH.GT.1.E-5) THRESH=0.5                                 
C                                                                        
C  FOR WORST PIVOT RATIO IN FACTOR                                       
      TGROW=1.D-6                                                        
CSP       TGROW=1.E-6                                                    
C                                                                        
      IF(EPMACH.GT.1.D-12) TGROW=1.D-5                                   
CSP       IF(EPMACH.GT.1.E-12) TGROW=1.E-5                               
      IF(EPMACH.GT.1.D-8) TGROW=1.D-3                                    
CSP       IF(EPMACH.GT.1.E-8) TGROW=1.E-3                                
      IF(EPMACH.GT.1.D-5) TGROW=1.D-1                                    
CSP       IF(EPMACH.GT.1.E-5) TGROW=1.E-1                                
C                                                                        
C  COUNT NONZEROS IN LU DECOMPOSITION                                    
      NZAMAX=0                                                           
C                                                                        
C----------------------------------------------------------------------- 
C  INITIAL PREPARATIONS TO START THE COMPUTATION                         
C----------------------------------------------------------------------- 
C  
      IFLAG=0
      NZV=0                                                              
      IF(MODEL.GE.5) NZV=NSP                                             
      TEMPH=ZERO                                                         
      IF(MODEL.GT.2) TEMP=Y(NSP+1)                                       
      DO 1 I=1,5                                                         
1     COLLI(I)=ZERO                                                      
C  SET MOLAR ENTHALPIES,HEAT CAPACITIES AND ENTROPIES TO ZERO            
      DO 2 I=1,NHCP                                                      
2     HCP(I)=ZERO                                                        
      DO 3 I=1,NSP                                                       
3     SK(I)=ZERO                                                         
C  CHECK, IF MOLAR ENTROPIES ARE NEEDED IN THE COURSE OF THE INTEGRATION 
      LENTRO=.FALSE.                                                     
      DO 4 I=1,NCEQ                                                      
4     IF(KINEV(I).EQ.5) LENTRO=.TRUE.                                    
C  RECOMPUTATION OF MOLAR ENTH. AND HEAT CAP. FOR LHCP=.TRUE.            
      LHCP=.TRUE.                                                        
C                                                                        
C  CHOOSE DISCRETISATION                                                 
      LEUL=.FALSE.                                                       
      LMPR=.FALSE.                                                       
      IF(IDISC.EQ.0) LEUL=.TRUE.                                         
      IF(IDISC.EQ.1) LMPR=.TRUE.                                         
C                                                                        
C  SET BOOLIAN VALUES FOR INTEGRATOR OUTPUT                              
      PRERR=.FALSE.                                                      
      PRM=.FALSE.                                                        
      PRM1=.FALSE.                                                       
      PRM2=.FALSE.                                                       
      IF(KFLAG.GE.0) PRERR=.TRUE.                                        
      IF(KFLAG.EQ.1) PRM1=.TRUE.                                         
      IF(KFLAG.GE.2) PRM2=.TRUE.                                         
      IF(PRM1.OR.PRM2) PRM=.TRUE.                                        
      IF(.NOT.PRERR) LP=0                                                
C                                                                        
C  FOR PRINT OUTPUT OF ENTHALPY PER REACTION                             
      LDHPR=.TRUE.                                                       
      IF(IPRINT.EQ.4) LDHPR=.FALSE.                                      
      IF(IPRINT.LT.2) LDHPR=.FALSE.                                      
      IF(MODEL.LE.2) LDHPR=.FALSE.                                       
C                                                                        
      KFLAG=0                                                            
      IF(PRM) WRITE(MOUT,9000)                                           
      IF(PRM .AND. LEUL) WRITE(MOUT,8998)                                
      IF(PRM .AND. LMPR) WRITE(MOUT,8999)                                
C                                                                        
C  COMPUTE POINTERS                                                      
      ISTEP=0                                                            
      NP1=N+1                                                            
      N2=N+N                                                             
      NP2=N2+1                                                           
      N3=N2+N                                                            
      NP3=N3+1                                                           
      N4=N3+N                                                            
      NP4=N4+1                                                           
      N5=N4+N                                                            
      NP5=N5+1                                                           
C                                                                        
C  PREPARATIONS FOR SPARSE SOLVER                                        
      IELBOW=N                                                           
      XKAMIN=1.3D0                                                       
CSP       XKAMIN=1.3                                                     
      LIRN=LDIM+IELBOW                                                   
      LICN=NLEQ-LIRN                                                     
CI4       CONTINUE                                                       
      IF(LICN.GT.32000) LICN=32000                                       
      LICNR=NWTMP-N                                                      
      NFREE=LICN-LICNR                                                   
      IF(NFREE.GT.N2.AND.PRM2) WRITE(MOUT,9300) NFREE                    
      NFREE=LICNR-LICN                                                   
      IF(NFREE.GT.N2.AND.PRM2) WRITE(MOUT,9305) NFREE                    
      IF(LICNR.LT.LICN) LICN=LICNR                                       
C  CHECK FOR SUFFICIENT STORAGE PLACE                                    
      XKAACT=DBLE(LICN)/DBLE(LDIM) 
CSP       XKAACT=FLOAT(LICN)/FLOAT(LDIM)                                 
      IF(XKAACT.LT.XKAMIN) GOTO 299                                      
      LICNP1=LICN+1                                                      
      IDISP(1)=1                                                         
      BANFA=.TRUE.                                                       
      DEPS=TWO                                                           
C                                                                        
C----------------------------------------------------------------------- 
C  PREPARATIONS FOR ORDER AND STEPSIZE CONTROL                           
C----------------------------------------------------------------------- 
C                                                                        
C  ASSOCIATED MAXIMUM COLUMN NUMBER (1.LE.KM.LE.6)                       
      KM=4                                                               
      IF(IDISC.EQ.0) KM=6                                                
C                                                                        
C  ASSOCIATED MAXIMUM ROW NUMBER (2.LE.JM.LE.7)                          
      JM=KM+1                                                            
C                                                                        
C  STEPSIZE SEQUENCE A (DUE TO /2/)                                      
      DO 1000 J=1,JM                                                     
1000  NJ(J)=J                                                            
      IF(LEUL) GOTO 1020                                                 
      DO 1010 I=1,JM                                                     
      IF(I.EQ.1) J=2                                                     
      IF(I.EQ.2) J=6                                                     
      IF(I.EQ.3) J=10                                                    
      IF(I.EQ.4) J=14                                                    
      IF(I.EQ.5) J=22                                                    
      IF(I.EQ.6) J=34                                                    
      IF(I.EQ.7) J=50                                                    
1010  NJ(I)=J                                                            
C                                                                        
1020  FN=DBLE(N)                                                         
CSP 1020  FN=FLOAT(N)                                                    
      FJ1=DBLE(NJ(1))                                                    
CSP       FJ1=FLOAT(NJ(1))                                               
C                                                                        
C  INFORMATION COEFFICIENTS                                              
      EPH=RO*EPS                                                         
      AJ(1)=FJ1+ONE                                                      
      DO 1050 J=2,JM                                                     
      J1=J-1                                                             
      INCR(J1)=0                                                         
      NRED(J1)=0                                                         
      FJ=DBLE(NJ(J))                                                     
CSP       FJ=FLOAT(NJ(J))                                                
      V=AJ(J1)+FJ                                                        
      AJ(J)=V                                                            
      DO 1030 K=1,J1                                                     
      W=FJ/DBLE(NJ(K))                                                   
CSP       W=FJ/FLOAT(NJ(K))                                              
      IF(LMPR) W=W*W                                                     
1030  D(J,K)=W                                                           
      IF(J.EQ.2) GOTO 1050                                               
      W=V-FJ1                                                            
      DO 1040 KH1=2,J1                                                   
      K=KH1-1                                                            
      IF(LEUL) U=(AJ(KH1)-V)/(W*DBLE(KH1))                               
CSP       IF(LEUL) U=(AJ(KH1)-V)/(W*FLOAT(KH1))                          
      IF(LMPR) U=(AJ(KH1)-V)/(W*DBLE(K+KH1))                             
CSP       IF(LMPR) U=(AJ(KH1)-V)/(W*FLOAT(K+KH1))                        
      U=EPH**U                                                           
1040  AL(J1,K)=U                                                         
1050  CONTINUE                                                           
C                                                                        
C  WORK COEFFICIENTS                                                     
      COSTF=ONE                                                          
      COSTJ=10.*COSTF                                                    
      COSTS=ONE                                                          
      COSTLR=4.*COSTS                                                    
      IF((COSTS+COSTLR+COSTJ).EQ.ZERO) GOTO 1070                         
      AJ(1)=COSTJ+COSTLR+(COSTF+COSTS)*(FJ1+ONE)                         
      DO 1060 J=2,JM                                                     
1060  AJ(J)=AJ(J-1)+(COSTF+COSTS)*DBLE(NJ(J))+COSTS+COSTLR               
CSP 1060  AJ(J)=AJ(J-1)+(COSTF+COSTS)*FLOAT(NJ(J))+COSTS+COSTLR          
1070  KOH=1                                                              
      JOH=2                                                              
1080  IF(JOH.GE.JM) GOTO 1090                                            
      IF(AJ(JOH+1)*ONE1.GT.AJ(JOH)*AL(JOH,KOH)) GOTO 1090                
      KOH=JOH                                                            
      JOH=JOH+1                                                          
      GOTO 1080                                                          
1090  K=0                                                                
C                                                                        
C  ASSOCIATED OPTIMAL MAXIMUM COLUMN                                     
      KM=KOH                                                             
      JM=JOH                                                             
      INCR(JM)=-1                                                        
      OMJO=ZERO                                                          
      IF(PRM1) WRITE(MOUT,9100) EPS                                      
      IF(PRM2) WRITE(MOUT,9200) EPS                                      
      IF(PRM2) WRITE(MOUT,9201) KM                                       
      IF(PRM2) WRITE(MOUT,9202) CON1                                     
      IF(PRM2) WRITE(MOUT,9203) CON2                                     
      IF(HMAX.LT.TP(ITPM)-TP(1).AND.PRM) WRITE(MOUT,9001) HMAX           
C                                                                        
C  MINIMUM PERMITTED STEPSIZE                                            
      HMAXU=HMAX                                                         
      HMIN=TWO*(TP(ITPM)-TP(1))/DBLE(ISMAX)                              
CSP       HMIN=TWO*(TP(ITPM)-TP(1))/FLOAT(ISMAX)                         
      FCM=ONE                                                            
      IF(HMAXU.GT.HMIN) GOTO 1100                                        
      HMAXU=HMIN                                                         
      IF(PRERR) WRITE(EOUT,9002) HMAXU                                   
C                                                                        
C  MINIMUM PERMITTED RELATIVE PRECISION                                  
1100  EMIN=EPMACH*1.D3                                                   
CSP 1100  EMIN=EPMACH*1.E3                                               
      IF(EPS.GT.EMIN) GOTO 1110                                          
      EPS=EMIN                                                           
      IF(PRERR) WRITE(EOUT,9004) EPS                                     
C                                                                        
C  MINIMUM PERMITTED VALUES FOR SCALING PARAMETER                        
1110  IF(CON1.GT.SMALL) GOTO 1120                                        
      CON1=SMALL                                                         
      IF(PRERR) WRITE(EOUT,9006) CON1                                    
1120  IF(CON2.GT.EPMACH) GOTO 1130                                       
      CON2=EPMACH                                                        
      IF(PRERR) WRITE(EOUT,9008) CON2                                    
1130  CON1H=CON1                                                         
      CON2H=CON2                                                         
      IF(PRM1) WRITE(MOUT,9110)                                          
      IF(PRM2) WRITE(MOUT,9210)                                          
C                                                                        
C----------------------------------------------------------------------- 
C  INITIAL SCALING                                                       
C----------------------------------------------------------------------- 
      YMAX=ZERO                                                          
      IF(N.GT.360) GOTO 899                                              
      DO 1160 I=1,N                                                      
      ETA(I)=ETADIF                                                      
      DO 1150 J=1,JM                                                     
1150  DT(I,J)=ZERO                                                       
1160  IF(Y(I).GT.YMAX .AND. I.LE.NSP) YMAX=Y(I)                          
      IF(YMAX.EQ.ZERO) YMAX=ONE                                          
      YMAXH=YMAX*CON1                                                    
      SCALET=CON1*YMAX                                                   
      IF(SCALET.GT.SMALL) SCALET=SMALL                                   
      DO 1180 I=1,NSP                                                    
      U=Y(I)                                                             
      IF(U.LT.ZERO) GOTO 399                                             
      IF(U.GE.SCALET) GOTO 1170                                          
      SM(I)=YMAX                                                         
      SH(I)=YMAXH                                                        
      IS(I)=0                                                            
      GOTO 1180                                                          
1170  SM(I)=U                                                            
      SH(I)=U                                                            
      IS(I)=1                                                            
1180  CONTINUE                                                           
      IF(NSP.EQ.N) GOTO 1190                                             
      NSP1=NSP+1                                                         
      SM(NSP1)=Y(NSP1)                                                   
      SH(NSP1)=Y(NSP1)                                                   
      IS(NSP1)=1                                                         
      IF(N-NSP.LT.2) GOTO 1190                                           
      SM(N)=Y(N)                                                         
      SH(N)=Y(N)                                                         
      IS(N)=1                                                            
1190  CONTINUE                                                           
C                                                                        
C----------------------------------------------------------------------- 
C  INITIAL STEPSIZE ( IF HI.LE.ZERO, COMPUTE HI INTERNALLY)              
C----------------------------------------------------------------------- 
      IFCNJ=0                                                            
      IF(HI.GT.ZERO) GOTO 1270                                           
      IFCNJ=1                                                            
      DM=ZERO                                                            
      T=TP(1)                                                            
      KFAIL=1                                                            
C     CALL MONON(1)                                                      
      CALL FCNINT (N,NSP,NCEQ,NTB,NHCP,T,Y,DZ,HCP,SCCP,SUMC,SUMDC,GAS1,  
     & GAS2,COLLI,MODEL,BV0,EPS,NB1,NCTB,COE,ICOE,SK,LENTB,KFAIL,ICOLLI, 
     & KINEV,AEA,ENTH,SK0,NAE,NEN,NSK,NCO,IPTD,TEMPH,LENTRO,LHCP,DDOT)   
C     CALL MONOFF(1)                                                     
      NFCN=NFCN+1                                                        
      IF(KFAIL.NE.0) GOTO 699                                            
C     CALL MONON(2)                                                      
      CALL JACOB1 (N,NSP,NCEQ,NHCP,LDIM,TP(1),Y,A,DZ,HCP,SCCP,SUMC,      
     & SUMDC,GAS1,GAS2,COLLI,NCTB,COE,ICOE,LENTB,MODEL,ICOLLI,KINEV,JA,  
     & IA,IDA,AEA,NIA,NAE,DDOT)                                          
C     CALL MONOFF(2)                                                     
      IMAX=0                                                             
      DO 1260 I=1,N                                                      
      SHI=ONE/SH(I)                                                      
      IMIN=IMAX+1                                                        
      IMAX=IA(I+1)-1                                                     
      V=ZERO                                                             
      DO 1250 L=IMIN,IMAX                                                
      IH=JA(L)                                                           
      U=SH(IH)*A(L)*SHI                                                  
      V=V+DABS(U)                                                        
CSP       V=V+ABS(U)                                                     
1250  CONTINUE                                                           
1260  IF(V.GT.DM) DM=V                                                   
      IF(MODEL.NE.5) GOTO 1269                                           
      SB=ZERO                                                            
      DO 1261 I=1,NZV                                                    
      U=DABS(BV0(I)*SH(N)/SH(I))                                         
CSP       U=ABS(BV0(I)*SH(N)/SH(I))                                      
1261  IF(U.GT.SB) SB=U                                                   
      SB=SB+ONE                                                          
      SBY=ZERO                                                           
      BPR1=DZ(N)/Y(N)                                                    
      BPR2=SH(N)*BPR1/Y(N)                                               
      DO 1262 I=1,NZV                                                    
      U=DABS(BPR2*Y(I)/SH(I))                                            
CSP       U=ABS(BPR2*Y(I)/SH(I))                                         
1262  IF(U.GT.SBY) SBY=U                                                 
      SBY=SBY+BPR1                                                       
      IF(ONE.GT.SBY) SBY=ONE                                             
      DM=SB*(DM+SBY)                                                     
1269  IF(DM.EQ.ZERO) DM=ONE                                              
      HI=HALF/DM                                                         
C                                                                        
1270  H=HI                                                               
C                                                                        
      ITP=2                                                              
      T=TP(1)                                                            
      TEND=TP(2)                                                         
      TOLD=T                                                             
      H1=TEND-T                                                          
      WINDOW=.FALSE.                                                     
      LPRINT=.TRUE.                                                      
C                                                                        
C*********************************************************************** 
C                NEXT BASIC INTEGRATION STEP                             
C*********************************************************************** 
C                                                                        
2000  CONTINUE                                                           
C      LHCP=.TRUE.                                                       
C      LENTRO=.TRUE.                                                     
      IF(IFCNJ.EQ.1) GOTO 2030                                           
      KFAIL=1                                                            
C     CALL MONON(1)                                                      
      CALL FCNINT (N,NSP,NCEQ,NTB,NHCP,T,Y,DZ,HCP,SCCP,SUMC,SUMDC,GAS1,  
     & GAS2,COLLI,MODEL,BV0,EPS,NB1,NCTB,COE,ICOE,SK,LENTB,KFAIL,ICOLLI, 
     & KINEV,AEA,ENTH,SK0,NAE,NEN,NSK,NCO,IPTD,TEMPH,LENTRO,LHCP,DDOT)   
C     CALL MONOFF(1)                                                     
      NFCN=NFCN+1                                                        
2030  CONTINUE                                                           
C      LHCP=.FALSE.                                                      
C      LENTRO=.FALSE.                                                    
C                                                                        
C  COMPUTE ENTHALPY PER REACTION IF LDHPR=.TRUE.                         
C  (FOR LATER TEST PURPOSES ONLY.)                                       
      IF(.NOT.LDHPR) GOTO 2020                                           
      IR2=0                                                              
      DO 2015 L=1,NCEQ                                                   
      IL1=IR2+1                                                          
      IL2=PLR(L)                                                         
      IR2=PLR(NCEQ+L)                                                    
      DHR(L)=ZERO                                                        
      DO 2010 I=IL1,IR2                                                  
      LRI=LR(I)                                                          
      IF(I.GT.IL2) GOTO 2005                                             
      DHR(L)=DHR(L) - HCP(LRI)                                           
      GOTO 2010                                                          
2005  DHR(L)=DHR(L) + HCP(LRI)                                           
2010  CONTINUE                                                           
2015  CONTINUE                                                           
C                                                                        
2020  IF(PRM1) WRITE(MOUT,9120) ISTEP,T                                  
      IF(PRM2) WRITE(MOUT,9220) ISTEP,NFCN,T,H,K,KOH                     
      IF(H1.GT.T*EPMACH) GOTO 2090                                       
C  HANDLING OF OUTPUT-POINTS TP(ITP)                                     
      LPRINT=.TRUE.                                                      
      H=HR                                                               
      HMAX=HMAXU                                                         
      IF(ITP.GE.ITPM) GOTO 7000                                          
C                                                                        
2070  ITP=ITP+1                                                          
      TEND=TP(ITP)                                                       
      H1=TEND-T                                                          
      IF(H1.GT.ZERO) GOTO 2080                                           
C  PRESCRIBED OUTPUT POINTS INVALID                                      
      IF(ITP.LT.ITPM) GOTO 2070                                          
      GOTO 599                                                           
2080  WINDOW=.FALSE.                                                     
C                                                                        
2090  CONTINUE                                                           
      CALL SIMDAT (N,NSP,NCEQ,MODEL,Y,GAS1,T,IPRINT,LPRINT)              
C                                                                        
      LPRINT=.FALSE.                                                     
      IF(H1.GE.ONE2*H) GOTO 2100                                         
      HR=H                                                               
      H=H1                                                               
2100  JRED=0                                                             
      NSTC=0                                                             
      KINRED=0                                                           
      DMH=SAFEDM                                                         
      DO 2110 K=1,KM                                                     
2110  INCR(K)=INCR(K)+1                                                  
      HMAX=H1                                                            
      IF(HMAXU.LT.HMAX) HMAX=HMAXU                                       
C                                                                        
C  ANALYTIC JACOBIAN                                                     
      IF(IFCNJ.EQ.1) GOTO 2150                                           
C     CALL MONON(2)                                                      
      CALL JACOB1 (N,NSP,NCEQ,NHCP,LDIM,T,Y,A,DZ,HCP,SCCP,SUMC,SUMDC,    
     & GAS1,GAS2,COLLI,NCTB,COE,ICOE,LENTB,MODEL,ICOLLI,KINEV,JA,IA,IDA, 
     & AEA,NIA,NAE,DDOT)                                                 
C     CALL MONOFF(2)                                                     
2150  IFCNJ=0                                                            
C                                                                        
C  ADDITION OF MATRIX BPRIME*DY (FOR MODEL.GE.5)                         
      IF(MODEL.LE.4) GOTO 2140                                           
      BPR1=DZ(N)/Y(N)                                                    
      BPR2=BPR1/Y(N)                                                     
      DO 2120 I=1,NSP                                                    
      IH=IDA(I)                                                          
2120  A(IH)=A(IH) + BPR1                                                 
      DO 2130 I=1,NSP                                                    
      IH=IA(I+1)-1                                                       
2130  A(IH)=A(IH) - BPR2*Y(I)                                            
2140  CONTINUE                                                           
C                                                                        
C  SCALING OF THE JACOBIAN                                               
      IMAX=0                                                             
      DO 2160 I=1,N                                                      
      SHI=ONE/SH(I)                                                      
      IH=IDA(I)                                                          
      IMIN=IMAX+1                                                        
      IMAX=IA(I+1)-1                                                     
      DO 2160 L=IMIN,IMAX                                                
      IH=JA(L)                                                           
2160  A(L)=-SH(IH)*A(L)*SHI                                              
C                                                                        
C  SCALING OF THE REST OF THE LINEAR SYSTEM                              
      DO 2410 I=1,N                                                      
2410  DZ(I)=DZ(I)/SH(I)                                                  
      IF(NZV.EQ.0) GOTO 2430                                             
      DO 2420 I=1,NZV                                                    
2420  BV0(I)=SH(N)*BV0(I)/SH(I)                                          
2430  CONTINUE                                                           
C                                                                        
3000  TN=T+H                                                             
      IEST=1                                                             
      FCM=H/HMAX                                                         
      IF(FCM.LT.FMIN) FCM=FMIN                                           
      DM=DMH*TWO                                                         
      GH=ZERO                                                            
C                                                                        
      DO 300 J=1,JM                                                      
C                                                                        
      M=NJ(J)                                                            
      G=H/DBLE(M)                                                        
CSP       G=H/FLOAT(M)                                                   
      GD=GH                                                              
      GH=ONE/G                                                           
      GD=GH-GD                                                           
C                                                                        
C----------------------------------------------------------------------- 
C  SEMI-IMPLICIT EULER STARTING STEP                                     
C----------------------------------------------------------------------- 
      DO 3010 I=1,N                                                      
      DEL(I)=DZ(I)                                                       
3010  YM(I)=Y(I)                                                         
C                                                                        
C  SPARSE LU-DECOMPOSION OF (B/G - A)                                    
C------------------------------------                                    
      IF(.NOT.BANFA) GOTO 3060                                           
C                                                                        
C  ANALYSE/FACTOR SPARSE MATRIX                                          
      IDISP2=LICN-LDIM                                                   
      IDISP(2)=IDISP2+1                                                  
      IMAX=0                                                             
      DO 3050 I=1,N                                                      
      IMIN=IMAX+1                                                        
      IMAX=IA(I+1)                                                       
C  SET ROW LENGTHS                                                       
      IKEEP(I)=IMAX-IMIN                                                 
      IMAX=IMAX-1                                                        
C  SET PERMUTATION TO IDENTITY                                           
      IKEEP(N+I)=I                                                       
      IKEEP(N2+I)=I                                                      
      DO 3040 II=IMIN,IMAX                                               
      IH=IDISP2+II                                                       
      ILEQ(IH)=JA(II)                                                    
      WTMP(IH)=A(II)                                                     
      IF(JA(II).EQ.I) WTMP(IH)=WTMP(IH) + GH                             
3040  CONTINUE                                                           
3050  CONTINUE                                                           
      IF(NZV.EQ.0) GOTO 3052                                             
      DO 3051 I=1,NZV                                                    
      L=IA(I+1)-1                                                        
      IH=IDISP2+L                                                        
3051  WTMP(IH)=WTMP(IH) + GH*BV0(I)                                      
3052  IDISP2=IDISP2+1                                                    
C                                                                        
C     CALL MONON(3)                                                      
      CALL MA30LA (N,ILEQ,WTMP,LICN,IKEEP,IKEEP(NP3),IDISP,IKEEP(NP1),   
     &             IKEEP(NP2),ILEQ(LICNP1),LIRN,IW1,IW1(NP1),IW1(NP2),    
     &             IW1(NP3),IW1(NP4),IW1(NP5),IWO,IWO(NP1),THRESH,IFLAG)  
C     CALL MONOFF(3)                                                     
C                                                                        
C  EMERGENCY EXIT FOR SINGULAR MATRIX AND ARRAY OVERFLOW                 
      IF(IFLAG.LT.0) GOTO 6070                                           
C                                                                        
      IELBO2=LICN-IDISP(2)                                               
C                                                                        
C  COMPUTE MAPPING                                                       
      CALL MA30LM (N,NP1,IDISP(2),ILEQ,IA,JA,LDIM,ILEQ(LICNP1),IKEEP,    
     &             IKEEP(NP1),IKEEP(NP2),IWO)                             
C                                                                        
      NANFA=NANFA+1                                                      
      IF(PRM2) WRITE(MOUT,9230)                                          
      IDISP2=IDISP(2)                                                    
      IF(IDISP2.GT.NZAMAX) NZAMAX=IDISP2                                 
C                                                                        
C  EMERGENCY EXIT FOR SINGULAR MATRIX AND ARRAY OVERFLOW 
      IF(IFLAG.LT.0) GOTO 6070                                           
C                                                                        
      BANFA=.FALSE.                                                      
      GOTO 3110                                                          
C                                                                        
C  FACTOR SPARSE MATRIX                                                  
C----------------------                                                  
3060  CONTINUE                                                           
      DO 3070 I=1,IDISP2                                                 
3070  WTMP(I)=ZERO                                                       
      DO 3080 I=1,LDIM                                                   
      II=LICN+I                                                          
      IH=ILEQ(II)                                                        
3080  WTMP(IH)=A(I)                                                      
      DO 3081 I=1,N                                                      
      IDI=IDA(I)                                                         
      IH=ILEQ(LICN+IDI)                                                  
3081  WTMP(IH)=WTMP(IH) + GH                                             
      IF(NZV.EQ.0) GOTO 3083                                             
      DO 3082 I=1,NZV                                                    
      L=IA(I+1)-1                                                        
      IH=ILEQ(LICN+L)                                                    
3082  WTMP(IH)=WTMP(IH) + GH*BV0(I)                                      
3083  IF(J.EQ.1) DEPS=HALF                                               
C                                                                        
C     IF(DEPS.EQ.TWO) THEN                                               
C       CALL MONON(5)                                                    
C     ELSE                                                               
C       CALL MONON(4)                                                    
C     ENDIF                                                              
      CALL MA30LB (N,ILEQ,WTMP,LICN,IKEEP,IKEEP(NP3),IDISP,              
     &             WTMP(LICNP1),IWO,IFLAG)                                
      IF(DEPS.EQ.TWO) THEN                                               
C       CALL MONOFF(5)                                                   
        NFAC=NFAC+1                                                      
C     ELSE                                                               
C       CALL MONOFF(4)                                                   
      ENDIF                                                              
      IF(DEPS.NE.TWO) NTFAC=NTFAC+1                                      
      DEPS=TWO
C                                                                        
C  MATRIX SINGULAR                                                       
      IF(IFLAG.LT.0) GOTO 3090                                           
C                                                                        
C  CHECK PIVOT RATIO (FOR J.EQ.1)                                        
      IF(J.NE.1.OR.RMIN.GE.TGROW) GOTO 3110                              
C                                                                        
C  NEW ANALYSE/FACTOR NECESSARY                                          
3090  BANFA=.TRUE.                                                       
      GOTO 3000                                                          
C                                                                        
C  SOLUTION OF LINEAR SYSTEM (B/G - A)*DEL = DEL                         
3110  CONTINUE                                                           
C                                                                        
C     CALL MONON(6)                                                      
      CALL MA30LC (N,ILEQ,WTMP,LICN,IKEEP,IKEEP(NP3),IKEEP(NP1),         
     &             IKEEP(NP2),DEL,WTMP(IDISP2+1))                         
C     CALL MONOFF(6)                                                     
      NSOL=NSOL+1                                                        
C                                                                        
      M=M-1                                                              
C                                                                        
400   IF(LMPR) GOTO 500                                                  
C                                                                        
C----------------------------------------------------------------------- 
C  SEMI-IMPLICIT EULER STEPS                                             
C----------------------------------------------------------------------- 
      DO 4010 I=1,N                                                      
4010  YM(I)=YM(I)+DEL(I)*SH(I)                                           
      IF(M.EQ.0) GOTO 3500                                               
      DO 4000 K=1,M                                                      
      TH=T+DBLE(K)*G                                                     
CSP       TH=T+FLOAT(K)*G                                                
      KFAIL=0                                                            
C     CALL MONON(1)                                                      
      CALL FCNINT (N,NSP,NCEQ,NTB,NHCP,TH,YM,DEL,HCP,SCCP,SUMC,SUMDC,    
     & GAS1,GAS2,COLLI,MODEL,BVK,EPS,NB1,NCTB,COE,ICOE,SK,LENTB,KFAIL,   
     & ICOLLI,KINEV,AEA,ENTH,SK0,NAE,NEN,NSK,NCO,IPTD,TEMPH,LENTRO,      
     & LHCP,DDOT)                                                        
C     CALL MONOFF(1)                                                     
      NFCN=NFCN+1                                                        
      IF(KFAIL.NE.0) GOTO 6060                                           
      DO 4020 I=1,N                                                      
4020  DEL(I)=DEL(I)/SH(I)                                                
      IF(MODEL.LE.4) GOTO 4060                                           
C  SCALING OF MATRIX BVK                                                 
      DO 4040 I=1,NZV                                                    
4040  BVK(I)=SH(N)*BVK(I)/SH(I)                                          
      DO 4050 I=1,NSP                                                    
4050  DEL(I)=DEL(I) + (BV0(I)-BVK(I))*DEL(N)                             
4060  CONTINUE                                                           
C                                                                        
C  SOLUTION OF LINEAR SYSTEM (B/G - A)*DEL = DEL                         
C     CALL MONON(6)                                                      
      CALL MA30LC (N,ILEQ,WTMP,LICN,IKEEP,IKEEP(NP3),                    
     &                 IKEEP(NP1), IKEEP(NP2),DEL,WTMP(IDISP2+1))         
C     CALL MONOFF(6)                                                     
      NSOL=NSOL+1                                                        
C                                                                        
      DO 4120 I=1,N                                                      
4120  YM(I)=YM(I)+DEL(I)*SH(I)                                           
4000   CONTINUE                                                          
C----------------------------------------------------------------------- 
C  END OF SEMI-IMPLICIT EULER STEPS                                      
C----------------------------------------------------------------------- 
C                                                                        
      GOTO 3500                                                          
C                                                                        
C----------------------------------------------------------------------- 
C  SEMI-IMPLICIT MID-POINT STEPS                                         
C----------------------------------------------------------------------- 
500   DO 5000 K=1,M                                                      
      DO 5010 I=1,N                                                      
5010  YM(I)=YM(I)+DEL(I)*SH(I)                                           
      TH=T+DBLE(K)*G                                                     
CSP       TH=T+FLOAT(K)*G                                                
      KFAIL=0                                                            
C     CALL MONON(1)                                                      
      CALL FCNINT (N,NSP,NCEQ,NTB,NHCP,TH,YM,DELQ,HCP,SCCP,SUMC,SUMDC,   
     & GAS1,GAS2,COLLI,MODEL,BVK,EPS,NB1,NCTB,COE,ICOE,SK,LENTB,KFAIL,   
     & ICOLLI,KINEV,AEA,ENTH,SK0,NAE,NEN,NSK,NCO,IPTD,TEMPH,LENTRO,      
     & LHCP,DDOT)                                                        
C     CALL MONOFF(1)                                                     
      NFCN=NFCN+1                                                        
      IF(KFAIL.NE.0) GOTO 6060                                           
      IF(MODEL.GE.5) GOTO 5030                                           
      DO 5020 I=1,N                                                      
5020  DELQ(I)=DELQ(I)/SH(I) - DEL(I)*GH                                  
      GOTO 5040                                                          
5030  CONTINUE                                                           
C  SCALING OF MATRIX BVK                                                 
      DO 5031 I=1,NZV                                                    
5031  BVK(I)=SH(N)*BVK(I)/SH(I)                                          
      DELQN=DELQ(N)/SH(N)                                                
      DELN=DEL(N)                                                        
      DO 5032 I=1,NZV                                                    
5032  DELQ(I)=DELQ(I)/SH(I) - (DEL(I)+BV0(I)*DELN)*GH                    
     &       + (BV0(I)-BVK(I))*DELQN                                     
      DELQ(NSP1)=DELQ(NSP1)/SH(NSP1) - DEL(NSP1)*GH                      
      DELQ(N)=DELQ(N)/SH(N) - DEL(N)*GH                                  
5040  CONTINUE                                                           
C                                                                        
C  SOLUTION OF LINEAR SYSTEM (B/G - A)*DELQ = DELQ                       
C     CALL MONON(6)                                                      
      CALL MA30LC (N,ILEQ,WTMP,LICN,IKEEP,IKEEP(NP3),                    
     &                 IKEEP(NP1),IKEEP(NP2),DELQ,WTMP(IDISP2+1))         
C     CALL MONOFF(6)                                                     
      NSOL=NSOL+1                                                        
C                                                                        
      DO 5110 I=1,N                                                      
5110  DEL(I)=TWO*DELQ(I)+DEL(I)                                          
5000  CONTINUE                                                           
C                                                                        
C----------------------------------------------------------------------- 
C  SMOOTHING FINAL STEP (WITH STABILITY CHECK)                           
C----------------------------------------------------------------------- 
      DO 5160 I=1,N                                                      
5160  YM(I)=YM(I)+DEL(I)*SH(I)                                           
      KFAIL=0                                                            
C     CALL MONON(1)                                                      
      CALL FCNINT (N,NSP,NCEQ,NTB,NHCP,TN,YM,DELQ,HCP,SCCP,SUMC,SUMDC,   
     & GAS1,GAS2,COLLI,MODEL,BVK,EPS,NB1,NCTB,COE,ICOE,SK,LENTB,KFAIL,   
     & ICOLLI,KINEV,AEA,ENTH,SK0,NAE,NEN,NSK,NCO,IPTD,TEMPH,LENTRO,      
     & LHCP,DDOT)                                                        
C     CALL MONOFF(1)                                                     
      NFCN=NFCN+1                                                        
      IF(KFAIL.NE.0) GOTO 6060                                           
      IF(MODEL.GE.5) GOTO 5180                                           
      DO 5170 I=1,N                                                      
5170  DEL(I)=DELQ(I)/SH(I) - DEL(I)*GH                                   
      GOTO 5190                                                          
5180  CONTINUE                                                           
C  SCALING OF MATRIX BVK                                                 
      DO 5181 I=1,NZV                                                    
5181  BVK(I)=SH(N)*BVK(I)/SH(I)                                          
      DELQN=DELQ(N)/SH(N)                                                
      DELN=DEL(N)                                                        
      DO 5182 I=1,NZV                                                    
5182  DEL(I)=DELQ(I)/SH(I) - (DEL(I)+BV0(I)*DELN)*GH                     
     &       + (BV0(I)-BVK(I))*DELQN                                     
      DEL(NSP1)=DELQ(NSP1)/SH(NSP1) - DEL(NSP1)*GH                       
      DEL(N)=DELQ(N)/SH(N) - DEL(N)*GH                                   
5190  CONTINUE                                                           
C                                                                        
C  SOLUTION OF LINEAR SYSTEM (B/G - A)*DELQ = DELQ                       
C     CALL MONON(6)                                                      
      CALL MA30LC (N,ILEQ,WTMP,LICN,IKEEP,IKEEP(NP3),                    
     &            IKEEP(NP1),IKEEP(NP2),DEL,WTMP(IDISP2+1))               
C     CALL MONOFF(6)                                                     
      NSOL=NSOL+1                                                        
C                                                                        
      DMA=DM                                                             
      DM=ZERO                                                            
      DO 5290 I=1,N                                                      
      U=DEL(I)*SH(I)                                                     
      YM(I)=YM(I)+U                                                      
      IF(J.GT.3) GOTO 5290                                               
C                                                                        
      U=DABS(U)                                                          
CSP       U=ABS(U)                                                       
      V=DABS(YM(I))                                                      
CSP       V=ABS(YM(I))                                                   
      IF(V.LT.SH(I)) V=SH(I)                                             
      U=U/V                                                              
      IF(U.GT.DM) IERRMX=I                                               
      IF(U.GT.DM) DM=U                                                   
5290  CONTINUE                                                           
C                                                                        
C STABILITY CHECK                                                        
      IF(J.GT.3) GOTO 3500                                               
      IF(DM.LT.DMA*HALF) GOTO 5300                                       
      IF(DM.LE.TENTH.AND.J.GT.1) GOTO 3500                               
C                                                                        
C  EMERGENCY EXIT IF STABILITY CHECK IS ACTIVATED                        
      NSTC=NSTC+1                                                        
      GOTO 6040                                                          
C                                                                        
5300  IF(J.EQ.1) DMH=DM*HALF                                             
      IF(DM.LT.HALF) GOTO 3500                                           
C  PREVENTION OF POSSIBLE ORDER INCREASE                                 
      DO 5310 L=JOH,JM                                                   
      IF(INCR(L).GT.0) INCR(L)=0                                         
      INCR(L)=INCR(L)-2                                                  
5310  CONTINUE                                                           
C                                                                        
C----------------------------------------------------------------------- 
C  END OF SEMI-IMPLICIT MID-POINT RULE                                   
C----------------------------------------------------------------------- 
C                                                                        
3500  CONTINUE                                                           
C                                                                        
C----------------------------------------------------------------------- 
C  POLYNOMIAL EXTRAPOLATION                                              
C----------------------------------------------------------------------- 
      ERR=ZERO                                                           
      DO 3540 I=1,N                                                      
      C=YM(I)                                                            
      V=DT(I,1)                                                          
      DT(I,1)=C                                                          
      IF(J.EQ.1) GOTO 3540                                               
      TA=C                                                               
      DO 3530 K=2,J                                                      
      JK=J-K+1                                                           
      B1=D(J,JK)                                                         
      W=C-V                                                              
      U=W/(B1-ONE)                                                       
      C=B1*U                                                             
      V=DT(I,K)                                                          
      DT(I,K)=U                                                          
3530  TA=U+TA                                                            
      YM(I)=TA                                                           
      TA=DABS(TA)                                                        
CSP       TA=ABS(TA)                                                     
      CONTRA=CON1                                                        
      IF(IS(I).GE.1) CONTRA=CON2H                                        
      SHH=SM(I)*CONTRA                                                   
      IF(TA.LT.SHH) TA=SHH                                               
      IF(TA.LT.SMALL) TA=SMALL                                           
      U=U/TA                                                             
      ERR=ERR+U*U                                                        
3540  CONTINUE                                                           
C                                                                        
      IF(J.EQ.1) GOTO 300                                                
C ERROR ( SCALED ROOT MEAN SQUARE )                                      
      ERR=DSQRT(ERR/FN)                                                  
CSP       ERR=SQRT(ERR/FN)                                               
      KONV=0                                                             
      IF(ERR.LT.EPS) KONV=1                                              
      ERR=ERR/EPH                                                        
C                                                                        
C----------------------------------------------------------------------- 
C  ORDER AND STEPSIZE CONTROL                                            
C----------------------------------------------------------------------- 
      K=J-1                                                              
      IF(LEUL) ROOT=ONE/DBLE(J)                                          
CSP       IF(LEUL) ROOT=ONE/FLOAT(J)                                     
      IX=K+J                                                             
      IF(LMPR) ROOT=ONE/DBLE(K+J)                                        
CSP       IF(LMPR) ROOT=ONE/FLOAT(K+J)                                   
      FC=ERR**ROOT                                                       
      IF(FC.LT.FCM) FC=FCM                                               
C                                                                        
C----------------------------------------------------------------------- 
C  OPTIMAL ORDER DETERMINATION                                           
C----------------------------------------------------------------------- 
      OMJ=FC*AJ(J)                                                       
      IF(J.GT.2.AND.OMJ*ONE1.GT.OMJO.OR.K.GT.JOH) GOTO 3560              
      KO=K                                                               
      JO=J                                                               
      OMJO=OMJ                                                           
      FCO=FC                                                             
3560  IF(J.LT.KOH.AND.WINDOW) GOTO 300                                   
      IF(KONV.EQ.0) GOTO 3580                                            
      IF(KO.LT.K.OR.INCR(J).LT.0) GOTO 200                               
C                                                                        
C----------------------------------------------------------------------- 
C  POSSIBLE INCREASE OF ORDER                                            
C----------------------------------------------------------------------- 
      IF(NRED(KO).GT.0) NRED(KO)=NRED(KO)-1                              
      FC=FCO/AL(J,K)                                                     
      IF(FC.LT.FCM) FC=FCM                                               
      J1=J+1                                                             
      IF(AJ(J+1)*FC*ONE1.GT.OMJO) GOTO 200                               
      FCO=FC                                                             
      KO=JO                                                              
      JO=JO+1                                                            
      GOTO 200                                                           
C                                                                        
C----------------------------------------------------------------------- 
C  CONVERGENCE MONITOR                                                   
C----------------------------------------------------------------------- 
3580  RED=ONE/FCO                                                        
      JK=KM                                                              
      IF(JOH.LT.KM) JK=JOH                                               
      IF(K.GE.JK) GOTO 6000                                              
      IF(KO.LT.KOH) RED=AL(KOH,KO)*RED                                   
      IF(AL(JK,KO).LT.FCO) GOTO 6000                                     
300   CONTINUE                                                           
C                                                                        
C----------------------------------------------------------------------- 
C STEPSIZE REDUCTION (DUE TO EXTRAPOLATION TABLE)                        
C----------------------------------------------------------------------- 
6000  RED=RED*SAFE                                                       
      IF(RED.GT.RMAX) RED=RMAX                                           
      H=H*RED                                                            
6010  IF(.NOT.WINDOW) GOTO 6030                                          
      NRED(KOH)=NRED(KOH)+1                                              
      DO 6020 L=KOH,KM                                                   
6020  INCR(L)=-2-NRED(KOH)                                               
6030  JRED=JRED+1                                                        
      IF(PRM2) WRITE(MOUT,9015) JRED,RED,KOH                             
      IF(JRED.GT.JRMAX) GOTO 499                                         
      GOTO 3000                                                          
C                                                                        
C----------------------------------------------------------------------- 
C  STEPSIZE REDUCTION (DUE TO STABILITY CHECK)                           
C----------------------------------------------------------------------- 
6040  HMAX=G*FJ1*HALF                                                    
      RED=HMAX/H                                                         
      H=HMAX                                                             
      IF(PRM2) WRITE(MOUT,9010)                                          
C      IF(PRM2) WRITE(MOUT,12345) DM,IERRMX                              
C12345 FORMAT(' SCALED ERRMAX=',D21.12,' COMPONENT',I5)                  
      IF(JRED.GT.0.OR.NSTC.GT.0) GOTO 6030                               
      GOTO 6010                                                          
C                                                                        
C----------------------------------------------------------------------- 
C  STEPSIZE REDUCTION DUE TO OVERFLOW IN KIP1 OR NEGATIVE TEMPERATURE    
C  (EMPIRICAL DEVICE)                                                    
C----------------------------------------------------------------------- 
6060  KINRED=KINRED+1                                                    
C      IF(PRM2 .AND. KFAIL.EQ.-1) WRITE(MOUT,9016)                       
C      IF(PRM2 .AND. KFAIL.EQ.1) WRITE(MOUT,9013)                        
C      IF(PRM2 .AND. KFAIL.EQ.2) WRITE(MOUT,9014)                        
      RED=QUART                                                          
      HMAX=H*RED                                                         
      H=HMAX                                                             
      IF(JRED.GT.0 .OR. KINRED.GT.0) GOTO 6030                           
      GOTO 6010                                                          
C                                                                        
C----------------------------------------------------------------------- 
C  STEPSIZE REDUCTION (DUE TO SINGULAR MATRIX OR ARRAY OVERFLOW)         
C  (EMPIRICAL DEVICE)                                                    
C----------------------------------------------------------------------- 
6070  HMAX=G*FJ1*QUART                                                   
      RED=HMAX/H                                                         
      H=HMAX                                                             
      IF(PRM2) WRITE(MOUT,9011)                                          
      IF(JRED.LE.1) GOTO 6010                                            
      IF(IFLAG.GT.-3) GOTO 6010                                          
      IF(IFLAG.NE.-3) GOTO 799                                           
      IELBOW=MINIRN+2*N                                                  
      LIRN=LDIM+IELBOW                                                   
      LICN=NLEQ-LIRN                                                     
      LICNR=NWTMP-N                                                      
      IF(LICNR.LT.LICN) LICN=LICNR                                       
C  CHECK FOR SUFFICIENT STORAGE PLACE                                    
      XKAACT=DBLE(LICN)/DBLE(LDIM)                                       
CSP       XKAACT=FLOAT(LICN)/FLOAT(LDIM)
      IF(XKAACT.LT.XKAMIN) GOTO 799                                      
      LICNP1=LICN+1                                                      
      IDISP(1)=1                                                         
      BANFA=.TRUE.                                                       
      DEPS=TWO                                                           
      GOTO 6010                                                          
C                                                                        
C----------------------------------------------------------------------- 
C  PREPARATION FOR NEXT BASIC INTEGRATION STEP                           
C----------------------------------------------------------------------- 
200   TOLD=T                                                             
      T=TN                                                               
      H1=TEND-T                                                          
      ISTEP=ISTEP+1                                                      
      WINDOW=.TRUE.                                                      
C                                                                        
C  RESCALING                                                             
7112  DO 2350 I=1,N                                                      
      U=YM(I)                                                            
      Y(I)=U                                                             
C                                                                        
C  NEGATIVE CONCENTRATIONS ARE SET TO ZERO                               
      IF(U.LT.ZERO) Y(I)=ZERO                                              
C                                                                        
      U=DABS(U)                                                          
CSP       U=ABS(U)                                                       
      IF(U.GT.SH(I)) GOTO 2330                                           
      IF(IS(I).LE.0) GOTO 2350                                           
      IF(IS(I).EQ.2) GOTO 2310                                           
      IF(IS(I).EQ.3) GOTO 2320                                           
C COMPONENT WAS JUDGED DUE TO RELATIVE ERROR CRITERION AND FALLS NOW     
      IF(U.GE.SM(I)*CON2H) GOTO 2340                                     
C COMPONENT LOSES RELATIVE PRECISION                                     
      IS(I)=2                                                            
2310  IF(U.LT.SM(I)*CON2*EPS) IS(I)=3                                    
      GOTO 2350                                                          
2320  IF(U.GT.SM(I)/ONE2 .AND. PRM1) WRITE(MOUT,9020) I                  
      IS(I)=2                                                            
      GOTO 2350                                                          
C                                                                        
C  HANDLING OF RISING COMPONENTS                                         
2330  IF(IS(I).EQ.3 .AND. PRM1) WRITE(MOUT,9020) I                       
      IF(U.GT.SM(I).OR.IS(I).LE.0) SM(I)=U                               
2340  SH(I)=U                                                            
      IS(I)=1                                                            
2350  CONTINUE                                                           
C                                                                        
      IF(ISTEP.GT.ISMAX) GOTO 199                                        
C STEPSIZE PREDICTION                                                    
      H=H/FCO                                                            
      KOH=KO                                                             
      JOH=KOH+1                                                          
      IF(H.GT.T*EPMACH) GOTO 2000                                        
      GOTO 99                                                            
C                                                                        
C*********************************************************************** 
C                  END OF BASIC INTEGRATION STEP                         
C*********************************************************************** 
C                                                                        
C----------------------------------------------------------------------- 
C  SOLUTION EXIT                                                         
C----------------------------------------------------------------------- 
C                                                                        
7000  CONTINUE                                                           
      CALL SIMDAT (N,NSP,NCEQ,MODEL,Y,GAS1,T,IPRINT,LPRINT)              
C                                                                        
C  CHECK FOR POSSIBLE REDUCTION OF WORK SPACE                            
      NWTMP=NZAMAX                                                       
      XKAACT=DBLE(NZAMAX)/DBLE(LDIM)                                     
CSP       XKAACT=FLOAT(NZAMAX)/FLOAT(LDIM)                               
      MFILL=NZAMAX-LDIM                                                  
      FRATIO=(XKAACT-ONE)*1.D2                                           
CSP       FRATIO=(XKAACT-ONE)*1.E2                                       
      WRITE(MOUT,9350) THRESH,TGROW,LDIM,MFILL,FRATIO                    
      IF(ICNCP.GT.10.AND.PRM2) WRITE(MOUT,9310)                          
      IF(ICNCP.GT.3) GOTO 7050                                           
      NFREE=LICN-NZAMAX                                                  
      NFREEH=LICN-(LDIM*14)/10                                           
      IF(NFREE.GT.NFREEH) NFREE=NFREEH                                   
      NFREEH=LICN-(MINICN*12)/10                                         
      IF(NFREE.GT.NFREEH) NFREE=NFREEH                                   
      NFREE=(NFREE*9)/10                                                 
      IF(NFREE.LE.10) GOTO 7050                                          
      IF(PRM2) WRITE(MOUT,9320) NFREE                                    
7050  NSTEP=NSTEP+ISTEP                                                  
C                                                                        
C----------------------------------------------------------------------- 
C  NORMAL RETURN, NO ERROR                                               
C----------------------------------------------------------------------- 
C                                                                        
      RETURN                                                             
C----------------------------------------------------------------------- 
C  ERROR EXITS                                                           
C----------------------------------------------------------------------- 
C                                                                        
 99   IF(PRERR) WRITE(EOUT,9030) H,T                                     
      KFLAG=-1                                                           
      GOTO 999                                                           
199   IF(PRERR) WRITE(EOUT,9031) ISMAX                                   
      KFLAG=-2                                                           
      GOTO 999                                                           
299   IF(PRERR) WRITE(EOUT,9032)                                         
      KFLAG=-3                                                           
      T=TP(1)                                                            
      GOTO 999                                                           
399   IF(PRERR) WRITE(EOUT,9033)                                         
      KFLAG=-4                                                           
      GOTO 999                                                           
499   IF(PRERR) WRITE(EOUT,9034) JRMAX                                   
      KFLAG=-5                                                           
      GOTO 999                                                           
599   IF(PRERR) WRITE(EOUT,9035) (TP(I),I=1,ITPM)                        
      KFLAG=-6                                                           
      GOTO 999                                                           
699   IF(PRERR) WRITE(EOUT,9036)                                         
      KFLAG=-1                                                           
      GOTO 999                                                           
799   IF(PRERR) WRITE(EOUT,9032)                                         
      KFLAG=-3                                                           
      GOTO 999                                                           
899   IF(PRERR) WRITE(EOUT,9038)                                         
      KFLAG=-3                                                           
      GOTO 999                                                           
999   IF(PRERR) WRITE(EOUT,9039)                                         
      NSTEP=NSTEP+ISTEP                                                  
      TP(ITPM)=T                                                         
      HMAX=HMAXU                                                         
      RETURN                                                             
C                                                                        
8998  FORMAT(/,' SEMI-IMPLICIT    E U L E R  DISCRETISATION')            
8999  FORMAT(/,' SEMI-IMPLICIT    M I D - P O I N T  DISCRETISATION')    
9000  FORMAT(///,19X,'INTEGRATION MONITOR',/,19X,19('-'))                
9001  FORMAT(' USER PRESCRIBED MAXIMUM STEPSIZE: HMAX=',D10.3,/)         
CSP 9001  FORMAT(' USER PRESCRIBED MAXIMUM STEPSIZE: HMAX=',E10.3,/)     
9002  FORMAT(/,' PRESCRIBED STEPSIZE HMAX TOO SMALL; SET TO',D10.3)      
CSP 9002  FORMAT(/,' PRESCRIBED STEPSIZE HMAX TOO SMALL; SET TO',E10.3)  
9004  FORMAT(/,'    PRESCR. REL. PREC. EPS TOO SMALL; SET TO',D10.3)     
CSP 9004  FORMAT(/,'    PRESCR. REL. PREC. EPS TOO SMALL; SET TO',E10.3) 
9006  FORMAT(/,' CON1 TOO SMALL; REDUCED TO',D10.3)                      
CSP 9006  FORMAT(/,' CON1 TOO SMALL; REDUCED TO',E10.3)                  
9008  FORMAT(/,' CON2 TOO SMALL; REDUCED TO',D10.3)                      
CSP 9008  FORMAT(/,' CON1 TOO SMALL; REDUCED TO',E10.3)                  
9010  FORMAT(/,' STABILITY TEST FAILED')                                 
9011  FORMAT(/,' LU-DECOMPOSITION FAILED')                               
C9013  FORMAT(/,' OVERFLOW IN SUBROUTINE KIP1')                          
C9014  FORMAT(/,' TEMP NEGATIV IN SUBROUTINE KIP1')                      
9015  FORMAT(' STEPSIZE REDUCTION',I7,'.RED. FACTOR:',D10.3,I5)          
CSP 9015  FORMAT(' STEPSIZE REDUCTION',I7,'.RED. FACTOR:',E10.3,I5)      
C9016  FORMAT(' NEGATIVE TEMPERATURE')                                   
9020  FORMAT(/,' COMPONENT NO.:',I4,' WHICH HAS POSSIBLY NO RIGHT',      
     &          ' DIGIT BEGINNS TO ARISE',/,                              
     &          ' IF THE O.D.E. IS NOT STABLE, RESULTS MAY BE   WRONG',/, 
     &          ' PROBLEM SEEMS TO BE OSCILLARITY, SET CON2:=EPMACH')     
9030  FORMAT(/,' STEPSIZE TOO SMALL; H=',D18.10,' AT T=',D18.10)         
CSP 9030  FORMAT(/,' STEPSIZE TOO SMALL; H=',E15.7,' AT T=',E15.7)       
9031  FORMAT(/,' MORE THAN ISMAX=',I5,' BASIC STEPS')                    
9032  FORMAT(/,' NOT ENOUGH SPACE FOR LU-DECOMPOSITION',/,               
     &          ' ENLARGE INTEGER AND REAL WORK SPACE')                   
9033  FORMAT(/,' NEGATIVE INITIAL CONCENTRATION GIVEN:',D21.12)          
CSP 9033  FORMAT(/,' NEGATIVE INITIAL CONCENTRATION GIVEN:',E21.12)      
9034  FORMAT(/,' MORE THAN JRMAX=',I2,' STEPSIZE REDUCTIONS',            
     &          ' PER BASIC STEP')                                        
9035  FORMAT(/,' PRESCRIBED OUTPUT POINTS INVALID:',17(/,3D21.12))       
CSP 9035  FORMAT(/,' PRESCRIBED OUTPUT POINTS   INVALID:',17(/,3E21.12)) 
9036  FORMAT(/,' INTERNAL COMPUTATION OF INITIAL STEPSIZE GUESS NOT',    
     &  ' POSSIBLE DUE TO OVERFLOW IN KIP1.')                             
9038  FORMAT(/,' NUMBER OF EQUATIONS EXCEEDS 360. ADJUST ARRAYS ',       
     &  'ETA AND ASQU IN SUBROUTINE METAS1.')                             
9039  FORMAT(///,' INTEGRATION FAILED',/)                                
9100  FORMAT(/,' PRESCRIBED RELATIVE PRECISION: EPS=',D10.3,/)           
CSP 9100  FORMAT(/,' PRESCRIBED RELATIVE PRECISION: EPS=',E10.3,/)       
9110  FORMAT(//,8X,'STEP',13X,'T')                                       
9120  FORMAT(/,I10,7X,'    T=',D18.10)                                   
CSP 9120  FORMAT(/,I10,7X,'    T=',E18.10)                               
9200  FORMAT(/,' PRESCRIBED RELATIVE PRECISION: EPS=',D10.3)             
CSP 9200  FORMAT(/,' PRESCRIBED RELATIVE PRECISION: EPS=',E10.3)         
9201  FORMAT(' MAXIMUM COLUMN NUMBER:           KM=',I2)                 
9202  FORMAT(' SCALING PARAMETERS:              CON1=',D10.3)            
CSP 9202  FORMAT(' SCALING PARAMETERS:              CON1=',E10.3)        
9203  FORMAT(32X,'CON2=',D10.3,/)                                        
CSP 9203  FORMAT(32X,'CON2=',E10.3,/)                                    
9210  FORMAT(//,3X,'STEP',4X,'FCN',9X,'TIME',12X,'PREDICTED H',          
     &           9X,'K',4X,'KOH')                                         
9220  FORMAT(/,2I7,'   T=',D16.10,' H=',D16.10,2I6)                      
CSP 9220  FORMAT(/,2I7,'   T=',E16.10,' H=',E16.10,2I6)                  
9230  FORMAT(/,' ** NEW ANALYSE FACTOR **')                              
9300  FORMAT(/,' REAL AND INTEGER WORK SPACES RW AND IW NOT   BALANCED', 
     & /,' IW MAY BE REDUCED BY',I6,' STORAGE LOCATIONS')                 
9305  FORMAT(/,' REAL AND INTEGER WORK SPACES RW AND IW NOT   BALANCED', 
     & /,' RW MAY BE REDUCED BY',I6,' STORAGE LOCATIONS')                 
9310  FORMAT(/,' REAL AND INTEGER WORK SPACES RW AND IW MAY BE',         
     & ' ENLARGED TO SPEED UP',/,' LINEAR ALGEBRA')                       
9320  FORMAT(/,' REAL AND INTEGER WORK SPACES RW AND IW MAY BE',         
     & ' REDUCED BY',I6,/,' STORAGE LOCATIONS')                           
9350  FORMAT(/,1X,17('*'),/,' FILL-IN MONITOR,'/,1X,17('*'),/,           
     &' THRESHOLD VALUE            :',D12.5,/,                           
CSP      &' THRESHOLD VALUE            :',E12.5,/,                       
     &' WORST PIVOT RATIO          :',D12.5,/,                           
CSP      &' WORST PIVOT RATIO          :',E12.5,/,                       
     &' ACTUAL NNZ OF MATRIX A     :',I6,/,                              
     &' MAXIMUM ADDITIONAL FILL-IN :',I6,/,                              
     &' PERCENTAGE OF FILL-IN      :',F12.4)                             
C                                                                        
C  END METAS1                                                            
C                                                                        
      END 
      SUBROUTINE FCNINT (N,NSP,NCEQ,NTB,NHCP,T,Y,DY,HCP,SCCP,SUMC,       
     &   SUMDC,GAS1,GAS2,COLLI,MODEL,B,EPS,NB,NCTB,COE,ICOE,SK,LENTB,    
     &   KFAIL,ICOLLI,KINEV,AEA,ENTH,SK0,NAE,NEN,NSK,NCO,IPTD,TEMOLD,    
     &   LENTRO,LHCP,DDOT)                                               
C                                                                        
C*********************************************************************** 
C                                                                      * 
C DATE OF LATEST CHANGE: MARCH 1, '86                                  * 
C                                                                      * 
C*********************************************************************** 
C  FCNINT COMPUTES THE RIGHT HAND SIDE OF THE SYSTEM OF ODE'S          * 
C  ACCORDING TO THE CHOSEN MODEL.                                      * 
C*********************************************************************** 
C                                                                      * 
C  INPUT PARAMETERS:                                                   * 
C-------------------                                                   * 
C  N         NUMBER OF ODES                                            * 
C  NSP       NUMBER OF CHEMICAL SPECIES                                * 
C  NCEQ      NUMBER OF CHEMICAL REACTIONS                              * 
C  NTB       NUMBER OF THIRD BODIES                                    * 
C  T         ACTUAL TIME POINT                                         * 
C  Y(N)      Y(1),...,Y(NSP) : CONCENTRATIONS                          * 
C            Y(NSP+1)        : TEMPERATURE (FOR MODEL.GT.2)            * 
C            Y(N)            : PRESSURE    (FOR MODEL.EQ.3,4)          * 
C                            : DENSITY     (FOR MODEL.EQ.5,6)          * 
C  NB        DIMENSION OF B (= NSP FOR MODEL.EQ.5, = 1 ELSE)           * 
C  TEMOLD    TEMPERATURE VALUE OF THE PREVIOUS CALL                    * 
C  LENTRO    = .TRUE.  : RECOMPUTATION OF MOLAR ENTROPIES OF SPECIES   * 
C            = .FALSE. : NO RECOMPUTATION OF MOLAR ENTROPIES OF SPECIES* 
C  LHCP      = .TRUE.  : RECOMPUTATION OF MOLAR ENTHALPIES OF SPECIES  * 
C            = .FALSE. : NO RECOMPUTATION OF MOLAR ENTHALPIES OF SPECIE* 
C                                                                      * 
C  OUTPUT PARAMETERS:                                                  * 
C--------------------                                                  * 
C  DY(N)     TIME DERIVATIVES OF Y                                     * 
C  SCCP      AVERAGE HEAT CAPACITY                                     * 
C  SUMC      SUM OF MOLAR CONCENTRATIONS                               * 
C  SUMDC      = DY(1) + ... + DY(NSP)                                  * 
C  B(NB)     VARIABLE PART OF LEFT HAND SIDE MATRIX (FOR MODE.GE.5)    * 
C  KFAIL     =  0 SUCCESSFUL CALL OF FCN                               * 
C            = -1 FAIL CALL (NEGATIVE TEMPERATURE OR OVERFLOW IN KIP1) * 
C                                                                      * 
C  FOR A MORE DETAILED DESCRIPTION SEE SUBROUTINES 'METAS1', 'SIMULA'  * 
C  'KIP1', 'THERMO', 'ENTROP'.                                         * 
C                                                                      * 
C*********************************************************************** 
C                                                                        
      DOUBLE PRECISION Y(N),DY(N),HCP(NHCP),B(NB),SK(NSP),COE(NCTB)      
CSP       REAL Y(N),DY(N),HCP(NHCP),B(NB),SK(NSP),COE(NCTB)              
      DOUBLE PRECISION AEA(NAE),ENTH(NEN),SK0(NSK),COLLI(5)              
CSP       REAL AEA(NAE),ENTH(NEN),SK0(NSK),COLLI(5)                      
      DOUBLE PRECISION DENSY,DHSYS,DELTMP,DDOT,EPS,GAS1,GAS2,HEAT        
CSP       REAL DENSY,DHSYS,DELTMP,DDOT,EPS,GAS1,GAS2,HEAT                
      DOUBLE PRECISION PDOT,PRESS,RK,RT,S,SUMC,SUMDC                     
CSP       REAL PDOT,PRESS,RK,RT,S,SUMC,SUMDC                             
      DOUBLE PRECISION SCCP,SCCPH,T,TEMP                                 
CSP       REAL SCCP,SCCPH,T,TEMP                                         
      DOUBLE PRECISION TEMOLD,U,VOLUME,XMASS,ZERO                        
CSP       REAL TEMOLD,U,VOLUME,XMASS,ZERO                                
C                                                                        
CI4       INTEGER ICOLLI(NCEQ),KINEV(NCEQ),ICOE(NCTB),IPTD(NCO),LR,PLR   
      INTEGER*2 ICOLLI(NCEQ),KINEV(NCEQ),ICOE(NCTB),IPTD(NCO),LR,PLR     
      INTEGER LENTB(5)                                                   
      INTEGER I,ICOLK,IFLAG,IH,IL1,IL2,IR2,ITB,JJ,JJH,K,LRI,NSP1         
      LOGICAL LENTRO,LHCP                                                
C                                                                        
      REAL R                                                             
C                                                                        
      EXTERNAL THERMO,KIP1,ENTROP                                        
C                                                                        
      EXTERNAL HEAT,DPROF,PPROF                                          
C                                                                        
      COMMON / LARK1  / LR(5000)                                         
      COMMON / LARK2  / PLR(2000)                                        
      COMMON / KINPAR / RK(1000)                                         
      COMMON / RATE1  / R(1000)                                          
      COMMON / RECOND / TEMP,DENSY,PRESS,VOLUME,XMASS 
      SAVE /LARK1/, /LARK2/, /KINPAR/, /RATE1/, /RECOND/
C                                                                        
      DATA ZERO/0.D0/                                                    
CSP       DATA ZERO/0.E0/                                                
C                                                                        
      IFLAG=KFAIL                                                        
      KFAIL=0                                                            
      NSP1=NSP+1                                                         
      PDOT=ZERO                                                          
      DDOT=ZERO                                                          
      U=ZERO                                                             
C                                                                        
      IF(NTB.EQ.0) GOTO 10                                               
      JJ=0                                                               
      DO 9 ITB=1,NTB                                                     
      COLLI(ITB)=ZERO                                                    
      IH=LENTB(ITB)                                                      
      IF(IH.LT.0) GOTO 11                                                
      DO 8 I=1,IH                                                        
      JJ=JJ+1                                                            
      JJH=ICOE(JJ)                                                       
8     COLLI(ITB)=COLLI(ITB) + COE(JJ)*Y(JJH)                             
      GOTO 9                                                             
11    S=ZERO                                                             
      DO 12 I=1,NSP                                                      
12    S=S + Y(I)                                                         
      COLLI(ITB)=S                                                       
9     CONTINUE                                                           
10    CONTINUE                                                           
C                                                                        
      IF(MODEL.LE.2) GOTO 25                                             
      TEMP=Y(NSP1)                                                       
      DELTMP=DABS(TEMOLD-TEMP)                                           
CSP       DELTMP=ABS(TEMOLD-TEMP)                                        
      IF(DELTMP.LE.0.D0) GOTO 24                                         
      TEMOLD=TEMP                                                        
C  RECALCULATION OF MOLAR ENTHALPIES AND HEAT CAPACITIES(IF LHCP=.TRUE.) 
      IF(LHCP) CALL THERMO (NSP,NHCP,HCP,NEN,ENTH,NCO,IPTD,TEMP,         
     &                      GAS1,IFLAG)                                  
      KFAIL=IFLAG                                                        
      IF(KFAIL.NE.0) RETURN                                              
C  RECALCULATION OF MOLAR ENTROPIES (IF LENTR=.TRUE.)                    
      IF(LENTRO) CALL ENTROP (NSP,TEMP,GAS1,SK,NEN,NSK,ENTH,SK0,         
     &                        NCO,IPTD,KFAIL)                            
      IF(KFAIL.NE.0) RETURN                                              
24    IF(DELTMP.LE.ZERO) GOTO 25                                         
      CALL KIP1 (NSP,NCEQ,NHCP,HCP,SK,NAE,AEA,KINEV,                     
     &           TEMP,RK,GAS1,GAS2,KFAIL)                                
      IF(KFAIL.NE.0) RETURN                                              
25    CONTINUE                                                           
C                                                                        
C  TIME DERIVATIVES OF SPECIES-CONCENTRATIONS (MOL/M**3-SEC)             
C  ---------------------------------------------------------             
C                                                                        
C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
C  ONLY FOR MODEL=2:                                                     
C  -----------------                                                     
C  IF A  T I M E  -  DEPENDENT TEMPERATURE PROFIL IS GIVEN,              
C  LARKIN CALLS THE USER PROVIDED FUNCTION                               
C                                                                        
C            DOUBLE PRECISION FUNCTION HEAT (T)                          
C                                                                        
C  WHICH RETURNS THE TEMPERATURE TEMP=HEAT(T) AT ACTUAL TIME             
C  POINT T. OTHERWISE LARKIN CALLS A DUMMY FUNCTION.                     
C                                                                        
      IF(MODEL.EQ.2) TEMP=HEAT(T)                                        
      IF(MODEL.EQ.2) CALL KIP1 (NSP,NCEQ,NHCP,HCP,SK,NAE,AEA,KINEV,      
     &                          TEMP,RK,GAS1,GAS2,KFAIL)                 
C                                                                        
      DO 30 I=1,N                                                        
30    DY(I)=ZERO                                                         
      IR2=0                                                              
      DO 80 K=1,NCEQ                                                     
      RT=RK(K)                                                           
      ICOLK=ICOLLI(K)                                                    
      IF(ICOLK.NE.0) RT=RT*COLLI(ICOLK)                                  
      IL1=IR2+1                                                          
      IL2=PLR(K)                                                         
      IR2=PLR(NCEQ+K)                                                    
      IF(IL2.LT.IL1) GOTO 50                                             
      DO 40 I=IL1,IL2                                                    
      LRI=LR(I)                                                          
40    RT=RT*Y(LRI)                                                       
50    CONTINUE                                                           
      DO 70 I=IL1,IR2                                                    
      LRI=LR(I)                                                          
      IF(I.GT.IL2) GOTO 60                                               
      DY(LRI)=DY(LRI) - RT                                               
      GOTO 70                                                            
60    DY(LRI)=DY(LRI) + RT                                               
70    CONTINUE                                                           
C  R(K) STORES THE REACTION RATE OF EQUATION K (FOR EXTERNAL USE)        
80    R(K)=RT                                                            
C                                                                        
      IF(MODEL.LE.2) RETURN                                              
C                                                                        
      DHSYS=ZERO                                                         
      SUMDC=ZERO                                                         
      SCCP=ZERO                                                          
      DO 100 I=1,NSP                                                     
      DHSYS=DHSYS + HCP(I)*DY(I)                                         
      SUMDC=SUMDC + DY(I)                                                
100   SCCP=SCCP + HCP(NSP+I)*Y(I)                                        
C                                                                        
      IF(MODEL.EQ.5) GOTO 125                                            
C                                                                        
C  ADDITIONAL DERIVATIVES FOR MODEL=3,4                                  
C--------------------------------------                                  
C                                                                        
C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
C  ONLY FOR MODEL=4:                                                     
C-------------------                                                     
C  IF A  T I M E  -  DEPENDENT DENSITY-PROFILE IS ASSUMED TO BE GIVEN,   
C  LARKIN CALLS THE USER-PROVIDED SUBROUTINE                             
C                                                                        
C            SUBROUTINE DPROF (T,DENSY,DDOT)                             
C                                                                        
C  WHERE     T     :  ACTUAL TIME POINT                                  
C            DENSY :  DENSITY AT TIME T                                  
C            DDOT  :  TIME DERIVATIVE OF DENSITY AT TIME T               
C  OTHERWISE LARKIN CALLS A DUMMY-SUBROUTINE                             
C                                                                        
      IF(MODEL.EQ.4) CALL DPROF (T,DENSY,DDOT)                           
      IF(MODEL.EQ.4) U=DDOT/DENSY                                        
C                                                                        
      SUMC=ZERO                                                          
      DO 210 I=1,NSP                                                     
210   SUMC=SUMC+Y(I)                                                     
      PRESS=GAS1*TEMP*SUMC                                               
      IF(DENSY.NE.ZERO) VOLUME=XMASS/DENSY                               
      SCCPH=GAS1*SUMC - SCCP                                             
C                                                                        
C  TIME DERIVATIVE OF TEMPERATURE                                        
      S=DHSYS - (U*PRESS + GAS1*TEMP*SUMDC)                              
      DY(NSP1)=S/SCCPH                                                   
C                                                                        
      IF(DDOT.EQ.ZERO) RETURN                                            
      DO 221 I=1,NSP                                                     
221   DY(I)=DY(I) + Y(I)*U                                               
C                                                                        
      RETURN                                                             
C                                                                        
125   CONTINUE                                                           
C  ADDITIONAL DERIVATIVES FOR MODEL=5,6                                  
C--------------------------------------                                  
C                                                                        
C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
C  ONLY FOR MODEL=6:                                                     
C-------------------                                                     
C  IF A  T I M E  -  DEPENDENT PRESSURE-PROFILE IS ASSUMED TO BE GIVEN,  
C  LARKIN CALLS THE USER-PROVIDED SUBROUTINE                             
C                                                                        
C            SUBROUTINE PPROF (T,PRESS,PDOT)                             
C                                                                        
C  WHERE     T     :  ACTUAL TIME POINT                                  
C            PRESS :  PRESSURE AT TIME T                                 
C            PDOT  :  TIME DERIVATIVE OF PRESSURE AT TIME T              
C  OTHERWISE LARKIN CALLS A DUMMY-SUBROUTINE                             
C                                                                        
      IF(MODEL.EQ.6) CALL PPROF (T,PRESS,PDOT)                           
      IF(MODEL.EQ.6) U=PDOT/PRESS                                        
C                                                                        
      DENSY=Y(N)                                                         
      VOLUME=XMASS/DENSY                                                 
      SUMC=PRESS/GAS1/TEMP                                               
C  TIME DERIVATIVE OF TEMPERATURE                                        
      DY(NSP1)=-( (DHSYS-PDOT)/SCCP)                                     
C                                                                        
C  TIME DERIVATIVE OF DENSITY                                            
      DY(N)=DENSY * (U - ( SUMDC/SUMC + DY(NSP1)/TEMP ))                 
C                                                                        
C  UPDATE VARIABLE PART OF MATRIX B                                      
      DO 150 I=1,NSP                                                     
150   B(I)=-Y(I)/DENSY                                                   
C                                                                        
      RETURN                                                             
C                                                                        
C----------------------------------------------------------------------- 
C               END OF SUBROUTINE FCN                                    
C----------------------------------------------------------------------- 
C                                                                        
      END  
      SUBROUTINE JACOB1(N,NSP,NCEQ,NHCP,LDIM,T,Y,A,DY,HCP,SCCP,SUMC,
     & SUMDC,GAS1,GAS2,COLLI,NCTB,COE,ICOE,LENTB,MODEL,ICOLLI,KINEV,JA,  
     & IA,IDA,AEA,NIA,NAE,DDOT)                                          
C----------------------------------------------------------------------- 
C  JACOBI COMPUTES THE JACOBIANS OF THE FCN'S ACCORDING TO THE CHOSEN    
C  MODEL                                                                 
C----------------------------------------------------------------------- 
C                                                                        
C  DATE OF LATEST CHANGE: 28 - 1 -86                                     
C                                                                        
      DOUBLE PRECISION Y(N),DY(N),HCP(NHCP),A(LDIM)                      
CSP       REAL Y(N),DY(N),HCP(NHCP),A(LDIM)                              
      DOUBLE PRECISION COE(NCTB),AEA(NAE)                                
CSP       REAL COE(NCTB),AEA(NAE)                                        
      DOUBLE PRECISION COLLI(5),DCKDT,DENSY,DYT                          
CSP       REAL COLLI(5),DCKDT,DENSY,DYT                                  
      DOUBLE PRECISION GAS1,GAS2,PRESS,RK,RT                             
CSP       REAL GAS1,GAS2,PRESS,RK,RT                                     
      DOUBLE PRECISION RKK,RATEK,SCCP,SCCPH                              
CSP       REAL RKK,RATEK,SCCP,SCCPH                                      
      DOUBLE PRECISION SUMC,SUMDC,SHPC,T,TH,TEMP                         
CSP       REAL SUMC,SUMDC,SHPC,T,TH,TEMP                                 
      DOUBLE PRECISION TERM,TERM1,U,VOLUME,XMASS,ZERO,DDOT               
CSP       REAL TERM,TERM1,U,VOLUME,XMASS,ZERO,DDOT                       
C                                                                        
CI4       INTEGER LR,PLR,ICOLLI(NCEQ),ICOE(NCTB),KINEV(NCEQ)             
      INTEGER*2 LR,PLR,ICOLLI(NCEQ),ICOE(NCTB),KINEV(NCEQ)               
CI4       INTEGER JA(LDIM),IA(NIA),IDA(N)                                
      INTEGER*2 JA(LDIM),IA(NIA),IDA(N)                                  
      INTEGER IN(10),I1(10),LENTB(5)                                     
      INTEGER I,I3,I4,ICOLK,ID,IDL,IDT,IK,IL1,IL2,IM,IR2,IT,ITL,         
     & J,JH1,JH2,JJ,K,K1,K2,KINEVK,L,LCK1,LD,LL,LRI,LT,M,NSP1,NSP2       
C                                                                        
      COMMON / LARK1  / LR(5000)                                         
      COMMON / LARK2  / PLR(2000)                                        
      COMMON / KINPAR / RK(1000)                                         
      COMMON / RECOND / TEMP,DENSY,PRESS,VOLUME,XMASS
      SAVE /LARK1/, /LARK2/, /KINPAR/, /RECOND/
C                                                                        
      DATA ZERO /0.D0/                                                   
CSP       DATA ZERO /0.E0/                                               
C                                                                        
      DO 10 L=1,LDIM                                                     
10    A(L)=ZERO                                                          
      NSP1=NSP+1                                                         
C                                                                        
C----------------------------------------------------------------------- 
C  SPECIES DERIVATIVES OF RATES                                          
C----------------------------------------------------------------------- 
C                                                                        
      RT=GAS2*TEMP                                                       
      IR2=0                                                              
      DO 180 K=1,NCEQ                                                    
      IDT=0                                                              
      ICOLK=ICOLLI(K)                                                    
      KINEVK=KINEV(K)                                                    
      K1=NCEQ+K                                                          
      K2=NCEQ+K1                                                         
      IL1=IR2+1                                                          
      IL2=PLR(K)                                                         
      IK=IL2-IR2                                                         
      IR2=PLR(NCEQ+K)                                                    
      IF(IK.LE.0) GOTO 180                                               
      IM=0                                                               
      DO 20 M=IL1,IR2                                                    
      IM=IM+1                                                            
      LRI=LR(M)                                                          
      I1(IM)=IDA(LRI)                                                    
20    IN(IM)=LRI                                                         
      RKK=RK(K)                                                          
      DO 90 J=1,IK                                                       
      IDT=IDT+1                                                          
      TERM=RKK                                                           
      JJ=IN(J)                                                           
      DO 30 I=1,IK                                                       
      L=IN(I)                                                            
30    IF(J.NE.I) TERM=TERM*Y(L)                                          
      RATEK=TERM*Y(JJ)                                                   
      IF(ICOLK.GT.0) TERM=TERM*COLLI(ICOLK)                              
C                                                                        
      IF(MODEL.GT.2.AND.IDT.EQ.1) GOTO 700                               
C                                                                        
C  TERM IS A SPECIES DERIVATIVE OF THE RATE OF REACTION K                
      DO 80 I=1,IM                                                       
      L=I1(I)                                                            
      IF(JA(L)-JJ) 40,60,50                                              
40    L=L+1                                                              
      IF(JA(L).NE.JJ) GOTO 40                                            
      GOTO 60                                                            
50    L=L-1                                                              
      IF(JA(L).NE.JJ) GOTO 50                                            
60    IF(I.GT.IK) GOTO 70                                                
      A(L)=A(L) - TERM                                                   
      GOTO 80                                                            
70    A(L)=A(L) + TERM                                                   
80    CONTINUE                                                           
      GOTO 90                                                            
C                                                                        
700   CONTINUE                                                           
C  TERM1 IS THE TEMPERATURE DERIVATIVE OF THE RATE OF REACTION K         
      TERM1=ZERO                                                         
      IF(KINEVK.GT.1) TERM1=RATEK*(AEA(K2)+AEA(K1)/RT)/TEMP              
      IF(ICOLK.GT.0) TERM1=TERM1*COLLI(ICOLK)                            
C                                                                        
      DO 420 I=1,IM                                                      
      L=I1(I)                                                            
      LL=JA(L)                                                           
      LL=IA(LL+1)-2                                                      
      IF(MODEL.LE.4) LL=LL+1                                             
      IF(JA(L)-JJ) 380,400,390                                           
380   L=L+1                                                              
      IF(JA(L).NE.JJ) GOTO 380                                           
      GOTO 400                                                           
390   L=L-1                                                              
      IF(JA(L).NE.JJ) GOTO 390                                           
400   IF(I.GT.IK) GOTO 410                                               
      A(L)=A(L) - TERM                                                   
      A(LL)=A(LL) - TERM1                                                
      GOTO 420                                                           
410   A(L)=A(L) + TERM                                                   
      A(LL)=A(LL) + TERM1                                                
420   CONTINUE                                                           
C                                                                        
90    CONTINUE                                                           
      IF(ICOLK.EQ.0) GOTO 180                                            
      LCK1=LENTB(ICOLK)                                                  
      IF(LCK1.GE.0) GOTO 95                                              
      JH1=1                                                              
      JH2=NSP                                                            
      GOTO 110                                                           
95    JH2=0                                                              
      DO 100 I=1,ICOLK                                                   
100   JH2=JH2 + LENTB(I)                                                 
      JH1=1                                                              
      IF(ICOLK.GT.1) JH1=JH2-LENTB(ICOLK)+1                              
110   DO 170 J=JH1,JH2                                                   
      TERM=RATEK                                                         
      JJ=J                                                               
      IF(LCK1.LT.0) GOTO 115                                             
      TERM=RATEK * COE(J)                                                
      JJ=ICOE(J)                                                         
115   CONTINUE                                                           
      DO 160 I=1,IM                                                      
      L=I1(I)                                                            
      IF(JA(L)-JJ) 120,140,130                                           
120   L=L+1                                                              
      IF(JA(L).NE.JJ) GOTO 120                                           
      GOTO 140                                                           
130   L=L-1                                                              
      IF(JA(L).NE.JJ) GOTO 130                                           
140   IF(I.GT.IK) GOTO 150                                               
      A(L)=A(L) - TERM                                                   
      GOTO 160                                                           
150   A(L)=A(L) + TERM                                                   
160   CONTINUE                                                           
170   CONTINUE                                                           
180   CONTINUE                                                           
C                                                                        
      IF(MODEL.LE.2) RETURN                                              
C                                                                        
      IF(MODEL.EQ.3) GOTO 600                                            
C----------------------------------------------------------------------- 
C  ADDITIONAL DERIVATIVES FOR MODEL=5,6                                  
C----------------------------------------------------------------------- 
C                                                                        
      DENSY=Y(N)                                                         
C  SPECIES-DERIVATIVES OF T-DOT AND OF D-DOT                             
      DYT=DY(NSP1)                                                       
      IT=IA(NSP1)-1                                                      
      ID=IA(N)-1                                                         
      TH=SUMC*SUMC                                                       
      DO 530 I=1,NSP                                                     
      L=IT+I                                                             
530   A(L)=HCP(NSP+I)*DYT                                                
      DO 550 K=1,NSP                                                     
      I3=IA(K)                                                           
      I4=IA(K+1)-3                                                       
      DO 540 I=I3,I4                                                     
      L=JA(I)                                                            
      ITL=IT+L                                                           
      IDL=ID+L                                                           
      A(IDL)=A(IDL) + A(I)                                               
540   A(ITL)=A(ITL) + A(I)*HCP(K)                                        
550   CONTINUE                                                           
      DO 560 I=1,NSP                                                     
      LT=IT+I                                                            
      LD=ID+I                                                            
      A(LT)=-A(LT)/SCCP                                                  
560   A(LD)=-DENSY*( A(LT)/TEMP + A(LD)/SUMC )                           
C                                                                        
C TEMPERATURE-DERIVATIVE OF T-DOT                                        
      TERM1=ZERO                                                         
      DCKDT=ZERO                                                         
      NSP2=2*NSP                                                         
      DO 570 I=1,NSP                                                     
      L=IA(I+1)-2                                                        
      DCKDT=DCKDT + A(L)                                                 
570   TERM1=TERM1+DYT*Y(I)*HCP(NSP2+I) + DY(I)*HCP(NSP+I) + HCP(I)*A(L)  
      LT=IA(NSP+2)-2                                                     
      A(LT)=-TERM1/SCCP                                                  
C                                                                        
C  TEMPERATURE-DERIVATIVE OF D-DOT                                       
      LD=IA(NSP+3)-2                                                     
      A(LD)=-DENSY*( (TEMP*A(LT)-DY(NSP1))/(TEMP*TEMP) +                 
     &              GAS1*(TEMP*DCKDT+SUMDC)/PRESS  )                     
C                                                                        
C  DENSITY-DERIVATIVE OF D-DOT                                           
      LD=IA(NSP+3)-1                                                     
      A(LD)=DY(N)/DENSY                                                  
C                                                                        
      RETURN                                                             
C                                                                        
600   CONTINUE                                                           
C----------------------------------------------------------------------- 
C  ADDITIONAL DERIVATIVES FOR MODEL=3,4                                  
C----------------------------------------------------------------------- 
C                                                                        
C     U=DDOT/DENSY                                                       
C      U=-( (DY(NSP1)/Y(NSP1)+SUMDC/SUMC) - DY(N)/Y(N) )                 
      U=ZERO                                                             
C                                                                        
      SCCPH=GAS1*SUMC - SCCP                                             
C  SPECIES-DERIVATIVES OF T-DOT                                          
      DYT=DY(NSP1)                                                       
      IT=IA(NSP1)-1                                                      
      SHPC=GAS1*TEMP                                                     
      DO 610 I=1,NSP                                                     
      L=IT+I                                                             
610   A(L)=(HCP(NSP+I)-GAS1)*DYT-SHPC*U                                  
      DO 630 K=1,NSP                                                     
      I3=IA(K)                                                           
      I4=IA(K+1)-2                                                       
      DO 620 I=I3,I4                                                     
      L=JA(I)                                                            
      ITL=IT+L                                                           
620   A(ITL)=A(ITL) + A(I)*(HCP(K)-SHPC)                                 
630   CONTINUE                                                           
      DO 640 I=1,NSP                                                     
      LT=IT+I                                                            
640   A(LT)=A(LT)/SCCPH                                                  
C                                                                        
C  TEMPERATURE-DERIVATIVES OF T-DOT AND P-DOT                            
      NSP2=2*NSP                                                         
      TERM1=ZERO                                                         
      DO 650 I=1,NSP                                                     
      L=IA(I+1)-1                                                        
650   TERM1=TERM1 + DYT*Y(I)*HCP(NSP2+I) + DY(I)*(HCP(NSP+I)-GAS1) +     
     &              (HCP(I)-SHPC)*A(L)                                   
C  T-DERIVATIVE OF T-DOT                                                 
      LT=IA(NSP+2)-1                                                     
      A(LT)=TERM1/SCCPH                                                  
C                                                                        
C  UPDATE PARTIAL SPECIES DERIVATIVES                                    
      IF(U.EQ.ZERO) RETURN                                               
      DO 660 I=1,NSP                                                     
      L=IDA(I)                                                           
660   A(L)=A(L) + U                                                      
C                                                                        
      RETURN                                                             
C                                                                        
C----------------------------------------------------------------------- 
C  END OF SUBROUTINE JACOB1                                              
C----------------------------------------------------------------------- 
      END 
      SUBROUTINE MA30LA (N,ICN,A,LICN,LENR,LENRL,IDISP,IP,IQ,IRN,        
     &LIRN,LENC,IFIRST,LASTR,NEXTR,LASTC,NEXTC,IPTR,IPC,UIN,IFLAG)       
C                                                                        
C  ANALYSE AND FACTOR SPARSE MATRIX                                      
C                                                                        
C                                                                        
C DATE OF LATEST CHANGE: APRIL 27, '82                                   
C                                                                        
      DOUBLE PRECISION A(LICN),U,UIN,AU,UMAX,AMAX,ZERO                   
CSP       REAL A(LICN)                                                   
      DOUBLE PRECISION DABS,DMAX1,DMIN1                                  
CSP       REAL ABS,AMAX1,AMIN1                                           
      INTEGER PIVOT,PIVEND,DISPC,OLDPIV,OLDEND,PIVROW                    
      INTEGER ROWI                                                       
      INTEGER I, I1, I2, IACTIV, IBEG, IDISPC, IDUMMY, IEND, IFILL,      
     & IFIR, II, III, IJFIR, IJP1, IJPOS, INDROW, IOP, IPIV, IPOS,       
     & ISING, ISW, ISW1, ITOP, J, J1, J2, JBEG, JCOST, JCOUNT, JDIFF,    
     & JEND, JJ, JMORE, JNPOS, JOLD, JPIV, JPOS, JROOM, JVAL, JZER,      
     & JZERO, K, KCOST, L, LC, LENPIV, LL, LR, MOREI, NC, NM1, NR, NZ,   
     & NZ2, NZCOL, NZMIN, NZPC, NZROW                                    
      INTEGER IDISP(2)                                                   
      LOGICAL ABORT1,ABORT2,ABORT3                                       
CI4       INTEGER ICN(LICN),LENR(N),LENRL(N),IP(N),IQ(N),LENC(N)         
      INTEGER*2 ICN(LICN),LENR(N),LENRL(N),IP(N),IQ(N),LENC(N)           
C                                                                        
C FOR HANLING OF MORE THAN 32000 NONZEROS WITH INTEGER*2 VERSION CHANGE  
C ARRAYS IPC AND IPATR TO STANDARD INTEGER                               
C                                                                        
CI4       INTEGER IRN(LIRN),IPC(N),IPTR(N)                               
      INTEGER*2 IRN(LIRN),IPC(N),IPTR(N)                                 
CI4       INTEGER IFIRST(N),LASTR(N),NEXTR(N),LASTC(N),NEXTC(N)          
      INTEGER*2 IFIRST(N),LASTR(N),NEXTR(N),LASTC(N),NEXTC(N)            
C                                                                        
      COMMON /MA30LE/ LP,ABORT1,ABORT2,ABORT3                            
      COMMON /MA30LF/ IRNCP,ICNCP,IRANK,MINIRN,MINICN 
      SAVE /MA30LE/, /MA30LF/
C                                                                        
      DATA UMAX/.999999999D0/                                            
CSP       DATA UMAX/.9999/                                               
      DATA ZERO/0.0D0/                                                   
CSP       DATA ZERO/0.0/                                                 
C                                                                        
      MINIRN=0                                                           
      MINICN=IDISP(1)-1                                                  
      MOREI=0                                                            
      IRANK=N                                                            
      IRNCP=0                                                            
      ICNCP=0                                                            
      IFLAG=0                                                            
C RESET U IF NECESSARY.                                                  
      U=DMIN1(UIN,UMAX)                                                  
CSP       U=AMIN1(UIN,UMAX)                                              
      U=DMAX1(U,ZERO)                                                    
CSP       U=AMAX1(U,ZERO)                                                
C IBEG IS THE POSITION OF THE NEXT PIVOT ROW AFTER ELIMINATION STEP      
C     USING IT.                                                          
      IBEG=IDISP(1)                                                      
C IACTIV IS THE POSITION OF THE FIRST ENTRY IN THE ACTIVE PART OF A/ICN. 
      IACTIV=IDISP(2)                                                    
C NZROW IS CURRENT NUMBER OF NON-ZEROS IN ACTIVE AND UNPROCESSED PART    
C     OF ROW FILE ICN.                                                   
      NZROW=LICN-IACTIV+1                                                
      MINICN=NZROW+MINICN                                                
C                                                                        
C     SET UP POINTERS TO THE                                             
C     BEGINNINGS OF THE ROWS.                                            
      IPTR(1)=IACTIV                                                     
      IF (N.EQ.1) GO TO 20                                               
      NM1=N-1                                                            
      DO 10 I=1,NM1                                                      
   10 IPTR(I+1)=IPTR(I)+LENR(I)                                          
C                                                                        
C ***********************************************                        
C ****    LU DECOMPOSITION OF MATRIX         ****                        
C ***********************************************                        
C                                                                        
C IACTIV IS THE POSITION OF THE FIRST ELEMENT IN THE MATRIX              
C ITOP IS THE POSITION OF THE LAST ELEMENT IN THE MATRIX.                
  20  ITOP=LICN                                                          
C                                                                        
C SET UP COLUMN ORIENTED STORAGE.                                        
      DO 110 I=1,N                                                       
      LENRL(I)=0                                                         
  110 LENC(I)=0                                                          
      IF (ITOP-IACTIV.LT.LIRN) GO TO 120                                 
      MINIRN=ITOP-IACTIV+1                                               
      PIVOT=0                                                            
      GO TO 1050                                                         
C                                                                        
C CALCULATE COLUMN COUNTS.                                               
  120 DO 130 II=IACTIV,ITOP                                              
      I=ICN(II)                                                          
  130 LENC(I)=LENC(I)+1                                                  
C SET UP COLUMN POINTERS SO THAT IPC(J) POINTS TO POSITION AFTER END     
C     OF COLUMN J IN COLUMN FILE.                                        
      IPC(N)=LIRN+1                                                      
      IF(N.EQ.1) GOTO 142                                                
      DO 140 JJ=2,N                                                      
      J=N-JJ+1                                                           
  140 IPC(J)=IPC(J+1)-LENC(J+1)                                          
  142 DO 160 INDROW=1,N                                                  
      J1=IPTR(INDROW)                                                    
      J2=J1+LENR(INDROW)-1                                               
      IF (J1.GT.J2) GO TO 160                                            
      DO 150 JJ=J1,J2                                                    
      J=ICN(JJ)                                                          
      IPOS=IPC(J)-1                                                      
      IRN(IPOS)=INDROW                                                   
      IPC(J)=IPOS                                                        
  150 CONTINUE                                                           
  160 CONTINUE                                                           
C DISPC IS THE LOWEST INDEXED ACTIVE LOCATION IN THE COLUMN FILE.        
      DISPC=IPC(1)                                                       
      NZCOL=LIRN-DISPC+1                                                 
      MINIRN=MAX0(NZCOL,MINIRN)                                          
      NZMIN=1                                                            
C                                                                        
C INITIALIZE ARRAY IFIRST.  IFIRST(I) = +/- K INDICATES THAT ROW/COL     
C     K HAS I NON-ZEROS.  IF IFIRST(I) = 0, THERE IS NO ROW OR COLUMN    
C     WITH I NON ZEROS.                                                  
      DO 170 I=1,N                                                       
  170 IFIRST(I)=0                                                        
C                                                                        
C COMPUTE ORDERING OF ROW AND COLUMN COUNTS.                             
C FIRST RUN THROUGH COLUMNS (FROM COLUMN N TO COLUMN 1).                 
      DO 190 JJ=1,N                                                      
      J=N-JJ+1                                                           
      NZ=LENC(J)                                                         
      IF (NZ.NE.0) GO TO 180                                             
      IPC(J)=0                                                           
      LASTC(J)=0                                                         
      GO TO 190                                                          
  180 ISW=IFIRST(NZ)                                                     
      IFIRST(NZ)=-J                                                      
      LASTC(J)=0                                                         
      NEXTC(J)=-ISW                                                      
      ISW1=IABS(ISW)                                                     
      IF (ISW.NE.0) LASTC(ISW1)=J                                        
  190 CONTINUE                                                           
C NOW RUN THROUGH ROWS (AGAIN FROM N TO 1).                              
      DO 210 II=1,N                                                      
      I=N-II+1                                                           
      NZ=LENR(I)                                                         
      IF (NZ.NE.0) GO TO 200                                             
      IPTR(I)=0                                                          
      LASTR(I)=0                                                         
      GO TO 210                                                          
  200 ISW=IFIRST(NZ)                                                     
      IFIRST(NZ)=I                                                       
      IF (ISW.GT.0) GO TO 205                                            
      NEXTR(I)=0                                                         
      LASTR(I)=ISW                                                       
      GO TO 210                                                          
 205  NEXTR(I)=ISW                                                       
      LASTR(I)=LASTR(ISW)                                                
      LASTR(ISW)=I                                                       
  210 CONTINUE                                                           
C                                                                        
C                                                                        
C **********************************************                         
C ****    START OF MAIN ELIMINATION LOOP    ****                         
C **********************************************                         
      DO 930 PIVOT=1,N                                                   
C                                                                        
C FIRST FIND THE PIVOT USING MARKOWITZ CRITERION WITH STABILITY          
C     CONTROL.                                                           
C JCOST IS THE MARKOWITZ COST OF THE BEST PIVOT SO FAR,.. THIS           
C     PIVOT IS IN ROW IPIV AND COLUMN JPIV.                              
      NZ2=NZMIN                                                          
      JCOST=N*N                                                          
C                                                                        
C EXAMINE ROWS/COLUMNS IN ORDER OF ASCENDING COUNT.                      
      DO 290 L=1,2                                                       
      LL=L                                                               
C A PASS WITH L EQUAL TO 2 IS ONLY PERFORMED IN THE CASE OF SINGULARITY. 
      DO 280 NZ=NZ2,N                                                    
      IF (JCOST.LE.(NZ-1)**2) GO TO 380                                  
      IJFIR=IFIRST(NZ)                                                   
      IF (IJFIR) 212,211,215                                             
 211  IF (LL.EQ.1) NZMIN=NZ+1                                            
      GO TO 280                                                          
 212  LL=2                                                               
      IJFIR=-IJFIR                                                       
      GO TO 245                                                          
 215  LL=2                                                               
C SCAN ROWS WITH NZ NON-ZEROS.                                           
      DO 235 IDUMMY=1,N                                                  
      IF (IJFIR.EQ.0) GO TO 240                                          
C ROW IJFIR IS NOW EXAMINED.                                             
      I=IJFIR                                                            
      IJFIR=NEXTR(I)                                                     
C FIRST CALCULATE MULTIPLIER THRESHOLD LEVEL.                            
      AMAX=ZERO                                                          
      J1=IPTR(I)+LENRL(I)                                                
      J2=IPTR(I)+LENR(I)-1                                               
      DO 220 JJ=J1,J2                                                    
 220  AMAX=DMAX1(AMAX,DABS(A(JJ)))                                       
CSP  220  AMAX=AMAX1(AMAX,ABS(A(JJ)))                                    
      AU=AMAX*U                                                          
C SCAN ROW FOR POSSIBLE PIVOTS                                           
      DO 230 JJ=J1,J2                                                    
      IF (DABS(A(JJ)).LE.AU.AND.L.EQ.1) GO TO 230                        
CSP       IF (ABS(A(JJ)).LE.AU.AND.L.EQ.1) GO TO 230                     
      J=ICN(JJ)                                                          
      KCOST=(NZ-1)*(LENC(J)-1)                                           
      IF (KCOST.GE.JCOST) GO TO 230                                      
C BEST PIVOT SO FAR IS FOUND.                                            
      JCOST=KCOST                                                        
      IJPOS=JJ                                                           
      IPIV=I                                                             
      JPIV=J                                                             
      IF (JCOST.LE.(NZ-1)**2) GO TO 380                                  
  230 CONTINUE                                                           
 235  CONTINUE                                                           
C                                                                        
C COLUMNS WITH NZ NON-ZEROS NOW EXAMINED.                                
 240  IJFIR=IFIRST(NZ)                                                   
      IJFIR=-LASTR(IJFIR)                                                
 245  IF (JCOST.LE.NZ*(NZ-1)) GO TO 380                                  
      DO 270 IDUMMY=1,N                                                  
      IF (IJFIR.EQ.0) GO TO 280                                          
      J=IJFIR                                                            
      IJFIR=NEXTC(IJFIR)                                                 
      I1=IPC(J)                                                          
      I2=I1+NZ-1                                                         
C SCAN COLUMN J.                                                         
      DO 260 II=I1,I2                                                    
      I=IRN(II)                                                          
      KCOST=(NZ-1)*(LENR(I)-LENRL(I)-1)                                  
      IF (KCOST.GE.JCOST) GO TO 260                                      
C PIVOT HAS BEST MARKOWITZ COUNT SO FAR ... NOW CHECK ITS                
C     SUITABILITY ON NUMERIC GROUNDS BY EXAMINING THE OTHER NON-ZEROS    
C     IN ITS ROW.                                                        
      J1=IPTR(I)+LENRL(I)                                                
      J2=IPTR(I)+LENR(I)-1                                               
C WE NEED A STABILITY CHECK ON SINGLETON COLUMNS BECAUSE OF POSSIBLE     
C     PROBLEMS WITH UNDERDETERMINED SYSTEMS.                             
      AMAX=ZERO                                                          
      DO 250 JJ=J1,J2                                                    
      AMAX=DMAX1(AMAX,DABS(A(JJ)))                                       
CSP       AMAX=AMAX1(AMAX,ABS(A(JJ)))                                    
  250 IF (ICN(JJ).EQ.J) JPOS=JJ                                          
      IF (DABS(A(JPOS)).LE.AMAX*U.AND.L.EQ.1) GO TO 260                  
CSP       IF (ABS(A(JPOS)).LE.AMAX*U.AND.L.EQ.1) GO TO 260               
      JCOST=KCOST                                                        
      IPIV=I                                                             
      JPIV=J                                                             
      IJPOS=JPOS                                                         
      IF (JCOST.LE.NZ*(NZ-1)) GO TO 380                                  
  260 CONTINUE                                                           
C                                                                        
  270 CONTINUE                                                           
C                                                                        
  280 CONTINUE                                                           
C                                                                        
C MATRIX IS NUMERICALLY OR STRUCTURALLY SINGULAR  ... WHICH IT IS WILL   
C     BE DIAGNOSED LATER.                                                
      IRANK=IRANK-1                                                      
  290 CONTINUE                                                           
C ASSIGN REST OF ROWS AND COLUMNS TO ORDERING ARRAY.                     
C MATRIX IS STRUCTURALLY SINGULAR.                                       
      IF (IFLAG.NE.2.AND.IFLAG.NE.-5) IFLAG=1                            
      IRANK=IRANK-N+PIVOT+1                                              
      IF (.NOT.ABORT1) GO TO 300                                         
      IDISP(2)=IACTIV                                                    
      IFLAG=-1                                                           
      IF (LP.NE.0) WRITE(LP,50)                                          
 50   FORMAT(' ERROR RETURN FROM MA30LA BECAUSE MATRIX IS STRUCTURALLY', 
     &' SINGULAR')                                                       
      GO TO 1110                                                         
  300 K=PIVOT-1                                                          
      DO 350 I=1,N                                                       
      IF (LASTR(I).NE.0) GO TO 350                                       
      K=K+1                                                              
      LASTR(I)=K                                                         
      IF (LENRL(I).EQ.0) GO TO 340                                       
      MINICN=MAX0(MINICN,NZROW+IBEG-1+MOREI+LENRL(I))                    
      IF (IACTIV-IBEG.GE.LENRL(I)) GO TO 320                             
      CALL MA30LD(A,ICN,IPTR(1),N,IACTIV,ITOP,.TRUE.)                    
C CHECK NOW TO SEE IF MA30LD HAS CREATED ENOUGH AVAILABLE SPACE.         
      IF (IACTIV-IBEG.GE.LENRL(I)) GO TO 320                             
C CREATE MORE SPACE BY DESTROYING PREVIOUSLY CREATED LU FACTORS.         
      MOREI=MOREI+IBEG-IDISP(1)                                          
      IBEG=IDISP(1)                                                      
      IF (LP.NE.0) WRITE(LP,310)                                         
      IFLAG=-5                                                           
  310 FORMAT(' LU DECOMPOSITION DESTROYED TO CREATE MORE SPACE')         
      IF (ABORT3) GO TO 1030                                             
  320 J1=IPTR(I)                                                         
      J2=J1+LENRL(I)-1                                                   
      IPTR(I)=0                                                          
      DO 330 JJ=J1,J2                                                    
      A(IBEG)=A(JJ)                                                      
      ICN(IBEG)=ICN(JJ)                                                  
      ICN(JJ)=0                                                          
  330 IBEG=IBEG+1                                                        
      NZROW=NZROW-LENRL(I)                                               
  340 IF (K.EQ.N) GO TO 360                                              
  350 CONTINUE                                                           
  360 K=PIVOT-1                                                          
      DO 370 I=1,N                                                       
      IF (LASTC(I).NE.0) GO TO 370                                       
      K=K+1                                                              
      LASTC(I)=-K                                                        
      IF (K.EQ.N) GO TO 940                                              
  370 CONTINUE                                                           
C                                                                        
C THE PIVOT HAS NOW BEEN FOUND IN POSITION (IPIV,JPIV) IN LOCATION       
C     IJPOS IN ROW FILE.                                                 
C UPDATE COLUMN AND ROW ORDERING ARRAYS TO CORRESPOND WITH REMOVAL       
C     OF THE ACTIVE PART OF THE MATRIX.                                  
  380 ISING=PIVOT                                                        
      IF (A(IJPOS).NE.ZERO) GO TO 390                                    
C NUMERICAL SINGULARITY IS RECORDED HERE.                                
      ISING=-ISING                                                       
      IF (IFLAG.NE.-5) IFLAG=2                                           
      IF (.NOT.ABORT2) GO TO 390                                         
      IDISP(2)=IACTIV                                                    
      IFLAG=-2                                                           
      IF (LP.NE.0) WRITE(LP,70)                                          
 70   FORMAT(' ERROR RETURN FROM MA30LA BECAUSE MATRIX IS NUMERICALLY',  
     &' SINGULAR')                                                       
      GO TO 1110                                                         
  390 OLDPIV=IPTR(IPIV)+LENRL(IPIV)                                      
      OLDEND=IPTR(IPIV)+LENR(IPIV)-1                                     
C CHANGES TO COLUMN ORDERING.                                            
      DO 460 JJ=OLDPIV,OLDEND                                            
      J=ICN(JJ)                                                          
      LC=LASTC(J)                                                        
      NC=NEXTC(J)                                                        
      IF (NC.NE.0) LASTC(NC)=LC                                          
      IF (LC.EQ.0) GO TO 440                                             
      NEXTC(LC)=NC                                                       
      GO TO 460                                                          
 440  NZ=LENC(J)                                                         
      ISW=IFIRST(NZ)                                                     
      IF (ISW.GT.0) LASTR(ISW)=-NC                                       
      IF (ISW.LT.0) IFIRST(NZ)=-NC                                       
  460 CONTINUE                                                           
C CHANGES TO ROW ORDERING.                                               
      I1=IPC(JPIV)                                                       
      I2=I1+LENC(JPIV)-1                                                 
      DO 530 II=I1,I2                                                    
      I=IRN(II)                                                          
      LR=LASTR(I)                                                        
      NR=NEXTR(I)                                                        
      IF (NR.NE.0) LASTR(NR)=LR                                          
      IF (LR.LE.0) GO TO 500                                             
      NEXTR(LR)=NR                                                       
      GO TO 530                                                          
 500  NZ=LENR(I)-LENRL(I)                                                
      IF (NR.NE.0) IFIRST(NZ)=NR                                         
      IF (NR.EQ.0) IFIRST(NZ)=LR                                         
  530 CONTINUE                                                           
C     RECORD THE COLUMN PERMUTATION IN LASTC(JPIV) AND THE ROW           
C     PERMUTATION IN LASTR(IPIV).                                        
      LASTC(JPIV)=ISING                                                  
      LASTR(IPIV)=PIVOT                                                  
C                                                                        
C MOVE PIVOT TO POSITION LENRL+1 IN PIVOT ROW AND MOVE PIVOT ROW         
C     TO THE BEGINNING OF THE AVAILABLE STORAGE.                         
C THE L PART AND THE PIVOT IN THE OLD COPY OF THE PIVOT ROW IS           
C     NULLIFIED WHILE, IN THE STRICTLY UPPER TRIANGULAR PART, THE        
C     COLUMN INDICES, J SAY, ARE OVERWRITTEN BY THE CORRESPONDING        
C     ELEMENT OF IQ (IQ(J)) AND IQ(J) IS SET TO THE NEGATIVE OF THE      
C     DISPLACEMENT OF THE COLUMN INDEX FROM THE PIVOT ELEMENT.           
      IF (OLDPIV.EQ.IJPOS) GO TO 540                                     
      AU=A(OLDPIV)                                                       
      A(OLDPIV)=A(IJPOS)                                                 
      A(IJPOS)=AU                                                        
      ICN(IJPOS)=ICN(OLDPIV)                                             
      ICN(OLDPIV)=JPIV                                                   
C CHECK TO SEE IF THERE IS SPACE IMMEDIATELY AVAILABLE IN A/ICN TO       
C     HOLD NEW COPY OF PIVOT ROW.                                        
  540 MINICN=MAX0(MINICN,NZROW+IBEG-1+MOREI+LENR(IPIV))                  
      IF (IACTIV-IBEG.GE.LENR(IPIV)) GO TO 550                           
      CALL MA30LD(A,ICN,IPTR(1),N,IACTIV,ITOP,.TRUE.)                    
      OLDPIV=IPTR(IPIV)+LENRL(IPIV)                                      
      OLDEND=IPTR(IPIV)+LENR(IPIV)-1                                     
C CHECK NOW TO SEE IF MA30LD HAS CREATED ENOUGH AVAILABLE SPACE.         
      IF (IACTIV-IBEG.GE.LENR(IPIV)) GO TO 550                           
C CREATE MORE SPACE BY DESTROYING PREVIOUSLY CREATED LU FACTORS.         
      MOREI=MOREI+IBEG-IDISP(1)                                          
      IBEG=IDISP(1)                                                      
      IF (LP.NE.0) WRITE(LP,310)                                         
      IFLAG=-5                                                           
      IF (ABORT3) GO TO 1030                                             
      IF (IACTIV-IBEG.GE.LENR(IPIV)) GO TO 550                           
C THERE IS STILL NOT ENOUGH ROOM IN A/ICN.                               
      IFLAG=-4                                                           
      GO TO 1030                                                         
C COPY PIVOT ROW AND SET UP IQ ARRAY.                                    
  550 IJPOS=0                                                            
      J1=IPTR(IPIV)                                                      
C                                                                        
      DO 570 JJ=J1,OLDEND                                                
      A(IBEG)=A(JJ)                                                      
      ICN(IBEG)=ICN(JJ)                                                  
      IF (IJPOS.NE.0) GO TO 560                                          
      IF (ICN(JJ).EQ.JPIV) IJPOS=IBEG                                    
      ICN(JJ)=0                                                          
      GO TO 570                                                          
  560 K=IBEG-IJPOS                                                       
      J=ICN(JJ)                                                          
      ICN(JJ)=IQ(J)                                                      
      IQ(J)=-K                                                           
  570 IBEG=IBEG+1                                                        
C                                                                        
      IJP1=IJPOS+1                                                       
      PIVEND=IBEG-1                                                      
      LENPIV=PIVEND-IJPOS                                                
      NZROW=NZROW-LENRL(IPIV)-1                                          
      IPTR(IPIV)=OLDPIV+1                                                
      IF (LENPIV.EQ.0) IPTR(IPIV)=0                                      
C                                                                        
C REMOVE PIVOT ROW (INCLUDING PIVOT) FROM COLUMN ORIENTED FILE.          
      DO 600 JJ=IJPOS,PIVEND                                             
      J=ICN(JJ)                                                          
      I1=IPC(J)                                                          
      LENC(J)=LENC(J)-1                                                  
C I2 IS LAST POSITION IN NEW COLUMN.                                     
      I2=IPC(J)+LENC(J)-1                                                
      IF (I2.LT.I1) GO TO 590                                            
      DO 580 II=I1,I2                                                    
      IF (IRN(II).NE.IPIV) GO TO 580                                     
      IRN(II)=IRN(I2+1)                                                  
      GO TO 590                                                          
  580 CONTINUE                                                           
  590 IRN(I2+1)=0                                                        
  600 CONTINUE                                                           
      NZCOL=NZCOL-LENPIV-1                                               
C                                                                        
C GO DOWN THE PIVOT COLUMN AND FOR EACH ROW WITH A NON-ZERO ADD          
C     THE APPROPRIATE MULTIPLE OF THE PIVOT ROW TO IT.                   
C WE LOOP ON THE NUMBER OF NON-ZEROS IN THE PIVOT COLUMN SINCE           
C     MA30L MAY CHANGE ITS ACTUAL POSITION.                              
C                                                                        
      NZPC=LENC(JPIV)                                                    
      IF (NZPC.EQ.0) GO TO 870                                           
      DO 820 III=1,NZPC                                                  
      II=IPC(JPIV)+III-1                                                 
      I=IRN(II)                                                          
C SEARCH ROW I FOR NON-ZERO TO BE ELIMINATED, CALCULATE MULTIPLIER,      
C     AND PLACE IT IN POSITION LENRL+1 IN ITS ROW.                       
      J1=IPTR(I)+LENRL(I)                                                
      IEND=IPTR(I)+LENR(I)-1                                             
      DO 610 JJ=J1,IEND                                                  
      IF (ICN(JJ).NE.JPIV) GO TO 610                                     
C IF PIVOT IS ZERO, REST OF COLUMN IS AND SO MULTIPLIER IS ZERO.         
      AU=ZERO                                                            
      IF (A(IJPOS).NE.ZERO) AU=-A(JJ)/A(IJPOS)                           
      A(JJ)=A(J1)                                                        
      A(J1)=AU                                                           
      ICN(JJ)=ICN(J1)                                                    
      ICN(J1)=JPIV                                                       
      LENRL(I)=LENRL(I)+1                                                
      GO TO 620                                                          
  610 CONTINUE                                                           
C GO TO 870 IF PIVOT ROW IS A SINGLETON.                                 
  620 IF (LENPIV.EQ.0) GO TO 820                                         
C NOW PERFORM NECESSARY OPERATIONS ON REST OF NON-PIVOT ROW I.           
      ROWI=J1+1                                                          
      IOP=0                                                              
C IF ALL THE PIVOT ROW CAUSES FILL-IN GO TO 640                          
      IF (ROWI.GT.IEND) GO TO 640                                        
C PERFORM OPERATIONS ON CURRENT NON-ZEROS IN ROW I.                      
C INNERMOST LOOP.                                                        
      DO 630 JJ=ROWI,IEND                                                
      J=ICN(JJ)                                                          
      IF (IQ(J).GT.0) GO TO 630                                          
      IOP=IOP+1                                                          
      PIVROW=IJPOS-IQ(J)                                                 
      A(JJ)=A(JJ)+AU*A(PIVROW)                                           
      ICN(PIVROW)=-ICN(PIVROW)                                           
  630 CONTINUE                                                           
  640 IFILL=LENPIV-IOP                                                   
C IF THERE IS NO FILL-IN GO TO 740.                                      
      IF (IFILL.EQ.0) GO TO 740                                          
C NOW FOR THE FILL-IN.                                                   
      MINICN=MAX0(MINICN,MOREI+IBEG-1+NZROW+IFILL+LENR(I))               
C SEE IF THERE IS ROOM FOR FILL-IN.                                      
C GET MAXIMUM SPACE FOR ROW I IN SITU.                                   
      DO 650 JDIFF=1,IFILL                                               
      JNPOS=IEND+JDIFF                                                   
      IF (JNPOS.GT.LICN) GO TO 660                                       
      IF (ICN(JNPOS).NE.0) GO TO 660                                     
  650 CONTINUE                                                           
C THERE IS ROOM FOR ALL THE FILL-IN AFTER THE END OF THE ROW SO IT       
C     CAN BE LEFT IN SITU.                                               
C NEXT AVAILABLE SPACE FOR FILL-IN.                                      
      IEND=IEND+1                                                        
      GO TO 740                                                          
C JMORE SPACES FOR FILL-IN ARE REQUIRED IN FRONT OF ROW.                 
  660 JMORE=IFILL-JDIFF+1                                                
      I1=IPTR(I)                                                         
C WE NOW LOOK IN FRONT OF THE ROW TO SEE IF THERE IS SPACE FOR           
C     THE REST OF THE FILL-IN.                                           
      DO 670 JDIFF=1,JMORE                                               
      JNPOS=I1-JDIFF                                                     
      IF (JNPOS.LT.IACTIV) GO TO 680                                     
      IF (ICN(JNPOS).NE.0) GO TO 690                                     
  670 CONTINUE                                                           
  680 JNPOS=I1-JMORE                                                     
      GO TO 700                                                          
C WHOLE ROW MUST BE MOVED TO THE BEGINNING OF AVAILABLE STORAGE.         
  690 JNPOS=IACTIV-LENR(I)-IFILL                                         
C IF THERE IS SPACE IMMEDIATELY AVAILABLE FOR THE SHIFTED ROW GO TO 720. 
  700 IF (JNPOS.GE.IBEG) GO TO 720                                       
      CALL MA30LD(A,ICN,IPTR(1),N,IACTIV,ITOP,.TRUE.)                    
      I1=IPTR(I)                                                         
      IEND=I1+LENR(I)-1                                                  
      JNPOS=IACTIV-LENR(I)-IFILL                                         
      IF (JNPOS.GE.IBEG) GO TO 720                                       
C NO SPACE AVAILABLE SO TRY TO CREATE SOME BY THROWING AWAY PREVIOUS     
C     LU DECOMPOSITION.                                                  
      MOREI=MOREI+IBEG-IDISP(1)-LENPIV-1                                 
      IF (LP.NE.0) WRITE(LP,310)                                         
      IFLAG=-5                                                           
      IF (ABORT3) GO TO 1030                                             
C KEEP RECORD OF CURRENT PIVOT ROW.                                      
      IBEG=IDISP(1)                                                      
      ICN(IBEG)=JPIV                                                     
      A(IBEG)=A(IJPOS)                                                   
      IJPOS=IBEG                                                         
      DO 710 JJ=IJP1,PIVEND                                              
      IBEG=IBEG+1                                                        
      A(IBEG)=A(JJ)                                                      
  710 ICN(IBEG)=ICN(JJ)                                                  
      IJP1=IJPOS+1                                                       
      PIVEND=IBEG                                                        
      IBEG=IBEG+1                                                        
      IF (JNPOS.GE.IBEG) GO TO 720                                       
C THIS STILL DOES NOT GIVE ENOUGH ROOM.                                  
      IFLAG=-4                                                           
      GO TO 1030                                                         
  720 IACTIV=MIN0(IACTIV,JNPOS)                                          
C MOVE NON-PIVOT ROW I.                                                  
      IPTR(I)=JNPOS                                                      
      DO 730 JJ=I1,IEND                                                  
      A(JNPOS)=A(JJ)                                                     
      ICN(JNPOS)=ICN(JJ)                                                 
      JNPOS=JNPOS+1                                                      
  730 ICN(JJ)=0                                                          
C FIRST NEW AVAILABLE SPACE.                                             
      IEND=JNPOS                                                         
  740 NZROW=NZROW+IFILL                                                  
C INNERMOST FILL-IN LOOP WHICH ALSO RESETS ICN.                          
      DO 810 JJ=IJP1,PIVEND                                              
      J=ICN(JJ)                                                          
      IF (J.LT.0) GO TO 800                                              
      A(IEND)=AU*A(JJ)                                                   
      ICN(IEND)=J                                                        
      IEND=IEND+1                                                        
C                                                                        
C PUT NEW ENTRY IN COLUMN FILE.                                          
      MINIRN=MAX0(MINIRN,NZCOL+LENC(J)+1)                                
      JEND=IPC(J)+LENC(J)                                                
      JROOM=NZPC-III+1+LENC(J)                                           
      IF (JEND.GT.LIRN) GO TO 750                                        
      IF (IRN(JEND).EQ.0) GO TO 790                                      
  750 IF (JROOM.LT.DISPC) GO TO 760                                      
C COMPRESS COLUMN FILE TO OBTAIN SPACE FOR NEW COPY OF COLUMN.           
      CALL MA30LD(A,IRN,IPC(1),N,DISPC,LIRN,.FALSE.)                     
      IF (JROOM.LT.DISPC) GO TO 760                                      
      JROOM=DISPC-1                                                      
      IF (JROOM.GE.LENC(J)+1) GO TO 760                                  
C COLUMN FILE IS NOT LARGE ENOUGH.                                       
      GO TO 1050                                                         
C COPY COLUMN TO BEGINNING OF FILE.                                      
  760 JBEG=IPC(J)                                                        
      JEND=IPC(J)+LENC(J)-1                                              
      JZERO=DISPC-1                                                      
      DISPC=DISPC-JROOM                                                  
      IDISPC=DISPC                                                       
      DO 770 II=JBEG,JEND                                                
      IRN(IDISPC)=IRN(II)                                                
      IRN(II)=0                                                          
  770 IDISPC=IDISPC+1                                                    
      IPC(J)=DISPC                                                       
      JEND=IDISPC                                                        
      DO 780 II=JEND,JZERO                                               
  780 IRN(II)=0                                                          
  790 IRN(JEND)=I                                                        
      NZCOL=NZCOL+1                                                      
      LENC(J)=LENC(J)+1                                                  
C END OF ADJUSTMENT TO COLUMN FILE.                                      
      GO TO 810                                                          
C                                                                        
  800 ICN(JJ)=-J                                                         
  810 CONTINUE                                                           
      LENR(I)=LENR(I)+IFILL                                              
C END OF SCAN OF PIVOT COLUMN.                                           
  820 CONTINUE                                                           
C                                                                        
C                                                                        
C REMOVE PIVOT COLUMN FROM COLUMN ORIENTED STORAGE AND UPDATE ROW        
C     ORDERING ARRAYS.                                                   
      I1=IPC(JPIV)                                                       
      I2=IPC(JPIV)+LENC(JPIV)-1                                          
      NZCOL=NZCOL-LENC(JPIV)                                             
      DO 860 II=I1,I2                                                    
      I=IRN(II)                                                          
      IRN(II)=0                                                          
      NZ=LENR(I)-LENRL(I)                                                
      IF (NZ.NE.0) GO TO 830                                             
      LASTR(I)=0                                                         
      GO TO 860                                                          
  830 IFIR=IFIRST(NZ)                                                    
      IFIRST(NZ)=I                                                       
      IF (IFIR) 840,855,850                                              
 840  LASTR(I)=IFIR                                                      
      NEXTR(I)=0                                                         
      GO TO 860                                                          
 850  LASTR(I)=LASTR(IFIR)                                               
      NEXTR(I)=IFIR                                                      
      LASTR(IFIR)=I                                                      
      GO TO 860                                                          
 855  LASTR(I)=0                                                         
      NEXTR(I)=0                                                         
      NZMIN=MIN0(NZMIN,NZ)                                               
 860  CONTINUE                                                           
C RESTORE IQ AND NULLIFY U PART OF OLD PIVOT ROW.                        
  870 IPC(JPIV)=0                                                        
      IF (LENPIV.EQ.0) GO TO 930                                         
      NZROW=NZROW-LENPIV                                                 
      JVAL=IJP1                                                          
      JZER=IPTR(IPIV)                                                    
      IPTR(IPIV)=0                                                       
      DO 880 JCOUNT=1,LENPIV                                             
      J=ICN(JVAL)                                                        
      IQ(J)=ICN(JZER)                                                    
      ICN(JZER)=0                                                        
      JVAL=JVAL+1                                                        
  880 JZER=JZER+1                                                        
C ADJUST COLUMN ORDERING ARRAYS.                                         
      DO 920 JJ=IJP1,PIVEND                                              
      J=ICN(JJ)                                                          
      NZ=LENC(J)                                                         
      IF (NZ.NE.0) GO TO 890                                             
      LASTC(J)=0                                                         
      GO TO 920                                                          
  890 IFIR=IFIRST(NZ)                                                    
      LASTC(J)=0                                                         
      IF (IFIR) 900,910,915                                              
 900  IFIRST(NZ)=-J                                                      
      IFIR=-IFIR                                                         
      LASTC(IFIR)=J                                                      
      NEXTC(J)=IFIR                                                      
      GO TO 920                                                          
 910  IFIRST(NZ)=-J                                                      
      NEXTC(J)=0                                                         
      NZMIN=MIN0(NZMIN,NZ)                                               
      GO TO 920                                                          
 915  LC=-LASTR(IFIR)                                                    
      LASTR(IFIR)=-J                                                     
      NEXTC(J)=LC                                                        
      IF (LC.NE.0) LASTC(LC)=J                                           
  920 CONTINUE                                                           
  930 CONTINUE                                                           
C ********************************************                           
C ****    END OF MAIN ELIMINATION LOOP    ****                           
C ********************************************                           
C                                                                        
C RESET IACTIV TO POINT TO THE BEGINNING OF THE NEXT BLOCK.              
  940 IF (N.NE.N) IACTIV=IPTR(N+1)                                       
  950 CONTINUE                                                           
C                                                                        
C ********************************************                           
C ****    END OF DEOMPOSITION OF MATRIX   ****                           
C ********************************************                           
C                                                                        
C                                                                        
C RUN THROUGH LU DECOMPOSITION CHANGING COLUMN INDICES TO THAT OF NEW    
C     ORDER AND PERMUTING LENR AND LENRL ARRAYS ACCORDING TO PIVOT       
C     PERMUTATIONS.                                                      
      IEND=IBEG-1                                                        
      DO 980 JJ=1,IEND                                                   
      JOLD=ICN(JJ)                                                       
  980 ICN(JJ)=LASTC(JOLD)                                                
      DO 990 II=1,N                                                      
      I=LASTR(II)                                                        
      NEXTR(I)=LENR(II)                                                  
  990 NEXTC(I)=LENRL(II)                                                 
      DO 1000 I=1,N                                                      
      LENRL(I)=NEXTC(I)                                                  
 1000 LENR(I)=NEXTR(I)                                                   
C                                                                        
C UPDATE PERMUTATION ARRAYS IP AND IQ.                                   
      DO 1010 II=1,N                                                     
      I=LASTR(II)                                                        
      J=LASTC(II)                                                        
      NEXTR(I)=IP(II)                                                    
 1010 NEXTC(J)=IQ(II)                                                    
      DO 1020 I=1,N                                                      
      IP(I)=NEXTR(I)                                                     
 1020 IQ(I)=NEXTC(I)                                                     
      IDISP(2)=IEND                                                      
      GO TO 1110                                                         
C                                                                        
C   ***    ERROR RETURNS    ***                                          
 1030 IDISP(2)=IACTIV                                                    
      IF (LP.EQ.0) GO TO 1110                                            
      WRITE(LP,1040)                                                     
 1040 FORMAT(' ERROR RETURN FROM MA30LA BECAUSE LICN NOT BIG ENOUGH')    
      GO TO 1080                                                         
 1050 IF (IFLAG.EQ.-5) IFLAG=-6                                          
      IF (IFLAG.NE.-6) IFLAG=-3                                          
      IDISP(2)=IACTIV                                                    
      IF (LP.EQ.0) GO TO 1110                                            
      IF (IFLAG.EQ.-3) WRITE(LP,1060)                                    
      IF (IFLAG.EQ.-6) WRITE(LP,1070)                                    
 1060 FORMAT(' ERROR RETURN FROM MA30LA BECAUSE LIRN NOT BIG ENOUGH')    
 1070 FORMAT(' ERROR RETURN FROM MA30LA LIRN AND LICN TOO SMALL')        
 1080 WRITE(LP,1090) PIVOT                                               
 1090 FORMAT(' AT STAGE ',I5)                                            
      IF (PIVOT.EQ.0) WRITE(LP,1100) MINIRN                              
 1100 FORMAT(' TO CONTINUE SET LIRN TO AT LEAST ',I8)                    
C                                                                        
C                                                                        
 1110 RETURN                                                             
C                                                                        
C  END MA30LA                                                            
C                                                                        
      END                                                                
      SUBROUTINE MA30LD(A,ICN,IPTR,N,IACTIV,ITOP,REALS)                  
C                                                                        
C  COMPRESS ARRAYS                                                       
C                                                                        
      DOUBLE PRECISION A(ITOP)                                           
CSP       REAL A(ITOP)                                                   
      LOGICAL REALS                                                      
C                                                                        
C FOR HANLING OF MORE THAN 32000 NONZEROS WITH INTEGER*2 VERSION CHANGE  
C ARRAY IPTR TO STANDARD INTEGER                                         
C                                                                        
CI4       INTEGER ICN(ITOP),IPTR(N)                                      
      INTEGER*2 ICN(ITOP),IPTR(N)                                        
      INTEGER J, JPOS, K, KL, KN                                         
      COMMON /MA30LF/ IRNCP,ICNCP,IRANK,MINIRN,MINICN 
      SAVE /MA30LF/
C IACTIV IS THE FIRST POSITION IN ARRAYS A/ICN FROM WHICH THE            
C     COMPRESS STARTS.                                                   
C ON EXIT IACTIV EQUALS THE POSITION OF THE FIRST ELEMENT IN THE         
C     COMPRESSED PART OF A/ICN                                           
      IF (REALS) ICNCP=ICNCP+1                                           
      IF (.NOT.REALS) IRNCP=IRNCP+1                                      
C SET THE FIRST NON-ZERO ELEMENT IN EACH ROW TO THE NEGATIVE OF THE      
C     ROW/COL NUMBER AND HOLD THIS ROW/COL INDEX IN THE ROW/COL          
C     POINTER.  THIS IS SO THAT THE BEGINNING OF EACH ROW/COL CAN        
C     BE RECOGNIZED IN THE SUBSEQUENT SCAN.                              
      DO 10 J=1,N                                                        
      K=IPTR(J)                                                          
      IF (K.LT.IACTIV) GO TO 10                                          
      IPTR(J)=ICN(K)                                                     
      ICN(K)=-J                                                          
   10 CONTINUE                                                           
      KN=ITOP+1                                                          
      KL=ITOP-IACTIV+1                                                   
C GO THROUGH ARRAYS IN REVERSE ORDER COMPRESSING TO THE BACK SO          
C     THAT THERE ARE NO ZEROS HELD IN POSITIONS IACTIV TO ITOP IN ICN.   
C     RESET FIRST ELEMENT OF EACH ROW/COL AND POINTER ARRAY IPTR.        
      DO 30 K=1,KL                                                       
      JPOS=ITOP-K+1                                                      
      IF (ICN(JPOS).EQ.0) GO TO 30                                       
      KN=KN-1                                                            
      IF (REALS) A(KN)=A(JPOS)                                           
      IF (ICN(JPOS).GE.0) GO TO 20                                       
C FIRST NON-ZERO OF ROW/COL HAS BEEN LOCATED                             
      J=-ICN(JPOS)                                                       
      ICN(JPOS)=IPTR(J)                                                  
      IPTR(J)=KN                                                         
   20 ICN(KN)=ICN(JPOS)                                                  
   30 CONTINUE                                                           
      IACTIV=KN                                                          
      RETURN                                                             
      END                                                                
      SUBROUTINE MA30LB(N,ICN,A,LICN,LENR,LENRL,IDISP,W,IW,              
     &IFLAG)                                                             
C                                                                        
C  FACTOR SPARSE MATRIX                                                  
C                                                                        
C                                                                        
C DATE OF LATEST CHANGE: APRIL 27, '82                                   
C                                                                        
      DOUBLE PRECISION A(LICN),W(N),AU,EPS,ROWMAX,ZERO,ONE,RMIN          
CSP       REAL A(LICN),W(N)                                              
      DOUBLE PRECISION DABS,DMAX1                                        
CSP       REAL ABS,AMAX1                                                 
      INTEGER IDISP(2),PIVPOS                                            
      INTEGER I, IFIN, ILEND, IPIVJ, ISING, ISTART, J, JAY, JAYJAY,      
     & JFIN, JJ                                                          
      LOGICAL ABORT1,ABORT2,ABORT3,STAB                                  
C                                                                        
C FOR HANLING OF MORE THAN 32000 NONZEROS WITH INTEGER*2 VERSION CHANGE  
C ARRAY IW TO STANDARD INTEGER                                           
C                                                                        
CI4       INTEGER ICN(LICN),LENR(N),LENRL(N),IW(N)                       
      INTEGER*2 ICN(LICN),LENR(N),LENRL(N),IW(N)                         
C                                                                        
      COMMON /MA30LE/ LP,ABORT1,ABORT2,ABORT3                            
      COMMON /MA30LG/ EPS,RMIN 
      SAVE /MA30LE/, /MA30LG/
C                                                                        
      DATA ZERO/0.0D0/,ONE/1.0D0/                                        
CSP       DATA ZERO/0.0/,ONE/1.0/                                        
C                                                                        
      STAB=EPS.LE.ONE                                                    
      RMIN=EPS                                                           
      ISING=0                                                            
      IFLAG=0                                                            
C 170 = RETURN...                                                        
      IF (N.EQ.1) GO TO 170                                              
      DO 10 I=1,N                                                        
   10 W(I)=ZERO                                                          
C SET UP POINTERS TO THE BEGINNING OF THE ROWS.                          
      IW(1)=IDISP(1)                                                     
      DO 20 I=2,N                                                        
   20 IW(I)=IW(I-1)+LENR(I-1)                                            
C                                                                        
C   ****   START  OF MAIN LOOP    ****                                   
C AT STEP I, ROW I OF A IS TRANSFORMED TO ROW I OF L/U BY ADDING         
C     APPROPRIATE MULTIPLES OF ROWS 1 TO I-1.                            
C     .... USING ROW-GAUSS ELIMINATION.                                  
      DO 140 I=1,N                                                       
C ISTART IS BEGINNING OF ROW I OF A AND ROW I OF L.                      
      ISTART=IW(I)                                                       
C IFIN IS END OF ROW I OF A AND ROW I OF U.                              
      IFIN=ISTART+LENR(I)-1                                              
C ILEND IS END OF ROW I OF L.                                            
      ILEND=ISTART+LENRL(I)-1                                            
      IF (ISTART.GT.ILEND) GO TO 70                                      
C LOAD ROW I OF A INTO VECTOR W.                                         
      DO 30 JJ=ISTART,IFIN                                               
      J=ICN(JJ)                                                          
   30 W(J)=A(JJ)                                                         
C                                                                        
C ADD MULTIPLES OF APPROPRIATE ROWS OF  I TO I-1  TO ROW I.              
      DO 50 JJ=ISTART,ILEND                                              
      J=ICN(JJ)                                                          
C IPIVJ IS POSITION OF PIVOT IN ROW J.                                   
      IPIVJ=IW(J)+LENRL(J)                                               
C FORM MULTIPLIER AU.                                                    
      AU=-W(J)/A(IPIVJ)                                                  
      W(J)=AU                                                            
C AU * ROW J (U PART) IS ADDED TO ROW I.                                 
      IPIVJ=IPIVJ+1                                                      
      JFIN=IW(J)+LENR(J)-1                                               
      IF (IPIVJ.GT.JFIN) GO TO 50                                        
C INNERMOST LOOP.                                                        
      DO 40 JAYJAY=IPIVJ,JFIN                                            
      JAY=ICN(JAYJAY)                                                    
   40 W(JAY)=W(JAY)+AU*A(JAYJAY)                                         
C                                                                        
   50 CONTINUE                                                           
C RELOAD W BACK INTO A (NOW L/U)                                         
      DO 60 JJ=ISTART,IFIN                                               
      J=ICN(JJ)                                                          
      A(JJ)=W(J)                                                         
   60 W(J)=ZERO                                                          
C WE NOW PERFORM THE STABILITY CHECKS.                                   
   70 PIVPOS=ILEND+1                                                     
C MATRIX HAD NON-ZERO PIVOT IN MA30LA AT THIS STAGE.                     
      IF (PIVPOS.GT.IFIN) GO TO 150                                      
      IF (A(PIVPOS).EQ.ZERO) GO TO 150                                   
      IF (.NOT.STAB) GO TO 140                                           
      ROWMAX=ZERO                                                        
      DO 130 JJ=PIVPOS,IFIN                                              
  130 ROWMAX=DMAX1(ROWMAX,DABS(A(JJ)))                                   
CSP   130 ROWMAX=AMAX1(ROWMAX,ABS(A(JJ)))                                
      IF (DABS(A(PIVPOS))/ROWMAX.GE.RMIN) GO TO 140                      
CSP       IF (ABS(A(PIVPOS))/ROWMAX.GE.RMIN) GO TO 140                   
      IFLAG=I                                                            
      RMIN=DABS(A(PIVPOS))/ROWMAX                                        
CSP       RMIN=ABS(A(PIVPOS))/ROWMAX                                     
C   ****    END OF MAIN LOOP    ****                                     
  140 CONTINUE                                                           
C                                                                        
      GO TO 170                                                          
C   ***   ERROR RETURN   ***                                             
  150 IF (LP.NE.0) WRITE(LP,160) I                                       
  160 FORMAT(' ERROR RETURN FROM MA30LB SINGULARITY DETECTED IN ROW',    
     1I8)                                                                
      IFLAG=-I                                                           
C                                                                        
  170 RETURN                                                             
C                                                                        
C  END MA30LB                                                            
C                                                                        
      END                                                                
      SUBROUTINE MA30LC (N,ICN,A,LICN,LENR,LENRL,IP,IQ,X,W)              
C                                                                        
C  PERFORMS FOR/BACKWARD SUBSTITUTION OF LINEAR SYSTEM                   
C                                                                        
C                                                                        
C DATE OF LATEST CHANGE: APRIL 27, '82                                   
C                                                                        
      DOUBLE PRECISION A(LICN),X(N),W(N),WII,WI,RESID,ZERO               
CSP       REAL A(LICN),X(N),W(N)                                         
CI4       INTEGER ICN(LICN),LENR(N),LENRL(N),IP(N),IQ(N)                 
      INTEGER*2 ICN(LICN),LENR(N),LENRL(N),IP(N),IQ(N)                   
      integer I, IBLOCK, IEND, II, III, J, J1, J2, JJ, JPIV, JPIVP1      
C                                                                        
      COMMON /MA30LH/ RESID 
      SAVE /MA30LH/
C                                                                        
      DATA ZERO/0.0D0/                                                   
CSP       DATA ZERO/0.0/                                                 
C THE FINAL VALUE OF RESID IS THE MAXIMUM RESIDUAL FOR AN INCONSISTENT   
C     SET OF EQUATIONS.                                                  
      RESID=ZERO                                                         
C                                                                        
C WE NOW SOLVE   A * X = B.                                              
C PREORDER VECTOR ... W(I) = X(IP(I))                                    
      DO 10 II=1,N                                                       
   10 W(II)=X(IP(II))                                                    
C IBLOCK HOLDS THE POSITION OF THE FIRST NON-ZERO IN THE                 
C     LU DECOMPOSITION (OF THE DIAGONAL BLOCKS).                         
      IBLOCK=1                                                           
C     A PASS THROUGH THIS LOOP                                           
C     PERFORMS FORWARD ELIMINATION USING ROW I OF THE LU                 
C     DECOMPOSITION.                                                     
      DO 120 I=1,N                                                       
      IF (LENRL(I).EQ.0) GO TO 50                                        
      WI=W(I)                                                            
C FORWARD ELIMINATION PHASE.                                             
C IEND IS THE END OF THE L PART OF ROW I IN THE LU DECOMPOSITION.        
      IEND=IBLOCK+LENRL(I)-1                                             
      DO 40 JJ=IBLOCK,IEND                                               
      J=ICN(JJ)                                                          
   40 WI=WI+A(JJ)*W(J)                                                   
      W(I)=WI                                                            
C IBLOCK IS ADJUSTED TO POINT TO THE START OF THE NEXT ROW.              
   50 IBLOCK=IBLOCK+LENR(I)                                              
 120  CONTINUE                                                           
C BACK SUBSTITUTION PHASE.                                               
C J1 IS POSITION IN A/ICN AFTER END OF A                                 
      J1=IBLOCK                                                          
C EACH PASS THROUGH THIS LOOP PERFORMS THE BACK-SUBSTITUTION             
C     OPERATIONS FOR A SINGLE ROW,                                       
C     WORKING THROUGH A IN REVERSE ORDER.                                
      DO 100 III=1,N                                                     
      II=N-III+1                                                         
C J2 IS END OF ROW II.                                                   
      J2=J1-1                                                            
C J1 IS BEGINNING OF ROW II.                                             
      J1=J1-LENR(II)                                                     
C JPIV IS THE POSITION OF THE PIVOT IN ROW II.                           
      JPIV=J1+LENRL(II)                                                  
      JPIVP1=JPIV+1                                                      
C IF ROW  II OF U HAS NO NON-ZEROS GO TO 90.                             
      IF (J2.LT.JPIVP1) GO TO 90                                         
      WII=W(II)                                                          
      DO 80 JJ=JPIVP1,J2                                                 
      J=ICN(JJ)                                                          
   80 WII=WII-A(JJ)*W(J)                                                 
      W(II)=WII                                                          
   90 W(II)=W(II)/A(JPIV)                                                
  100 CONTINUE                                                           
C                                                                        
C REORDER SOLUTION VECTOR ... X(I) = W(IQINVERSE(I))                     
      DO 130 II=1,N                                                      
      I=IQ(II)                                                           
  130 X(I)=W(II)                                                         
      RETURN                                                             
C                                                                        
C  END MA30LC                                                            
C                                                                        
      END                                                                
      SUBROUTINE MA30LM (N,NP1,NZLU,ICN,IA,JA,NZ,IMAP,LENR,IP,IQ,IW)     
C                                                                        
C  THIS ROUTINE GENERATES A MAPPING FROM THE ORIGINAL MATRIX TO THE      
C  FACTORS                                                               
C                                                                        
C                                                                        
C DATE OF LATEST CHANGE: APRIL 27, '82                                   
C                                                                        
C                                                                        
C FOR HANLING OF MORE THAN 32000 NONZEROS WITH INTEGER*2 VERSION CHANGE  
C ARRAY IW TO STANDARD INTEGER                                           
C                                                                        
CI4       INTEGER IW(N),IQ(N)                                            
      INTEGER*2 IW(N),IQ(N)                                              
CI4       INTEGER ICN(NZLU),IA(NP1),JA(NZ),IMAP(NZ),LENR(N),IP(N)        
      INTEGER*2 ICN(NZLU),IA(NP1),JA(NZ),IMAP(NZ),LENR(N),IP(N)          
      INTEGER INEW, IOLD, J, J1, J2, JAY1, JAY2, JJ, JNEW, JOLD          
      J1=1                                                               
      DO 100 INEW=1,N                                                    
      J2=J1+LENR(INEW)-1                                                 
      DO 30 JJ=J1,J2                                                     
      JNEW=ICN(JJ)                                                       
      JOLD=IQ(JNEW)                                                      
30    IW(JOLD)=JJ                                                        
      IOLD=IP(INEW)                                                      
      JAY1=IA(IOLD)                                                      
      JAY2=IA(IOLD+1)-1                                                  
      DO 40 JJ=JAY1,JAY2                                                 
      J=JA(JJ)                                                           
      IMAP(JJ)=IW(J)                                                     
40    CONTINUE                                                           
      J1=J2+1                                                            
100   CONTINUE                                                           
      RETURN                                                             
C                                                                        
C  END MA30LM                                                            
C                                                                        
      END                                                                
      DOUBLE PRECISION FUNCTION HEAT (T)                                 
CSP       FUNCTION HEAT (T)                                              
      DOUBLE PRECISION T
      HEAT=1.0D0                                                            
      RETURN                                                             
      END                                                                
      SUBROUTINE ENTROP (NSP,TEMP,GAS,SK,NEN,NSK,ENTH,SK0,NCO,IPTD,INFO) 
C*********************************************************************** 
C                                                                      * 
C DATE OF LATEST CHANGE: JANUARY 3, '86                                * 
C                                                                      * 
C*********************************************************************** 
C  ENTROP COMPUTES THE MOLAR ENTROPIES OF THE SPECIES                  * 
C  THE NASA-FITS OF THE THERMODYNAMICAL POLYNOMIALS, WHICH ARE READ    * 
C  AND PREPROCESSED IN SUBBLOCK 'COMPIL', ARE USED.                    * 
C*********************************************************************** 
C                                                                      * 
C  FOR DESCRIPTION OF PARAMETERS SEE SUBROUTINE 'METAS1'.              * 
C                                                                      * 
C*********************************************************************** 
C                                                                        
CI4       INTEGER IPTD(NCO)                                              
      INTEGER*2 IPTD(NCO)                                                
      DOUBLE PRECISION SK(NSP),ENTH(NEN),SK0(NSK)                        
CSP       REAL SK(NSP),ENTH(NEN),SK0(NSK)                                
      DOUBLE PRECISION A,B,C,C1,C2,C3,GAS,TLOG,TEMP,TEN3,TWO,ZERO        
CSP       REAL A,B,C,C1,C2,C3,GAS,TLOG,TEMP,TEN3,TWO,ZERO                
      INTEGER I,IPO,M1,N1                                                
C                                                                        
      DATA TEN3/1.D3/ , C1/1.25D0/ , C2/1.333334D0/ , C3/1.5D0/          
CSP       DATA TEN3/1.E3/ , C1/1.25E0/ , C2/1.333334E0/ , C3/1.5E0/      
      DATA TWO/2.D0/ , ZERO/0.D0/                                        
CSP       DATA TWO/2.E0/ , ZERO/0.E0/                                    
C                                                                        
      INFO=0                                                             
      IF(TEMP.LE.ZERO) GOTO 999                                          
      TLOG=DLOG(TEMP)                                                    
CSP       TLOG=ALOG(TEMP)                                                
      DO 3 I=1,NCO                                                       
      IPO=IPTD(I)                                                        
      N1=12*I                                                            
      M1=2*I                                                             
      IF(TEMP.LE.TEN3) GOTO 2                                            
      N1=N1-6                                                            
      M1=M1-1                                                            
2     A=ENTH(N1) * C1                                                    
      B=ENTH(N1-1)*C2                                                    
      C=ENTH(N1-2)*C3                                                    
      A=B + A*TEMP                                                       
      A=C + A*TEMP                                                       
      A=ENTH(N1-3)*TWO + A*TEMP                                          
      SK(IPO)=( ENTH(N1-4)*TLOG + SK0(M1) + A*TEMP )*GAS                 
  3   CONTINUE                                                           
      RETURN                                                             
C                                                                        
999   INFO=-1                                                            
      RETURN                                                             
C                                                                        
C----------------------------------------------------------------------- 
C                 END OF SUBROUTINE ENTROP                               
C----------------------------------------------------------------------- 
      END 
      SUBROUTINE DPROF (T,DENSY,DDOT)
      DOUBLE PRECISION T,DENSY,DDOT
      RETURN
      END
      SUBROUTINE PPROF (T,PRESS,PDOT)
      DOUBLE PRECISION T,PRESS,PDOT
      RETURN
      END
