C       PROGRAM MASS2
C
C  (INPUT,OUTPUT,TAPE5=INPUT,TAPE6=OUTPUT)
C
C  COMPUTE EXACT (NON-EQUILIBRIUM) RATES FOR INTERMEDIATE
C  AND (IF REQUIRED) LOW-PRESSURE REGIME BY
C  SOLUTION OF THE MASTER EQUATION (LARGEST EIGENVALUE ONLY).
C  THE GAS/GAS COLLISION PROBABILITY IS ASSUMED TO BE A FUNCTION
C  OF THE ENERGY DIFFERENCE OF THE INITIAL AND FINAL STATES,WITH A
C  PARAMETER WHICH MAY BE A FUNCTION OF THE INITIAL ENERGY.
C
C  ANGULAR MOMENTUM CONSERVATION IN THE FALL-OFF AND LOW-PRESSURE
C  REGIMES IS TREATED USING THE METHOD OF SMITH AND GILBERT.
C
C  THIS PROGRAM WILL HANDLE UP TO THREE REACTION CHANNELS.
C
C  OPTION PARAMETER IOPTHT ENABLES ONE TO DO A FULL FALL-OFF
C  CALCULATION AND, IF REQUIRED, THE LOW-PRESSURE LIMIT.
C  DETERMINE EIGENVALUE BY NESBET METHOD.
C
C  THE COMPLETE INPUT FILE FOR THIS PROGRAM IS PREPARED BY RRKM
C
C 
C*************************************************************************
C  COSMETIC IMPROVEMENTS TO PROGRAM: COMMENT LINES, IF BLOCKS AND DO LOOPS 
C  PUTTING LOGICAL, INTEGER AND DOUBLE PRECISION VARIABLES INTO SEPARATE
C  COMMON BLOCKS
C  REMOVE SOME OF THE ARGUMENTS IN CALLS INTO COMMON BLOCKS
C  PUTTING EXTRA CHECKS TO AVOID TAKING LOGARITHM OF ZERO 
C  PRINT OUT K0*[M] RATHER THAN INCORRECT K0*OMEGA
C  DON'T PRINT ANYTHING INTO COLUMN 1
C**************************************************************************
C
C  ALTERATIONS BY Gary P Knight 9th September 1994.
C
C  SINCE RRKM NEEDS AN INC VALUE OF LESS THAN 100 cm-1 FOR SOME ACCURATE
C  CALCULATIONS AT LOW TEMPERATUES. THIS PROGRAM WILL NOW WORK WITH
C  INC=100 cm-1 OR 10 cm-1.
C  THE VARIABLE TESTINC TAKES THE VALUES 1 AND 10 RESPECTIVELY TO 
C  ACHIEVE THIS.
C  THIS NEW VERSION OF MASNEW (MAS55) IS DESIGNED TO WORK WITH ARRAYS
C  SET TO 10x BIGGER. (NOW PARAMETERISED)
C  NOTE THE RESULTING .EXE FILE WILL BE VERY LARGE.
C
C  THE PROGRAM LOOPS HAVE BEEN TRUNCATED FOR INC=100 cm-1, USING THE 
C  VARIABLE TESTINC TO SPEED UP THE CALCULATION.
C
       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
       DOUBLE PRECISION KT
       PARAMETER (NMAX100=10000, NMAX50=5000, NMAX25=2500)
       COMMON /AMIN/ ID,NCUT,NILOC
       COMMON /BLOCK1/ NRATES,NI,NE0,NEFF
       COMMON /BLOCK2/ DELT,CORRF
       COMMON /BLOCK3/ NREACT,NBAND,NRPL1
       COMMON /BLOCK6/ Q(NMAX100),H(NMAX100),RTRHB(NMAX100)
       COMMON /BLOCK7/ D(NMAX100)
       COMMON /BLK10/ AVRATE(NMAX50)
       COMMON /BLK12/ AVR2(NMAX50,2)
       COMMON /BLK13/ AROW(NMAX100),OMEGA,ANORM(NMAX100)
       COMMON /BLK14/ NRCTH
       COMMON /BLK15/ NCH2
       COMMON /B2/ PROB(NMAX25)
       COMMON /DEOC/ EI,WE0
       COMMON /DENS/ RHO(NMAX100)
       COMMON /EXP/ ERR1,ERR2,ERR3
       COMMON /F1/ IFIRST,NSPEC,IMIN
       COMMON /F2/ ALMT,ALMT1,DET
       COMMON /JLP1/ JAVX
       COMMON /JLP2/ CRF2(NMAX100),EXL2(NMAX100)
       COMMON /JLP3/ CRF(NMAX100),EXL(NMAX100)
       COMMON /LOW/ LOWP
       COMMON /OPT/ IOPTPR
       COMMON /ROTTX/ ROTHR(NMAX100),ROTKT(NMAX100)
       COMMON /TRANS1/ RATTOT,RATE1,RATE2,RATEA(2)
       COMMON /TRANS2/ ITR
       COMMON /VERSCH/ IXV,JXV
       COMMON /NEWKIDS/ TESTINC
C
       DIMENSION PR(20),R1(NMAX100),R2(NMAX100,2),R1S(NMAX100),
     *           ALPHAV(10),TITLE(20),TEMPV(10),CORRAT(10),RATHP(3)
C
       LOGICAL LOWP,JAVX,BRW
C
       DATA R1/NMAX100*0./,R2/20000*0./,R1S/NMAX100*0./,
     *      WKA/349.689/,RATHP/3*0./
C
C  READ IN AND PRINT OUT TITLE.
C
       READ(5,1) TITLE
C
   1   FORMAT(1X,20A4)
C
       WRITE(6,2) TITLE
C
   2   FORMAT(1H1,20X,20A4,/,
     * 10X,'MASTER EQUATION SOLUTION FOR FALL-OFF AND LOW PRESSURE',
     * ' REGIMES',/,
     * 10X,'LATEST UPDATE: Sept 1994')
C
C  READ IN INITIAL GRAIN SIZE, NUMBER OF TIMES (+1) THAT
C  ONE WANTS TO HALVE THE GRAIN SIZE, STARTING AT A GRAIN SIZE (CM-1) WHICH
C  MUST BE 100,200,400,800,1600,...  CM-1
C  OR 10,20,40,80,160,... CM-1 (GPK 9/9/94)
C
C  B IS A PARAMETER, NORMALLY BETWEEN 0.01 AND 0.9, WHICH HELPS
C  THE EIGENVALUE ROUTINE FROM GOING WILD; B GT 0.1 MEANS SLOWER
C  CONVERGENCE BUT LESS LIKELIHOOD OF NUMERICAL INSTABILITY.
C  A VALUE OF 0.1 GIVES BEST CONVERGENCE PROPERTIES, BUT THIS
C  MAY BE TESTED FOR EACH CASE BY VARYING B BETWEEN 0.01 AND 1.
C
       READ(5,*) INC,NCHAN,INCCHK
       IDM=2
       NCH2=MAX0(NCHAN,2)-1
C
C  CHANGED TO DELTS=DFLT(INCCHK) INSTEAD OF INC WHICH IS ALWAYS 400
C
       DELTS=DFLT(INCCHK)
       TESTINC=1.0
       IF(INCCHK .EQ. 10) TESTINC=10.
       IF(INCCHK .EQ. 100) GO TO 100
       IF(INCCHK .EQ. 10) GO TO 100       
       WRITE(6,3)
C
   3   FORMAT(' RRKM PROGRAM WAS NOT CALCULATED 
     * WITH INC=10 OR 100 - ABORT')
C
       STOP
 100   CONTINUE
C
C  THE FOLLOWING IS THE FIXED VALUE OF THE CONVERGENCE-GOVERNING PARAMETER.
C
       B=0.999
C
C  READ IN ERRORS
C
       READ(5,*) ERR1,ERR2,ERR3
       WRITE(6,4) ERR1,ERR2,ERR3
C
   4   FORMAT(1P,/,10X,'ERROR FOR TRUNCATION OF MATRIX',T50,E10.1,
     * /,10X,'EIGENVALUE TOLERANCE',T50,E10.1,/,
     * 10X,'ERROR FOR TRUNCATION OF P(E,E',1H',')',T50,E10.1)
C
       WRITE(6,5) NCHAN
C
   5   FORMAT(/,10X,'NUMBER OF REACTION CHANNELS',T58,I2)
C
       ITMAX=800
C
C  ITMAX IS THE MAXIMUM NUMBER OF NESBET ITERATIONS PERMITTED.
C
       READ(5,*) E0
       E0T=E0
       WE0=E0T*WKA
C
C  E0 IS THE CRITICAL ENERGY; THE MAXIMUM ENERGY CONSIDERED IS COMPUTED
C  BELOW AS THE ENERGY OF THE HIGHEST INPUT MICROSCOPIC RATE.
C  NALPHA IS THE NUMBER OF ALPHA VALUES,AND ALPHAV THE VECTOR OF THEIR
C  VALUES.
C
       READ(5,*) NALPHA
       READ(5,*) (ALPHAV(I),I=1,NALPHA)
C
C  READ IN PARAMETERS GIVING THE
C  FUNCTIONAL FORM OF P(E,E').
C
       READ(5,*) IXV,JXV
C
C  IF THE FIRST VALUE OF ALPHA IS GT 0, AN EXPONENTIAL MODEL
C  FOR P(E,E') IS ASSUMED. IF FIRST ALPHA VALUE LT 0, THEN ADDITIONAL
C  PARAMETERS MUST BE READ IN TO SPECIFIY FUNCTIONAL FORM OF P.
C  READ IN PRESSURES (TORR) AT WHICH CALCS TO BE PERFORMED
C
       READ(5,*) NP
       READ(5,*) (PR(I),I=1,NP)
C
C  IF JAV=0, SINGLE FILE OF K(E)'S WITH AVERAGE (I+/I) J CORRECTION
C  FACTOR ASSUMED.
C  IF JAV.NE.0, K(E)= ASSUMED.(SEPARATE K(E) FILE FOR EACH
C  PRESSURE AND TEMPERATURE).
C
       READ(5,*) PALMT
C
C  PALMT IS THE PARAMETER SUCH THAT THE POPULATION VECTOR FOR
C  ALL ENERGIES LT E0*PALMT ARE GIVEN THEIR EQUILIBRIUM VALUES.
C
       READ(5,*)JAV
       JAVX = JAV.NE.0
       IF(JAVX) WRITE(6,7)
C
   7   FORMAT(10X,'FULL ANGULAR MOMENTUM CONSERVATION USED (SMITH-',
     * 'GILBERT METHOD)')
C
       IF(.NOT.JAVX) WRITE(6,8)
C
   8   FORMAT(10X,'THIS CALCULATION DOES NOT USE FULL ANGULAR MOMENTUM',
     * ' CONSERVATION;',/,10X,'ANGULAR MOMENTUM CONSERVATION ONLY',
     * ' IN HIGH PRESSURE LIMIT')
C
       IF((NALPHA.EQ.1).OR.(NP.EQ.1).OR.(.NOT.JAVX))GO TO 110
       WRITE(6,9)
C
   9   FORMAT(/,' DUE TO THE STRUCTURE OF THE INPUT FILE,',/,
     * ' YOU MAY NOT LOOP SEVERAL ALPHAS AND SEVERAL PRESSURES ',
     * ' IN THE SAME RUN.',/,' MAKE EITHER NALPHA OR NP =1',/,' ABORT.')
C
       STOP
 110   CONTINUE
       READ(5,*)NTEMP,(TEMPV(I),I=1,NTEMP)
       IF(.NOT.JAVX)GO TO 120
       IF(PR(1).GE.0)GO TO 120
       WRITE(6,10)
C
  10   FORMAT(/,' OPTION TO GENERATE FALL-OFF AUTOMATICALLY IS NOT',
     * ' VALID WHEN USING J-AVERAGED K(E)S. ABORT')
C
       STOP
 120   CONTINUE
C
C  READ IN COLLISION DIAMETER (ANGSTROM),COLLISION MASSES
C
       READ(5,*) SGMA,WT1,WT2,EPS
       WRITE(6,11) SGMA,WT1,WT2
C
  11   FORMAT(/,10X,'COLLISION DIAMETER =',T50,F10.2,' ANGSTROMS',
     * /,10X,'MASS OF REACTANT=',T50,F10.2,' A.M.U.'
     *,/,10X,'MASS OF BATH GAS  ='
     * ,T50,F10.2,' A.M.U.',/)
C
       IF(EPS.LE.0.) WRITE(6,12)
C
  12   FORMAT(10X,'HARD SPHERE COLLISION FREQUENCY USED',/)
C
       IF(EPS.GT.0.) WRITE(6,13) EPS
C
  13   FORMAT(10X,'LENNARD JONES COLLISION FREQUENCY USED',/,
     * 10X,'WELL DEPTH =',T50,F10.1,' K',/)
C
       AMASS=1./((1./WT1)+(1./WT2))
C
C  READ IN OPTIONS FOR DOING   CALCULATION
C  THE OPTION PARAMETERS ARE AS FOLLOWS:
C  IOPTHT = 0, DOES LOW PRESSURE CALCULATION IN ADDITION TO
C  FALL-OFF CALCULATION, IOPT.NE.0, DOES FALL-OFF CALC. ONLY.
C  N.B., IF LOW PRESSURE CALCULATION DONE, LOWEST GRAIN SIZE IS MADE 100 CM-1.
C
C  IOPTPR     = 1 OR -1  WEIGHTS CONVERGENCE TESTS BY OVERALL RATE
C                    (SUITABLE FOR 1-CHANNEL REACTION)
C         = 2 OR -2  WEIGHTS CONVEGENCE TESTS FOR BOTH CHANNELS EQUALLY
C                    (SUITABLE FOR 2-CHANNEL SYSTEM).
C
       READ(5,*) IOPTHT,IOPTPR
       IOPTPR=IABS(IOPTPR)
C
C  READ IN NUMBER OF DENSITIES OF STATES THEN THE LIST OF DENSITIES OF STATES;
C  IT IS ASSUMED THAT THESE HAVE BEEN COMPUTED WITH A SPACING OF 100 CM-1.
C
       WRITE(6,14) NP
C
  14   FORMAT(10X,'NO. OF PRESSURES =',T50,I10,/)
C
C  FOR LOW PRESSURE CALCULATION IDM=2 DELTS = 200 OR 20
C
       IDM=2
       DELTS=200./TESTINC
       INC=DELTS
       WRITE(6,15) INC
C
  15   FORMAT(10X,'INITIAL GRAIN SIZE =',T50,I10,' CM-1')
C
       WRITE(6,16) NALPHA,IDM
C
  16   FORMAT(/,10X,'NO. OF VALUES OF ALPHA =',T50,I10,/,
     * 10X,'NO. OF TIMES,+1,GRAIN SIZE HALVED =',T50,I10)
C
       BRW=IXV.LE.0.AND.JXV.LE.0
C
C  BRW IS LOGICAL VARIABLE SPECIFIYING THAT BIASED RANDOM WALK MODEL
C  USED FOR P(E,E').
C
       WRITE(6,17) PALMT
C
  17   FORMAT(10X,'CUTOFF PARAMETER FOR H(E)=1',T50,F10.2)
C
       I=IABS(IOPTPR)
       WRITE(6,18) I
C
  18   FORMAT(/,10X,'NUMBER OF REACTION CHANNELS =',T50,I10)
C
       READ(5,*) NDEGS
       IF(NDEGS.LE.NMAX100) GO TO 140
       WRITE(6,20) NMAX100
C
  20   FORMAT(' MORE THAN ',I6,' ELEMENTS IN RHO INPUT,ABORT')
C
       STOP
 140   CONTINUE
       READ(5,*) (RHO(I),I=1,NDEGS)
C
C  RENORMALIZE RHO
C
       DO I=1,NDEGS
         IF(RHO(I).GT.0.) GO TO 150
       ENDDO
 150   DNM=RHO(I)
       DO I=1,NDEGS
         RHO(I)=RHO(I)/DNM
       ENDDO
       IF(JAVX)GO TO 160
       DO I=1,NTEMP
         READ(5,*)T,CORRAT(I)
       ENDDO
       NE0=INT(WE0*TESTINC/100.)
       GO TO 170
 160   CONTINUE
C
C  READ IN ARRAY OF ROTATIONAL THRESHOLD ENERGIES.
C
       READ(5,*)NTHR
       READ(5,*)(ROTHR(I),I=1,NTHR)
       NE0=NTHR
       NEON2=(NE0/2)-1
 170  CONTINUE
C
C  INCREASE NUMBER OF PRESSURES CONSIDERED BY 1 SO THAT THE LOW PRESSURE
C  CALCULATION CAN BE PERFORMED THE FIRST TIME THE PRESSURE LOOP IS 
C  ACCESSED FOR EACH TEMP AND EACH ALPHA
C
       NP=NP+1
       DO 400 ITEMPR=1,NTEMP
         T=TEMPV(ITEMPR)
         WRITE(6,23) T
C
  23   FORMAT(/////,10X,'TEMPERATURE =',T50,F10.1,' KELVINS',//)
C
         KT=T/1.43879
         IF(.NOT.JAVX)GO TO 180
         DO I=NEON2,NE0
           ROTKT(I) = ROTHR(I)/KT
           CRF(I)=0.
           IF(ROTKT(I).LE.30.) CRF(I)=DEXP(-ROTKT(I))
         ENDDO
         DO I=NE0+1,NDEGS
           CRF(I)=1.
           ROTKT(I) = 0.
         ENDDO
 180     CONTINUE
         ICHP=0
         P1=PR(1)
C
C  LOOP THROUGH ALL THE COLLISIONAL ENERGY TRANSFER PARAMETERS CONSIDERED
C
         DO 399 IALPHA=1,NALPHA
           L=-1-(NP/2)
           BETA=1.
           ALPHA=ALPHAV(IALPHA)
           IF(.NOT.BRW) WRITE(6,24) ALPHA,IXV,JXV
C
  24   FORMAT(///,T30,' ALPHA =',F8.0,' CM-1',
     * //,10X,'PROBABILITY FUNCTION IS R(E-E',1H',')/N(E',1H',
     * '), WHERE',/,8X,12HR(E-E')=(X**,I2,
     * 10H) EXP(-X**,I2,21H), WHERE X=E-E'/ALPHA)
C
           IPRT=0
           DO 398 N=1,NP
C
C  DO LOW PRESSURE CALCULATION FIRST TIME THROUGH  PRESSURE LOOP
C
             LOWP=N.EQ.1
             IF(JAVX)GO TO 190
             IF((N.NE.1).OR.(IALPHA.NE.1).OR.(ITEMPR.NE.1))GO TO 220
             GO TO 200
 190         CONTINUE
             IF((IALPHA.NE.1).OR.(N.EQ.2))GO TO 220
             READ(5,*)(RATHP(I),I=1,NCHAN),CORAV,CORPF
             READ(5,*)STLPJ
             READ(5,*) E0EFF
             E0=E0EFF
 200         CONTINUE
C
C  READ IN NUMBER OF MICROSCOPIC RATES TO BE READ IN, AND THEN THE ACTUAL
C  RATES; IT IS ASSUMED THAT THESE HAVE BEEN COMPUTED WITH A SPACING OF
C  EITHER 100 OR 10 CM-1.
C
C  MICROSCOPIC RATES AND DENSITIES OF STATES,FOR EITHER 10 CM-1 OR
C  100 CM-1 SPACING,ARE READ IN R1 (K1), R2 (K2) AND RHO (RHO).
C
             READ(5,*) NRATES
             READ(5,*) (R1(I),I=1,NRATES)
             DO I=1,NCH2
               READ(5,*) (R2(II,I),II=1,NRATES)
             ENDDO
             WRITE(6,25) NRATES,NDEGS
C
  25   FORMAT(/,10X,'NUMBER OF INPUT MICROSCOPIC RATES =',T50,I10,/,
     * 10X,'NUMBER OF INPUT DENSITIES OF STATES =',T50,I10)
C
             IF(NRATES.LE.NMAX50) GO TO 210
             WRITE(6,26) NMAX50,NMAX50
C
  26   FORMAT(' MAXIMUM OF ',I10,' INPUT RATES ONLY ALLOWED, ABORT.',
     * /,' INCREASE DIMENSION FROM',I10,' TO RUN.')
C
             STOP
 210         CONTINUE
C
C  MAXIMUM ENERGY CONSIDERED CORRESPONDS WITH THE HIGHEST INPUT RATE
C
             NRTSTM=NRATES
             EMAX=E0+(0.2859*DFLT(NRATES)/TESTINC)
             NIORIG=INT((EMAX*200./DELTS)/0.5719)-1
C
C  NI IS NUMBER OF LEVELS ALTOGETHER
C  NREACT IS NUMBER NOT REACTING
C  DELT IS ENERGY DIFFERENCE BETWEEN LEVELS IN CM-1
C  NBAND IS THE BANDWIDTH OF THE MATRIX
C
             EJM=EMAX*4.184
             ECM=E0*4.184
             WRITE(6,27) EMAX,EJM,E0,ECM
C
  27   FORMAT(//,10X,'MAXIMUM ENERGY CONSIDERED',T50,F10.1,' KCAL MOL-1'
     * ,/,T45,'=',T50,F10.1,' KJ MOL-1'
     * ,//,10X,'CRITICAL ENERGY',T50,F10.2,' KCAL MOL-1',/,
     * 45X,'=',T50,F10.2,' KJ MOL-1')
C
 220         CONTINUE
             NI=NIORIG
             PN=PR(MAX0(1,N-1))
             NEFF=INT(E0EFF*WKA*TESTINC/100.)
             IF(P1.LE.0.) PN=1.
             L=L+1
             IF(.NOT.LOWP.AND.P1.LE.0.) PN=CSET*(10.**L)
             OMEGA=(4.412E+7*SGMA*SGMA/DSQRT(AMASS))*PN/(DSQRT(T))
C
C  CALCULATE LENNARD-JONES COLLISION FREQUENCY IF REQUIRED.
C
             IF(EPS.GT.0.) OMEGA=OMEGA/(0.636+0.567*DLOG10(T/EPS))
C
C  THE FOLLOWING SET OF STATEMENTS IS FOR THE SPECIAL CASE OF THE
C  LOW PRESSURE LIMIT CALCULATION, IE IF IOPTHT.NE.0 THEN THE LOW
C  PRESSURE CALCULATION IS NOT PERFORMED.
C
             IF(IOPTHT.NE.0.AND.LOWP) GO TO 398
C
C  CALL SUBROUTINE TO CARRY OUT PREPARATION FOR LOW-PRESSURE CALCULATION.
C
             IF(JAVX.AND.(.NOT.LOWP))GO TO 230
             E0=E0T
             GO TO 240
 230         E0=E0EFF
 240         CONTINUE
             IF (LOWP) CALL SETUPL(E0,T,R1S,DELTS,NDEGS,ALPHA,IPRT,
     *                                                      IOPTEM)
             DELT=DELTS*2.
C
C  END OF SEPARATE SECTION FOR LOW PRESSURE LIMIT.
C  ALL INITIALIZATION HAS NOW BEEN COMPLETED.
C  --------------------------------------------------------------------
C  START LOOP OVER GRAIN SIZES DELTA E
C
             NCUT=2
             IF(LOWP) NCUT=1
C
C  LOOP OVER THE NUMBER OF TIMES THAT THE GRAIN SIZE, DELT, IS HALVED
C
             DO 397 IDT=1,IDM
               ID=IDT
               DELT=DELT/2.
               IF((.NOT.LOWP).OR.(ID.EQ.1))GO TO 250
               NE01=NE0*INT(101./(DELT*TESTINC))
               NI=NE01
 250           CONTINUE
               NREACT=INT(E0/(DELT*2.859E-3))
               IF(LOWP)NREACT=NE0*100/(TESTINC*2*INT(DELT+0.1))
C
C  PALMT IS THE PARAMETER SUCH THAT THE POPULATION VECTOR FOR
C  ALL ENERGIES LT E0*PALMT ARE GIVEN THEIR EQUILIBRIUM VALUES.
C
               ALMT1=PALMT*E0/(DELT*2.859E-3)
               EXPON=1.43879*DELT/T
               IF(NI.LE.(1000*TESTINC)) GO TO 260
               WRITE(6,28)
C
  28   FORMAT(' NI TRUNCATED ')
C
               NI=1000*TESTINC
 260           CONTINUE
               SUMSCT=0.
               SUMSC2=0.
C
C  GENDEG AND RATES SORT THE DENSITIES OF STATES (RHO) AND THE MICROSCOPIC
C  RATE COEFFICIENTS (R1,R2) - FOR REACTING LEVELS - INTO NEW ARRAYS D AND 
C  AVRATE WHERE ENERGY SPACING BETWEEN ELEMENTS IS 100 CM-1.
C  BOTH SUBROUTINES INTERPOLATE OR TABLE SELECT
C
               DO I=1,(500*TESTINC)
                 AVRATE(I)=0.0
               ENDDO
               IF(.NOT.LOWP) GO TO 270
               CALL RATES(R1S,R2)
C
C  GENCF SORTS ROTATIONAL THRESHOLDS, DEXP(-R0/KT), INTO EITHER
C  10 OR 100 CM-1 SPACINGS DEPENDING ON VALUE OF TESTINC
C
               IF(JAVX)CALL GENCF
               GO TO 280
 270           CALL RATES(R1,R2)
 280           CONTINUE
               CALL GENDEG(NDEGS)
C
C  COMPUTE STRONG COLLISION RATES (SUMSC1 AND SUMSC2) AND (FOR THE FIRST
C  VALUE OF DELTA E) THE STRONG COLLISION VECTOR AS INITIAL GUESS TO
C  EIGENVECTOR Q(I).
C
               ALMT=0.
               DENST=0.
               IF(LOWP) BETA=1.
               DO I=1,NI
                 CI2=0. 
                 IF(D(I).GT.0.) CI2=DEXP(DLOG(D(I))-DFLT(I)*EXPON)
                 DENST=DENST+CI2
                 RTRHB(I)=DSQRT(CI2)
C
C  RTRHB IS BOLTZMANN VECTOR (SYMMETRIZED),Q IS EIGENVECTOR.
C  FOR ENERGIES BELOW NREACT Q(I) IS ASSIGNED THE VALUE RTRHB(I)
C
                 IF(ID.NE.1)GO TO 290
                 Q(I)=RTRHB(I)
C
C  EXTRA SYMMETRIZATION IS REQUIRED FOR LOW PRESSURE J-AVERAGED
C  CALCULATION.
C
                 IF(JAVX.AND.LOWP)Q(I)=Q(I)/DSQRT(1.-CRF2(I))
 290             CONTINUE
C
C  STRONG COLLISION VECTOR IS FIRST GUESS
C
                 AA=1.
                 TT=1.
                 IF(I.LE.NREACT) GO TO 310
                 AK=AVRATE(I-NREACT)
                 IF(JAVX.AND.LOWP)AK=AK+EXL2(I)
                 AA=OMEGA/(OMEGA+AK)
                 IF(.NOT.LOWP) TT=AA
C
C  IN THE FOLLOWING, THE FIRST GUESS AT THE EIGENVECTOR, Q(I), FOR ENERGIES
C  ABOVE NREACT IS THE STRONG COLLISION FORM.
C
C                    BETA*OMEGA
C     QSC(E)=F(E)-------------------
C                 BETA*OMEGA + K(E)
C
C
                 IF(ID.NE.1)GO TO 300
                 IF(JAVX.AND.LOWP)Q(I)=Q(I)*(1.-CRF2(I))
                 Q(I)=Q(I)*OMEGA*BETA/(OMEGA*BETA+AK)
 300             CONTINUE

C  THE STRONG COLLISION LOW PRESSURE RATE IS COMPUTED IN SUBROUTINE STRONG.
C  N.B., STRONG ALWAYS USES A GRAIN SIZE OF 100 CM-1.
C
C           __           OMEGA
C           \  F(E) --------------
C           /_       OMEGA + K(E)
C    KSC = --------------------------
C                __
C                \  F(E)
C                /_
C
C
                 STERM=CI2*AK*AA
                 SUMSCT=SUMSCT+STERM
                 SUMCH=0.
                 DO ISUM=1,NCH2
                   SUMCH=SUMCH+AVR2(I-NREACT,ISUM)
                 ENDDO
                 SUMSC2=SUMSC2+CI2*SUMCH*AA
 310             CONTINUE
C
C  ALMT IS USED AS A CORRETED NUMERATOR IN SUBROUTINE NESBET
C  WHEN KUNI IS CALCULATED, IE:
C
C
C         NREACT     NI
C          __        __          OMEGA  
C   ALMT = \ F(E) +  \ F(E) --------------
C          /_        /_      OMEGA + K(E)
C         E=0       NREACT
C
C
                 ALMT=ALMT+TT*CI2
               ENDDO 
               SUMSCT=SUMSCT/DENST
               SUMSC2=SUMSC2/DENST
               SUMSC1=SUMSCT-SUMSC2
               IF(ID.EQ.1) GO TO 320
C
C  FOR SUBSEQUENT GUESS TO Q(I), USE LINEARLY-INTERPOLATED RESULT FROM
C  PREVIOUS GRAIN SIZE.
C
               AROW(1)=0.5*Q(1)
               DO I=2,NILST*2,2
                 ION2=I/2
                 AROW(I)=Q(ION2)
                 AROW(I+1)=0.5*(Q(ION2)+Q(ION2+1))
               ENDDO
               DO I=NILST*2+1,NI
                 AROW(I)=AROW(I-1)*0.5
               ENDDO
               DO I=1,NI
                 Q(I)=AROW(I)
                 IF(D(I).LE.0.) Q(I)=0.
                 IF(I.LT.INT(ALMT1)) Q(I)=RTRHB(I)
               ENDDO
C
C  INTERPOLATED EIGENVECTOR GUESS HAS NOW BEEN FOUND.
C
 320           CONTINUE
C
C  STORE D(I)*B(I) IN AROW.
C
               NILOC=NI
               IF(LOWP) NILOC=INT(DFLT(NDEGS)*100./DELT)
               if(lowp) niloc=niloc/testinc
               DO 330 J=1,NILOC
                 IF(N.NE.1) AROW(J)=RTRHB(J)**2
                 IF(.NOT.LOWP) GO TO 330
                 AROW(J)=0.
C
C  AROW IS NOW IN (100/TESTINC) CM-1 INTERVALS
C
                 IF(D(J).GT.0.) AROW(J)=
     *           DEXP(DLOG(D(J))-DFLT(J)*EXPON)
 330           CONTINUE
               NRPL1=NREACT+1
               IPRT=IPRT+1
               CALL ENEXPT(DELT,ALPHA,IPRT)
C
C  THIS SUBROUTINE GENERATES PROBABILITY MATRIX
C  THE RESULTING MATRIX IS STORED IN ARRAY PROB.
C
               IF(NBAND.LE.0) GO TO 397
C
C  NESBET IS SUBROUTINE WHICH CARRIES OUT COMPUTATION OF EIGENVALUE.
C  SET STARTING VALUE AT LOWEST NON-ZERO POPULATION.
C
               DO I=1,(200*TESTINC)
                 IF(AROW(I).GT.0.) GO TO 340
               ENDDO
 340           IMIN=I
C
C  THIS LAST SET OF INSTRUCTIONS IS TO ALLOW FOR SOME RHO(E) ( D(I) ) BEING 0.
C
               CALL NESBET(B,ITMAX)
C
C  FIRST TIME THROUGH BOTH DELTAE HALVING LOOP AND PRESSURE LOOP PRINT OUT
C  LOW PRESSURE RESULTS:
C
               IF(ID.EQ.1)THEN
                 IF(LOWP)THEN
                    WRITE(6,29)
C
  29   FORMAT(//,10X,'THE FOLLOWING IS FOR THE LOW PRESSURE',
     *' LIMIT')
C
                   WRITE(6,30)
C
  30   FORMAT(/,' NO. NOT',T12,'DELTA',T20,'BAND-',T30,'MATRIX',
     * T44,'TOTAL',T54,'NO.OF',/, ' REACTING',T13,'E',
     * T20,'WIDTH',T31,'SIZE',T40,' RATE (K0)',T52,'ITERATIONS')
C
                 ELSE
                   PPAS=PN*133.3
                   WRITE(6,31) PN,PPAS,OMEGA
C
  31   FORMAT(1P,6X,'PRESSURE  =',T20,E12.2,' TORR (=',E10.2,' PASCAL)'
     * ,/,10X,'OMEGA  =',T20,E12.2,' SEC-1',///,
     * ' EXACT (EIGENVALUE) RATE   STRONG COLLISION ',
     * 3X,'TOTAL','   DELTA ',' BAND-',' NO.NOT',4X,'NUMBER',/,
     * 5X,'K1',4X,'K2,3',4X,'TOTAL',6X,'K1',4X,'K2,3',
     * 5X,'DIMENSION',3X,'E',3X,'WIDTH',' REACTING',1X,'ITERATIONS')
C
                 ENDIF
               ENDIF
               IF(JAVX)THEN
                 IF(LOWP)THEN
                   CORRF=CORPF
                 ELSE
                   CORRF=CORAV
                 ENDIF
               ELSE
                 CORRF=CORRAT(ITEMPR)
               ENDIF 
               RATE1=RATE1*CORRF
               RATE2=RATE2*CORRF
               RATTOT=RATTOT*CORRF
               SUMSC1=SUMSC1*CORRF
               SUMSC2=SUMSC2*CORRF
               SUMSCT=SUMSCT*CORRF
               IF(.NOT.LOWP) WRITE(6,32) RATE1,RATE2,
     *                 RATTOT,SUMSC1,SUMSC2,NI,DELT,NBAND,
     *                 NREACT,ITR
C
  32   FORMAT(1P,3E9.2,1X,E7.1,1X,E7.1,0P,I7,F10.0,I5,
     * I7,3X,I6)
C
C  FOR THE LOW PRESSURE CASE NESBET CALCULATES [M]*K0 IN S-1
C  NEED TO CONVERT TO K0 WITH UNITS CM3 MOL S-1
C
C  USE IDEAL GAS RELATION PV=NRT TO CALCULATE V/N:
C
               IF(LOWP) RSEC=RATTOT*T/(1.603E-5*PN)
               IF(LOWP) WRITE(6,33) NREACT,DELT,NBAND,NI,RSEC,ITR
C
  33   FORMAT(1X,I4,5X,F5.0,2X,I5,5X,I5,6X,1P,E10.2,0P,I9)
C
               NILST=NI
               NI=NI*2
 397         CONTINUE
             IF(N.NE.1.OR.IOPTHT.NE.0) WRITE(6,34) PN,OMEGA,RATE1
C
  34   FORMAT(1P,//,10X,25(' *'),/,
     * 10X,'FINAL RESULT: AT PRESSURE =',E10.2,
     * ' TORR, OMEGA =',E10.2,' S-1',/,
     * 20X,' Kuni (S-1) =',T40,E10.2,' (CHANNEL 1)')
C
             IF(.NOT.LOWP.AND.NCHAN.GE.2)
     *       WRITE(6,35)(RATEA(I23)*CORRF,I23=1,NCH2)
C
  35   FORMAT(1P,T40,E10.2, ' (CHANNEL 2)',/,T40,E10.2,' (CHANNEL 3)')
C
             IF((.NOT.LOWP).AND.RATTOT.GT.1.05*SUMSCT) WRITE(6,36)
C
  36   FORMAT(5X,'WARNING,WEAK COLL. RESULT GREATER THAN STRONG'
     * ,', INDICATES NUMERICAL INSTABILITY,'
     * /,5X, ' REDUCE GRAIN SIZE AND/OR TOLERANCES')
C
C  CARRY OUT HIGH PRESSURE AND STRONG COLLISION LOW PRESSURE CALCULATIONS.
C
             NRCTH=INT(E0*WKA*TESTINC/100.)
             IF(ICHP.EQ.0) CALL STRONG(HP,ALPT,NDEGS,T,R1,R2)
             IF(.NOT.JAVX)THEN
               IF(ICHP.EQ.0)THEN
                 HP=HP*CORRF
                 ALPT=ALPT*CORRF
               ENDIF
             ELSE
               HP=RATHP(1)+RATHP(2)+RATHP(3)
               ALPT=STLPJ*PN/(T*6.237E+4)
               ALPT=ALPT/OMEGA
             ENDIF
             ICHP=ICHP+1
C
C  CALL SUBROUTINE TO CALCULATE TROE QUANTITIES FLH, FSC, FSC.
C
             IF(.NOT.LOWP) CALL TROEF(HP,ALPT,SUMSCT,IOPTHT,
     *                              RATOTP,T,BETA,NDEGS,R1,R2)
             WRITE(6,37)
C
  37   FORMAT(/,25X,'*******************************',///)
C
             NI=NI/2
             IF(.NOT.LOWP) GO TO 398
C
C  CALCULATE AND PRINT OUT LOW PRESSURE RATE COEFFICIENT, RSEC
C  AND COLLISION EFFICIENCY, BETA
C
             RATOTP=RATTOT/OMEGA
             CSET=HP/RATTOT
             NRATES=NRTSTM
             WRITE(6,38) RSEC
C
  38   FORMAT(1P,10X,'LOW PRESSURE K0 =',E10.2,' CM3 MOL-1 S-1')
C
             BETA=RATOTP/ALPT
             WRITE(6,39) BETA
C
  39   FORMAT(1P,/,20X,'BETA =',T56,E10.2)
C
             NEFF=INT(E0EFF*WKA/DELT)
             WRITE(6,*)
             IF(NCHAN.GT.1) CALL MULTCH(R1,R2,E0,NDEGS,
     *                                           RSEC,T,PN)
             WRITE(6,40)
C
  40   FORMAT(////,19X,'-------FALL-OFF CALCULATION-------',///)
C
C  LOW PRESSURE CALC. NOW DONE, SO START FALL-OFF CALCS.
C
             IOPTPR=IOPTEM
 398       CONTINUE
 399     CONTINUE
 400   CONTINUE
       STOP
       END
C
C****************************************************************************
C
      SUBROUTINE SETUPL(E0,T,R1S,DELTS,NDEGS,ALPHA,IPRT,IOPTEM)
C
C  PREPARATION FOR LOW-PRESSURE CALCULATION
C
       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
       PARAMETER (NMAX100=10000, NMAX50=5000, NMAX25=2500)
       COMMON /AMIN/ ID,NCUT,NILOC
       COMMON /BLOCK1/ NRATES,NI,NE0,NEFF
       COMMON /BLOCK2/ DELT,CORRF
       COMMON /BLOCK3/ NREACT,NBAND,NRPL1
       COMMON /BLK10/ AVRATE(NMAX50)
       COMMON /BLK13/ AROW(NMAX100),OMEGA,ANORM(NMAX100)
       COMMON /B2/ PROB(NMAX25)
       COMMON /DENS/ RHO(NMAX100)
       COMMON /EXP/ ERR1,ERR2,ERR3
       COMMON /F1/ IFIRST,NSPEC,IMIN
       COMMON /F2/ ALMT,ALMT1,DET
       COMMON /JLP1/ JAVX
       COMMON /OPT/ IOPTPR
       COMMON /TRANS1/ RATTOT,RATE1,RATE2,RATEA(2)
       COMMON /VERSCH/ IXV,JXV
       COMMON /NEWKIDS/ TESTINC
C
       DIMENSION R1S(NMAX100)
C
       LOGICAL JAVX
C
C  SET UP RATE VECTOR AS COLLISION RATES
C
       DELT=(100./TESTINC)
       NI=NE0
       EXPON=1.43879*DELT/T
       ITEST=TESTINC*1000
       NI=MIN0(ITEST,NI)
       NREACT=NI/2
       NRPL1=NREACT+1
       ID=5
       IOPTEM=IOPTPR
       IOPTPR=-1
C
C  AROW IS THE BOLTZMANN POPULATION VECTOR AT ENERGY I*DELT
C
C  IE  AROW(E) = F(E) = RHO(E).EXP(-E/KT)
C
C  NB: THE NORMALIZING FACTOR IS NOT INCLUDED IN AROW
C
       DO I=1,NDEGS
         AROW(I)=0.
         IF(RHO(I).GT.0.) AROW(I)=DEXP(DLOG(RHO(I))-DFLT(I)*EXPON)
       ENDDO
       DO I=1,(20*TESTINC)
         IF(RHO(I).GT.0.) GO TO 1000
       ENDDO
 1000  IMIN=I
       NCUT=1
       NILOC=NDEGS
C
C  THE MAXIMUM NUMBER OF LEVELS, FOR THE PURPOSE OF SETTING UP
C  THE DIAGONAL COLLISION TERMS IN THE LOW PRESSURE CALCULATION,
C  MUST BE THE MAXIMUM NUMBER CONSIDERED (NDEGS) FOR CALCULATING
C  THE NORMALIZING FACTOR ANORM.
C
C  ENEXPT CALCULATES THE PROBABILITY MATRIX FOR COLLISIONAL ENERGY 
C  TRANSFER, IE MATRIX J0, WHERE ALL MICROSCOPIC REACTION RATES ARE
C  ZERO
C
C  IN A LOW PRESSURE CALCULATION WE CONSIDER THE MATRIX K OF RATE
C  COEFFICIENTS AS BEING PERTURBED BY THE MATRIX J0 - BUT ALL K'S
C  ARE ZERO BELOW THE CRITICAL ENERGY E0 SO CONSIDERING THE LARGEST
C  EIGENVALUE ONLY WE CAN TRUNCATE THE MATRIX AT E0
C
       CALL ENEXPT(DELT,ALPHA,IPRT)
C
       NIP1=NI+1
       INDI=0
       DO 1020 I=NRPL1,NI
         INDI=INDI+1
         R1S(INDI)=0.
         IF(RHO(I).LE.0.) GO TO 1020
         SUM=0.
C
C  ONLY COUNT TERMS WITHIN THE BANDWIDTH OF PROB MATRIX
C
         IF(NIP1-I.GT.NBAND) GO TO 1020
         JMAX=MIN0(NDEGS,NIP1+NBAND)
         IF(NIP1.GE.JMAX) GO TO 1020
C
C  IN LOW PRESSURE LIMIT THE EIGENVALUE RELATIONSHIP IS:
C
C            EO                               EMAX
C             __                              __
C  -K0 G(I) = \ PROB(J-I)G(J)-PROB(I-J)G(I) - \ PROB(J-I)G(I)
C             /_                              /_
C            J=0                              J=E0
C
C  WHERE G(I) IS THE EIGENVECTOR OF THE MATRIX J0, = F(E)
C  AND THE NORMALISATION IS ASSUMED
C
         DO 1010 J=NIP1,JMAX
           II=J-I+1
           IF(II.LE.0.OR.II.GT.NBAND) GO TO 1010
C
C  TERM IS THE SECOND TERM IN THE ABOVE SUM DIVIDED THROUGH BY G(I)
C  THIS TERM REPRESENTS THE 'LEAKAGE' TO ENERGY LEVELS ABOVE THE
C  CRITICAL ENERGY AND IS EQUIVALENT TO THE RATE TERM IN THE NORMAL
C  FALLOFF PROCEDURE
C
           TERM=(RHO(J)/RHO(I))*DEXP(DFLT(I-J)*EXPON)*
     *                                    PROB(II)*ANORM(J)
C
C  FOR THE LOW PRESSURE LIMIT, ONE DOES NOT USE THE CORRECT
C  TRAPEZOIDAL RULE (0.5*FIRST AND LAST TERMS) AS PER THE
C  FOLLOWING COMMENT CARD, SINCE IT CAN BE SHOWN BY APPROPRIATE
C  MATRIX PERTURBATION THEORY THAT, FOR THE LOW PRESSURE LIMIT
C  ONLY, CONVERGENCE IN GRAIN SIZE IS ACCELERATED BY OMITTING IT.
C
C     IF(J.EQ.NIP1.OR.J.EQ.JMAX) TERM=TERM*0.5
C
           SUM=SUM+TERM
 1010    CONTINUE
C
C  IE:       E0 EMAX
C             __ __  
C             \  \  DELTAE.PROB(J-I).ANORM(J).G(I)
C             /_ /_
C            I=0 J=E0
C  [M].K0 = ----------------------------------------
C                        E0
C                         __
C                         \  G(I)
C                         /_
C                        I=0
C
         R1S(INDI)=SUM*OMEGA*DELT
 1020  CONTINUE
C
C  COLLF GENERATES EXPRESSION FOR EXL:
C
       IF(JAVX)CALL COLLF(DELT)
       ITEST=TESTINC*1000
       NRATES=MIN0(ITEST,INDI)
       NI=INT(E0/(DELTS*2.859E-3))
       RETURN
       END
C
C****************************************************************************
C
       SUBROUTINE RATES(R1,R2)
C
       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
       PARAMETER (NMAX100=10000, NMAX50=5000, NMAX25=2500)
       COMMON /LOW/ LOWP
       COMMON /BLOCK1/ NRATES,NI,NE0,NEFF
       COMMON /BLOCK2/ DELT,CORRF
       COMMON /BLK10/ AVRATE(NMAX50)
       COMMON /BLK12/ AVR2(NMAX50,2)
       COMMON /BLK15/ NCH2
       COMMON /DEOC/ EI,WE0
       COMMON /NEWKIDS/ TESTINC
C
       DIMENSION R1(1),R2(NMAX100,2),ALAST2(2),ANEXT2(2),SP2(2)
C
       LOGICAL LOWP
C
C  FORM VECTORS OF MICROSCOPIC RATES REQUIRED FOR GIVEN ENERGY SPACING
C  FROM THE INPUT MICRSCOPIC RATES FOR 100 CM-1 SPACING (R1,R2).IF
C  SPACING GE 100 CM-1, DO BY TABLE SELECTION, OTHERWISE BY LINEAR
C  INTERPOLATION.
C  BRANCH DEPENDING ON WHETHER DELT LT OR GE (100. OR 10.) 
C
       IF(DELT.LT.(99.99/TESTINC)) GO TO 2030
C
C  FIND INTEGER FOR TABLE SELECTION, EXTRA 1. IS TO ACCOUNT FOR ANY ROUNDOFF
C
       N=INT((DELT+1.)*(0.01*TESTINC))
       J=0
       IF(DELT.GT.(101./TESTINC))GO TO 2000
       N3=1
       GO TO 2010
C
C  DELT GT (100/TESTINC) CM-1:
C
 2000  CONTINUE
       INTE0=INT(WE0*TESTINC/100.)
       IF(LOWP)INTE0=INT(WE0*TESTINC/200.)
       IF(MOD(INTE0,2).EQ.1)THEN
         N3=1
       ELSE
         N3=2
       ENDIF
C
C  DELT GE 100/TESTINC CM-1
C  WANT TO END UP WITH 100/TESTINC CM-1 SPACING REGARDLESS
C
 2010  CONTINUE
       DO 2020 II=N3,NRATES,N
         I=MIN0(II,NRATES)
         J=J+1
         DO ISUM=1,NCH2
           AVR2(J,ISUM)=0.
         ENDDO
         AVRATE(J)=R1(I)
         IF(LOWP) GO TO 2020
         SUMCH=0.
         DO ISUM=1,NCH2
           AVR2(J,ISUM)=R2(I,ISUM)
           SUMCH=SUMCH+AVR2(J,ISUM)
         ENDDO
         AVRATE(J)=R1(I)+SUMCH
 2020  CONTINUE
       GO TO 2055
C
C  DELT LT (100/TESTINC) CM-1
C  FORM RATE VECTOR USING LINEAR INTERPOLATION
C
 2030  N=INT(101./DELT)
       JJ=0
       ALAST1=(R1(1))
       DO ISUM=1,NCH2
         ALAST2(ISUM)=0.
       ENDDO
       DO 2050 I=1,NRATES
         IF(I.LT.NRATES) ANEXT1=R1(I+1)
         DO ISUM=1,NCH2
           ANEXT2(ISUM)=0.
         ENDDO
         IF(I.EQ.NRATES) GO TO 2040
         IF(LOWP.OR.R2(I+1,1).LE.0.) GO TO 2040
         DO ISUM=1,NCH2
           ANEXT2(ISUM)=R2(I+1,ISUM)
         ENDDO
 2040    CONTINUE
         SP1=(ANEXT1-ALAST1)/DFLT(N)
         DO ISUM=1,NCH2
           SP2(ISUM)=(ANEXT2(ISUM)-ALAST2(ISUM))/DFLT(N)
         ENDDO
         DO J=1,N
           JJ=JJ+1
           AVRATE(JJ)=(ALAST1+SP1*DFLT(J-1))
           SUMCH=0.
           DO ISUM=1,NCH2
             AVR2(JJ,ISUM)=ALAST2(ISUM)+SP2(ISUM)*DFLT(J-1)
             SUMCH=SUMCH+AVR2(JJ,ISUM)
           ENDDO
           AVRATE(JJ)=AVRATE(JJ)+SUMCH
         ENDDO
         ALAST1=ANEXT1
         DO ISUM=1,NCH2
           ALAST2(ISUM)=ANEXT2(ISUM)
         ENDDO
 2050  CONTINUE
 2055  RETURN
       END
C
C***********************************************************************
C
       SUBROUTINE GENCF
C
C  SELECTS DEXP(-R0/KT) BY GRAIN SIZE, ONLY UTILISED IF JAVX.
C
       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
       PARAMETER (NMAX100=10000, NMAX50=5000, NMAX25=2500)
       COMMON /BLOCK1/ NRATES,NI,NE0,NEFF
       COMMON /BLOCK2/ DELT,CORRF
       COMMON /JLP2/ CRF2(NMAX100),EXL2(NMAX100)
       COMMON /JLP3/ CRF(NMAX100),EXL(NMAX100)
       COMMON /NEWKIDS/ TESTINC
C
       IF(DELT.LT.(99.99/TESTINC))GO TO 3000
       N=INT((DELT+0.01)*(0.01*TESTINC))
       J=0
C
C  NT IS THE MAXIMUM ENERGY FOR WHICH RATES WERE ENTERED
C  CRF IS EXPRESSION FOR EXP(-ROTHR/KT)
C  EXL IS CALCULATED IN SUBROUTINE COLLF CALLED FROM SETUPL
C
       NT=NE0+NRATES
       DO II=N,NT,N
         I=II
         IF(I.GT.NT)I=NT
         J=J+1
         CRF2(J)=CRF(I)
         EXL2(J)=EXL(I)
       ENDDO
       GO TO 3001
C
C  SPLIT FOR DELTA LT (100/TESTINC) CM-1, 
C  FORM VECTORS CRF2 AND EXL2 BY LINEAR INTERPOLATION.
C
 3000  CONTINUE
       N=INT(101./(DELT*TESTINC))
       JJ=0
       ALAST2=CRF(1)
       ALAST3=EXL(1)
       DO I=1,NE0
         IF(I.LT.NE0)ANEXT2=CRF(I+1)
         IF(I.LT.NE0)ANEXT3=EXL(I+1)
         SP2=(ANEXT2-ALAST2)/DFLT(N)
         SP3=(ANEXT3-ALAST3)/DFLT(N)
         DO J=1,N
           JJ=JJ+1
           CRF2(JJ)=ALAST2+SP2*DFLT(J-1)
           EXL2(JJ)=ALAST3+SP3*DFLT(J-1)
         ENDDO
         ALAST2=ANEXT2
         ALAST3=ANEXT3
       ENDDO
 3001  RETURN
       END
C
C**************************************************************************
C
      SUBROUTINE GENDEG(NDEGS)
C
C  FINDS VECTOR OF DEGENERACIES FOR ANY GRAIN SIZE BY INTERPOLATION
C  IN THIS CASE VECTOR IS THE DENSITY OF STATES, RHO
C
       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
       PARAMETER (NMAX100=10000, NMAX50=5000, NMAX25=2500)
       COMMON /BLOCK1/ NRATES,NI,NE0,NEFF
       COMMON /BLOCK2/ DELT,CORRF
       COMMON /BLOCK7/ D(NMAX100)
       COMMON /DENS/ RHO(NMAX100)
       COMMON /NEWKIDS/ TESTINC
C
C  BRANCH DEPENDING ON WHETHER DELT IS GT OR LT 100.
C
       IF(DELT.LT.(99.99/TESTINC)) GO TO 4000
C
C  FIND INTEGER FOR TABLE SELECTION
C
       N=INT((DELT+0.01)*(0.01*TESTINC))
       J=0
       DO II=N,NDEGS,N
         I=II
         IF(I.GT.NDEGS) I=NDEGS
         J=J+1
         D(J)=RHO(I)
       ENDDO
       GO TO 4001
 4000  N=INT(101./(DELT*TESTINC))
       JJ=0
       ALAST=DLOG(RHO(1))
       DO I=1,NDEGS
         IF(I.LT.NDEGS) ANEXT=DLOG(RHO(I+1))
         SP=(ANEXT-ALAST)/DFLT(N)
         DO J=1,N
           JJ=JJ+1
           D(JJ)=DEXP(ALAST+SP*DFLT(J-1))
         ENDDO
         ALAST=ANEXT
       ENDDO
 4001  RETURN
       END
C
C***********************************************************************
C
      SUBROUTINE ENEXPT(DELT,ALPHA,IPRT)
C
C  COMPUTES COLLISION PROBABILITY ARRAY PROB.
C  TO CHANGE PROBABILITY FUNCTION,  CHANGE FUNCTIONAL FORM
C
       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
       PARAMETER (NMAX100=10000, NMAX50=5000, NMAX25=2500)
       COMMON /BLOCK1/ NRATES,NI,NE0,NEFF
       COMMON /VERSCH/ IXV,JXV
       COMMON /LOW/ LOWP
       COMMON /BLOCK3/ NREACT,NBAND,NRPL1
       COMMON /AMIN/ ID,NCUT,NILOC
       COMMON /BLK13/ AROW(NMAX100),OMEGA,ANORM(NMAX100)
       COMMON /B2/ PROB(NMAX25)
       COMMON /EXP/ ERR1,ERR2,ERR3
       COMMON /NEWKIDS/ TESTINC
C
       LOGICAL LOWP,BRW
C
       NBAND=0
       BRW=IXV.LE.0.AND.JXV.LE.0
       NE00=NREACT
       IF(LOWP) NE00=NI
       Z=DLOG(AROW(NE00+1)/AROW(NE00))*ALPHA*ALPHA/DELT
       FS2=4.*ALPHA*ALPHA
       NTEST=INT(249.9*TESTINC)
       NTOPX=MIN0(NTEST,NILOC)
       DO I=1,NTOPX
         PROB(I)=0.
       ENDDO
C
C  GENERATE UN-NORMALIZED PROBABILITY MATRIX, J0, WITH ELEMENTS:
C
C                J0(I,J)=R(I,J)
C                          __
C                J0(I,I)=- \  R(J,I)
C                          /_
C                        I.NE.J
C
C  WHERE R(I,J), THE RATE OF ENERGY TRANSFER FROM LEVEL J TO I,
C  DEPENDS ONLY ON THE ENERGY DIFFERENCE J-I AND IS GIVEN BY
C  PROB(J-I) - WHICH NEEDS TO BE NORMALISED BY ANORM(J)
C
       BOT=0.
       DEDOWN=0.
       DO I=1,NTOPX
         DE=DFLT(I-1)*DELT
         X=DE/ALPHA
         IF(BRW) GO TO 5000
C
C  USE EXPONENTIAL DOWN MODEL FOR PROB(J-I)
C
         XX=0.
         IF(X.GT.0.) XX=X**JXV
         PWR=1.
         IF((X.GT.0.).AND.IXV.NE.0) PWR=X**IXV
         IF((X.LE.0.).AND.IXV.NE.0) PWR=0.
         IF(XX.LT.20.) PROB(I)=PWR*DEXP(-XX)
         GO TO 5010
C
C  USE BIASED RANDOM WALK MODEL FOR PROB(J-I)
C
 5000    TOP=-(Z+DE)
         PROB(I)=DEXP(-TOP*TOP/FS2)
C
C  FOLLOWING ADDED BY RGG JUN 29 1989
C
         IF(PROB(I).LT.ERR3) GO TO 5020
 5010    CONTINUE
         IF(PROB(I).LE.ERR3.AND.X.GT.0.5) GO TO 5020
C
C  INTEGRATE USING THE TRAPEZOIDAL RULE, OMITTING THE FACTOR OF 0.5
C  AT THE HIGH ENERGY END
C
         IF(I.GT.1)THEN
           DEDOWN=DEDOWN+DE*PROB(I)
           BOT=BOT+PROB(I)
         ELSE
           DEDOWN=0.5*DE*PROB(1)
           BOT=0.5*PROB(1)
         ENDIF
C
C  NBAND IS THE BANDWIDTH OF THE MATRIX PROB, IE HOW FAR YOU MUST GO 
C  FROM THE DIAGONAL BEFORE THE OFF-DIAGONAL ELEMENTS ARE LESS THAN
C  ERR3
C
         NBAND=MAX0(NBAND,I)
       ENDDO
C
C                                          E'
C                                          __
C                                          \  (E'-E)PROB(E'-E)
C                                          /_
C                                          E=0
C  SET PROB(250*TESTINC) =  = ------------------------
C                                          E'
C                                          __
C                                          \  PROB(E'-E)
C                                          /_
C                                          E=0
C
 5020  PROB(I)=0.
       PROB(TESTINC*250)=DEDOWN/BOT
       NBAND=NBAND-1
       IF(NBAND.EQ.1) NBAND=0
       IF(NBAND.EQ.0) GO TO 5125
       IF(.NOT.LOWP) NILOC=NI
C
C  COMMENCE FINITE DIFFERENCE SOLUTION OF INTEGRAL EQUATION FOR C.
C
       JJ=NILOC+1
C
C  GENERATE ANORM, THE VECTOR OF NORMALIZERS OF THE ENERGY TRANSFER
C  RATE, R(I,J), INITIALLY AS C, ITS RECIPROCAL.
C
       DO 5060 J=1,NILOC
         JJ=JJ-1
         ANORM(JJ)=0.
         IF(AROW(JJ).LE.0.) GO TO 5060
         JP=JJ+1
         JM=JJ-1
         IF(JJ.GE.(NREACT/NCUT)) GO TO 5030
C
C  TO AVOID OCCASIONAL PROBLEMS WITH NORMALIZATION ALGORITHM FOR STATES
C  FOR STATES OF LOW ENERGY (BELOW NREACT/NCUT), ALL ARE GIVEN THE
C  SAME NORMALIZATION. THIS HAS NO PHYSICAL EFFECT, SINCE THESE LEVELS
C  ALWAYS HAVE THEIR EQUILIBIRIUM POPULATIONS.
C
         ANORM(JJ)=ANORM(JJ+1)
         GO TO 5060
 5030    CONTINUE
         DE=PROB(1)
         IMIN=MAX0(1,JJ-NBAND)
         IF(JM.LT.IMIN) GO TO 5040
C
C  I.E. FOR THE FIRST LOOP THROUGH FIND NORMALISATION AT MAX ENERGY:
C
C                                NMAX
C                                 __
C          ANORM(NMAX) = DELTAE * \  PROB(NMAX-I)
C                                 /_
C                                I=1
C
         DO I=IMIN,JM
           DE=DE+PROB(JP-I)
         ENDDO
 5040    ANORM(JJ)=DE*DELT
         IF(JJ.EQ.NILOC) GO TO 5060
         IUP=MIN0(NILOC,JJ+NBAND)
         DE=0.
         IF(IUP.LT.JP) GO TO 5060
C
C  SUBSEQUENTLY AT ENERGY JJ:
C
C
C                                JJ
C                                __
C                        DELTAE  \ PROB(JJ-I)
C                                /_
C                               I=1
C  ANORM(JJ) = -----------------------------------------------
C                            NMAX
C                   DELTAE   __
C              1 - --------  \  AROW(I)*PROB(I-JJ)/ANORM(I)
C                  AROW(JJ)  /_
C                           I=JP
C
C
         DO 5050 I=JP,IUP
           IF(ANORM(I).LE.0.) GO TO 5050
           DE=DE+(PROB(I-JM)*AROW(I)/ANORM(I))
 5050    CONTINUE
         ANORM(JJ)=ANORM(JJ)/(1.-(DELT*DE/AROW(JJ)))
 5060  CONTINUE
C
C  INVERT ANORM TO OBTAIN NORMALISATION VECTOR
C
       JS=NILOC+1
       DO J=1,NILOC
         JS=JS-1
         IF(ANORM(JS).GT.0.) ANORM(JS)=1./ANORM(JS)
         IF(ANORM(JS).LE.0..AND.JS.LT.NILOC) ANORM(JS)=ANORM(JS+1)
         IF(AROW(J).LE.0.) ANORM(J)=0.
       ENDDO
       IF(LOWP.AND.ID.EQ.5) GO TO 5070
       IF(IPRT.NE.1) GO TO 5125
       IF(BRW) WRITE(6,5001) DABS(ALPHA)
C
 5001  FORMAT(//,10X,'BIASED RANDOM WALK MODEL FOR P(E,E',1H',')',/,
     * 20X,'S =',F10.1,' CM-1')
C
       IF (TESTINC .EQ. 1)WRITE(6,5002) (PROB(I),I=1,NBAND)
C
C   Print out only every 10 R(X) values in inc=10
C
       IF (TESTINC .EQ. 10)WRITE(6,5005) (PROB(I),I=1,NBAND,10)
C
 5002  FORMAT(/,'  VALUES OF R(X) (FOR LARGEST GRAIN SIZE):',
     * /,1P,10(3X,6E12.2,/))
 5005  FORMAT(/,'  VALUES OF R(X) (EVERY TEN )
     *    (FOR LARGEST GRAIN SIZE):',
     * /,1P,10(3X,6E12.2,/))
C
       GO TO 5125
C
C  COMPUTE AVERAGE ENERGY TRANSFER QUANTITIES
 5070  WRITE(6,5003)
C
 5003  FORMAT(//,'  ENERGY (KJ/MOL) DELTA E TOTAL (CM-1)',
     * ' SQRT(DE**2)',T56,'EQUILIBRIUM POPULATION')
C
       NZ=NI/2
       SUMP=0.
       DO J=1,NILOC
         SUMP=SUMP+AROW(J)
       ENDDO                           
       ATEST=1.0
       DO 5120 II=NZ,NILOC,10
         E=DFLT(II)*DELT*1.195E-2
         ANM=ANORM(II)
         IF(ANM.LE.0.) GO TO 5120
         JP=II+1
         JM=II-1
         DE=0.
         DE2=0.
         IF(JM.LT.1) GO TO 5090
         DO 5080 I=1,JM
           IF(II-I.GE.NBAND) GO TO 5080
           FAC=ANM*PROB(JP-I)
           DIF=DFLT(I-II)
           DE=DE+DIF*FAC
           DE2=DE2+DIF*DIF*FAC
 5080    CONTINUE
 5090    IF(JP.GT.NILOC.OR.AROW(II).LE.0.) GO TO 5110
         RBJ=1./AROW(II)
         DO 5100 J=JP,NILOC
           IF((J-II.GE.NBAND).OR.(ANORM(J).LE.0.)) GO TO 5100
           FAC=PROB(J-JM)*ANORM(J)*AROW(J)*RBJ
           DIF=DFLT(J-II)
           DE=DE+DIF*FAC
           DE2=DE2+DIF*DIF*FAC
 5100    CONTINUE
 5110    DE=DE*DELT*DELT
         DE2=DSQRT(DE2*DELT)*DELT
         AROWN=AROW(II)/(DELT*SUMP)
C
C  If INC = 10cm-1 then print out only very ten sets of calculations.
C  i.e. if testinc=10
C
         IF (TESTINC .EQ. 1) WRITE(6,5004) E,DE,DE2,AROWN 
         IF (TESTINC .EQ. 10 .AND. (ATEST/10) .EQ. 
     *    INT(ATEST/10)) WRITE(6,5004) E,DE,DE2,AROWN
C
 5004  FORMAT(1X,F9.1,10X,F10.1,6X,F10.1,T62,1P,E10.2)
C
       ATEST=ATEST+1
 5120  CONTINUE  
 5125  RETURN
       END
C
C************************************************************************
C
      SUBROUTINE NESBET(B,ITMAX)
C
C  SUBROUTINE TO  SOLVE EIGENVALUE PROBLEM.
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NMAX100=10000, NMAX50=5000, NMAX25=2500)
      COMMON /AMIN/ ID,NCUT,NILOC
      COMMON /BLOCK1/ NRATES,NI,NE0,NEFF
      COMMON /BLOCK2/ DELT,CORRF
      COMMON /BLOCK3/ NREACT,NBAND,NRPL1
      COMMON /BLOCK6/ Q(NMAX100),H(NMAX100),RTRHB(NMAX100)
      COMMON /BLK10/ AVRATE(NMAX50)
      COMMON /BLK12/ AVR2(NMAX50,2)
      COMMON /BLK13/ AROW(NMAX100),OMEGA,ANORM(NMAX100)
      COMMON /BLK15/ NCH2
      COMMON /BLOK17/ AVEC(NMAX100)
      COMMON /B2/ PROB(NMAX25)
      COMMON /EXP/ ERR1,ERR2,ERR3
      COMMON /F1/ IFIRST,NSPEC,IMIN
      COMMON /F2/ ALMT,ALMT1,DET
      COMMON /JLP1/ JAVX
      COMMON /JLP2/ CRF2(NMAX100),EXL2(NMAX100)
      COMMON /LOW/ LOWP
      COMMON /OPT/ IOPTPR
      COMMON /TRANS1/ RATTOT,RATE1,RATE2,RATEA(2)
      COMMON /TRANS2/ ITR
      COMMON /NEWKIDS/ TESTINC
C
      DIMENSION RATK(NMAX50,2)
C
      LOGICAL DONE,LOWP,JAVX
C
C  INTEGRATION IS DONE USING COMPLETE TRAPEZOIDAL RULE, FOR INNER LOOPS
C  (J LE I OR J GE I). FOR OUTER LOOPS, THE FIRST AND LAST VALUES OF THE
C  INTEGRAND CONTRIBUTE NEGLIGIBLY TO THE INTEGRAL, SO THE 0.5*INTEGRAND
C  CORRECTION USED IN PROPER TRAPEZOIDAL INTEGRATION IS NOT INCLUDED.
C  AROW(I) IS RHO(I)*DEXP(-E/KT) = F(I)
C  GAMMA IS Q(I)/RTRHB(I) = Q(I)/DSQRT(F(I))
C  Q(I) IS THE EIGENVECTOR AND INITIALLY WE START WITH:
C
C  Q(I) = SQRT(F(I))   FOR I < ALMT1
C       = F(I)         FOR ALMT1 < I < NI
C
       RATEL=0.
       IFIRST=INT(ALMT1)
       IF(IFIRST.LT.IMIN) IFIRST=IMIN
       DO I=1,2
         DO II=1,(500*TESTINC)
           RATK(II,I)=0.
         ENDDO
       ENDDO
       RATEL2=0.
       DEOMGA=OMEGA*DELT
       DET=0.
       DET2=0.
       SS=0.
C
C  THE FOLLOWING LOOP ALSO GENERATES FIRST GUESS AT RATE COEFFICIENTS
C  LOOP FROM THE FIRST NON-ZERO ELEMENT OF AROW
C
       DO 6000 I=IMIN,NI
         H(I)=1.
         AVEC(I)=0.
         IF(AROW(I).LE.0.) GO TO 6000
         AVEC(I)=AROW(I)*ANORM(I)
         CI2=Q(I)*RTRHB(I)
         IF(JAVX.AND.LOWP)CI2=CI2*DSQRT(1.-CRF2(I))
         SS=SS+CI2
         H(I)=Q(I)/RTRHB(I)
         IF(I.LT.IFIRST) H(I)=1.
         IF(I.LE.NREACT) GO TO 6000
         EFR=AVRATE(I-NREACT)
         IF(JAVX.AND.LOWP)EFR=EFR+EXL2(I)
         DET=DET+EFR*CI2
         SUMCH=0.
         DO ISUM=1,NCH2
           SUMCH=SUMCH+AVR2(I-NREACT,ISUM)
         ENDDO
         DET2=DET2+SUMCH*CI2
 6000  CONTINUE
       IF(IOPTPR.LT.0) SS=ALMT
C
C  FIRST GUESS AT KUNI IS:
C
C           __
C           \  Q(I)*K(I)*SQRT(F(I))
C           /_
C  KUNI =  -------------------------
C             __
C             \  Q(I)*SQRT(F(I))
C             /_
C
C
C  WHERE Q(I) IS THE FIRST GUESS AT THE EIGENVECTOR
C
       RATE=DET/SS
       IF(ID.NE.1) RATE=RATTOT
       RATE2=DET2/SS
C
C  INITIALIZATION PROCESS NOW COMPLETE.
C
C  COMMENCE NESBET ITERATIVE PROCESS.
C
C  LOOP THROUGH THE NUMBER OF NESBET ITERATIONS:
C
       DO 6060 ITRX=1,ITMAX
         ITR=ITRX
C
C  THE FOLLOWING CONVERGENCE TEST GIVES EQUAL WEIGHTING TO THE ACCURACY
C  OF THE OVERALL RATE AND THE RATE FROM THE SECOND CHANNEL.
C  ALTERNATIVELY,THE WEIGHTING CAN BE FOR OVERALL RATE,DEPENDING ON INPUT
C  OPTION.
C
         DONE=ABS(1.-(RATEL/RATE)).LE.ERR2
         IF(IABS(IOPTPR).EQ.2) DONE=DONE.AND.ABS(1.-(RATEL2/RATE2)).LE.
     *       ERR2
         IF(DONE) GO TO 6070
         RATEL2=RATE2
         RATEL=RATE
C
C  LOOP THROUGH THE ENERGY:
C
         DO 6050 I=IFIRST,NI
           IF(AROW(I).LE.0.) GO TO 6050
           QI=Q(I)
           CRFI=(1.-CRF2(I))
           RTI=RTRHB(I)
           IP=I+1
           HI=H(I)
           IL1=I-1
           JMIN=MAX0(IMIN,I-NBAND+1)
           S=0.
           IF(I.EQ.NI) GO TO 6010
           JMAX=MIN0(NI,I+NBAND-1)
C
C  COMMENCE LOOP OVER J > I.
C
C                 NMAX
C                 __
C       1         \  F(J) * PROB(J-I) * ANORM(J) * (H(I)-H(J))  
C   ----------    /_
C   SQRT(F(I))    J=I
C
C
           IF(JMAX-IL1.LE.0) GO TO 6010
           HT=H(JMAX)
           IF(JAVX.AND.LOWP)HT=HT*DSQRT(CRFI*(1.-CRF2(JMAX)))
           S=0.5*(HI-HT)*AVEC(JMAX)*PROB(JMAX-IL1)
           JM=JMAX-1
           IF(JM.LT.IP) GO TO 6010
           DO J=IP,JM
             HT=H(J)
             IF(JAVX.AND.LOWP)HT=HT*DSQRT(CRFI*(1.-CRF2(J)))
             S=S+(HI-HT)*AVEC(J)*PROB(J-IL1)
           ENDDO
 6010      CONTINUE
           S=S/RTI
           S1=0.
           IF(I.LE.JMIN) GO TO 6020
C
C  COMMENCE LOOP OVER J < I.
C
C                             I
C                             __
C    SQRT(F(I)) * ANORM(I) *  \  PROB(I-J) * (H(I)-H(J))  
C                             /_
C                             J=1
C
C
           IF(IP-JMIN.LE.0) GO TO 6020
           HL=H(JMIN)
           IF(JAVX.AND.LOWP)HL=HL*DSQRT(CRFI*(1.-CRF2(JMIN)))
           S1=0.5*PROB(IP-JMIN)*(HI-HL)
           JM=JMIN+1
           IM=I-1
           IF(IM.LT.JM) GO TO 6020
           DO J=JM,IM
             HL=H(J)
             IF(JAVX.AND.LOWP)HL=HL*DSQRT(CRFI*(1.-CRF2(J)))
             S1=S1+PROB(IP-J)*(HI-HL)
           ENDDO
 6020      AA=-RATE
           S=S+S1*RTI*ANORM(I)
           IF(JAVX.AND.LOWP)S=S+PROB(1)*QI*ANORM(I)*CRF2(I)
           IF(I.LE.NREACT) GO TO 6030
           ILNR=I-NREACT
           AKE=AVRATE(ILNR)
           AA=AA+AKE
 6030      CONTINUE
C
C                                             (  __    __  )
C                                             (  \  +  \   )
C  S(I) = (K(I) - KUNI)*Q(I) + OMEGA*DELTAE * (  /_    /_  )
C                                             (  J>I   J (BY NUMERICAL INTEGRATION) =',F10.1,
     * ' CM-1')
C
       RATIOT=RATTOT/HP
       APOM=ALPT*OMEGA
       WRITE(6,8002) HP,RATIOT,APOM
C
 8002  FORMAT(1P,20X,'TOTAL K(INF) =',T56,E10.2,' S-1'
     * ,/,20X,'K(TOT)/K(INF) ='
     * ,T56,E10.2,/,20X,'STRONG COLLISION K0*[M] =',T56
     * ,E10.2,' S-1')
C
       IF(APOM.LT.SUMSCT) WRITE(6,8003)
C
 8003  FORMAT(' WARNING, STRONG COLLISION RATE GT STRONG COLLISION',
     * /,' VALUE FOR K0*[M]          USE SMALLER GRAIN SIZE')
C
       IF(IOPTHT.NE.0) RETURN
C
C  COMPUTE LINDEMANN-HINSHELWOOD AND STRONG COLLISION BROADENING FACTORS.
C
       WCK0=RATOTP*OMEGA
       IF(BETA.GT.1.) WRITE(6,8004)
C
 8004  FORMAT(' WARNING, BETA GREATER THAN 1, INDICATES'
     * ,' TOLERANCES SHOULD BE ALTERED')
C
       WCRAT=WCK0/HP
       WRITE(6,8005) WCK0,WCRAT
C
 8005  FORMAT(1P,20X,'WEAK COLLISION K0*[M] =',T56,E10.2,' S-1',
     * /,20X,'WEAK COLLISION K0*[M]/K(INF)=',T56,E10.2)
C
       WRITE(6,8006) BETA
C
 8006  FORMAT(1P,/,20X,'BETA =',T56,E10.2)
C
       FLH=WCRAT/(1.+WCRAT)
       OMSC=RATOTP*OMEGA/ALPT
       AKSC=0.
       EXPON=143.879/(T*TESTINC)
       BOT=0.
       ICSC=0.
       IF(JAVX)NN11=NEFF
       IF(.NOT.JAVX)NN11=NE0
       DO 8000 I=1,NDEGS
         TERM=0.
         IF(RHO(I).GT.0.) TERM=DEXP(DLOG(RHO(I))-DFLT(I)*EXPON)
         BOT=BOT+TERM
         IF(I.LE.(NN11)) GO TO 8000
         ICSC=ICSC+1
         AKE=R1(ICSC)+(R2(ICSC,1)+R2(ICSC,2))
         AKSC=AKSC+(TERM*OMSC*AKE/(AKE+OMSC))
 8000  CONTINUE
       AKSC=AKSC*CORRF/BOT
       FSC=AKSC/(FLH*HP)
       FWC=RATTOT/(HP*FLH*FSC)
       WRITE(6,8007) FLH,FSC,FWC
C
 8007  FORMAT(1P,//,20X,'FLH =',E10.2,',FSC =',E10.2,',FWC =',E10.2)
C
       IF(FWC.GT.1.5) WRITE(6,8008)
C
 8008  FORMAT(' WARNING,FWC SIGNIFICANTLY GREATER THAN 1 :'
     * ,/,' TOLERANCES SHOULD BE CHANGED',
     * ' OR THAT GRAIN SIZE SHOULD BE'
     * ,/,' REDUCED TO 100 CM-1 (OR 10 CM-1 AT LOW TEMPERATURES)')
C
       RETURN
       END
C
C****************************************************************************
C
      SUBROUTINE MULTCH(R1,R2,E0,NDEGS,RSEC,T,PN)
C
C  CALCULATES INDIVIDUAL LOW-PRESSURE RATES FOR MULTIPLE CHANNEL REACTIONS.
C
       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
       PARAMETER (NMAX100=10000, NMAX50=5000, NMAX25=2500)
       COMMON /AMIN/ ID,NCUT,NILOC
       COMMON /BLOCK1/ NRATES,NI,NE0,NEFF
       COMMON /BLOCK2/ DELT,CORRF
       COMMON /BLOCK3/ NREACT,NBAND,NRPL1
       COMMON /BLOCK6/ Q(NMAX100),H(NMAX100),RTRHB(NMAX100)
       COMMON /BLK13/ AROW(NMAX100),OMEGA,ANORM(NMAX100)
       COMMON /BLK15/ NCH2
       COMMON /B2/ PROB(NMAX25)
       COMMON /DENS/ RHO(NMAX100)
       COMMON /F1/ IFIRST,NSPEC,IMIN
       COMMON /F2/ ALMT,ALMT1,DET
       COMMON /JLP1/ JAVX
       COMMON /JLP2/ CRF2(NMAX100),EXL2(NMAX100)
       COMMON /OPT/ IOPTPR
       COMMON /NEWKIDS/ TESTINC
C
       DIMENSION R1(NMAX100),R2(NMAX100,2),RTE2(2),SUM(2),
     * SUMR(NMAX100,2),CHMRT(2)
C
       LOGICAL JAVX
C
       DO ISUM=1,NCH2
         RTE2(ISUM)=0.
         CHMRT(ISUM)=0.
       ENDDO
       RTOT=0.
       DELT=100./TESTINC
       EXPON = 1.43879*DELT/T
       NIT=INT(NE0*100/DELT)
       NIP1=NIT+1
       S=0.
       DO I=1,NRATES
         SUMX=R1(I)
         DO ISUM=1,NCH2
           SUMX=SUMX+R2(I,ISUM)
         ENDDO
         DO ISUM=1,NCH2
           SUMR(I,ISUM)=R2(I,ISUM)/SUMX
         ENDDO
       ENDDO
       IF(.NOT.JAVX)NEFF=NIT
       DO 9040 I=IMIN,NIT
         BBI=Q(I)*RTRHB(I)
         IF(JAVX )BBI=BBI*DSQRT(1.D0-CRF2(I))
         S=S+BBI
         DO ISUM=1,NCH2
           SUM(ISUM)=0.
         ENDDO
         SUMX=0.
         IF(I.LE.NREACT) GO TO 9040
         IF(NIP1-I.GT.NBAND) GO TO 9010
         JMAX=MIN0(NDEGS,NIP1+NBAND)
         IF(NIP1.GE.JMAX) GO TO 9040
         DO 9000 J=NIP1,JMAX
           II=J-I+1
           IF(II.LE.0) GO TO 9000
           TERM=(RHO(J)/RHO(I))*DEXP(DFLT(I-J)*EXPON)*
     *                                       PROB(II)*ANORM(J)
           SUMX=SUMX+TERM
           IF(R2(J-NEFF,1).LE.0.) GO TO 9000
           DO ISUM=1,NCH2
             SUM(ISUM)=SUM(ISUM)+TERM*SUMR(J-NEFF,ISUM)
           ENDDO
 9000    CONTINUE
 9010    CONTINUE
         IF(.NOT.JAVX)GO TO 9030
         IB=MAX0(NEFF+1,I-NBAND+1)
         IT=MIN0(NIT,I+NBAND-1)
         IF(I.GE.IT)GO TO 9020
         IF(I.LE.NEFF)GO TO 9030
         DO II=I+1,IT
           IK=II-I+1
           TERM=(RHO(II)/RHO(I))*DEXP(DFLT(I-II)*EXPON)*PROB(IK)
     *                                         *ANORM(II)*CRF2(II)
           DO ISUM=1,NCH2
             SUM(ISUM)=SUM(ISUM)+TERM*SUMR(II-NEFF,ISUM)
           ENDDO
           SUMX=SUMX+TERM
         ENDDO
 9020    CONTINUE
         DO II=IB,I
           IK=I+1-II
           TERM=PROB(IK)*ANORM(I)*CRF2(II)
           DO ISUM=1,NCH2
             SUM(ISUM)=SUM(ISUM)+TERM*SUMR(II-NEFF,ISUM)
           ENDDO
           SUMX=SUMX+TERM
         ENDDO
 9030    CONTINUE
         DO ISUM=1,NCH2
           RTE2(ISUM)=RTE2(ISUM)+BBI*SUM(ISUM)
         ENDDO
         RTOT=RTOT+BBI*SUMX
 9040  CONTINUE
C
C  RSEC IS TOTAL SECOND-ORDER LOW PRESSURE RATE COEFFT; HERE WE HAVE
C  CALCULATED RATE 2, RATE 1 IS FOUND BY DIFFERENCE.
C
       SUMCH=0.
       DO ISUM=1,NCH2
         CHMRT(ISUM)=RTE2(ISUM)/RTOT
         SUMCH=SUMCH+CHMRT(ISUM)
         RTE2(ISUM)=CHMRT(ISUM)*RSEC
       ENDDO
       RATE1=(1.-SUMCH)*RSEC
       WRITE(6,9001) RATE1,(RTE2(I7),I7=1,NCH2)
C
 9001  FORMAT(//,'  MULTIPLE CHANNEL LOW-PRESSURE RATE COEFFICIENTS',
     * ' (CM3 MOL-1 S-1):',/,15X,1P,3E10.2)
C
       IF(JAVX)WRITE(6,9002)
C
 9002  FORMAT(/,' NOTE: THIS RATIO WILL ONLY BE VALID IF THE ',
     * 'FIRST FILE OF MULTICHANNEL K(E) ',/,' INPUT CORRESPONDS ',
     * 'TO A SUFFICIENTLY LOW PRESSURE.(SEE PROGRAM NOTES).',/)
C
       RETURN
       END
C
C***********************************************************************
C
      SUBROUTINE COLLF(DELT)
C
C  FINDS TOTAL COLLISION RATE OUTSIDE ROTATIONAL THRESHOLD.
C
       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
       PARAMETER (NMAX100=10000, NMAX50=5000, NMAX25=2500)
       COMMON /BLOCK1/ NRATES,NI,NE0,NEFF
       COMMON /BLOCK3/ NREACT,NBAND,NRPL1
       COMMON /BLK13/ AROW(NMAX100),OMEGA,ANORM(NMAX100)
       COMMON /B2/ PROB(NMAX25)
       COMMON /JLP3/ CRF(NMAX100),EXL(NMAX100)
       COMMON /NEWKIDS/ TESTINC
C
       DO I=1,(TESTINC*1000)
         EXL(I)=0.
       ENDDO
       DO N=NREACT,NI
         IMAX=MIN0(NI,N+NBAND-1)
         IMIN=MAX0(NREACT,N-NBAND+1)
         IP=N+1
         IM=N-1
         S33=0.
         IF(N.EQ.NI)GO TO 1021
         DO I=IP,IMAX
           IK=I-IM
           S33=S33+PROB(IK)*AROW(I)*CRF(I)*ANORM(I)
         ENDDO
         S33=S33*DELT*OMEGA/AROW(N)
 1021    CONTINUE
         S44=0.
         IF(N.EQ.NREACT)GO TO 1022
         DO I=IMIN,IM
           IK=IP-I
           S44=S44+PROB(IK)*CRF(I)
         ENDDO
         S44=S44+PROB(1)*CRF(N)
         S44=S44*OMEGA*DELT*ANORM(N)
 1022    CONTINUE
         EXL(N)=S44+S33
       ENDDO
       RETURN
       END
C
C*************************************************************************
C
       BLOCK DATA
C
       IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
       PARAMETER (NMAX100=10000, NMAX50=5000, NMAX25=2500)
       COMMON /BLOCK6/ Q(NMAX100),H(NMAX100),RTRHB(NMAX100)
       COMMON /BLOCK7/ D(NMAX100)
       COMMON /BLK10/ AVRATE(NMAX50)
       COMMON /BLK12/ AVR2(NMAX50,2)
       COMMON /BLK13/ AROW(NMAX100),OMEGA,ANORM(NMAX100)
       COMMON /B2/ PROB(NMAX25)
       COMMON /DENS/ RHO(NMAX100)
       COMMON /JLP2/ CRF2(NMAX100),EXL2(NMAX100)
       COMMON /JLP3/ CRF(NMAX100),EXL(NMAX100)
       COMMON /ROTTX/ ROTHR(NMAX100),ROTKT(NMAX100)
C
       LOGICAL JAVX
C
       DATA Q/NMAX100*0./,RTRHB/NMAX100*0./,D/NMAX100*0./
       DATA AVRATE/NMAX50*0./,AVR2/NMAX100*0./,AROW/NMAX100*0./
       DATA ANORM/NMAX100*0./,PROB/NMAX25*0./,ROTHR/NMAX100*0./
       DATA CRF/NMAX100*0./,RHO/NMAX100*0./
       END
C
C**************************************************************************
C
       DOUBLE PRECISION FUNCTION DFLT(I)
C
       DFLT = DBLE(FLOAT(I))
       RETURN
       END