C     PROGRAM MASTER(INPUT,OUTPUT,TAPE5=INPUT,TAPE6=OUTPUT)
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 RRKM2
C
C  LATEST UPDATE OF THIS PROGRAM: Oct 23 1991
C
C
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      REAL*8 KT,KCAPT
      COMMON /F1/ IFIRST,NSPEC,ALMT,ALMT1,DET
       COMMON /ROTTX/ ROTHR(500),ROTKT(1000),KT
      COMMON /VERSCH/ IXV,JXV
      COMMON /LOW/ LOWP
      COMMON /IMNX/ IMIN
      COMMON/BLOCK1/ NRATES,NI,NE0,NEFF
       COMMON /AMIN/ ID,NCUT,NILOC
       COMMON /TRANS/ ITR,RATTOT,RATE1,RATE2,RATEA(2),FNE
      COMMON/BLOCK2/ DELT,CORRF
      COMMON /BLOCK3/ NREACT,NBAND,NRPL1
      COMMON/EXP/ ERR1,ERR2,ERR3
      COMMON /BLOCK6/ Q(1000),H(1000),RTRHB(1000)
      COMMON /OPT/ IOPTPR
      COMMON/BLOCK7/ D(1000)
      COMMON/BLK10/ AVRATE(500)
      COMMON /BLK12/ AVR2(500,2),NCH2
      COMMON /BLK13/ AROW(1000),OMEGA,ANORM(1000)
      COMMON /BLK14/ NRCTH
      COMMON /B2/ PROB(250)
      COMMON /JLP/ CRF2(1000),JAVX,EXL2(1000)
      COMMON /DEOC/ EI,WE0
      COMMON /JLP2/ CRF(1000),EXL(500)
      COMMON /WEAKJ/ GAMMA,DELTA,IWKJ
      DIMENSION PR(20),R1(1000),R2(1000,2),R1S(1000),RHO(1000),
     1 ALPHAV(10),TITLE(20),TEMPV(10),CORRAT(10),RATHP(3),GAMMAV(10)
      LOGICAL LOWP,JAVX,BRW,ION
      EXTERNAL ES,QU
      DATA R1/1000*0./,R2/2000*0./,R1S/1000*0./,RHO/1000*0./
     1 ,WKA/349.689/,RATHP/3*0./
C      OPEN(5,FILE='MASFIL',STATUS='OLD')
C      OPEN(6,FILE='MASOUT',STATUS='UNKNOWN')
C  READ IN AND PRINT OUT TITLE.
      READ(5,913) TITLE
 913   FORMAT(20A4)
1906  FORMAT(1P,6X,'PRESSURE  =',T20,E12.2,' TORR (=',E10.2,' PASCAL)'
     1 ,/,10X,'OMEGA  =',T20,E12.2,' SEC-1',///,
     2 ' EXACT (EIGENVALUE) RATE   STRONG COLLISION ',
     3 3X,'TOTAL','   DELTA ',' BAND-',' NO.NOT',4X,'NUMBER',/,
     4 5X,'K1',4X,'K2,3',4X,'TOTAL',6X,'K1',4X,'K2,3',
     5 5X,'DIMENSION',3X,'E',3X,'WIDTH',' REACTING',1X,'ITERATIONS')
 9905 FORMAT(1P,3E9.2,1X,E7.1,1X,E7.1,0P,I7,F10.0,I5,
     1 I7,3X,I6)
  5   FORMAT(/,' NO. NOT',T12,'DELTA',T20,'BAND-',T30,'MATRIX',
     1 T44,'TOTAL',T54,'NO.OF',/, ' REACTING',T13,'E',
     2 T20,'WIDTH',T31,'SIZE',T40,' RATE (K0)',T52,'ITERATIONS')
 5903  FORMAT(I5,5X,F5.0,2X,I5,5X,I5,6X,1P,E10.2,0P,I9)
      WRITE(6,3912) TITLE
 3912  FORMAT(1H1,20X,20A4,/,
     1 10X,'MASTER EQUATION SOLUTION FOR FALL-OFF AND LOW PRESSURE',
     2 ' REGIMES')
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  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.
      READ(5,*) INC,NCHAN,INCCHK
      IDM=2
      NCH2=MAX0(NCHAN,2)-1
      DELTS=DBLE(INC)
      IF(INCCHK.EQ.100) GO TO 3920
      WRITE(6,3921)
 3921 FORMAT(' RRKM PROGRAM WAS NOT CALCULATED WITH INC=100,ABORT')
      STOP
 3920 CONTINUE
C  THE FOLLOWING IS THE FIXED VALUE OF THE CONVERGENCE-GOVERNING PARAMETER.
      B=0.1
 5006  FORMAT(I1,I6,8X,2I5)
 7241  FORMAT(10X,'INITIAL GRAIN SIZE =',T50,I10,' CM-1')
C  READ IN ERRORS
      READ(5,*) ERR1,ERR2,ERR3
      WRITE(6,4784) ERR1,ERR2,ERR3
 4784  FORMAT(1P,/,10X,'ERROR FOR TRUNCATION OF MATRIX',T50,E10.1,
     1 /,10X,'EIGENVALUE TOLERANCE',T50,E10.1,/,
     2 10X,'ERROR FOR TRUNCATION OF P(E,E',1H',')',T50,E10.1)
      WRITE(6,1) NCHAN
  1   FORMAT(/,10X,'NUMBER OF REACTION CHANNELS',T58,I2)
      ITMAX=800
C  ITMAX IS THE MAXIMUM NUMBER OF NESBET ITERATIONS PERMITTED.
      READ(5,5002) E0
 5002  FORMAT(4F12.3)
      E0T=E0
      WE0=E0T*WKA
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.
      READ(5,*)JAV
      JAVX = JAV.NE.0
      READ(5,*) NALPHA
      READ(5,*) (ALPHAV(I),I=1,NALPHA)
      IF(NALPHA.GT.1)WRITE(6,4441)
 4441 FORMAT(//,' NALPHA IS GREATER THAN 1: ARE YOU ALSO CHANGING ',
     1 'GAMMA ?',/,' - IF SO, YOU WILL NEED TO RE-RUN THE RRKM ',
     2 'PROGRAM FOR EACH VALUE OF GAMMA.',//)
C  READ IN PARAMETERS GIVING THE
C  FUNCTIONAL FORM OF P(E,E').
      READ(5,*) IXV,JXV
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
      IF(JAVX)READ(5,*)(GAMMAV(I),I=1,NALPHA)
      READ(5,*) NP
      READ(5,*) (PR(I),I=1,NP)
C  IF FIRST PRESSURE IS LT 0, THEN THE INPUT PRESSURE VALUES ARE
C  IGNORED, AND THE PRESSURES ARE CHOSEN SO THAT K0*OMEGA/K(INF)
C  IS 10**N, WHERE N=0,1,-1,2,-2,... UP TO NP VALUES.
C  THIS PROCEDURE IS NOT VALID FOR JAVX=TRUE,WHERE K(E)=
C  ITSELF CHANGES WITH PRESSURE.
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).
      READ(5,*) PALMT
C  PALMT IS THE PARAMETER SUCH THAT THE POPULATION VECTOR FOR
C  ALL ENERGIES LT E0*PALMT ARE GIVEN THEIR EQUILIBRIUM VALUES.
      IF(JAVX) WRITE(6,314)
 314   FORMAT(10X,'FULL ANGULAR MOMENTUM CONSERVATION USED (SMITH-',
     1 'GILBERT METHOD)')
      IF(.NOT.JAVX) WRITE(6,315)
 315  FORMAT(10X,'THIS CALCULATION DOES NOT USE FULL ANGULAR MOMENTUM',
     1 ' CONSERVATION;',/,10X,'ANGULAR MOMENTUM CONSERVATION ONLY',
     2 ' IN HIGH PRESSURE LIMIT')
      IF((NALPHA.EQ.1).OR.(NP.EQ.1).OR.(.NOT.JAVX))GO TO 612
      WRITE(6,46)
   46 FORMAT(/,' DUE TO THE STRUCTURE OF THE INPUT FILE,',/,
     1 ' YOU MAY NOT LOOP SEVERAL ALPHAS AND SEVERAL PRESSURES ',
     2 ' IN THE SAME RUN.',/,' MAKE EITHER NALPHA OR NP =1',/,' ABORT.')
      STOP
  612 CONTINUE
      READ(5,*)NTEMP,(TEMPV(I),I=1,NTEMP)
      IF(.NOT.JAVX)GO TO 121
      IF(PR(1).GE.0)GO TO 121
      WRITE(6,123)
  123 FORMAT(/,'OPTION TO GENERATE FALL-OFF AUTOMATICALLY IS NOT',
     1 ' VALID WHEN USING J-AVERAGED K(E)S. ABORT')
      STOP
  121 CONTINUE
C  READ IN COLLISION DIAMETER (ANGSTROM),COLLISION MASSES
      READ(5,*) SGMA,WT1,WT2,EPS
      ION=.FALSE.
      IF(SGMA.LT.0)ION=.TRUE.
      IF(ION)GO TO 6711
      WRITE(6,4731) SGMA,WT1,WT2
      IF(EPS.LE.0.) WRITE(6,2210)
      IF(EPS.GT.0.) WRITE(6,2211) EPS
 2210 FORMAT(10X,'HARD SPHERE COLLISION FREQUENCY USED',/)
 2211 FORMAT(10X,'LENNARD JONES COLLISION FREQUENCY USED',/,
     1 10X,'WELL DEPTH =',T50,F10.1,' K',/)
 4731  FORMAT(/,10X,'COLLISION DIAMETER =',T50,F10.2,' ANGSTROMS',
     1 /,10X,'MASS OF REACTANT=',T50,F10.2,' A.M.U.'
     2,/,10X,'MASS OF BATH GAS  ='
     3 ,T50,F10.2,' A.M.U.',/)
      GO TO 6712
 6711 CONTINUE
      KCAPT=DABS(SGMA)
      WRITE(6,6715)KCAPT
 6715 FORMAT(/,' FOR ION-MOLECULE REACTIONS, PARENT ION/NEUTRAL BATH',
     1 ' GAS',/,' CAPTURE RATE IS THE COLLISION RATE OMEGA:',
     2 //,T30,'K(CAPTURE) =',1P,D11.4,' CM3 S-1')
 6712 CONTINUE
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).
      READ(5,*) IOPTHT,IOPTPR
      IOPTPR=IABS(IOPTPR)
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.
      WRITE(6,8816) NP
 8816  FORMAT(10X,'NO. OF PRESSURES =',T50,I10,/)
 2872  FORMAT(/,10X,'NO. OF VALUES OF ALPHA =',T50,I10,/,
     3 10X,'NO. OF TIMES,+1,GRAIN SIZE HALVED =',T50,I10)
      IF(IOPTHT.NE.0) GO TO 2862
      IDM=2
      DELTS=200.
 2862  CONTINUE
      INC=DELTS
      WRITE(6,7241) INC
      WRITE(6,2872) NALPHA,IDM
      BRW=IXV.LE.0.AND.JXV.LE.0
C  BRW IS LOGICAL VARIABLE SPECIFIYING THAT BIASED RANDOM WALK MODEL
C  USED FOR P(E,E').
      WRITE(6,8152) PALMT
 8152  FORMAT(10X,'CUTOFF PARAMETER FOR H(E)=1',T50,F10.2)
      I=IABS(IOPTPR)
      WRITE(6,4) I
  4   FORMAT(/,10X,'NUMBER OF REACTION CHANNELS =',T50,I10)
      READ(5,1905) NDEGS
      IF(NDEGS.LE.1000) GO TO 4929
      WRITE(6,3928)
 3928  FORMAT(' MORE THAN 1000 ELEMENTS IN RHO INPUT,ABORT')
      STOP
 4929   CONTINUE
      READ(5,1207) (RHO(I),I=1,NDEGS)
C  RENORMALIZE RHO
      DO 24 I=1,NDEGS
      IF(RHO(I).GT.0.) GO TO 25
  24  CONTINUE
  25  DNM=RHO(I)
      DO 26 I=1,NDEGS
  26  RHO(I)=RHO(I)/DNM
      IF(JAVX)GO TO 117
      DO 119 I119=1,NTEMP
      READ(5,*)T,CORRAT(I119)
  119 CONTINUE
      NE0=INT(WE0/100.)
      GO TO 411
  117 CONTINUE
C  READ IN ARRAY OF ROTATIONAL THRESHOLD ENERGIES.
      READ(5,1905)NTHR
      READ(5,299)(ROTHR(I),I=1,NTHR)
  299 FORMAT(6E11.4)
      NE0=NTHR
      NEON2=(NE0/2)-1
  411 CONTINUE
      NP=NP+1
      DO 720 ITEMPR=1,NTEMP
      T=TEMPV(ITEMPR)
      WRITE(6,8635) T
      KT=T/1.4387
      IF(.NOT.JAVX)GO TO 412
      DO 164 I164=NEON2,NE0
      ROTKT(I164) = ROTHR(I164)/KT
      ROTKT(I164-NEON2+1)=50.
      CRF(I164)=0.
      IF(ROTKT(I164).LE.30.) CRF(I164)=DEXP(-ROTKT(I164))
  164 CONTINUE
      DO 263 I263=NE0+1,NDEGS
      CRF(I263)=1.
      ROTKT(I263) = 0.
  263 CONTINUE
  412 CONTINUE
       ICHP=0
      P1=PR(1)
      DO 720 IALPHA=1,NALPHA
      L=-1-(NP/2)
      BETA=1.
      ALPHA=ALPHAV(IALPHA)
      IF(.NOT.JAVX)GO TO 2234
      GAMMA=GAMMAV(IALPHA)
      IWKJ=1
      IF(GAMMA.LT.11.)IWKJ=0
      IF(GAMMA.GT.600.)IWKJ=2
      DELTA=GAMMA*KT/(GAMMA+KT)
 2234 IF(.NOT.BRW) WRITE(6,974) ALPHA,IXV,JXV
      IF(JAVX.AND.(IWKJ.LT.2))WRITE(6,977)GAMMA,DELTA
  977 FORMAT(10X,'WEAK COLLISION (EXPONENTIAL DOWN) FORM FOR ',
     1 'P(R,R") :',/,15X,'GAMMA =',F8.2,10X,' DELTA =',F8.2,//)
      IF(JAVX.AND.(IWKJ.EQ.2))WRITE(6,978)
  978 FORMAT(10X,'STRONG COLLISION FORM FOR P(R,R") :',
     1 /,15X,'P(R,R") = EXP(-R/KT)',//)
      IPRT=0
      DO 700 N=1,NP
      IF(JAVX)GO TO 113
      IF((N.NE.1).OR.(IALPHA.NE.1).OR.(ITEMPR.NE.1))GO TO 297
      GO TO 115
  113 CONTINUE
      IF((IALPHA.NE.1).OR.(N.EQ.2))GO TO 297
      READ(5,*)(RATHP(I),I=1,NCHAN),CORAV,CORPF
      READ(5,*)STLPJ
      READ(5,*) E0EFF
      E0=E0EFF
  115 CONTINUE
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 100
C  CM-1.
C
C  MICROSCOPIC RATES AND DENSITIES OF STATES,FOR 100 CM-1 SPACING,ARE READ IN
C  IN R1 (K1), R2 (K2) AND RHO (RHO).
      READ(5,1905) NRATES
      READ(5,1207) (R1(I),I=1,NRATES)
 1207  FORMAT(8E10.3)
      DO 6592 I6592=1,NCH2
      READ(5,1207) (R2(I,I6592),I=1,NRATES)
 6592   CONTINUE
      WRITE(6,8265) NRATES,NDEGS
      IF(NRATES.LE.500) GO TO 5929
      WRITE(6,5928)
      STOP
 5929  CONTINUE
 5928  FORMAT(' MAXIMUM OF 500 INPUT RATES ONLY ALLOWED, ABORT.',
     1 /,' INCREASE DIMENSION FROM (500) TO RUN.')
 8265  FORMAT(/,10X,'NUMBER OF INPUT MICROSCOPIC RATES =',T50,I10,/,
     1 10X,'NUMBER OF INPUT DENSITIES OF STATES =',T50,I10)
      NRTSTM=NRATES
      EMAX=E0+0.2859*DBLE(NRATES)
      NIORIG=INT((EMAX*200./DELTS)/0.5718)-1
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
      NI=NIORIG
      EJM=EMAX*4.184
      ECM=E0*4.184
      WRITE(6,848) EMAX,EJM,E0,ECM
 848  FORMAT(//,10X,'MAXIMUM ENERGY CONSIDERED',T50,F10.1,' KCAL MOL-1'
     1 ,/,T45,'=',T50,F10.1,' KJ MOL-1'
     2 ,//,10X,'CRITICAL ENERGY',T50,F10.2,' KCAL MOL-1',/,
     3 45X,'=',T50,F10.2,' KJ MOL-1')
 8635 FORMAT(/////,10X,'TEMPERATURE =',T50,F10.1,' KELVINS',//)
  974 FORMAT(///,T30,' ALPHA =',F8.0,' CM-1',
     1 //,10X,'PROBABILITY FUNCTION IS R(E-E',1H',')/N(E',1H',
     2 '), WHERE',
     3 /,8X,12HR(E-E')=(X**,I2,
     4 11H) DEXP(-X**,I2,23H), WHERE X=(E-E')/ALPHA,//)
  297 CONTINUE
      LOWP=N.EQ.1
      NI=NIORIG
      PN=PR(MAX0(1,N-1))
      IF(JAVX)NEFF=INT(E0EFF*WKA/100.)
      IF(P1.LE.0.) PN=1.
      L=L+1
      IF(.NOT.LOWP.AND.P1.LE.0.) PN=CSET*(10.**L)
      IF((WT1.GT.0.).AND.(WT2.GT.0.))AMASS=1./((1./WT1)+(1./WT2))
      IF(ION)GO TO 6713
      OMEGA=(4.412D+7*SGMA*SGMA/DSQRT(AMASS))*PN/(SQRT(T))
C  CALCULATE LENNARD-JONES COLLISION FREQUENCY IF REQUIRED.
      IF(EPS.GT.0.) OMEGA=OMEGA/(0.636+0.567*DLOG10(T/EPS))
      GO TO 6714
C  CALCULATE PARENT ION/NEUTRAL BATH GAS CAPTURE RATE FOR
C  ION-MOLECULE REACTIONS.
 6713 TEM=T
      IF(EPS.LT.-.005)TEM=DABS(EPS)
      OMEGA=9.661D+18*PN*KCAPT/TEM
 6714 CONTINUE
C  THE FOLLOWING SET OF STATEMENTS IS FOR THE SPECIAL
CASE OF LOW PRESSURE LIMIT CALCULATION.
      IF(IOPTHT.NE.0.AND.N.EQ.1) GO TO 700
C  CALL SUBROUTINE TO CARRY OUT PREPARATION FOR LOW-PRESSURE CALCULATION.
      IF(JAVX.AND.(.NOT.LOWP))GO TO 244
      E0=E0T
      GO TO 246
  244 E0=E0EFF
  246 CONTINUE
      IF (LOWP) CALL SETUPL(E0,T,R1S,DELTS,NDEGS,ALPHA,IPRT,
     1 IOPTEM,RHO)
      DELT=DELTS*2.
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(N.EQ.1) NCUT=1
      DO 710 IDT=1,IDM
      ID=IDT
      DELT=DELT/2.
      IF((.NOT.LOWP).OR.(ID.EQ.1))GO TO 429
      NE01=NE0*INT(101./DELT)
      NI=NE01
  429 CONTINUE
      NREACT=INT(E0/(DELT*2.859E-3))
      IF(N.EQ.1)NREACT=NE0*100/(2*INT(DELT+0.1))
       ALMT1=PALMT*E0/(DELT*2.859E-3)
      EXPON=1.4387*DELT/T
      IF(NI.LE.1000) GO TO 73
      WRITE(6,72)
  72   FORMAT(' NI TRUNCATED ')
      NI=1000
  73   CONTINUE
      SUMSCT=0.
      SUMSC2=0.
C.....GENDEG IS SUBROUTINE WHICH CALCULATES DENSITIES OF STATES OF LEVELS
C.....RATES CALCULATES MICROSCOPIC RATE COEFFICIENTS FOR REACTING LEVELS
C  BOTH INTERPOLATE OR TABLE SELECT FROM INPUT RATES OR DENSITIES OF STATES
      DO 462 I462=1,500
      AVRATE(I462)=0.0
  462 CONTINUE
      IF(.NOT.LOWP) GO TO 15
      CALL RATES(R1S,R2)
      IF(JAVX)CALL GENCF
      GO TO 11
  15  CALL RATES(R1,R2)
  11  CONTINUE
      CALL GENDEG(NDEGS,RHO)
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.
       ALMT=0.
       DENST=0.
      IF(LOWP) BETA=1.
      DO 9047 I=1,NI
      CI2=0.
      IF(D(I).GT.0.) CI2=DEXP(DLOG(D(I))-DBLE(I)*EXPON)
      DENST=DENST+CI2
      RTRHB(I)=DSQRT(CI2)
C  RTRHB IS BOLTZMANN VECTOR (SYMMETRIZED),Q IS EIGENVECTOR.
      IF(ID.NE.1)GO TO 712
      Q(I)=RTRHB(I)
C  EXTRA SYMMETRIZATION IS REQUIRED FOR LOW PRESSURE J-AVERAGED
C   CALCULATION.
      IF(JAVX.AND.LOWP)Q(I)=Q(I)/DSQRT(1.-CRF2(I))
  712 CONTINUE
C  STRONG COLLISION VECTOR IS FIRST GUESS
      AA=1.
      TT=1.
      IF(I.LE.NREACT) GO TO 9046
      AK=AVRATE(I-NREACT)
      IF(JAVX.AND.LOWP)AK=AK+EXL2(I)
      AA=OMEGA/(OMEGA+AK)
      IF(.NOT.LOWP) TT=AA
C  IN THE FOLLOWING, THE FIRST GUESS AT THE EIGENVECTOR IS THE
C  STRONG COLLISION FORM.
      IF(ID.NE.1)GO TO 341
      IF(JAVX.AND.LOWP)Q(I)=Q(I)*(1.-CRF2(I))
      Q(I)=Q(I)*OMEGA*BETA/(OMEGA*BETA+AK)
  341 CONTINUE
C  THE STRONG COLLISION LOW PRESSURE RATE IS
C  COMPUTED IN SUBROUTINE STRONG. N.B., STRONG ALWAYS
C  USES A GRAIN SIZE OF 100 CM-1.
      STERM=CI2*AK*AA
      SUMSCT=SUMSCT+STERM
      SUMCH=0.
      DO 2931 ISUM=1,NCH2
      SUMCH=SUMCH+AVR2(I-NREACT,ISUM)
 2931 CONTINUE
      SUMSC2=SUMSC2+CI2*SUMCH*AA
 9046   CONTINUE
      ALMT=ALMT+TT*CI2
 9047 CONTINUE
      SUMSCT=SUMSCT/DENST
      SUMSC2=SUMSC2/DENST
      SUMSC1=SUMSCT-SUMSC2
      IF(ID.EQ.1) GO TO 735
C  FOR SUBSEQUENT GUESS TO Q(I), USE LINEARLY-INTERPOLATED RESULT FROM
C  PREVIOUS GRAIN SIZE.
      AROW(1)=0.5*Q(1)
      DO 740 I=2,NILST*2,2
      ION2=I/2
      AROW(I)=Q(ION2)
      AROW(I+1)=0.5*(Q(ION2)+Q(ION2+1))
 740  CONTINUE
      DO 761 I=NILST*2+1,NI
      AROW(I)=AROW(I-1)*0.5
  761 CONTINUE
      DO 745 I=1,NI
      Q(I)=AROW(I)
      IF(D(I).LE.0.) Q(I)=0.
      IF(I.LT.INT(ALMT1)) Q(I)=RTRHB(I)
  745 CONTINUE
C  INTERPOLATED EIGENVECTOR GUESS HAS NOW BEEN FOUND.
  735 CONTINUE
C  STORE D(I)*B(I) IN AROW.
      NILOC=NI
      IF(N.EQ.1) NILOC=INT(DBLE(NDEGS)*100./DELT)
      DO 7021 I=1,NILOC
      IF(N.NE.1) AROW(I)=RTRHB(I)**2
      IF(.NOT.LOWP) GO TO 7021
      AROW(I)=0.
      IF(D(I).GT.0.) AROW(I)=DEXP(DLOG(D(I))-DBLE(I)*EXPON)
 7021   CONTINUE
      NRPL1=NREACT+1
      IPRT=IPRT+1
      CALL ENEXPT(DELT,ALPHA,IPRT)
C  THIS SUBROUTINE GENERATES PROBABILITY MATRIX
C  THE RESULTING MATRIX IS STORED IN ARRAY PROB.
      IF(NBAND.LE.0) GO TO 710
C  NESBET IS SUBROUTINE WHICH CARRIES OUT COMPUTATION OF EIGENVALUE.
C  SET STARTING VALUE AT LOWEST NON-ZERO POPULATION.
      DO 10 I=1,20
      IF(AROW(I).GT.0.) GO TO 12
  10  CONTINUE
   12  IMIN=I
C  THIS LAST SET OF INSTRUCTIONS IS TO ALLOW FOR SOME RHO(E) ( D(I) ) BEING 0.
       CALL NESBET(B,ITMAX)
      IF(N.EQ.1.AND.ID.EQ.1) WRITE(6,1350)
 1350  FORMAT(//,10X,'THE FOLLOWING IS FOR THE LOW PRESSURE',
     1' LIMIT,WITH K=K0*OMEGA')
      PPAS=PN*133.3
      IF(ID.EQ.1.AND.(.NOT.LOWP)) WRITE(6,1906) PN,PPAS,OMEGA
      IF(LOWP.AND.ID.EQ.1) WRITE(6,5)
      IF(.NOT.JAVX)GO TO 131
      IF(LOWP)CORRF=CORPF
      IF(.NOT.LOWP)CORRF=CORAV
      GO TO 133
  131 CONTINUE
      CORRF=CORRAT(ITEMPR)
  133 CONTINUE
      RATE1=RATE1*CORRF
      RATE2=RATE2*CORRF
      RATTOT=RATTOT*CORRF
      SUMSC1=SUMSC1*CORRF
      SUMSC2=SUMSC2*CORRF
      SUMSCT=SUMSCT*CORRF
      IF(.NOT.LOWP) WRITE(6,9905) RATE1,RATE2,
     1 RATTOT,SUMSC1,SUMSC2,NI,DELT,NBAND,
     2 NREACT,ITR
      IF(LOWP)WRITE(*,4223)DELT,FNE
 4223 FORMAT(' LOW PRESSURE LIMIT, DELT =',F8.2,5X,'FNE =',F10.6)
      IF(.NOT.LOWP)WRITE(*,4224)PN,DELT,FNE
 4224 FORMAT(' PRESSURE =',D10.3,5X,'DELT =',F8.2,5X,'FNE =',F10.6)
      IF(LOWP) RSEC=RATTOT*TEM/(1.603E-5*PN)
      IF(LOWP) WRITE(6,5903) NREACT,DELT,NBAND,NI,RSEC,ITR
 1905 FORMAT(I4,5X,1P,2E10.2,0P,F6.0,2I4,1P,2E10.2,E10.3,0P,I4)
      NILST=NI
      NI=NI*2
  710    CONTINUE
      IF(N.NE.1.OR.IOPTHT.NE.0) WRITE(6,23) PN,OMEGA,RATE1
       IF(.NOT.LOWP.AND.NCHAN.GE.2)
     1 WRITE(6,2871)(RATEA(I23)*CORRF,I23=1,NCH2)
  23   FORMAT(1P,//,10X,25(' *'),/,
     1 10X,'FINAL RESULT: AT PRESSURE =',E10.2,
     2 ' TORR, OMEGA =',E10.2,' S-1',/,
     3 20X,' Kuni (S-1) =',T40,E10.2,' (CHANNEL 1)')
 2871 FORMAT(1P,T40,E10.2, ' (CHANNEL 2)',/,T40,E10.2,' (CHANNEL 3)')
      IF((.NOT.LOWP).AND.RATTOT.GT.1.05*SUMSCT) WRITE(*,28)
  28    FORMAT(5X,'WARNING,WEAK COLL. RESULT GREATER THAN STRONG'
     1 ,', INDICATES NUMERICAL INSTABILITY,'
     2 /,5X, ' REDUCE GRAIN SIZE AND/OR TOLERANCES')
C  CARRY OUT HIGH PRESSURE AND STRONG COLLISION LOW PRESSURE CALCULATIONS.
      NRCTH=INT(E0*WKA/100.)
      IF(ICHP.GT.0)GO TO 129
      IF(JAVX)GO TO 722
       CALL STRONG(HP,ALPT,NDEGS,T,R1,
     1 R2,RHO,NCH2)
       HP=HP*CORRF
       ALPT=ALPT*CORRF
      GO TO 129
  722 HP=RATHP(1)+RATHP(2)+RATHP(3)
      ALPT=STLPJ*PN/(TEM*6.237D+4)
      ALPT=ALPT/OMEGA
  129 CONTINUE
      ICHP=ICHP+1
C  CALL SUBROUTINE TO CALCULATE TROE QUANTITIES FLH, FSC, FSC.
      IF(.NOT.LOWP) CALL TROEF(NCH2,HP,ALPT,SUMSCT,IOPTHT,RATOTP,T,
     1 BETA,NDEGS,RHO,R1,R2)
      WRITE(6,112)
  112 FORMAT(
     3 /,25X,'*******************************',///)
      NI=NI/2
      IF(.NOT.LOWP) GO TO 700
      RATOTP=RATTOT/OMEGA
      CSET=HP/RATTOT
      NRATES=NRTSTM
      WRITE(6,1349) RSEC
      BETA=RATOTP/ALPT
 1349 FORMAT(1P,10X,'LOW PRESSURE K0 =',E10.2,' CM3 MOL-1 S-1')
      WRITE(6,4928) BETA
 4928  FORMAT(1P,/,20X,'BETA =',T56,E10.2)
      IF(JAVX)NEFF=INT(E0EFF*WKA/DELT)
      WRITE(6,*)
      IF(NCHAN.GT.1) CALL MULTCH(NCH2,R1,R2,E0,NDEGS,RSEC,RHO,T,PN)
       WRITE(6,5105)
 5105  FORMAT(////,19X,'-------FALL-OFF CALCULATION-------',///)
C  LOW PRESSURE CALC. NOW DONE, SO START FALL-OFF CALCS.
      IOPTPR=IOPTEM
 700   CONTINUE
720   CONTINUE
C      CLOSE(5,STATUS='KEEP')
C      CLOSE(6,STATUS='KEEP')
      STOP
      END
      SUBROUTINE NESBET(B,ITMAX)
C  SUBROUTINE TO  SOLVE EIGENVALUE PROBLEM.
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON/BLOCK1/ NRATES,NI,NE0,NEFF
      COMMON /AMIN/ ID,NCUT,NILOC
      COMMON /LOW/ LOWP
      COMMON /IMNX/ IMIN
      COMMON /F1/ IFIRST,NSPEC,ALMT,ALMT1,DET
      COMMON /OPT/ IOPTPR
      COMMON/BLOCK2/ DELT,CORRF
       COMMON /TRANS/ ITR,RATTOT,RATE1,RATE2,RATEA(2),FNE
      COMMON /BLOCK3/ NREACT,NBAND,NRPL1
      COMMON/EXP/ ERR1,ERR2,ERR3
      COMMON /BLOCK6/ Q(1000),H(1000),RTRHB(1000)
      COMMON/BLOK17/ AVEC(1000)
      COMMON/BLK10/ AVRATE(500)
      COMMON /BLK12/ AVR2(500,2),NCH2
      COMMON /BLK13/ AROW(1000),OMEGA,ANORM(1000)
      COMMON /B2/ PROB(250)
      COMMON /JLP/ CRF2(1000),JAVX,EXL2(1000)
      REAL*8 RATK(500,2),KT,S28(2),SS8(2)
      LOGICAL DONE,LOWP,JAVX,OUT
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 (DENSITY OF STATES)*DEXP(-E/KT)
C  GAMMA IS Q(I)/RTRHB(I)
      RATEL=0.
      IFIRST=INT(ALMT1)
      IF(IFIRST.LT.IMIN) IFIRST=IMIN
      DO 722 I722=1,2
      DO 724 I724=1,500
      RATK(I724,I722)=0.
  724 CONTINUE
  722 CONTINUE
      RATEL2=0.
      DEOMGA=OMEGA*DELT
       DET=0.
      DET2=0.
      SS=0.
C  THE FOLLOWING LOOP ALSO GENERATES FIRST GUESS AT RATE COEFFICIENTS
      DO 1041 I1041=IMIN,NI
      H(I1041)=1.
      AVEC(I1041)=0.
      IF(AROW(I1041).LE.0.) GO TO 1041
      AVEC(I1041)=AROW(I1041)*ANORM(I1041)
      CI2=Q(I1041)*RTRHB(I1041)
      IF(JAVX.AND.LOWP)CI2=CI2*DSQRT(1.-CRF2(I1041))
      SS=SS+CI2
      H(I1041)=Q(I1041)/RTRHB(I1041)
      IF(I1041.LT.IFIRST) H(I1041)=1.
      IF(I1041.LE.NREACT) GO TO 1041
      EFR=AVRATE(I1041-NREACT)
      IF(JAVX.AND.LOWP)EFR=EFR+EXL2(I1041)
      DET=DET+EFR*CI2
      SUMCH=0.
      DO 2931 ISUM=1,NCH2
      SUMCH=SUMCH+AVR2(I1041-NREACT,ISUM)
 2931 CONTINUE
      DET2=DET2+SUMCH*CI2
 1041   CONTINUE
      IF(IOPTPR.LT.0) SS=ALMT
      RATE=DET/SS
      IF(ID.NE.1) RATE=RATTOT
      RATE2=DET2/SS
C  INITIALIZATION PROCESS NOW COMPLETE.
C
C  COMMENCE NESBET ITERATIVE PROCESS.
       DO 864 ITRX=1,ITMAX
      ITR=ITRX
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.
       DONE=ABS(1.-(RATEL/RATE)).LE.ERR2
      IF(IABS(IOPTPR).EQ.2) DONE=DONE.AND.ABS(1.-(RATEL2/RATE2)).LE.
     1 ERR2
      IF(DONE) GO TO 865
      RATEL2=RATE2
      RATEL=RATE
      DO 15 I=IFIRST,NI
      IF(AROW(I).LE.0.) GO TO 15
      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 914
      JMAX=MIN0(NI,I+NBAND-1)
C  COMMENCE LOOP OVER J GE I.
      IF(JMAX-IL1.LE.0) GO TO 914
      HT=H(JMAX)
      IF(JAVX.AND.LOWP)HT=HT*DSQRT((1.-CRF2(JMAX))/CRFI)*
     1 QU(I,JMAX,DELT)
      S=0.5*(HI-HT)*AVEC(JMAX)*PROB(JMAX-IL1)
      JM=JMAX-1
      IF(JM.LT.IP) GO TO 914
      DO 913 J913=IP,JM
      HT=H(J913)
      IF(JAVX.AND.LOWP)HT=HT*DSQRT((1.-CRF2(J913))/CRFI)*
     1 QU(I,J913,DELT)
      S=S+(HI-HT)*AVEC(J913)*PROB(J913-IL1)
 913   CONTINUE
 914   CONTINUE
      S=S/RTI
      S1=0.
      IF(I.LE.JMIN) GO TO 38
C  COMMENCE LOOP OVER J LE I.
      IF(IP-JMIN.LE.0) GO TO 38
      HL=H(JMIN)
      IF(JAVX.AND.LOWP)HL=HL*DSQRT(CRFI/(1.-CRF2(JMIN)))*
     1 QU(JMIN,I,DELT)
      S1=0.5*PROB(IP-JMIN)*(HI-HL)
      JM=JMIN+1
      IM=I-1
      IF(IM.LT.JM) GO TO 38
      DO 917 J=JM,IM
      HL=H(J)
      IF(JAVX.AND.LOWP)HL=HL*DSQRT(CRFI/(1.-CRF2(J)))*
     1 QU(J,I,DELT)
      S1=S1+PROB(IP-J)*(HI-HL)
  917 CONTINUE
  38    AA=-RATE
      S=S+S1*RTI*ANORM(I)
      IF(JAVX.AND.LOWP)S=S+PROB(1)*QI*ANORM(I)*ES(I,I,DELT)
      IF(I.LE.NREACT) GO TO 1
      ILNR=I-NREACT
      AKE=AVRATE(ILNR)
      AA=AA+AKE
 1     CONTINUE
      SGMAI=AA*QI+DEOMGA*S
      IF(.NOT.LOWP)DQI=-SGMAI/(OMEGA+AA)
      IF(LOWP)DQI=-SGMAI/(OMEGA-RATE)
      IF(B* ABS(DQI).GT.QI) GO TO 15
      Q(I)=QI+DQI
      H(I)=Q(I)/RTI
      DET1=DQI*RTI
      IF(JAVX.AND.LOWP)DET1=DET1*DSQRT(CRFI)
      IF(I.LE.NREACT) GO TO 6
      EFR=AKE
      IF(JAVX.AND.LOWP)EFR=EFR+EXL2(I)
      DET=DET+EFR*DET1
      SUMCH=0.
      DO 4931 ISUM=1,NCH2
      SUMCH=SUMCH+AVR2(ILNR,ISUM)
 4931 CONTINUE
      DET2=DET2+SUMCH*DET1
  6   SS=SS+DET1
      IF(IOPTPR.LT.0) SS=ALMT
      RATE2=DET2/SS
       RATE=DET/SS
   15   CONTINUE
C  ITERATION FINISHED
  864   CONTINUE
      WRITE(6,866)
  866   FORMAT(' ++++++++++WARNING, DID NOT CONVERGE ++++++')
  865   CONTINUE
C  HAVING FOUND EIGENVECTOR,FIND EIGENVALUE FROM SIGMA(K(I)*G(I) FORMULA
C  COMPUTE RATE COEFFICIENTS AND TRUNCATE DIMENSIONS.
C  DET IS USED FOR TOTAL RATE, DET1 AND DET2 FOR PARTIAL RATES FROM THE TWO
C  PATHS.
      S=0.
      DET=0.
      SS=0.
      FNE1=0.
      FNE2=0.
      DO 2 I2=1,NCH2
      RATEA(I2)=0.
      SS8(I2)=0.
   2  CONTINUE
      DO 779 I779=IMIN,NI
      IF(Q(I779).GE.0.) GO TO 28
      Q(I779)=0.
      GO TO 779
  28  CONTINUE
      BBI=Q(I779)*RTRHB(I779)
      IF(JAVX.AND.LOWP)BBI=BBI*DSQRT(1.-CRF2(I779))
      S=S+BBI
      FNE1=FNE1+(BBI**2)/AROW(I779)
      FNE2=FNE2+AROW(I779)
C  AVRATE IS TOTAL,AND AVR2 CHANNEL 2, MICROSCOPIC RATES. RATE FOR
C  CHANNEL 1 FOUND BY DIFFERENCE.
      IF(I779.LE.NREACT) GO TO 779
      EFR=AVRATE(I779-NREACT)
      IF(JAVX.AND.LOWP)EFR=EFR+EXL2(I779)
      S2=BBI*EFR
      SS=DMAX1(S2,SS)
      DET=DET+S2
      DO 6931 ISUM=1,NCH2
      IF(EFR.GT.0.)RATK(I779-NREACT,ISUM)=AVR2(I779-NREACT,ISUM)/EFR
      S28(ISUM)=RATK(I779-NREACT,ISUM)*S2
      RATEA(ISUM)=RATEA(ISUM)+S28(ISUM)
      SS8(ISUM)=DMAX1(S28(ISUM),SS8(ISUM))
 6931 CONTINUE
      IF(I779-NREACT.LT.20) GO TO 779
      OUT=.TRUE.
      DO 444 ISUM=1,NCH2
      OUT=OUT.AND.(S28(ISUM).LT.ERR1*SS8(ISUM))
  444 CONTINUE
      IF(S2.LT.SS*ERR1.AND.(.NOT.LOWP).AND.OUT) GO TO 100
C THE FOLLOWING STATEMENT IS FOR COMPATABILITY WITH IBM AND CDC MACHINES
      IF(I779.EQ.NI) GO TO 100
  779 CONTINUE
C  TRUNCATE DIMENSIONS WHEN VALUE OF I GIVES WITHIN ERR1 OF INTEGRAND OF
C  RATE COEFFICIENT.
 100   NI=I779
      IF(IOPTPR.LT.0) S=ALMT
      RATTOT=DET/S
C Next lines added by RGG Oct 23 to guard for divide-by-zero
      IF(RATTOT.GT.0.) GO TO 5215
      WRITE(6,5216)
 5216 FORMAT(' INSTABILITY IN EIGENVALUE SUBROUTINE -',
     1 ' INCREASE EIGENVALUE TOLERANCE, ERR2',/,'  ABORT')
      STOP
 5215 CONTINUE
      FNE=(S**2)/(FNE1*FNE2)
      RATE2=0.
      DO 3 I3=1,NCH2
      RATEA(I3)=RATEA(I3)/S
      RATE2=RATE2+RATEA(I3)/RATTOT
  3   CONTINUE
      RATE1=(1.-RATE2)*RATTOT
      RATE2=RATE2*RATTOT
      RETURN
      END
      SUBROUTINE STRONG( HP,ALPT,NDEGS,T,
     1 R1,R2,RHO,NCH2)
C  CALCULATES HIGH PRESSURE RATE COEFFICIENT AND LOW PRESSURE
C  STRONG COLLISION K0.
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON /BLK14/ NRCTH
      DIMENSION R1(1000),R2(1000,2),RHO(1000)
      EXPON=143.87/T
      BOT=0.
      HP=0.
      ALPT=0.
      IN=0
      DO 1 I=1,NDEGS
      TERM=0.
      IF(RHO(I).GT.0.) TERM=DEXP(DLOG(RHO(I))-DBLE(I)*EXPON)
      BOT=BOT+TERM
      IF(I.LE.NRCTH) GO TO 1
      IN=IN+1
      SUMCH=0.
      DO 7067 I7067=1,NCH2
      SUMCH=SUMCH+R2(IN,I7067)
 7067 CONTINUE
      HP=HP+TERM*(R1(IN)+SUMCH)
  156 CONTINUE
      ALPT=ALPT+TERM
  1   CONTINUE
      HP=HP/BOT
      ALPT=ALPT/BOT
      RETURN
      END
      SUBROUTINE ENEXPT(DELT,ALPHA,IPRT)
C  COMPUTES COLLISION PROBABILITY ARRAY PROB.
C  TO CHANGE PROBABILITY FUNCTION,  CHANGE FUNCTIONAL FORM
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      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(1000),OMEGA,ANORM(1000)
      COMMON /B2/ PROB(250)
      COMMON /EXP/ ERR1,ERR2,ERR3
      LOGICAL LOWP,BRW
      NBAND=0
      BRW=IXV.LE.0.AND.JXV.LE.0
      NE00=NREACT
      IF(LOWP) NE00=NI
      Z=AROW(NE00+1)/AROW(NE00)
      Z=DLOG(Z)/DELT
      Z=Z*ALPHA*ALPHA
      FS2=4.*ALPHA*ALPHA
      NTOPX=MIN0(249,NILOC)
      DO 809 I809=1,NTOPX
      PROB(I809)=0.
  809  CONTINUE
C  GENERATE UN-NORMALIZED PROBABILITY MATRIX.
      BOT=0.
      DEDOWN=0.
      DO 1 I=1,NTOPX
      DE=DBLE(I-1)*DELT
      X=DE/ALPHA
      IF(BRW) GO TO 381
      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 380
  381 TOP=-(Z+DE)
      PROB(I)=DEXP(-TOP*TOP/FS2)
  380 CONTINUE
      IF(PROB(I).LE.ERR3.AND.X.GT.0.5) GO TO 2
            IF(I.GT.1) DEDOWN=DEDOWN+DE*PROB(I)
      IF(I.EQ.1) DEDOWN=0.5*DE*PROB(1)
      IF(I.GT.1) BOT=BOT+PROB(I)
      IF(I.EQ.1) BOT=0.5*PROB(1)
      NBAND=MAX0(NBAND,I)
  1    CONTINUE
  2   PROB(I)=0.
      PROB(250)=DEDOWN/BOT
C  GENERATE ANORM, THE VECTOR OF NORMALIZERS,INITIALLY AS C,ITS RECIPROCAL.
      NBAND=NBAND-1
      IF(NBAND.EQ.1) NBAND=0
      IF(NBAND.EQ.0) RETURN
      IF(.NOT.LOWP) NILOC=NI
COMMENCE FINITE DIFFERENCE SOLUTION OF INTEGRAL EQUATION FOR C.
      J=NILOC+1
      DO 5 JJ=1,NILOC
      J=J-1
      ANORM(J)=0.
      IF(AROW(J).LE.0.) GO TO 5
      JP=J+1
      JM=J-1
      IF(J.GE.(NREACT/NCUT)) GO TO 4
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
C  NO PHYSICAL EFFECT,SINCE THESE ALWAYS HAVE THEIR EQUILIBIRIUM POPULATIONS.
      ANORM(J)=ANORM(J+1)
      GO TO 5
  4   CONTINUE
      DE=PROB(1)
      IMIN=MAX0(1,J-NBAND)
      IF(JM.LT.IMIN) GO TO 3
      DO 7 I7=IMIN,JM
      DE=DE+PROB(JP-I7)
  7   CONTINUE
 3    ANORM(J)=DE*DELT
      IF(J.EQ.NILOC) GO TO 5
      IUP=MIN0(NILOC,J+NBAND)
      DE=0.
      IF(IUP.LT.JP) GO TO 5
      DO 6 I6=JP,IUP
      IF(ANORM(I6).LE.0.) GO TO 6
      DE=DE+(PROB(I6-JM)*AROW(I6)/ANORM(I6))
   6  CONTINUE
       ANORM(J)=ANORM(J)/(1.-(DELT*DE/AROW(J)))
  5   CONTINUE
      JS=NILOC+1
      DO 13 J13=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)
 13     CONTINUE
      DO 7791 I7791=1,NILOC
      IF(AROW(I7791).LE.0.) ANORM(I7791)=0.
 7791  CONTINUE
      IF(LOWP.AND.ID.EQ.5) GO TO 239
      IF(IPRT.NE.1) RETURN
      IF(BRW) WRITE(6,382) ALPHA
  382 FORMAT(//,10X,'BIASED RANDOM WALK MODEL FOR P(E,E',1H',')',/,
     1 20X,'S =',F10.1,' CM-1')
       WRITE(6,78) (PROB(I),I=1,NBAND)
  78   FORMAT(/,'  VALUES OF R(X) (FOR LARGEST GRAIN SIZE):',
     1 /,1P,10(3X,6E12.2,/))
      RETURN
COMPUTE AVERAGE ENERGY TRANSFER QUANTITIES
 239  WRITE(6,100)
  100 FORMAT(//,'  ENERGY (KJ/MOL) DELTA E TOTAL (CM-1)',
     1 ' DSQRT(DE**2)',T56,'EQUILIBRIUM POPULATION')
      NZ=NI/2
      SUMP=0.
      DO 9981 J9981=1,NILOC
      SUMP=SUMP+AROW(J9981)
 9981 CONTINUE
      DO 101 J101=NZ,NILOC,10
      E=DBLE(J101)*DELT*1.195E-2
      ANM=ANORM(J101)
      IF(ANM.LE.0.) GO TO 101
      JP=J101+1
      JM=J101-1
      DE=0.
      DE2=0.
      IF(JM.LT.1) GO TO 102
      DO 103 I103=1,JM
      IF(J101-I103.GE.NBAND) GO TO 103
      FAC=ANM*PROB(JP-I103)
      DIF=DBLE(I103-J101)
      DE=DE+DIF*FAC
      DE2=DE2+DIF*DIF*FAC
  103 CONTINUE
  102 IF(JP.GT.NILOC.OR.AROW(J101).LE.0.) GO TO 105
      RBJ=1./AROW(J101)
      DO 104 I104=JP,NILOC
      IF(I104-J101.GE.NBAND) GO TO 104
      IF(ANORM(I104).LE.0.) GO TO 104
      FAC=PROB(I104-JM)*ANORM(I104)*AROW(I104)*RBJ
      DIF=DBLE(I104-J101)
      DE=DE+DIF*FAC
      DE2=DE2+DIF*DIF*FAC
  104 CONTINUE
  105 DE=DE*DELT*DELT
      DE2=DSQRT(DE2*DELT)*DELT
      AROWN=AROW(J101)/(DELT*SUMP)
       WRITE(6,106) E,DE,DE2,AROWN
  106 FORMAT(F10.1,10X,F10.1,6X,F10.1,T62,1P,E10.2)
  101 CONTINUE
      RETURN
      END
      SUBROUTINE GENDEG(NDEGS,DSTORE)
C  FINDS VECTOR OF DEGENERACIES FOR ANY GRAIN SIZE BY INTERPOLATION
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON /BLOCK1/ NRATES,NI,NE0,NEFF
      COMMON /BLOCK2/ DELT,CORRF
      COMMON /BLOCK7/ D(1000)
      DIMENSION DSTORE(1)
C  BRANCH DEPENDING ON WHETHER DELT IS GT OR LT 100.
      IF(DELT.LT.99.99) GO TO 2
C  FIND INTEGER FOR TABLE SELECTION
      N=INT((DELT+0.01)*0.01)
      J=0
      DO 13 II=N,NDEGS,N
      I=II
      IF(I.GT.NDEGS) I=NDEGS
      J=J+1
      D(J)=DSTORE(I)
   13   CONTINUE
      RETURN
   2   N=INT(101./DELT)
      JJ=0
      ALAST=DLOG(DSTORE(1))
      DO 400 I=1,NDEGS
      IF(I.LT.NDEGS) ANEXT=DLOG(DSTORE(I+1))
      SP=(ANEXT-ALAST)/DBLE(N)
      DO 1000 J=1,N
      JJ=JJ+1
      D(JJ)=DEXP(ALAST+SP*DBLE(J-1))
 1000   CONTINUE
      ALAST=ANEXT
  400   CONTINUE
      RETURN
      END
      SUBROUTINE RATES(R1,R2)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON /LOW/ LOWP
      COMMON/BLOCK1/ NRATES,NI,NE0,NEFF
      COMMON/BLOCK2/ DELT,CORRF
      COMMON/BLK10/ AVRATE(500)
      COMMON /BLK12/ AVR2(500,2),NCH2
      COMMON /DEOC/ EI,WE0
      DIMENSION R1(1),R2(1000,2),ALAST2(2),ANEXT2(2),SP2(2)
      LOGICAL LOWP
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.
      IF(DELT.LT.99.99) GO TO 2
C  FIND INTEGER FOR TABLE SELECTION, EXTRA 1. IS TO ACCOUNT FOR ANY ROUNDOFF
      N=INT((DELT+1.)/100.)
      J=0
      IF(DELT.GT.101.)GO TO 463
      N3=1
      GO TO 462
  463 CONTINUE
      INTE0=INT(WE0/100.)
      IF(LOWP)INTE0=INT(WE0/200.)
      IF(MOD(INTE0,2).EQ.1)N3=1
      IF(MOD(INTE0,2).EQ.0)N3=2
  462 CONTINUE
      DO 13 II=N3,NRATES,N
      I=MIN0(II,NRATES)
      J=J+1
      DO 1 I1=1,NCH2
      AVR2(J,I1)=0.
   1   CONTINUE
      AVRATE(J)=R1(I)
      IF(LOWP) GO TO 13
      SUMCH=0.
      DO 3 I3=1,NCH2
      AVR2(J,I3)=R2(I,I3)
      SUMCH=SUMCH+AVR2(J,I3)
   3   CONTINUE
      AVRATE(J)=R1(I)+SUMCH
 13   CONTINUE
      RETURN
 2    N=INT(101./DELT)
      JJ=0
      ALAST1=(R1(1))
      DO 10 I10=1,NCH2
      ALAST2(I10)=0.
  10   CONTINUE
      DO 400 I=1,NRATES
      IF(I.LT.NRATES) ANEXT1=R1(I+1)
      DO 93 I93=1,NCH2
      ANEXT2(I93)=0.
  93   CONTINUE
      IF(I.EQ.NRATES) GO TO 100
      IF(LOWP.OR.R2(I+1,1).LE.0.) GO TO 14
      DO 15 I15=1,NCH2
      ANEXT2(I15)=R2(I+1,I15)
  15  CONTINUE
  14  CONTINUE
  100   CONTINUE
      SP1=(ANEXT1-ALAST1)/DBLE(N)
      DO 11 I11=1,NCH2
      SP2(I11)=(ANEXT2(I11)-ALAST2(I11))/DBLE(N)
  11    CONTINUE
      DO 1000 J=1,N
      JJ=JJ+1
      AVRATE(JJ)=(ALAST1+SP1*DBLE(J-1))
      SUMCH=0.
      DO 12 I12=1,NCH2
      AVR2(JJ,I12)=ALAST2(I12)+SP2(I12)*DBLE(J-1)
      SUMCH=SUMCH+AVR2(JJ,I12)
  12   CONTINUE
      AVRATE(JJ)=AVRATE(JJ)+SUMCH
 1000 CONTINUE
      ALAST1=ANEXT1
      DO 16 I16=1,NCH2
      ALAST2(I16)=ANEXT2(I16)
  16   CONTINUE
  400 CONTINUE
      RETURN
      END
      SUBROUTINE SETUPL(E0,T,R1S,DELTS,NDEGS,ALPHA,IPRT,IOPTEM,
     1 RHO)
C  PREPARATION FOR LOW-PRESSURE CALCULATION
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON /F1/ IFIRST,NSPEC,ALMT,ALMT1,DET
      COMMON /VERSCH/ IXV,JXV
      COMMON /IMNX/ IMIN
      COMMON/BLOCK1/ NRATES,NI,NE0,NEFF
       COMMON /AMIN/ ID,NCUT,NILOC
       COMMON /TRANS/ ITR,RATTOT,RATE1,RATE2,RATEA(2),FNE
      COMMON/BLOCK2/ DELT,CORRF
      COMMON /BLOCK3/ NREACT,NBAND,NRPL1
      COMMON/EXP/ ERR1,ERR2,ERR3
      COMMON /OPT/ IOPTPR
      COMMON/BLOCK7/ D(1000)
      COMMON/BLK10/ AVRATE(500)
      COMMON /BLK13/ AROW(1000),OMEGA,ANORM(1000)
      COMMON /B2/ PROB(250)
      COMMON /JLP2/ CRF(1000),EXL(500)
      COMMON /JLP/ CRF2(1000),JAVX,EXL2(1000)
      DIMENSION R1S(1000),RHO(1000)
      REAL*8 KT
      LOGICAL JAVX
C  SET UP RATE VECTOR AS COLLISION RATES
      DELT=100.
      NI=NE0
      EXPON=1.4387*DELT/T
      NI=MIN0(1000,NI)
      NREACT=NI/2
      NRPL1=NREACT+1
      ID=5
       IOPTEM=IOPTPR
      IOPTPR=-1
      DO 2391 I2391=1,NDEGS
      AROW(I2391)=0.
      IF(RHO(I2391).GT.0.) AROW(I2391)=
     1 DEXP(DLOG(RHO(I2391))-DBLE(I2391)*EXPON)
 2391  CONTINUE
      DO 3010 I3010=1,20
      IF(RHO(I3010).GT.0.) GO TO 3011
3010  CONTINUE
3011  IMIN=I
      NCUT=1
      NILOC=NDEGS
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.
      CALL ENEXPT(DELT,ALPHA,IPRT)
      NIP1=NI+1
      INDI=0
      DO 4016 I=NRPL1,NI
      INDI=INDI+1
      R1S(INDI)=0.
      IF(RHO(I).LE.0.) GO TO 4016
      SUM=0.
      IF(NIP1-I.GT.NBAND) GO TO 4016
      JMAX=MIN0(NDEGS,NIP1+NBAND)
      IF(NIP1.GE.JMAX) GO TO 4016
      DO 4015 J=NIP1,JMAX
      II=J-I+1
      IF(II.LE.0.OR.II.GT.NBAND) GO TO 4015
      TERM=(RHO(J)/RHO(I))*DEXP(DBLE(I-J)*EXPON)*
     1 PROB(II)*ANORM(J)
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     IF(J.EQ.NIP1.OR.J.EQ.JMAX) TERM=TERM*0.5
      SUM=SUM+TERM
 4015 CONTINUE
      R1S(INDI)=SUM*OMEGA*DELT
 4016 CONTINUE
      IF(.NOT.JAVX)GO TO 303
      CALL COLLF(CRF,EXL,DELT)
  303 CONTINUE
      NRATES=MIN0(1000,INDI)
      NI=INT(E0/(DELTS*2.859E-3))
      RETURN
      END
      SUBROUTINE TROEF(NCH2,HP,ALPT,SUMSCT,IOPTHT,RATOTP,T,
     1 BETA,NDEGS,RHO,R1,R2)
C  COMPUTES FUNCTIONS FLH, FSC, FWC DEFINED BY TROE.
C
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON /F1/ IFIRST,NSPEC,ALMT,ALMT1,DET
      COMMON /IMNX/ IMIN
      COMMON/BLOCK1/ NRATES,NI,NE0,NEFF
       COMMON /AMIN/ ID,NCUT,NILOC
       COMMON /TRANS/ ITR,RATTOT,RATE1,RATE2,RATEA(2),FNE
      COMMON/BLOCK2/ DELT,CORRF
      COMMON /BLOCK3/ NREACT,NBAND,NRPL1
      COMMON /OPT/ IOPTPR
      COMMON/BLOCK7/ D(1000)
      COMMON/BLK10/ AVRATE(500)
      COMMON /BLK13/ AROW(1000),OMEGA,ANORM(1000)
      COMMON /B2/ PROB(250)
      COMMON /BLK14/ NRCTH
      DIMENSION RHO(1000),R1(1000),R2(1000,2)
C
      WRITE(6,9201) PROB(250)
 9201 FORMAT(/,5X,'DELTA E DOWN (BY NUMERICAL INTEGRATION) =',T56,F10.1,
     1 ' CM-1')
      RATIOT=RATTOT/HP
      APOM=ALPT*OMEGA
      WRITE(6,113) HP,RATIOT,APOM
  113 FORMAT(1P,20X,'TOTAL K(INF) =',T56,E10.2,' S-1'
     1 ,/,20X,'K(TOT)/K(INF) ='
     2 ,T56,E10.2,/,20X,'STRONG COLLISION K0*OMEGA =',T56
     3 ,E10.2,' S-1')
      IF(APOM.LT.SUMSCT) WRITE(6,1898)
 1898 FORMAT(' WARNING, STRONG COLLISION RATE GT STRONG COLLISION',
     1 /,' VALUE FOR K0*OMEGA          USE SMALLER GRAIN SIZE')
      IF(IOPTHT.NE.0) RETURN
C  COMPUTE LINDEMANN-HINSHELWOOD AND STRONG COLLISION BROADENING FACTORS.
      WCK0=RATOTP*OMEGA
      IF(BETA.GT.1.) WRITE(6,2928)
 2928 FORMAT(' WARNING, BETA GREATER THAN 1, INDICATES'
     1 ,' TOLERANCES SHOULD BE ALTERED')
      WCRAT=WCK0/HP
      WRITE(6,114) WCK0,WCRAT
      WRITE(6,4928) BETA
 4928 FORMAT(1P,/,20X,'BETA =',T56,E10.2)
  114 FORMAT(1P,20X,'WEAK COLLISION K0*OMEGA =',T56,E10.2,' S-1',
     1 /,20X,'WEAK COLLISION K0*OMEGA/K(INF)=',T56,E10.2)
      FLH=WCRAT/(1.+WCRAT)
      OMSC=RATOTP*OMEGA/ALPT
      AKSC=0.
      EXPON=143.87/T
      BOT=0.
      ICSC=0.
      DO 630 I=1,NDEGS
      TERM=0.
      IF(RHO(I).GT.0.) TERM=DEXP(DLOG(RHO(I))-DBLE(I)*EXPON)
      BOT=BOT+TERM
      IF(I.LE.(NDEGS-NRATES)) GO TO 630
      ICSC=ICSC+1
      AKE=R1(ICSC)+(R2(ICSC,1)+R2(ICSC,2))
      AKSC=AKSC+(TERM*OMSC*AKE/(AKE+OMSC))
  630 CONTINUE
      AKSC=AKSC/BOT
      FSC=AKSC/(FLH*HP)
      FWC=RATTOT/(HP*FLH*FSC)
      WRITE(6,631) FLH,FSC,FWC
 631  FORMAT(1P,//,20X,'FLH =',E10.2,',FSC =',E10.2,',FWC =',E10.2)
      IF(FWC.GT.1.05) WRITE(6,2927)
 2927 FORMAT(' WARNING,FWC GREATER THAN 1, INDICATES THAT
     1 ',/,' TOLERANCES SHOULD BE CHANGED',
     2 ' OR THAT GRAIN SIZE SHOULD BE',/,' REDUCED TO 100 CM-1')
      RETURN
      END
      SUBROUTINE MULTCH(NCH2,R1,R2,E0,NDEGS,RSEC,RHO,T,PN)
C  CALCULATES INDIVIDUAL LOW-PRESSURE RATES FOR MULTIPLE CHANNEL REACTIONS.
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON /F1/ IFIRST,NSPEC,ALMT,ALMT1,DET
      COMMON /IMNX/ IMIN
      COMMON/BLOCK1/ NRATES,NI,NE0,NEFF
       COMMON /AMIN/ ID,NCUT,NILOC
      COMMON/BLOCK2/ DELT,CORRF
      COMMON /BLOCK3/ NREACT,NBAND,NRPL1
      COMMON /BLOCK6/ Q(1000),H(1000),RTRHB(1000)
      COMMON /OPT/ IOPTPR
      COMMON /BLK13/ AROW(1000),OMEGA,ANORM(1000)
      COMMON /B2/ PROB(250)
      COMMON /JLP/ CRF2(1000),JAVX,EXL2(1000)
      DIMENSION R1(1000),R2(1000,2),RHO(1000),RTE2(2),SUM(2),SUMR(500,2)
     1 ,CHMRT(2)
      REAL*8 KT
      LOGICAL JAVX
       DO 2 I2=1,NCH2
      RTE2(I2)=0.
      CHMRT(I2)=0.
  2    CONTINUE
      RTOT=0.
      DELT=100.
      EXPON = 1.4387*DELT/T
      NIT=INT(NE0*100/DELT)
      NIP1=NIT+1
      S=0.
      DO 199 I199=1,NRATES
      SUMX=R1(I199)
      DO 197 I197=1,NCH2
      SUMX=SUMX+R2(I199,I197)
  197 CONTINUE
      DO 195 I195=1,NCH2
      SUMR(I199,I195)=R2(I199,I195)/SUMX
  195 CONTINUE
  199 CONTINUE
      IF(.NOT.JAVX)NEFF=NIT
      DO 779 I=IMIN,NIT
      BBI=Q(I)*RTRHB(I)
      IF(JAVX)BBI=BBI*DSQRT(1.-CRF2(I))
      S=S+BBI
      DO 17 I17=1,NCH2
      SUM(I17)=0.
  17  CONTINUE
      SUMX=0.
      IF(I.LE.NREACT) GO TO 779
      IF(NIP1-I.GT.NBAND) GO TO 82
      JMAX=MIN0(NDEGS,NIP1+NBAND)
      IF(NIP1.GE.JMAX) GO TO 779
      DO 4015 J=NIP1,JMAX
      II=J-I+1
      IF(II.LE.0) GO TO 4015
      TERM=(RHO(J)/RHO(I))*DEXP(DBLE(I-J)*EXPON)*
     1 PROB(II)*ANORM(J)
      SUMX=SUMX+TERM
      IF(R2(J-NEFF,1).LE.0.) GO TO 4015
      DO 4 I4=1,NCH2
      SUM(I4)=SUM(I4)+TERM*SUMR(J-NEFF,I4)
  4   CONTINUE
 4015 CONTINUE
   82 IF(.NOT.JAVX)GO TO 81
      IB=MAX0(NEFF+1,I-NBAND+1)
      IT=MIN0(NIT,I+NBAND-1)
      IF(I.GE.IT)GO TO 69
      IF(I.LE.NEFF)GO TO 81
      DO 74 I74=I+1,IT
      IK=I74-I+1
      TERM=(RHO(I74)/RHO(I))*DEXP(DBLE(I-I74)*EXPON)*PROB(IK)
     1 *ANORM(I74)*CRF2(I74)
      DO 68 I68=1,NCH2
      SUM(I68)=SUM(I68)+TERM*SUMR(I74-NEFF,I68)
   68 CONTINUE
      SUMX=SUMX+TERM
   74 CONTINUE
   69 CONTINUE
      DO 67 I67=IB,I
      IK=I+1-I67
      TERM=PROB(IK)*ANORM(I)*CRF2(I67)
      DO 64 I64=1,NCH2
      SUM(I64)=SUM(I64)+TERM*SUMR(I67-NEFF,I64)
   64 CONTINUE
      SUMX=SUMX+TERM
   67 CONTINUE
   81 CONTINUE
      DO 5 I5=1,NCH2
      RTE2(I5)=RTE2(I5)+BBI*SUM(I5)
  5    CONTINUE
      RTOT=RTOT+BBI*SUMX
  779 CONTINUE
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 6 I6=1,NCH2
      CHMRT(I6)=RTE2(I6)/RTOT
      SUMCH=SUMCH+CHMRT(I6)
      RTE2(I6)=CHMRT(I6)*RSEC
   6   CONTINUE
      RATE1=(1.-SUMCH)*RSEC
       WRITE(6,1) RATE1,(RTE2(I7),I7=1,NCH2)
  1    FORMAT(//,'  MULTIPLE CHANNEL LOW-PRESSURE RATE COEFFICIENTS',
     1 ' (CM3 MOL-1 S-1):',/,15X,3D10.3)
      IF(JAVX)WRITE(6,772)
  772 FORMAT(/,'NOTE: THIS RATIO WILL ONLY BE VALID IF THE ',
     1 'FIRST FILE OF MULTICHANNEL K(E) ',/,'INPUT CORRESPONDS ',
     2 'TO A SUFFICIENTLY LOW PRESSURE.(SEE PROGRAMME NOTES).',/)
      RETURN
      END
      SUBROUTINE GENCF
C  SELECTS DEXP(-R0/KT) BY GRAIN SIZE.
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON /BLOCK1/ NRATES,NI,NE0,NEFF
      COMMON /JLP/ CRF2(1000),JAVX,EXL2(1000)
      COMMON /BLOCK2/ DELT,CORRF
      COMMON /JLP2/ CRF(1000),EXL(500)
      LOGICAL JAVX
      REAL*8 KT
      IF(DELT.LT.99.99)GO TO 222
      N=INT((DELT+0.01)*0.01)
      J=0
      NT=NE0+NRATES
      DO 433 II=N,NT,N
      I=II
      IF(I.GT.NT)I=NT
      J=J+1
      CRF2(J)=CRF(I)
      EXL2(J)=EXL(I)
  433 CONTINUE
      RETURN
  222 CONTINUE
      N=INT(101./DELT)
      JJ=0
      ALAST2=CRF(1)
      ALAST3=EXL(1)
      DO 437 I=1,NE0
      IF(I.LT.NE0)ANEXT2=CRF(I+1)
      IF(I.LT.NE0)ANEXT3=EXL(I+1)
      SP2=(ANEXT2-ALAST2)/DBLE(N)
      SP3=(ANEXT3-ALAST3)/DBLE(N)
      DO 1001 J=1,N
      JJ=JJ+1
      CRF2(JJ)=ALAST2+SP2*DBLE(J-1)
      EXL2(JJ)=ALAST3+SP3*DBLE(J-1)
 1001 CONTINUE
      ALAST2=ANEXT2
      ALAST3=ANEXT3
  437 CONTINUE
      RETURN
      END
      SUBROUTINE COLLF(CRF,EXL,DELT)
C  FINDS TOTAL COLLISION RATE OUTSIDE ROTATIONAL THRESHOLD.
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON /BLOCK1/ NRATES,NI,NE0,NEFF
      COMMON /BLK13/ AROW(1000),OMEGA,ANORM(1000)
      COMMON /ROTTX/ ROTHR(500),ROTKT(1000),KT
      COMMON /TEST/ I249
      COMMON /B2/ PROB(250)
      COMMON /BLOCK3/ NREACT,NBAND,NRPL1
      REAL*8 KT
      DIMENSION EXL(500),CRF(1000)
      DO 229 I229=1,500
      EXL(I229)=0.
  229 CONTINUE
      DO 249 I249=NREACT,NI
      IMAX=MIN0(NI,I249+NBAND-1)
      IMIN=MAX0(NREACT,I249-NBAND+1)
      IP=I249+1
      IM=I249-1
      S33=0.
      IF(I249.EQ.NI)GO TO 418
      DO 251 I251=IP,IMAX
      IK=I251-IM
      S33=S33+PROB(IK)*AROW(I251)*ANORM(I251)*ES(I251,I249,DELT)
  251 CONTINUE
      S33=S33*DELT*OMEGA/AROW(I249)
  418 CONTINUE
      S44=0.
      IF(I249.EQ.NREACT)GO TO 417
      DO 253 I253=IMIN,IM
      IK=IP-I253
      S44=S44+PROB(IK)*ES(I253,I249,DELT)
  253 CONTINUE
      S44=S44+PROB(1)*ES(I249,I249,DELT)
      S44=S44*OMEGA*DELT*ANORM(I249)
  417 CONTINUE
      EXL(I249)=S44+S33
  249 CONTINUE
      RETURN
      END
      REAL*8 FUNCTION ES(IE,IEP,DELT)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON /ROTTX/ ROTHR(500),ROTKT(1000),KT
      COMMON /WEAKJ/ GAMMA,DELTA,IWKJ
      REAL*8 KT
      NIE=0
      NIEP=0
      ES=0.
      N=NINT(DELT/100.)
      NIE=N*IE
      NIEP=N*IEP
      IF(ROTKT(NIEP).GT.30.)GO TO 24
      ZEP=1.-DEXP(-ROTKT(NIEP))
      ZE=1.-DEXP(-ROTKT(NIE))
      IF(IWKJ.NE.2)GO TO 61
      ES=1.-ZE
      RETURN
   61 CONTINUE
      IF(IEP.GT.IE)GO TO 21
      IF(IWKJ.NE.0)GO TO 71
      ES=1.-ZE/ZEP
      RETURN
   71 CONTINUE
      TM1=DEXP(ROTHR(NIE)/GAMMA)-1.
      TM1=TM1*DEXP(-ROTHR(NIEP)/DELTA)*GAMMA*DELTA
      TM1=TM1/(KT*(GAMMA+DELTA))
      ES=(TM1+ZEP-ZE)/ZEP
      RETURN
   21 CONTINUE
      IF(IWKJ.NE.0)GO TO 72
      ES=0.
      RETURN
   72 CONTINUE
      TM1=DEXP(ROTHR(NIEP)/GAMMA-ROTHR(NIE)/DELTA)
      TM1=TM1-DEXP(-ROTHR(NIEP)/DELTA)
      TM1=TM1*GAMMA*DELTA
      ES=TM1/(KT*(GAMMA+DELTA)*ZEP)
   24 RETURN
      END
      REAL*8 FUNCTION QU(IE,IEP,DELT)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      QU=1.-ES(IE,IEP,DELT)
      RETURN
      END
      BLOCK DATA
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
       COMMON /ROTTX/ ROTHR(500),ROTKT(1000),KT
      COMMON /BLOCK6/ Q(1000),H(1000),RTRHB(1000)
      COMMON/BLOCK7/ D(1000)
      COMMON/BLK10/ AVRATE(500)
      COMMON /BLK12/ AVR2(500,2),NCH2
      COMMON /BLK13/ AROW(1000),OMEGA,ANORM(1000)
      COMMON /B2/ PROB(250)
      COMMON /JLP/ CRF2(1000),JAVX,EXL2(1000)
      COMMON /JLP2/ CRF(1000),EXL(500)
      LOGICAL JAVX
      REAL*8 KT
      DATA Q/1000*0./,RTRHB/1000*0./,D/1000*0./,
     1 AVRATE/500*0./,AVR2/1000*0./,AROW/1000*0./,ANORM/1000*0./,
     2 JAVX/.TRUE./,PROB/250*0./,ROTHR/500*0./,CRF/1000*0./
      END