|
unimol
|
|
README,
brw.dat,
brw.for,
brw.out,
c2h5clgeom.dat,
c2h5cltsgeom.dat,
c2h6geom.dat,
c2h6tsgeom.dat,
c4h2geom.dat,
c4h2rrkm.dat,
c5h12geom.dat,
c5h12tsgeom.dat,
cf3geom.dat,
cf4geom.dat,
cf4tsgeom.dat,
ch3geom.dat,
ch3hcnrrkm.dat,
ch3ohgeom.dat,
ch4geom.dat,
ch4tsgeom.dat,
doc,
geom.dat,
geom.for,
geom.out,
h2ogeo.dat,
mas55.for,
mas77.for,
readme,
rrkm.for,
smpl1mas.dat,
smpl1mas.out,
smpl1out.out,
smpl1rrkm.dat,
smpl2mas.dat,
smpl2mas.out,
smpl2out.out,
smpl2rrkm.dat,
smpl3mas.dat,
smpl3mas.out,
smpl3out.out,
smpl3rrkm.dat,
smpl4mas.dat,
smpl4mas.out,
smpl4out.out,
smpl4rrkm.dat,
smpl5mas.dat,
smpl5mas.out,
smpl5out.out,
smpl5rrkm.dat
|
|
|
C
C PROGRAM RRKM(INPUT,OUTPUT,MASFIL,TAPE5=INPUT,TAPE6=OUTPUT,
C 1 TAPE10=MASFIL)
C
C RRKM MULTIPLE-REACTION PROGRAM, WITH PROVISION FOR OUTPUTTING MICROSCOPIC
C RATES,K(E), AND DENSITY OF STATES, RHOMOL(E), ONTO UNIT 10,
C SUITABLE FOR INPUT TO MASTER EQUATION PROGRAM.
C METHOD USED IS BEYER-SWINEHART ALGORITHM, WITH TROE MODIFICATION
C FOR COMPUTING EFFECTS OF FREE INTERNAL ROTORS.
C
C Authors: R.G. Gilbert, S.C. Smith, M.J.T. Jordan
C
C
C UPDATES:
C
C***********************************************************************
C AUGUST 1989 INCORPORATE SINUSOIDALLY HINDERED ROTOR OPTION
C (MJTJ) FOR ION-DIPOLE OR RADICAL-RADICAL REACTIONS.
C ACCESS WITH JAV.GT.1.
C ION.NE.0 SPECIFIES ION-DIPOLE.
C IHIND IS NO. 2-DIM SIN. HINDERED INTERNAL ROTORS
C (SYMMETRY ISYM1, ISYM2).
C IF IHIND.LT.0 THEN SINUSOIDAL ROTORS ALSO HAVE A
C STERIC CUTTOFF ANGLE (THETA1, THETA2).
C************************************************************************
C 10TH JUNE 1992 FIX LIN MOLECULE AND COMPLEX OPTIONS NINTR/N.LT.0.
C (MJTJ)
C*************************************************************************
C 2ND JULY 1992 FOR JAV.NE.0 AND WT1.LT.0 CAN ENTER NB
C (MJTJ) ROTATIONAL CONSTANTS AT SEPARATION RR.
C LINEAR INTERPOLATION ENABLES LOCATION OF CENT.
C BARRIER TO BE FOUND WHEN PIVOT ATOM IS NOT
C CENTRE OF MASS.
C*************************************************************************
C 29TH SEPTEMBER 1992 ENSURED NOTHING PRINTED IN COLUMN 1.
C (MJTJ)
C*************************************************************************
C 17TH NOVEMBER 1992 GROUPED VARIABLR TYPES TOGETHER IN COMMON BLOCKS.
C (MJTJ) PROGRAM ALL UPPER CASE.
C CONSTANTS ALL DOUBLE PRECISION.
C***********************************************************************
C 10TH DECEMBER 1992 CORRECTED EXPRESSION FOR RTHI.
C (MJTJ) CORRECTED CONVERGENCE FOR XKHP.
C**********************************************************************
C 5/8TH JANUARY 1993 CORRECTED STRONG COLLISION EXPRESSION FOR SCLX.
C (MJTJ) CORRECTED INEQUALITIES IN CALC OF RTHI.
C ENSURE EILO(I).GT.0.
C CONVERT ALL SUMS OVER E AND J TO INTEGRALS
C USING THE TRAPEZOIDAL RULE, IE HALF FIRST TERM.
C**********************************************************************
C 16TH MAY 1994 CORRECTED ION/DIPOLE FORMAT STATEMENT FOR DH0
C (MJTJ) CORRECTED LOOP IN JAVRGE CALCULATING SCLOW
C FOR JRECOM.EQ.1
C DISABLED CHECK FOR 2 2-DIM INTERNAL ROTORS IN
C ION/DIPOLE CASE
C**********************************************************************
C AUGUST 1994 ARRAY SIZES PARAMETERISED SO INC=10 CM-1 POSSIBLE.
C (Gary P Knight) ROTINC DEFAULT 50 CM-1 EXCEPT WHEN INC IS 10CM-1.
C WARNING THE .EXE FILE IS VERY LARGE.
C I(AMU ANGSTROM**2) TO B(CM-1) CONVERSION FACTOR
C NOW 16.8477.
C************************************************************************
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DOUBLE PRECISION K,KCAPT,KEQ(3),MORSE
PARAMETER (NMAX8=8000, NMAX5=5000, NMAX80=80000)
DIMENSION TITLE(20),
* IHIND(3),ION(3),
* ROTINC(3),ERR(3),
* THETA1(3),ISYM1(3),THETA2(3),ISYM2(3),GAMMA(3),
* NV(3),VCH(3,30),RVCH(3,30),
* IEXRTD(3),BVECM(20),SIGVCM(20),IRTDMM(20),NM(100),JM(100),
* WTA(3),WTB(3),NLTST(3),
* PRESS(30),TEMP(30),
* NCHK(3),NC1(100),JC1(100),BVEC1(20),SIGVC1(20),
* RATIO(3),ALNA(3),ALNAE(3),EAC(3),RATE(3),EJ(3),A(3),
* KKC(3),RATTH(3),AKG(3),RATHP(3),EKG(3),ISPARE(20),WRHP(3),
* CPC(3),SVC(3),SROTC(3),STC(3),STOTC(3),QROTC(3),QVC(3),STLP(3),
* AKHP(3),RHPN(3),REXQC(3),REXSRC(3),REXPC(3),
* CAPT(3),RATEMV(3),REXCPC(3),
* N(3,1),BVEC(3,1,20),SIGVC(3,1,20),IRTDMC(3,1,20),
* JF(3,1),NC(3,1,100),JC(3,1,100),V0(3,1),
* RR1(30),BB1(30)
C
LOGICAL IONX(3),IHINDX(3),STERIC(3),WKJ(3),
* TEST1D,LINM,LINC(3),CHECK,CHCKRC(3),
* BINPUT,JCBS,WKJ1,IHINDX1,IONX1,JAVX,STERICX
C
COMMON /ALL/ R,RT,H,WKA
COMMON /BFIT/ RR(3,30),BB(3,30),AA1(3),AA2(3)
COMMON /BINTEG/ NB(3)
COMMON /BLOGIC/ BINPUT
COMMON /BIG1/ K(3,NMAX8),AK(NMAX8),AK2(NMAX8),ROTHR(NMAX5)
COMMON /JCOLL/ DELTA,GAMMA1,OMEGA
COMMON /JHIGHP/ XKHP(NMAX8,3)
COMMON /JINT/ INCHAN,JRECOM,NJC
COMMON /JLIFE/ PLIFE,PPL1,PPL2,RADST,RSTAB,UNIR(3)
COMMON /JLOGIC/ JCBS(3),WKJ1
COMMON /JMICRO/ RHOMOL(NMAX8),ROTH(NMAX5),WE(NMAX8,3)
COMMON /JPARAM/ ROTINC1,ERR1
COMMON /JPOT/ BETA(3),DHO(3),DIP(3),EOK(3),POLZ(3)
COMMON /JRATES/ SCLE,SCLOW(2),RTHI(3),RTHI2(3),FWR(3)
COMMON /JROTC/ BCMPLX(3,1),RCPL(3,1),SRC(3,1)
COMMON /JROTM/ BMOLEC,BRATIO(3),REQ(3),SRM
COMMON /LOOP/ IN,NCHAN
COMMON /POT/ D,BE,REQ1
COMMON /ROTF/ GAMON2(40),ROTTM(3),V01,OMEGA1,OMEGA2
COMMON /ROTINT/ ISYMX1,ISYMX2,NH(3,1)
COMMON /STOR/ BVEC,SIGVC
COMMON /THML/ HRCORR,SQC
COMMON /LOGVAR/ IONX1,IHINDX1,JAVX,STERICX
COMMON /VFIT/ VCH1(30),RVCH1(30)
COMMON /INTVFIT/ NV1
C
EXTERNAL MORSE,DIFF
DATA THETA1/3*0./,THETA2/3*0./,
* BVECM/20*0./,SIGVCM/20*0./,IRTDMM/20*0/,NM/100*0/,JM/100*0/,
* IEXMD/3/,IEXRTD/3*3/,KEQ/3*0./,CAPT/3*0./,
* IBL/1H /,STC/3*0./,NMAX4/NMAX80/,N/3*0/,V0/3*0./
100 FORMAT(20A4)
105 FORMAT(' RRKM CALCULATION :',10X,20A4,/,
1 10X,'LATEST UPDATE August 1994')
115 FORMAT(//,20X,'NUMBER OF TERMS IN INTEGRATION',T55,I5,
1 /,20X,'ENERGY SPACING (CM-1)',T50,I10)
120 FORMAT(20X,'NUMBER OF INPUT PRESSURES:',T55,I5,
1 /,20X,'NUMBER OF INPUT TEMPERATURES:',T55,I5,/,
2 20X,'NUMBER OF CHANNELS:',T55,I5)
125 FORMAT(/,' FOR CHANNEL NUMBER',I4,':',/,
1' SINUSOIDALLY HINDERED NEUTRAL-NEUTRAL REACTION')
130 FORMAT(/,' FOR CHANNEL NUMBER',I4,':',/,
1' SINUSOIDALLY HINDERED ION-DIPOLE REACTION')
135 FORMAT(/,20X,'ROTATIONAL ENERGY GRAINSIZE:',T55,F10.4,
1/,20X,'CONVERGENCE PARAMETER:',T55,F10.4)
140 FORMAT(/,' CHANNEL NUMBER',I4,/,' HAS',I2,
1' SINUSOIDALLY HINDERED ROTORS')
145 FORMAT(/,' HARD-SPHERE STERIC HINDRANCE IS INCORPORATED',
1' INTO THE',/,' SINUSOIDALLY HINDERED ROTOR MODEL FOR',
2' CHANNEL',I4)
150 FORMAT(/,' ROTOR',T30,'SYMMETRY NUMBER',T50,'HINDRANCE ANGLE',
1' (DEGREES)',/,' 1',T35,I4,T60,F6.2)
155 FORMAT(' 2',T35,I4,T60,F6.2)
160 FORMAT(/,' CALCULATIONS USE ANGULAR MOMENTUM-CONSERVING',
1' TREATMENT: FULL K(E,J)')
C
C FORMAT STATEMENT 165 CORRECTED 23RD NOVEMBER 1992
C
165 FORMAT(/,' CALCULATIONS USE ANGULAR MOMENTUM-CONSERVING',
1' TREATMENT: FULL K(E,J):',/,T15,'WITH GAMMA =',F8.2,' CM-1')
170 FORMAT(/,' USING THE ION-DIPOLE OPTION FOR ONE CHANNEL',
1/, ' REQUIRES ITS USE FOR ALL OTHER CHANNELS. ABORT')
175 FORMAT(/,' CALCULATIONS USE ANGULAR MOMENTUM-CONSERVING',
1' TREATMENT: FULL K(E,J);'/,T15,'WITH STRONG COLLISIONS ',
2'FOR ROTATIONAL RELAXATION. ')
180 FORMAT(/,' FOR J-CONSERVING TREATMENT:',/,
1 ' R DAGGER (ANGSTROM)',T30,'R EQUILIBRIUM (ANGSTROM)',T60,
2 'CHANNEL',/,
3 3(F15.3,T35,F10.3,T60,I4,/))
185 FORMAT(/,' FOR J-CONSERVING TREATMENT:',/,
1' R DAGGER (ANGSTROM)',T60,'CHANNEL',/,
2 3(F15.3,T60,I4,/))
190 FORMAT(/,' FOR CHANNEL NUMBER',I4,':')
195 FORMAT(' NUMBER OF INPUT MORSE POTENTIAL POINTS =',I5)
200 FORMAT(/,' DISTANCE (ANGSTROM)',T30,'POTENTIAL (KCAL MOL-1)',
1 /,30(F15.3,T35,F10.3,/))
205 FORMAT(' CHEMICAL BARRIER AT R DAGGER =',F10.3,' ANGSTROMS',
1 /,' BARRIER HEIGHT =',F10.3,' KCAL MOL-1')
210 FORMAT(/,' CALCULATIONS USE MODEL THAT DOES NOT CONSERVE ANGULAR',
1 ' MOMENTUM IN FALL-OFF REGIME')
215 FORMAT(//,38X,'MOLECULE',T49,'COMPLEX(ES)',I2,2I15)
220 FORMAT(' CRITICAL ENERGY(KCAL/MOL)',T50,F12.3,T65,F12.3,F15.3)
225 FORMAT(/,' THE RATIO (BCMPLX/BMOLEC) FOR CHANNEL ',I2,
1 ' IS GREATER THAN 0.9',/,' THIS IMPLIES A TIGHT TRANSITION',
2 ' STATE :',/,' INPUT NV(',I1,') AS A NEGATIVE,WITH NO POTENTIAL '
3 ,'POINTS FOR THAT CHANNEL.',/,' ABORT.')
230 FORMAT(/,' USING THE K(E)= OPTION REQUIRES THAT THE',
1/, ' EXTERNAL ROTATIONS BE DECLARED AS INACTIVE, (EXTERNAL ',
2/, 'B VALUES SHOULD BE INPUT AS NEGATIVES). ABORT')
235 FORMAT(' EXTERNAL SYMMETRY NUMBERS',T35,F12.3,T50,F12.3,T65,
1 F12.3,F15.3)
240 FORMAT(' WARNING: NO 1-DIMENSIONAL ROTORS FOR REACTANT.',
1 /,' INCONSISTENT WITH INPUTTING EXTERNAL B VALUES AS NEGATIVE.',
2 /,' UNLESS THE MOLECULE IS LINEAR, THERE MUST BE AT LEAST',
3 /,' AN ACTIVE 1-D EXTERNAL ROTOR. TO SPECIFY A LINEAR MOLECULE,',
4 /,' PLACE A "-" SIGN IN FRONT OF NINTR.',/,' PROGRAM ABORTED.')
245 FORMAT(//,' INCORRECT NUMBER OF FREQUENCIES OR INTERNAL ROTORS',
1 ' FOR CHANNEL #',I3,/,'. ABORT.')
250 FORMAT(/,20X,'COLLISION DIAMETER (ANGSTROM)',T50,F10.2,/,
120X,'WEIGHT OF REACTANT (A.M.U.)',T50,F10.2,/,20X,
2'WEIGHT OF BATH GAS (A.M.U.)',T50,F10.2)
255 FORMAT(20X,'HARD SPHERE COLLISION FREQUENCY USED')
260 FORMAT(20X,'LENNARD JONES COLLISION FREQUENCY USED',/,
1 20X,'WELL DEPTH =',T50,F10.1,' K')
265 FORMAT(//' NO. OF INTERNAL ROTATIONAL DEGREES OF FREEDOM FOR'
1 ,' THE LOOSE',/,' TRANSITION STATE(S) OF CHANNEL ',I4,
2 ' IS LESS THAN 4.',/,' THERE MUST BE AT LEAST TWO ',
3 '2-DIMENSIONAL INTERNAL ROTORS.',/,' PROGRAM ABORTED.')
270 FORMAT(' NO 1-DIMENSIONAL ROTORS FOR CHANNEL #',I3,
1 /,' UNLESS THE TRANSITION STATE IS LINEAR,',/,
2 ' THERE MUST BE AT LEAST AN ACTIVE EXTERNAL 1-D ROTOR.',
3 /,' TO SPECIFY A LINEAR TRANSITION STATE, PLACE A "-" SIGN',
4 /,' IN FRONT OF N(IN,1). PROGRAM ABORTED.')
275 FORMAT(//,' ION/DIPOLE INTERACTION WITH DIPOLE MOMENT =',
1 F8.2,' DEBYE.')
280 FORMAT(10X,'(KJ/MOL)',T50,F12.3,T65,F12.3,F15.3)
285 FORMAT(' OVERALL ROTATION: B (CM-1)',T35,1P,E12.2,T50,E12.2,
1 T65,E12.2,E15.2)
290 FORMAT(' CORRESPONDING MOMENTS OF INERTIA (AMU A**2)',/,T35,F12.3,
1 T50,F12.3,T65,F12.3,F15.3)
295 FORMAT(' DIMENSIONS OF ADIABATIC ROTATIONS',T45,I2,T50,I12,
1 T65,I12,I15)
300 FORMAT(' EXTERNAL (2-DIMENSIONAL) ROTATION TREATED AS ADIABATIC.',
1 /,' 1-DIMENSIONAL EXTERNAL ROTOR ABOUT UNIQUE AXIS TREATED AS ',
2 /,' FULLY ACTIVE.')
305 FORMAT(/,' NO. OF FREQUENCIES',T35,I12,T50,I12,T65,I12,I15)
310 FORMAT(' FREQUENCIES AND DEGENERACIES')
315 FORMAT(T35,2I6,T50,2I6,T65,2I6,3X,2I6)
320 FORMAT(/,' FREE INTERNAL ROTATIONS:')
325 FORMAT(T25,4(A1,' B SYMMETRY DIMEN-'))
330 FORMAT(T25,4(A1,'(CM-1) NO. SION '))
335 FORMAT(23X,2(1P,E8.2,0P,F7.3,I5,2X),2X,2(1P,E8.2,0P,F7.3,I5,2X))
340 FORMAT(/,' HINDERED DIPOLE ROTATIONS (COMPLEX) :')
345 FORMAT(/,4(A1,' B SYMM.NO. DIM. V0',4X))
350 FORMAT(4(A1,'(CM-1)',13X,'(CM-1)'),2X)
355 FORMAT(3(2F7.3,I5,2X,F5.0,2X))
360 FORMAT(/,' CAPTURE RATE COEFFICIENT FOR ION/BATH GAS COLLISIONS:',
1 //,20X,'K(CAPTURE) = ',1P,D11.4,' CM3 S-1',/)
365 FORMAT(' CRITICAL ENERGIES MUST BE IN INCREASING ORDER,ABORT')
370 FORMAT(' IF USING J-AVERAGING, TEMPERATURES MUST BE IN',/,
1 ' DESCENDING ORDER (2000,1000,...); ABORT')
375 FORMAT(' ARRAY EXCEEDED, REDUCE NN OR INCREASE INC')
380 FORMAT(/,10X,'BETA VALUE OBTAINED FROM MORSE FITTING:',/)
385 FORMAT(/,' MINIM DID NOT CONVERGE. ABORT.')
390 FORMAT(10X,'FOR CHANNEL ',I3,' BETA = ',F8.4//)
395 FORMAT(10X,' FOR HINDERED ROTOR TREATMENT:',/
1,10X,'FOR CHANNEL ',I3,' V0 = ',F10.4,' CM-1'//)
400 FORMAT(1H1,1X,20A4,/,10X,'TEMPERATURE',T40,F6.1,' KELVIN')
405 FORMAT(33X,'CHANNEL(S)',I2,T56,I2,I8)
410 FORMAT(' CRITICAL ENERGIES'
2,' (KCAL/MOL)',T36,F8.2,T51,2F8.2)
415 FORMAT(' HIGH PRES. ACTVN. ENERGY (KCAL/MOL)'
4,T38,F6.2,T51,2F8.2)
420 FORMAT(10X, ' (KJ/MOL)',
5 T36,F8.2,T51,2F8.2)
425 FORMAT(' HIGH PRESSURE LOG A',T36,F8.2,T51,2F8.2)
430 FORMAT(' HIGH PRESSURE RATE COEFFICIENT (FROM NUMERICAL INTEG',
1 'RATION)',/,25X,'(S-1)',T34,1P,E10.2,T50,
1 2E10.2)
435 FORMAT(/,' HIGH PRESSURE RATE COEFFICIENT USING /I FACTOR ',
1 '(CALCULATED NUMERICALLY)',/,26X,'(S-1)',1X,'=',1P,T34,E10.2,
2 ' (CHANNEL 1)',/,T34,E10.2,' (CHANNEL 2)',/,T34,E10.2,
3 ' (CHANNEL 3)')
440 FORMAT(/,' HIGH PRESSURE /I FACTOR =',
1 3(F8.2,3X),/)
445 FORMAT(/,' WARNING : THE TWO HIGH PRESSURE RATES FOR '
1 ,'CHANNEL ',I1,' DIFFER',/,' BY MORE THAN 20%',/,
2 ' MOVE THIS TRANSITION STATE TO OBTAIN A BETTER AGREEMENT.')
450 FORMAT(/,' WARNING: DIFFERENCE BETWEEN EXACT AND NUMERICALLY',
1 /,' INTEGRATED LOG A GT 0.05; INCREASE NUMBER OF TERMS IN',/,
2 ' INTEGRATION AND/OR DECREASE GRAIN SIZE.')
455 FORMAT(/,' HIGH PRESSURE RATE COEFFICIENT (FROM THERMODYNAMICS)',
1 /,25X,'(S-1)',T34,1P,E10.2,T50,2E10.2)
460 FORMAT(' NUMERICAL AND EXACT RATES DIFFER;',
1 /,' THESE RATES AND RATES FROM MASTER WILL BE',/,
2 ' CORRECTED BY PARTITION FUNCTION RATIO (AT THIS TEMPERATURE)'
3 ,' =',F10.3)
465 FORMAT(' NUMERICALLY INTEGRATED HIGH PRESSURE RATES WITH',
1 ' THIS CORRECTION ARE',/
1 ,T20,'(S-1)',1P,T34,E10.2,T50,2E10.2)
470 FORMAT(/,' RATIO OF NUMERICAL TO EXACT PARTITION',
1 ' FUNCTION =',F6.3,/,' CORRECTED HIGH PRESSURE RATE:',
2 1P,3(E10.2,2X))
475 FORMAT(/,10X,'THERMODYNAMICS',/,30X,
1 'MOLECULE COMPLEX(ES):',I1,2I10)
C480 FORMAT(' SPEC. HEAT (KJ/MOL K)',T30,F10.2,5X,3F10.2)
485 FORMAT(' VIBRATION ENTROPY (J/MOL K)',T30,F10.2,5X,3F10.2)
490 FORMAT(' ROTATION ENTROPY (J/MOL K)',T30,F10.2,5X,3F10.2)
495 FORMAT(' (INTERNAL AND EXTERNAL ROTATIONS INCLUDED TOGETHER)')
500 FORMAT(' TRANSLATIONAL ENTROPY (J/MOL K)',T33,F7.2,5X,3F10.2)
505 FORMAT(' TOTAL ENTROPY (J/MOL K)',T30,F10.2,5X,3F10.2)
510 FORMAT(/,' LOG A FROM THERMODYNAMICS',T45,3F10.3)
511 FORMAT(/,' APPROX. LOG A FROM THERMODYNAMICS',T45,3F10.3)
515 FORMAT(' VIB. PARTITION FUNCTION',T30,1P,E10.2,5X,3E10.2)
520 FORMAT(' ROTNL PARTITION FUNCTION',T30,1P,E10.2,5X,3E10.2)
521 FORMAT(' APPROX. HIGH PRES. ACTVN. ENERGY (KCAL/MOL)',/,
1 T38,F6.2,T51,2F8.2)
525 FORMAT(//,20X,'STRONG COLLISION CALCULATION:',/,
1 5X,'PRESSURE',T24,'OMEGA',T35,'K1',T43,'K1/K1 INF',
2 T58,'K2',T65,'K2/K2 INF',/,6X,'(TORR)',T24,'(S-1)',T33,'(S-1)')
530 FORMAT(1P,E12.2,T20,E10.2,T31,3(E10.2,1X,E10.2,1X))
535 FORMAT(//,' STRONG COLLISION LOW-PRESSURE RATE COEFFICIENT',/,
1 ' (CM3 MOL-1 S-1)',T34,1P,E10.2,T50,2E10.2)
540 FORMAT(//,' TOTAL STRONG COLLISION (E,J) LOW PRESSURE RATE ',
1 'COEFFICIENT',/,' (CM3 MOL-1 S-1)',T34,1P,E10.2)
545 FORMAT(//,' STRONG COLLISION (E,J=0) LOW PRESSURE RATE',
1 ' COEFFICIENT',/,' (CM3 MOL-1 S-1)',T34,1P,E10.2)
550 FORMAT(1X,20A4,/,1X,I5,8X,2I5)
555 FORMAT(1X,7E10.3)
560 FORMAT(1X,F12.3,/,1X,' 1',/,1X,' 500.',/,1X,' 0 1',/,1X,I4)
565 FORMAT(1X,' 0.5')
570 FORMAT(1X,I4)
575 FORMAT(1X,I4,/,1X,15(F8.2,2X))
580 FORMAT(1X,4F12.3)
585 FORMAT(1X,2I2)
590 FORMAT(1X,F8.2,2X,F8.3)
595 FORMAT(1X,6E11.4)
600 FORMAT(1X,5E10.3)
605 FORMAT(1X,E10.2)
610 FORMAT(//,' TOTAL NUMBER OF INCREMENTS (NN) IS NOT LARGE',
1 ' ENOUGH TO ALLOW CONVERGENCE',/,' OF HIGH PRESSURE RATE',/,
2 ' INCREASE NN BY 50. (PROGRAM ABORTED).')
615 FORMAT(/,' EQUILIBRIUM CONSTANT FOR REACTANT CHANNEL :',10X,
1 1P,E10.2,3X,'(CM3 MOLEC-1)')
620 FORMAT(/,' CAPTURE RATE COEFFICIENT :',10X,1P,E10.2,3X,
1 ' CM3 MOLEC-1 S-1')
625 FORMAT(/,' LOW DENSITY LIMITING BRANCHING RATIOS :',/)
630 FORMAT(' CHANNEL ',I4,10X,' BRANCHING RATIO =',1P,E10.2)
635 FORMAT(/,' RADIATIVE STABILISATION BRANCHING RATIO :',10X,
1 1P,E10.2)
640 FORMAT(/,' COLLISION COMPLEX LIFETIME IN THE LOW DENSITY LIMIT',
1 ' :',5X,1P,E10.2,' MICROSEC')
645 FORMAT(/,' UNIMOLECULAR RATE COEFFICIENT(S) FOR DISSOCIATION OF',
1 ' COLLISION COMPLEX IN THE LOW DENSITY LIMIT :',/,' (S-1)',
2 T20,1P,3(E10.2,5X))
650 FORMAT(/,' UNIMOLECULAR RADIATIVE STABILISATION RATE',
1 ' COEFFICIENT :',/,10X,'(S-1)',T30,1P,E10.2)
660 FORMAT(/,' CONTRIBUTIONS FROM TWO DIMENSIONAL EXTERNAL'
1 ,' ROTATION :')
665 FORMAT(/,' FOR RECOMBINATION, TOTAL ENTROPY OF REACTANTS ',
1 '(J/MOL K) :',10X,1P,E10.2)
670 FORMAT(6X,'(TORR)',T24,'(S-1)',T33,'(S-1)',T57,'(S-1)',T77,
1 '(S-1)')
675 FORMAT(6X,'(TORR)',T24,'(S-1)',T32,'(CM3 S-1)',T56,'(CM3 S-1)'
1 ,T75,'(CM3 S-1)')
680 FORMAT(/,' STRONG COLLISION LOW PRESSURE ASSOCIATION RATE',
1 ' COEFFICIENT :',/,10X,'(CM6 MOLEC-2 S-1)',10X,1P,E10.2)
685 FORMAT(1X,F12.3,/,' 1',/,' 1',/,' 500.',/,' 0 1')
690 FORMAT(1X,F5.0)
695 FORMAT(1P,E11.4,' 0 0 0')
C
C OPEN(5,FILE='rrkm.dat',STATUS='OLD')
C OPEN(6,FILE='rrkm.out',STATUS='UNKNOWN')
C OPEN(10,FILE='mas.dat',STATUS='UNKNOWN')
READ(5,100) TITLE
READ(5,*)NN,INC,NP,NT,NCHAN,JAV
WRITE(6,105)TITLE
WRITE(6,115)NN,INC
WRITE(6,120)NP,NT,NCHAN
NCHP=MAX0(NCHAN,2)
JAVX=JAV.NE.0
C
C THE SINUSOIDALLY-HINDERED AND ION-DIPOLE OPTIONS ARE ONLY AVAILABLE
C IF THE J-AVERAGING REGIME IS USED.
C TO ACCESS THESE OPTIONS IT IS NECESSARY TO SPECIFY
C JAV GREATER THAN 1.
C
IF(JAV.GT.1)READ(5,*)(IHIND(IN),ION(IN),IN=1,NCHAN)
DO 900 IN=1,NCHAN
WKJ(IN)=.FALSE.
IHINDX(IN)=.FALSE.
IONX(IN)=.FALSE.
STERIC(IN)=.FALSE.
IF(.NOT.JAVX) GO TO 900
IHINDX(IN)=IHIND(IN).NE.0
IONX(IN)=ION(IN).NE.0
C
C INFORMATION READ IN FOR ION-DIPOLE AND NEUTRAL-NEUTRAL SINUSOIDALLY
C HINDERED ROTORS:
C
IF(.NOT.IHINDX(IN)) GO TO 900
IF(.NOT.IONX(IN))WRITE(6,125)IN
IF(IONX(IN))WRITE(6,130)IN
READ(5,*)ROTINC(IN),ERR(IN)
WRITE(6,135)ROTINC(IN),ERR(IN)
C
C IF IHIND(1) IS READ IN AS A NEGATIVE NUMBER THEN STERIC HINDRANCE
C IS USED EXPLICITLY IN THE CALCULATION OF THE SINUSOIDALY HINDERED
C PARTITION FUNCTION.
C THETA1 AND THETA2 ARE THE HINDRANCES IN THE THETA EULER ANGLE FOR
C ONE MOIETY AND - IF USED - THE OTHER.
C
STERIC(IN)=IHIND(IN).LT.0
THETA1(IN)=180.0D0
THETA2(IN)=180.0D0
ISYM1(IN)=1
ISYM2(IN)=1
IHIND(IN)=ABS(IHIND(IN))
WRITE(6,140)IN,IHIND(IN)
IF(.NOT.STERIC(IN))GO TO 800
WRITE(6,145)IN
C
C ISYM1&2 ARE THE SYMMETRY NUMBERS OF THE SIN HINDERED ROTORS
C EG PLANAR CH3 IS SYMMETRIC AND HAS SYM NO.=2 FOR THE 2-DIM
C ROTATION.
C
READ(5,*)THETA1(IN),ISYM1(IN)
WRITE(6,150)ISYM1(IN),THETA1(IN)
IF(THETA1(IN).EQ.180.0D0)ISYM1(IN)=1
IF(IHIND(IN).EQ.1)GO TO 800
C
C IF IHIND(IN)=-2 THEN WE HAVE 2 HINDRANCE ANGLES FOR THIS CHANNEL.
C
READ(5,*)THETA2(IN),ISYM2(IN)
WRITE(6,155)ISYM2(IN),THETA2(IN)
IF(THETA2(IN).EQ.180.0D0)ISYM2(IN)=1
800 CONTINUE
READ(5,*)GAMMA(IN)
C
C IF GAMMA > 600: STRONG j RELAXATION EVALUATION OF j-AVERAGED K(e)
C IF GAMMA < 600: WEAK COLLISIONAL EVALUATION OF j-AVERAGED K(e)
C
WKJ(IN)=GAMMA(IN).LE.600.0D0
IF(.NOT.WKJ(IN))WRITE(6,175)
IF(WKJ(IN))WRITE(6,165)GAMMA(IN)
900 CONTINUE
IF(.NOT.JAVX)GO TO 1030
C
C OBVIOUSLY IF WE START OF WITH AN ION-DIPOLE SYSTEM WE KEEP IT
C REGARDLESS OF WHAT THE TRANSITION STATE LOOKS LIKE, SO IT IS
C IMPLICITLY ASSUMED THAT IF IONX(1) IS TRUE THEN IONX(2) AND
C IONX(3) MUST ALSO BE TRUE:
C
IONX1=IONX(1)
DO 1000 IN=1,NCHAN
IF((IONX1).AND.(IONX(IN)))GO TO 1000
IF((.NOT.IONX1).AND.(.NOT.IONX(IN)))GO TO 1000
WRITE(6,170)
STOP
1000 CONTINUE
IF(IONX1) GO TO 1999
READ(5,*)(RCPL(IN,1),REQ(IN),IN=1,NCHAN)
DO 1010 I10=1,3
JCBS(I10)=.TRUE.
1010 CONTINUE
DO 1020 IN=1,NCHAN
READ(5,*) NV(IN)
JCBS(IN)=NV(IN).LE.0
IF(NV(IN).GT.0)READ(5,*)(RVCH(IN,II),VCH(IN,II),
1 II=1,NV(IN))
1020 CONTINUE
1030 CONTINUE
READ(5,*) (JF(IN,1),IN=1,NCHAN),NF
READ(5,*) (EOK(IN),IN=1,NCHAN)
IF(.NOT.JAVX) GO TO 1060
IF(JAV.EQ.1)WRITE(6,160)
WRITE(6,180) (RCPL(IN,1),REQ(IN),IN,IN=1,NCHAN)
DO 1050 IN=1,NCHAN
WRITE(6,190)IN
I50 = NV(IN)
IF(I50.LE.0)GO TO 1040
WRITE(6,195) I50
WRITE(6,200) (RVCH(IN,J),VCH(IN,J),J=1,I50)
GO TO 1050
1040 WRITE(6,205)RCPL(IN,1),EOK(IN)
1050 CONTINUE
GO TO 1070
1060 WRITE(6,210)
1070 CONTINUE
WRITE(6,215) (IN,IN=1,NCHAN)
WRITE(6,220) (EOK(IN),IN=1,NCHAN)
DO 1080 IN=1,NCHAN
EJ(IN)=4.184D0*EOK(IN)
1080 CONTINUE
WRITE(6,280) (EJ(IN),IN=1,NCHAN)
READ(5,*) (SRC(IN,1),IN=1,NCHAN),SRM
READ(5,*) (BCMPLX(IN,1),IN=1,NCHAN),BMOLEC
DO 1090 IN=1,NCHAN
IF((BCMPLX(IN,1)/BMOLEC).GT.0.9D0.AND.(NV(IN).GT.0))GO TO 1100
1090 CONTINUE
GO TO 1110
1100 WRITE(6,225)IN,IN
STOP
1110 CONTINUE
IF(BMOLEC.LT.0.0D0) IEXMD=2
TEST1D=BMOLEC.LT.0.0D0
BMOLEC=DABS(BMOLEC)
IF(.NOT.JAVX)GO TO 1120
IF(TEST1D)GO TO 1120
WRITE(6,230)
STOP
1120 CONTINUE
READ(5,*) (N(IN,1),IN=1,NCHAN),NINTR
C
C LINEAR MOLECULE OR COMPLEX SPECIFIED BY NEGATIVE
C N(IN,1) AND/OR NINTR.
C
LINM=NINTR.LE.0
NINTR=MAX0(NINTR,0)
DO 1140 IN=1,NCHAN
IF(BCMPLX(IN,1).GT.0.0D0) GO TO 1130
TEST1D=.TRUE.
IEXRTD(IN)=2
BCMPLX(IN,1)=DABS(BCMPLX(IN,1))
1130 CONTINUE
LINC(IN)=N(IN,1).LE.0
N(IN,1)=MAX0(N(IN,1),0)
J=N(IN,1)
IF(.NOT.LINC(IN))READ(5,*) (BVEC(IN,1,II),SIGVC(IN,1,II),
* IRTDMC(IN,1,II),II=1,J)
1140 CONTINUE
IF(.NOT.LINM) READ(5,*) (BVECM(I),SIGVCM(I),IRTDMM(I),I=1,NINTR)
READ(5,*) SGMA,WT1,WT2,EPS
C
C IF JAVX AND WT1.LT.0. THEN THE CENTRES OF MASS OF THE FRAGMENTS ARE
C NOT BOTH THE PIVOT ATOMS OF THE BREAKING BOND. IN THIS CASE AN ARRAY
C OF BREAKING BOND LENGTHS AND ROT. CONSTANTS IS READ IN FOR EACH
C CHANNEL AND A STRAIGHT LINE FIT IS PERFORMED.
C
C MAX NO. POINTS = 30
C
BINPUT=(JAVX.AND.(WT1.LT.0.0D0))
WT1=ABS(WT1)
IF(.NOT.BINPUT)GO TO 1146
DO 1145 IN=1,NCHAN
READ(5,*) NB(IN)
READ(5,*)(RR(IN,II),BB(IN,II),II=1,NB(IN))
NB1=NB(IN)
DO 1144 I1145=1,NB1
RR1(I1145)=RR(IN,I1145)
BB1(I1145)=BB(IN,I1145)
1144 CONTINUE
CALL BESTFIT(NB1,RR1,BB1,A1,A2)
AA1(IN)=A1
AA2(IN)=A2
1145 CONTINUE
1146 CONTINUE
DO 1150 IN=1,NCHAN
JF1=JF(IN,1)
READ(5,*) (NC(IN,1,I),JC(IN,1,I),I=1,JF1)
1150 CONTINUE
WRITE(6,235) SRM, (SRC(IN,1),IN=1,NCHAN)
READ(5,*) (NM(I),JM(I),I=1,NF)
C
C CHECK ON NUMBER OF DEGREES OF FREEDOM OF REACTANT
C AND (IF ACTIVE EXTERNAL ROTOR) THAT THERE IS A
C 1-DIMENSIONAL ACTIVE ROTOR.
C
CHECK=.FALSE.
NCHKM=0
NCHKM=IEXMD
IF(LINM)GO TO 1210
DO 1200 J1200=1,NINTR
CHECK=CHECK.OR.(IRTDMM(J1200).EQ.1)
NCHKM=NCHKM+IRTDMM(J1200)
1200 CONTINUE
1210 CONTINUE
IF(.NOT.TEST1D)GO TO 1220
IF(CHECK.OR.LINM)GO TO 1220
WRITE(6,240)
STOP
1220 CONTINUE
DO 1230 J1230=1,NF
NCHKM=NCHKM+JM(J1230)
1230 CONTINUE
DO 1300 IN=1,NCHAN
C
C CHECK ON NUMBER OF DEGREES OF FREEDOM OF COMPLEX(ES)
C AND (IF ACTIVE EXTERNAL ROTOR) THT THERE IS A
C 1-DIMENSIONAL ACTIVE ROTOR
C
CHCKRC(IN)=.FALSE.
NCHK(IN)=0
NCHK(IN)=IEXRTD(IN)
IF(LINC(IN))GO TO 1250
DO 1240 J1240=1,N(IN,1)
CHCKRC(IN)=CHCKRC(IN).OR.(IRTDMC(IN,1,J1240).EQ.1)
NCHK(IN)=NCHK(IN)+IRTDMC(IN,1,J1240)
1240 CONTINUE
1250 CONTINUE
IF(.NOT.TEST1D)GO TO 1260
IF(CHCKRC(IN).OR.LINC(IN))GO TO 1260
WRITE(6,270)IN
STOP
1260 CONTINUE
NL=JF(IN,1)
DO 1270 J1270=1,NL
NCHK(IN)=NCHK(IN)+JC(IN,1,J1270)
1270 CONTINUE
IF(NCHK(IN).EQ.NCHKM-1)GO TO 1280
WRITE(6,245)IN
STOP
1280 CONTINUE
C
C END OF CHECK
C
1300 CONTINUE
WRITE(6,250) SGMA,WT1,WT2
IF(EPS.LE.0.0D0) WRITE(6,255)
IF(EPS.GT.0.0D0) WRITE(6,260) EPS
go to 2535
1999 CONTINUE
C
C INPUT STATEMENTS FOR THE ION-DIPOLE CASE (PARAMETERS ARE INPUT IN A
C SLIGHTLY DIFFERENT WAY)
C
READ(5,*)BMOLEC
TEST1D=BMOLEC.LT.0.0D0
BMOLEC=DABS(BMOLEC)
IEXMD=2
READ(5,*)NINTR
C
C IF NINTR IS READ IN AS A NEGATIVE, A LINEAR MOLECULE IS SPECIFIED
C
LINM=NINTR.LE.0
NINTR=MAX0(NINTR,0)
IF(LINM)GO TO 2000
READ(5,*)(BVECM(I),SIGVCM(I),IRTDMM(I),I=1,NINTR)
2000 CONTINUE
READ(5,*)NF
IF(NF.LE.0)GO TO 2010
READ(5,*)(NM(I),JM(I),I=1,NF)
2010 CONTINUE
C
C CHECK ON NUMBER OF DEGREES OF FREEDOM OF REACTANT
C AND (IF ACTIVE EXTERNAL ROTOR) THAT THERE IS A
C 1-DIMENSIONAL ACTIVE ROTOR.
C
CHECK=.FALSE.
NCHKM=0
NCHKM=IEXMD
IF(LINM) GO TO 2030
DO 2020 J2020=1,NINTR
CHECK=CHECK.OR.(IRTDMM(J2020).EQ.1)
NCHKM=NCHKM+IRTDMM(J2020)
2020 CONTINUE
2030 CONTINUE
IF(.NOT.TEST1D)GO TO 2040
IF(CHECK.OR.LINM)GO TO 2040
WRITE(6,240)
STOP
2040 CONTINUE
NL=NF
DO 2050 I2050=1,NF
NCHKM=NCHKM+JM(I2050)
2050 CONTINUE
READ(5,*)KCAPT,WT1
DO 2500 IN=1,NCHAN
READ(5,*)DHO(IN)
JCBS(IN)=.FALSE.
LINC(1)=.FALSE.
IEXRTD(IN)=2
READ(5,*)REQ(IN)
READ(5,*)RCPL(IN,1),BCMPLX(IN,1)
TEST1D=BCMPLX(IN,1).LT.0.0D0
BCMPLX(IN,1)=DABS(BCMPLX(IN,1))
READ(5,*)DIP(IN),POLZ(IN)
READ(5,*)WTA(IN),WTB(IN)
READ(5,*)N(IN,1)
LINC(1)=N(IN,1).LT.0
c corrected 16 May 1994
N(IN,1)=MAX0(N(IN,1),0)
IF(DIP(IN).GT.0.05D0)N(IN,1)=N(IN,1)-1
IF(DIP(IN).GT.0.05D0)NH(IN,1)=1
IF(.NOT.LINC(1))READ(5,*)(BVEC(IN,1,II),SIGVC(IN,1,II),
1 IRTDMC(IN,1,II),II=1,N(IN,1)+NH(IN,1))
READ(5,*)JF(IN,1)
IF(JF(IN,1).GT.0)READ(5,*)(NC(IN,1,II),JC(IN,1,II),
1 II=1,JF(IN,1))
POTD=69.1235D0*DIP(IN)/(RCPL(IN,1)**2)
V0(IN,1)=POTD*WKA*2.0D+0
C
C CHECK ON NUMBER OF DEGREES OF FREEDOM OF COMPLEX(ES)
C AND (IF ACTIVE EXTERNAL ROTOR) THAT THERE IS A
C 1-DIMENSIONAL ACTIVE ROTOR.
C
CHCKRC(IN)=.FALSE.
NCHK(IN)=0
NRT=N(IN,1)+NH(IN,1)
NCHK(IN)=IEXRTD(IN)
IF(NRT.EQ.0) GO TO 2180
DO 2170 J2170=1,NRT
CHCKRC(IN)=CHCKRC(IN).OR.(IRTDMC(IN,1,J2170).EQ.1)
NCHK(IN)=NCHK(IN)+IRTDMC(IN,1,J2170)
2170 CONTINUE
C END-CHECK COMMENTED OUT IE WE DO NOT NEED 2 2-DIM INT ROTORS
C EG ONE FRAGMENT IS AN ATOM (MJTJ 16.5.94)
C IF(NCHK(IN).GE.6)GO TO 2180
C WRITE(6,265)IN
C STOP
2180 CONTINUE
IF(.NOT.TEST1D)GO TO 2190
IF(CHCKRC(IN).OR.LINC(1))GO TO 2190
WRITE(6,270)IN
STOP
2190 CONTINUE
NL=JF(IN,1)
DO 2200 I2200=1,NL
NCHK(IN)=NCHK(IN)+JC(IN,1,I2200)
2200 CONTINUE
IF(NCHK(IN).EQ.NCHKM-1) GO TO 2210
WRITE(6,245)IN
STOP
2210 CONTINUE
C
C END OF CHECK.
C
EOK(IN)=DHO(IN)
2500 CONTINUE
READ(5,*)JRECOM
IF(JRECOM.EQ.1)READ(5,*)INCHAN
IF(JRECOM.EQ.1)READ(5,*)RADST
WRITE(6,185) (RCPL(IN,1),IN,IN=1,NCHAN)
DO 2520 IN=1,NCHAN
WRITE(6,190)IN
WRITE(6,275)DIP(IN)
2520 CONTINUE
WRITE(6,215) (IN,IN=1,NCHAN)
C CORRECTED FORMAT STATEMENT 16.5.1994 (MJTJ)
WRITE(6,220) (DHO(IN),IN=1,NCHAN)
DO 2530 IN=1,NCHAN
EJ(IN)=4.184D0*DHO(IN)
2530 CONTINUE
C
C MERGE NEUTRAL-NEUTRAL AND ION-DIPOLE BITS
C
2535 CONTINUE
WRITE(6,280) (EJ(IN),IN=1,NCHAN)
WRITE(6,285) BMOLEC,(BCMPLX(IN,1),IN=1,NCHAN)
DO 2540 IN=1,NCHAN
A(IN)=16.8477D0/BCMPLX(IN,1)
2540 CONTINUE
AM=16.8477D0/BMOLEC
WRITE(6,290) AM,(A(IN),IN=1,NCHAN)
WRITE(6,295) IEXMD,(IEXRTD(IN),IN=1,NCHAN)
IF(TEST1D)WRITE(6,300)
WRITE(6,305) NF,(JF(IN,1),IN=1,NCHAN)
WRITE(6,310)
NL=NF
DO 2550 IN=1,NCHAN
NL=MAX0(NL,JF(IN,1))
2550 CONTINUE
DO 2560 I2560=1,NL
WRITE(6,315) NM(I2560),JM(I2560),(NC(IN,1,I2560),
1 JC(IN,1,I2560),IN=1,NCHAN)
2560 CONTINUE
NQ=NINTR
DO 2570 IN=1,NCHAN
NQ=MAX0(NQ,N(IN,1))
2570 CONTINUE
IF(NQ.LE.0) GO TO 2700
DO 2580 IN=1,NCHAN
IF(NH(IN,1).EQ.0)GO TO 2580
NFR=N(IN,1)
IF(NFR.EQ.NQ)GO TO 2580
BVEC(IN,1,NQ+1)=BVEC(IN,1,NFR+1)
BVEC(IN,1,NFR+1)=0.0D0
SIGVC(IN,1,NQ+1)=SIGVC(IN,1,NFR+1)
SIGVC(IN,1,NFR+1)=0.0D0
IRTDMC(IN,1,NQ+1)=IRTDMC(IN,1,NFR+1)
IRTDMC(IN,1,NFR+1)=0
2580 CONTINUE
WRITE(6,320)
I=NCHAN+1
WRITE(6,325) (IBL,J=1,I)
WRITE(6,330) (IBL,J=1,I)
DO 2590 I590=1,NQ
WRITE(6,335) BVECM(I590),SIGVCM(I590),IRTDMM(I590),
1 (BVEC(IN,1,I590),SIGVC(IN,1,I590),IRTDMC(IN,1,I590),IN=1,NCHAN)
2590 CONTINUE
IF(.NOT.IONX1)GO TO 2700
WRITE(6,340)
I=NCHAN
WRITE(6,345)(IBL,J=1,I)
WRITE(6,350)(IBL,J=1,I)
WRITE(6,355)(BVEC(IN,1,NQ+1),SIGVC(IN,1,NQ+1),IRTDMC(IN,1,NQ+1),
1 V0(IN,1),IN=1,NCHAN)
DO 2600 IN=1,NCHAN
IF(N(IN,1).EQ.NQ)GO TO 2600
IF(NH(IN,1).EQ.0)GO TO 2600
NP1=N(IN,1)+1
BVEC(IN,1,NP1)=BVEC(IN,1,NQ+1)
SIGVC(IN,1,NP1)=SIGVC(IN,1,NQ+1)
IRTDMC(IN,1,NP1)=IRTDMC(IN,1,NQ+1)
BVEC(IN,1,NQ+1)=0.0D0
SIGVC(IN,1,NQ+1)=0.0D0
IRTDMC(IN,1,NQ+1)=0.0D0
2600 CONTINUE
2700 CONTINUE
IF((JAVX).AND.(IONX1))WRITE(6,360)KCAPT
2800 CONTINUE
READ(5,*) (PRESS(I),I=1,NP)
READ(5,*) (TEMP(I),I=1,NT)
DO 3000 IN=1,NCHAN
IF(IN.LE.1) GO TO 3005
C
C CHECK THAT CRITICAL ENERGIES ARE IN INCREASING ORDER
C
IF(EOK(IN).GT.EOK(IN-1)) GO TO 3005
WRITE(6,365)
STOP
3005 CONTINUE
3000 CONTINUE
C
C CHECK TEMPERATURES IN DESCENDING ORDER
C
IF(.NOT.JAVX) GO TO 3020
CHECK = .TRUE.
IF(NT.EQ.1) GO TO 3020
DO 3010 I3010 = 2,NT
CHECK = CHECK.AND.(TEMP(I3010).LT.TEMP(I3010-1))
3010 CONTINUE
IF(CHECK) GO TO 3020
WRITE(6,370)
STOP
3020 CONTINUE
C
C COMMENCE GENERATION OF MICROSCOPIC RATES, K(E).
C
C SET UP GAMMA FUNCTION; GAMON2(I)=GAMMA(I/2), I=1,2,3 ...
C
GAMON2(1)=1.772454D0
GAMON2(2)=1.0D0
DO 4000 I4000=2,20
NG2=2*I4000
GAMON2(NG2)=DFLT(I4000-1)*GAMON2(NG2-2)
N2L1=NG2-1
GAMON2(N2L1)=(DFLT(I4000)-1.5D0)*GAMON2(N2L1-2)
4000 CONTINUE
DELT=DFLT(INC)
DO 4999 I4999=1,NT
T=TEMP(I4999)
C
C HERE R IS GIVEN IN UNITS OF CM-1*K-1 IE WE ARE ALWAYS
C LOOKING AT ENERGY IN UNITS OF CM-1.
C
RT=R*T
IF(WKJ(1))DELTA=GAMMA(1)*RT/(GAMMA(1)+RT)
CORRAT=1.0D0
CORRAT1=1.0D0
IF((.NOT.JAVX).AND.(I4999.GT.1))GO TO 4030
NREACT=INT((EOK(1)*WKA/DELT))
NI=NREACT+NN
C
C TEST FOR EXCEEDING ARRAY SIZE
C
IF(NI.LE.(NMAX8-1)) GO TO 4010
WRITE(6,375)
STOP
4010 CONTINUE
C
C FIND DENSITY OF STATES OF REACTANT.
C
CALL BSWINE(NM,JM,RHOMOL,NI,NF,INC,BVECM,SIGVCM,IRTDMM,NINTR,
1 .FALSE.)
C
C RESET DEGENERACIES AFTER RE-ARRANGEMENT BY DIRECT COUNT.
C
DO 4020 I4020=1,NF
JM(I4020)=1
4020 CONTINUE
C
C THIS COMPUTES DENSITIES OF STATES OF ALL LEVELS OF REACTANT.
C
CP=DFLT(IEXMD)/2.0D0
4030 DO 4990 J=1,NP
C
C SEPARATE BITS FOR ION-DIPOLE VS NEUTRAL-NEUTRAL TRANSITION STATES
C
IF(IONX1) GO TO 5000
NJC=NINT(EOK(1)*WKA/(2.0D0*DFLT(INC)))
C
C CALCULATE OMEGA FOR THIS PRESSURE
C
PRS=PRESS(J)
OMI=DSQRT(T*WT1*WT2/(WT1+WT2))
OMEGA=(4.41313D+07)*PRS*(SGMA**2)/OMI
IF(EPS.GT.0.0D0)OM22=0.636D0+0.567D0*LOG10(T/EPS)
IF(EPS.GT.0.0D0)OMEGA=OMEGA/OM22
IF(JAVX)GO TO 4040
IF(J.NE.1)GO TO 4470
IF((J.EQ.1).AND.(I4999.GT.1))GO TO 4250
4040 CONTINUE
IF(J.EQ.1.AND.JAVX)WRITE(6,380)
DO 4160 IN=1,NCHAN
ROTINC1=ROTINC(IN)
ERR1=ERR(IN)
WKJ1=WKJ(IN)
IHINDX1=IHINDX(IN)
STERICX=STERIC(IN)
GAMMA1=GAMMA(IN)
OMEGA1=THETA1(IN)*3.1415926D0/180.0D0
OMEGA2=THETA2(IN)*3.1415926D0/180.0D0
ISYMX1=ISYM1(IN)
ISYMX2=ISYM2(IN)
JF1=JF(IN,1)
IF(JAVX)STC(IN)=1.0D0
IF(.NOT.JAVX)STC(IN)=(SRM/SRC(IN,1))*
1 ((BMOLEC/BCMPLX(IN,1))**CP)
DO 4050 I4050=1,JF1
NC1(I4050)=NC(IN,1,I4050)
JC1(I4050)=JC(IN,1,I4050)
4050 CONTINUE
N1=N(IN,1)
N11=MAX0(N1,1)
DO 4060 I4060=1,N11
BVEC1(I4060)=BVEC(IN,1,I4060)
SIGVC1(I4060)=SIGVC(IN,1,I4060)
ISPARE(I4060)=IRTDMC(IN,1,I4060)
4060 CONTINUE
KK=INT((EOK(IN)-EOK(1))*WKA/(DFLT(INC))+0.5D0)
KKC(IN)=KK
NN2=NN-KK
EOK1=EOK(IN)
BCPLX1=BCMPLX(IN,1)
IF(.NOT.JAVX)GO TO 4085
BETAX=2.0D0
IF(JCBS(IN))GO TO 4080
DBETA=0.2D+0
ACB=0.50D-1
ICOUNT=100
CHECK=.FALSE.
D=EOK1
REQ1=REQ(IN)
NV1=NV(IN)
DO 4070 I4070=1,NV1
VCH1(I4070)=VCH(IN,I4070)
RVCH1(I4070)=RVCH(IN,I4070)
4070 CONTINUE
CALL MINIM(BETAX,ICOUNT,DBETA,ACB,CHECK)
C
C MINIM CALCULATES THE BEST FIT MORSE POTENTIAL FOR THE NEUTRAL-NEUTRAL
C CASE
C
IF(CHECK)GO TO 4080
WRITE(6,385)
STOP
4080 CONTINUE
BETA(IN)=BETAX
C
C FOR THE SINUSOIDALLY HINDERED SYSTEM, V0(IN,1)-THE MAXIMUM POTENTIAL
C INTERACTION FOR CHANNEL (IN)-IS CALCULATED FROM THE BEST FIT MORSE
C POTENTIAL FUNCTION FOR THE DATA.
C
C IE VHIND=V0*(1-COS(THETA))
C
IF(.NOT.IHINDX1)GO TO 4085
NH(IN,1)=IHIND(IN)
C
C IE THERE CAN BE ONE OR TWO 2-DIM SINUSOIDALLY HINDERED ROTORS
C SO NH=1 OR NH=2.
C
BE=BETA(IN)
V01=0.0D0
V0(IN,1)=EOK(IN)-MORSE(RCPL(IN,1))
V0(IN,1)=349.7574286D0*V0(IN,1)
C
C FACTOR 349.757... CONVERTS KCAL/MOL TO CM-1,
C IE WE NEED TO MULTIPLY BY:
C
C 4.184*1000/(H*C*6.022E+23)
C
V01=V0(IN,1)
4085 CONTINUE
C RGG correction Jan 9 1993 to eliminate printout when JAV = 0:
IF(.NOT.JCBS(IN).AND.J.EQ.1.AND.JAV.NE.0)
1 WRITE(6,390)IN,BETA(IN)
IF((IHINDX1).AND.(.NOT.JCBS(IN).AND.J.EQ.1))
* WRITE(6,395)IN,V0(IN,1)
C
C NOW CALCULATE THE RO-VIBRATIONAL DENSITY AND SUM OF STATES ELEMENTS
C FOR THE TRANSITION STATE (IE SUM=.TRUE.):
C
CALL BSWINE(NC1,JC1,AK,NN2,JF1,INC,BVEC1,SIGVC1,
1 ISPARE,N1,.TRUE.)
IF(.NOT.JAVX)GO TO 4110
IF(IN.NE.NCHAN)GO TO 4110
DO 4090 I4090=1,NREACT+10
IF(ROTH(I4090).LE.0.0D0)GO TO 4100
ROTHR(I4090)=ROTH(I4090)
4090 CONTINUE
4100 NTHR=I4090-1
4110 CONTINUE
C
C RESET DEGENERACIES AFTER RE-ARRANGEMENT BY DIRECT COUNT.
C
JF(IN,1)=JF1
DO 4120 I4120=1,JF1
NC(IN,1,I4120)=NC1(I4120)
JC(IN,1,I4120)=1
4120 CONTINUE
IF(JAVX)GO TO 4140
DO 4130 I4130=1,NN2
K(IN,I4130)=AK(I4130)
4130 CONTINUE
4140 NMAX4=MIN0(NMAX4,NN2)
IF((.NOT.JAVX).OR.(IN.NE.NCHAN))GO TO 4160
DO 4150 I4150=1,NI
DO 4145 I=1,NCHAN
K(I,I4150)=WE(I4150,I)
4145 CONTINUE
4150 CONTINUE
4160 CONTINUE
IF(.NOT.JAVX)GO TO 4180
NREACT=NJC
DO 4170 I4170=1,NCHAN
KKC(I4170)=0
4170 CONTINUE
4180 CONTINUE
NI=NREACT+NMAX4
DO 4220 IN=1,NCHAN
STEM=STC(IN)/H
ICC=0
DO 4200 I4200=1,NMAX4
IF(I4200.GT.KKC(IN))GO TO 4190
AK(I4200)=0.0D0
AK2(I4200)=0.0D0
GO TO 4200
4190 CONTINUE
ICC=ICC+1
AK(I4200)=K(IN,ICC)*STEM/(RHOMOL(I4200+NREACT))
IF(JAVX.AND.(IN.EQ.1))AK2(I4200)=XKHP(I4200,1)*STEM/
1 (RHOMOL(I4200+NREACT))
4200 CONTINUE
DO 4210 I4210=1,NMAX4
K(IN,I4210)=AK(I4210)
IF(JAVX.AND.(IN.EQ.1))XKHP(I4210,1)=AK2(I4210)
4210 CONTINUE
IF(.NOT.JAVX)GO TO 4220
RTHI(IN)=RTHI(IN)*SRM/(SRC(IN,1)*H)
RTHI2(IN)=RTHI2(IN)*SRM/(SRC(IN,1)*H)
4220 CONTINUE
C
C COMPUTATION OF DENSITY OF STATES AND MICROSCOPIC RATES NOW COMPLETE.
C
C NOTE THE DIRECT COUNT SUBROUTINE UNSCRAMBLES DEGENERATE FREQUENCIES
C TO AVOID BUNCHING CAUSED BY DIRECT COUNT ALGORITHM.
C
C COMMENCE INTEGRATION TO FIND RATE COEFFICIENTS,ETC.
C
IF(JAVX.AND.(J.NE.1))GO TO 4500
4250 DO 4260 IN=1,NCHAN
STLP(IN)=0.0D0
RATE(IN)=0.0D0
EKG(IN)=0.0D0
AKG(IN)=0.0D0
AKHP(IN)=0.0D0
4260 CONTINUE
EGI=0.0D0
GI=0.0D0
FEXP=DELT/RT
C
C THE FOLLOWING LOOP FINDS HIGH-PRESSURE PARAMETERS.
C
DO 4290 I=1,NI
G=0.0D0
IF(RHOMOL(I).LE.0.0D0)GO TO 4290
G=EXP(LOG(RHOMOL(I))-DFLT(I)*FEXP)
C
C USING TRAPEZOIDAL RULE TAKE HALF OF FIRST VALUE IN SUM
C MJTJ 7.1.93
C
GI=GI+G
IF(I.EQ.1)GI=0.5D0*GI
IF(JAVX)GO TO 4280
E=DFLT(I*INC)/WKA
EGI=EGI+G*E
IF(I.EQ.1)EGI=0.5D0*EGI
IF(I.LE.NREACT) GO TO 4290
DO 4270 IN=1,NCHAN
EKG(IN)=EKG(IN)+G*E*K(IN,I-NREACT)
AKG(IN)=AKG(IN)+G*K(IN,I-NREACT)
STLP(IN)=STLP(IN)+G*(K(IN,I-NREACT)/(K(1,I-NREACT)
* +K(2,I-NREACT)+K(3,I-NREACT)))
C
C USING TRAPEZOIDAL RULE FOR INTEGRATION, TAKE HALF OF FIRST VALUE IN SUM.
C
IF(I.EQ.NREACT+1)THEN
EKG(IN)=0.5D0*EKG(IN)
AKG(IN)=0.5D0*AKG(IN)
STLP(IN)=0.5D0*STLP(IN)
ENDIF
4270 CONTINUE
GO TO 4290
4280 IF(I.LE.NREACT)GO TO 4290
AKHP(1)=AKHP(1)+G*XKHP(I-NREACT,1)
IF(I.EQ.NREACT+1)AKHP(1)=0.5D0*AKHP(1)
4290 CONTINUE
DO 4320 IN=1,NCHAN
IF(.NOT.JAVX)GO TO 4300
RATE(IN)=RTHI(IN)/(GI*DFLT(INC))
IF(IN.EQ.1)RHPN(IN)=AKHP(IN)/GI
WRHP(IN)=RTHI2(IN)/(GI*DFLT(INC))
C
C STLPJ IS THE EXACT TOTAL LOW PRESSURE STRONG COLLISION (E,J) RATE
C COEFFICIENT.STLP IS THE STRONG COLLISION (E,J=0) LOW PRESSURE
C RATE COEFFICENT.
C
STLPJ=SCLOW(1)*6.237D+4*T/(PRS*GI*DFLT(INC))
STLP(1)=SCLE*6.237D+4*T/(PRS*GI*DFLT(INC))
GO TO 4310
4300 CONTINUE
EAC(IN)=(EKG(IN)/AKG(IN))-(EGI/GI)
RATE(IN)=AKG(IN)/GI
ALNA(IN)=LOG10(RATE(IN))+(EAC(IN)/(T*2.303D-3*1.987D0))
EJ(IN)=EAC(IN)*4.184D0
STLP(IN)=STLP(IN)/GI
4310 RATHP(IN)=RATE(IN)
4320 CONTINUE
WRITE(6,400) TITLE,T
WRITE(6,405) (I,I=1,NCHAN)
WRITE(6,410) (EOK(I),I=1,NCHAN)
IF(JAVX)GO TO 4330
WRITE(6,415) (EAC(I),I=1,NCHAN)
WRITE(6,420) (EJ(I),I=1,NCHAN)
WRITE(6,425) (ALNA(I),I=1,NCHAN)
4330 CONTINUE
WRITE(6,430) (RATE(I),I=1,NCHAN)
IF(.NOT.JAVX)GO TO 4350
WRITE(6,435)(WRHP(I),I=1,NCHAN)
WRITE(6,440)(FWR(I),I=1,NCHAN)
DO 4340 I4340=1,NCHAN
IF(DABS(1.0D0-RATE(I4340)/WRHP(I4340)).GT.0.2D0)
* WRITE(6,445)I4340
4340 CONTINUE
4350 CONTINUE
C
C CALCULATE AND PRINT OUT THERMODYNAMIC PARAMETERS
C
DO 4360 IN=1,NCHAN
BCPLX1=BCMPLX(IN,1)
N1=N(IN,1)
BVEC(IN,1,N1+1)=BCPLX1
SIGVC(IN,1,N1+1)=SRC(IN,1)
IRTDMC(IN,1,N1+1)=IEXRTD(IN)
IF(IN.NE.1) GO TO 4360
4355 BVECM(NINTR+1)=BMOLEC
SIGVCM(NINTR+1)=SRM
IRTDMM(NINTR+1)=IEXMD
4360 CONTINUE
C
C INCLUDE EXTERNAL AND INTERNAL ROTATIONS TOGETHER FOR THERMODYNAMICS.
C
NRTOT=NINTR+1
CALL THERM(QV,QROT,RT,NM,JM,NF,CP,SV,SROT,ST,STOT,WT1,NRTOT,
1 BVECM,SIGVCM,IRTDMM,T,.TRUE.)
CHECK=.TRUE.
DO 4400 IN=1,NCHAN
JF1=JF(IN,1)
DO 4370 I4370=1,JF1
NC1(I4370)=NC(IN,1,I4370)
JC1(I4370)=JC(IN,1,I4370)
4370 CONTINUE
N1=N(IN,1)+1
DO 4380 I4380=1,N1
BVEC1(I4380)=BVEC(IN,1,I4380)
SIGVC1(I4380)=SIGVC(IN,1,I4380)
ISPARE(I4380)=IRTDMC(IN,1,I4380)
4380 CONTINUE
CALL THERM(V1,R1,RT,NC1,JC1,JF1,CP1,SV1,SROT1,ST1,STOT1,WT1,
1 N1,BVEC1,SIGVC1,ISPARE,T,.FALSE.)
DELTAS=STOT1-STOT
ALNAE(IN)=(2.71828D0*RT/H)*EXP(DELTAS/8.314D0)
ALNAE(IN)=LOG10(ALNAE(IN))
IF(JAVX)GO TO 4390
IF(DABS(ALNAE(IN)-ALNA(IN)).GT.0.05D0) WRITE(6,450)
4390 CONTINUE
CPC(IN)=CP1
QROTC(IN)=R1
QVC(IN)=V1
SVC(IN)=SV1
SROTC(IN)=SROT1
STC(IN)=ST1
STOTC(IN)=STOT1
IF(JAVX)GO TO 4400
RATTH(IN)=2.084D10*T*(QVC(IN)/QV)*(QROTC(IN)/QROT)
1 *EXP(-EOK(IN)*503.25D0/T)
IF(DABS(1.0D0-(RATTH(IN)/RATE(IN))).GT.1.0D-3) CHECK=.FALSE.
4400 CONTINUE
IF(JAVX)GO TO 4410
WRITE(6,455) (RATTH(IN),IN=1,NCHAN)
IF(CHECK) GO TO 4440
C
C RE-GENERATE THERMODYNAMICS WITHOUT INACTIVE ROTORS
C
4410 CALL THERM(QX,QROX,RT,NM,JM,NF,CX,SX,SROX,SX,STOX,WT1,NINTR,
1 BVECM,SIGVCM,IRTDMM,T,.TRUE.)
CORRAT = GI*DFLT(INC)/(QX*QROX)
DO 4420 IN=1,NCHAN
RATE(IN)=RATE(IN)*CORRAT
RATHP(IN)=RATE(IN)
4420 CONTINUE
IF(JAVX)GO TO 4430
WRITE(6,460) CORRAT
WRITE(6,465) (RATE(IN),IN=1,NCHAN)
GO TO 4440
4430 CONTINUE
C
C CORRAT HERE CORRECTS THE HIGH PRESSURE RATE FOR NUMERICAL
C ERROR IN THE MOLECULAR PARTITION FUNCTION.
C
WRITE(6,470)CORRAT,(RATE(IN),IN=1,NCHAN)
4440 CONTINUE
WRITE(6,475) (I,I=1,NCHAN)
C SPECIFIC HEAT CALCULATION CONTAINS ERROR IN CERTAIN CASES;
C OMIT PRINTOUT (RGG, OCT 25 1991)
C WRITE(6,480) CP,(CPC(IN),IN=1,NCHAN)
WRITE(6,485) SV, (SVC(IN),IN=1,NCHAN)
WRITE(6,490) SROT, (SROTC(IN),IN=1,NCHAN)
WRITE(6,495)
WRITE(6,500) ST,(STC(IN),IN=1,NCHAN)
WRITE(6,505) STOT,(STOTC(IN),IN=1,NCHAN)
IF(.NOT.JAVX)WRITE(6,510) (ALNAE(IN),IN=1,NCHAN)
IF(JAVX)WRITE(6,511) (ALNAE(IN),IN=1,NCHAN)
IF(.NOT.JAVX)GO TO 4460
C
C GENERATE ACTIVATION ENERGY FROM NUMERICAL RATE COEFFICIENT
C AND THERMODYNAMIC A FACTOR
C
DO 4450 I4450=1,NCHAN
EAC(I4450)=-1.987D-3*T*(LOG(RATE(I4450))
1 -2.303D0*ALNAE(I4450))
4450 CONTINUE
4460 CONTINUE
WRITE(6,515) QV,(QVC(IN),IN=1,NCHAN)
WRITE(6,520) QROT,(QROTC(IN),IN=1,NCHAN)
WRITE(6,495)
IF(JAVX)WRITE(6,405)(I,I=1,NCHAN)
IF(JAVX)WRITE(6,521)(EAC(I),I=1,NCHAN)
C
C HAVING FINISHED COMPUTATION OF HIGH PRESSURE PARAMETERS, FIND FALL-OFF
C REACTION RATES, USING STRONG COLLISION FORMULAE.
C INTEGRALS USE TRAPEZOIDAL RULE WITH HALF OF THE FIRST TERM. MJTJ 8.1.93
C
WRITE(6,525)
4470 CONTINUE
4500 CONTINUE
DO 4510 IN=1,NCHAN
RATE(IN)=0.0D0
4510 CONTINUE
BOT=0.0D0
DO 4550 II=1,NI
G=0.0D0
IF(RHOMOL(II).GT.0.0D0) G=EXP(LOG(RHOMOL(II))-DFLT(II)*FEXP)
IF(II.LE.NREACT) GO TO 4540
AKTOT=0.0D0
DO 4520 IN=1,NCHAN
AKTOT=AKTOT+K(IN,II-NREACT)
4520 CONTINUE
G=G*OMEGA/(OMEGA+AKTOT)
DO 4530 IN=1,NCHAN
RATE(IN)=RATE(IN)+G*K(IN,II-NREACT)
IF(II.EQ.NREACT+1)RATE(IN)=0.5D0*RATE(IN)
4530 CONTINUE
4540 BOT=BOT+G
IF(II.EQ.1)BOT=0.5D0*BOT
4550 CONTINUE
DO 4580 IN=1,NCHAN
C
C ADJUST CORRAT TO INCLUDE NUMERICAL ERROR INCURRED IN THE AVERAGE
C OVER K(E,J) AND ENSUING INTEGRATION OVER ENERGY E.(THE HIGH
C PRESSURE RATE IS CALCULATED DIRECTLY AND IS ESSENTIALLY EXACT
C EXCEPT FOR THE MOLECULAR PARTITION FUNCTION ERROR.)
C
IF(.NOT.JAVX)GO TO 4560
IF(IN.EQ.1)CORAT1=RATHP(1)/RHPN(IN)
RATE(IN)=RATE(IN)*CORAT1
GO TO 4570
4560 RATE(IN)=RATE(IN)*CORRAT
4570 RATE(IN)=RATE(IN)/BOT
RATIO(IN)=RATE(IN)/RATHP(IN)
4580 CONTINUE
WRITE(6,530) PRS,OMEGA,(RATE(IN),RATIO(IN),IN=1,NCHAN)
IF((.NOT.JAVX).AND.(J.LT.NP))GO TO 4610
IF(JAVX)GO TO 4600
DO 4590 I4590=1,NCHAN
STLP(I4590)=STLP(I4590)*OMEGA*6.237D4*T*CORRAT/PRS
4590 CONTINUE
WRITE(6,535) (STLP(I),I=1,NCHAN)
GO TO 4610
4600 IF(J.LT.NP)GO TO 4610
WRITE(6,540)STLPJ*CORRAT
WRITE(6,545)STLP(1)*CORRAT
4610 CONTINUE
C
C PREPARE INPUT FILE FOR MASTER EQUATION PROGRAM
C
JAV1=JAV
IF(JAV1.GT.1)JAV1=1
IF((J.EQ.1).AND.(I4999.EQ.1))GO TO 4700
IF(.NOT.JAVX)GO TO 4710
IF(JAVX)GO TO 4730
4700 CONTINUE
IDELT=400
REWIND 10
WRITE(10,550) TITLE,IDELT,NCHAN,INC
ERR1=1.0D-3
ERR2=1.0D-3
WRITE(10,555) ERR1,ERR2,ERR1
WRITE(10,560) EOK(1),NP
WRITE(10,555) (PRESS(IP),IP=1,NP)
WRITE(10,565)
WRITE(10,570) JAV1
WRITE(10,575)NT,(TEMP(I),I=1,NT)
WRITE(10,580) SGMA,WT1,WT2,EPS
IOPTHT=0
IOPTPR=NCHAN
WRITE(10,585) IOPTHT,IOPTPR
WRITE(10,570) NI
WRITE(10,555) (RHOMOL(I),I=1,NI)
IF(JAVX)GO TO 4720
4710 IF(J.NE.NP)GO TO 4990
WRITE(10,590)T,CORRAT
IF(I4999.NE.NT)GO TO 4990
GO TO 4740
4720 CONTINUE
WRITE(10,570)NTHR
WRITE(10,595)(ROTHR(I),I=1,NTHR)
4730 CONTINUE
EOEFF=(DFLT(NJC*INC)+5.0D0)/WKA
WRITE(10,600)(RATHP(I),I=1,NCHAN),CORAT1,CORRAT
WRITE(10,605)STLPJ*CORRAT
WRITE(10,590)EOEFF
4740 CONTINUE
WRITE(10,570)NMAX4
WRITE(10,555) (K(1,I),I=1,NMAX4)
DO 4750 I4750=2,NCHP
WRITE(10,555) (K(I4750,I),I=1,NMAX4)
4750 CONTINUE
GO TO 4990
C
C FOR ION-DIPOLE REACTION:
C
5000 CONTINUE
IF(.NOT.JAVX)GO TO 5300
C
C CALCULATE OMEGA FOR THIS PRESSURE
C
TEM=T
PRS=PRESS(J)
OMEGA=9.661D+18*PRS*KCAPT/TEM
DO 5100 IN=1,NCHAN
ROTINC1=ROTINC(IN)
ERR1=ERR(IN)
WKJ1=WKJ(IN)
IHINDX1=IHINDX(IN)
GAMMA1=GAMMA(IN)
KK=INT((DHO(IN)-DHO(1))*WKA/(DBLE(INC))+0.5D0)
KKC(IN)=KK
NN2=NN-KK
DHO1=DHO(IN)
JF1=JF(IN,1)
DO 5010 I5010=1,JF1
NC1(I5010)=NC(IN,1,I5010)
JC1(I5010)=JC(IN,1,I5010)
5010 CONTINUE
N1=N(IN,1)+NH(IN,1)
N11=MAX0(N1,1)
DO 5020 I5020=1,N11
BVEC1(I5020)=BVEC(IN,1,I5020)
SIGVC1(I5020)=SIGVC(IN,1,I5020)
ISPARE(I5020)=IRTDMC(IN,1,I5020)
5020 CONTINUE
V01=V0(IN,1)
CALL BSWINE(NC1,JC1,AK,NN2,JF1,INC,BVEC1,SIGVC1,
1 ISPARE,N1,.TRUE.)
DO 5030 I5030=1,NREACT+10
IF(ROTH(I5030).LE.0.0D0)GO TO 5040
ROTHR(I5030)=ROTH(I5030)
5030 CONTINUE
5040 NTHR=I5030-1
C
C RESET DEGENERACIES AFTER RE-ARRANGEMENT BY DIRECT COUNT.
C
JF(IN,1)=JF1
DO 5060 I=1,JF1
NC(IN,1,I)=NC1(I)
JC(IN,1,I)=1
5060 CONTINUE
IF(J.EQ.1)NMAX4=MIN0(NMAX4,NN2)
IF(J.GT.1)NMAX4=NI-NJC
DO 5080 I5080=1,NI
DO 5070 I5070=1,NCHAN
K(I5070,I5080)=WE(I5080,I5070)
5070 CONTINUE
5080 CONTINUE
5100 CONTINUE
NREACT=NJC
IF(J.EQ.1)NI=NREACT+NMAX4
DO 5130 IN=1,NCHAN
STEM=1.0D0/H
ICC=0
DO 5110 I5110=1,NMAX4
ICC=ICC+1
AK(I5110)=K(IN,ICC)*STEM/(RHOMOL(I5110+NREACT))
AK2(I5110)=XKHP(I5110,IN)*STEM/(RHOMOL(I5110+NREACT))
5110 CONTINUE
DO 5120 I5120=1,NMAX4
K(IN,I5120)=AK(I5120)
XKHP(I5120,IN)=AK2(I5120)
5120 CONTINUE
RTHI(IN)=RTHI(IN)/H
5130 CONTINUE
C
C COMPUTATION OF DENSITY OF STATES AND MICROSCOPIC RATES NOW COMPLETE.
C
C NOTE THE DIRECT COUNT SUBROUTINE UNSCRAMBLES DEGENERATE FREQUENCIES
C TO AVOID BUNCHING CAUSED BY DIRECT COUNT ALGORITHM.
C
C COMMENCE INTEGRATION TO FIND RATE COEFFICIENTS,ETC.
C INTEGRATION USES TRAPEZOIDAL RULE WITH HALF OF FIRST TERM. MJTJ 8.1.93
C
IF(J.NE.1)GO TO 5500
DO 5200 IN=1,NCHAN
STLP(IN)=0.0D0
AKHP(IN)=0.0D0
RATE(IN)=0.0D0
5200 CONTINUE
GI=0.0D0
FEXP=DELT/RT
C
C THE FOLLOWING LOOP FINDS HIGH-PRESSURE PARAMETERS.
C
DO 5230 I5230=1,NI
G=0.0D0
IF(RHOMOL(I5230).LE.0)GO TO 5230
G=EXP(LOG(RHOMOL(I5230))-DBLE(I5230)*FEXP)
GI=GI+G
IF(I5230.EQ.1)GI=0.5D0*GI
IF(I5230.LE.NREACT)GO TO 5230
DO 5210 IN=1,NCHAN
AKHP(IN)=AKHP(IN)+G*XKHP(I5230-NREACT,IN)
IF(I5230.EQ.NREACT+1)AKHP(IN)=0.5D0*AKHP(IN)
5210 CONTINUE
CHECK=.TRUE.
DO 5220 IN=1,NCHAN
CHECK=CHECK.AND.(G*XKHP(I5230-NREACT,IN).LT.ERR1*AKHP(IN))
5220 CONTINUE
IF(CHECK)GO TO 5240
5230 CONTINUE
WRITE(6,610)
STOP
5240 NI=I5230
NMAX4=NI-NREACT
DO 5250 IN=1,NCHAN
RATE(IN)=RTHI(IN)/(GI*DBLE(INC))
RATEMV(IN)=AKHP(IN)/GI
RATHP(IN)=RATE(IN)
5250 CONTINUE
C
C STLP(1) IS THE TOTAL LOW PRESSURE STRONG COLLISION (E,J) RATE
C COEFFICIENT. STLP(2) IS THAT FOR THE REACTANT CHANNEL IF
C JRECOM=1
C
STLP(1)=SCLOW(1)*6.237D+4*TEM/(PRS*GI*DBLE(INC))
IF(JRECOM.EQ.1)STLP(2)=SCLOW(2)*6.237D+4*TEM/(PRS*GI*DBLE(INC))
C
C CALCULATE AND PRINT OUT THERMODYNAMIC PARAMETERS
C
5300 DO 5310 IN=1,NCHAN
NPNH=N(IN,1)+NH(IN,1)
BCPLX1=BCMPLX(IN,1)
BVEC(IN,1,NPNH+1)=BCPLX1
SIGVC(IN,1,NPNH+1)=1
IRTDMC(IN,1,NPNH+1)=IEXRTD(IN)
IF(IN.NE.1) GO TO 5310
BVECM(NINTR+1)=BMOLEC
SIGVCM(NINTR+1)=1
IRTDMM(NINTR+1)=IEXMD
5310 CONTINUE
C
C INCLUDE EXTERNAL AND INTERNAL ROTATIONS TOGETHER FOR THERMODYNAMICS.
C
NRTOT=NINTR+1
CALL THERM(QV,QROT,RT,NM,JM,NF,CP,SV,SROT,ST,STOT,WT1,NRTOT,
1 BVECM,SIGVCM,IRTDMM,T,.TRUE.)
C
C CALCULATE CONTRIBUTION FROM EXTERNAL ROTATION
C
REXQM=RT/BMOLEC
REXSRM=8.314D0*(LOG(REXQM)+1.0D0)
REXCPM=1.0D0
CHECK=.TRUE.
DO 5350 IN=1,NCHAN
NPNH=N(IN,1)+NH(IN,1)
JF1=JF(IN,1)
V01=V0(IN,1)
STERICX=STERIC(IN)
OMEGA1=THETA1(IN)*3.1415926D0/180.0D0
OMEGA2=THETA2(IN)*3.1415926D0/180.0D0
ISYMX1=ISYM1(IN)
ISYMX2=ISYM2(IN)
DO 5320 I5320=1,JF1
NC1(I5320)=NC(IN,1,I5320)
JC1(I5320)=JC(IN,1,I5320)
5320 CONTINUE
N1=NPNH+1
DO 5340 I5340=1,N1
BVEC1(I5340)=BVEC(IN,1,I5340)
SIGVC1(I5340)=SIGVC(IN,1,I5340)
ISPARE(I5340)=IRTDMC(IN,1,I5340)
5340 CONTINUE
HRCORR=1.0D0
SQC=0.0D0
CALL THERM(V1,R1,RT,NC1,JC1,JF1,CP1,SV1,SROT1,ST1,
1 STOT1,WT1,N1,BVEC1,SIGVC1,ISPARE,T,.FALSE.)
DELTAS=STOT1-STOT
ALNAE(IN)=(2.71828D0*RT/H)*EXP(DELTAS/8.314D0)
ALNAE(IN)=LOG10(ALNAE(IN))
CPC(IN)=CP1
QROTC(IN)=R1
QVC(IN)=V1
SVC(IN)=SV1
SROTC(IN)=SROT1
STC(IN)=ST1
STOTC(IN)=STOT1
REXQC(IN)=RT/BCMPLX(IN,1)
REXSRC(IN)=8.314D0*(LOG(REXQC(IN))+1.0D0)
REXCPC(IN)=1.0D0
IF(IN.NE.INCHAN)GO TO 5350
KEQ(IN)=5.33D-21*QROT*QV*REXQC(IN)*HRCORR/(QROTC(IN)*QVC(IN))
RED=WTA(IN)*WTB(IN)/(WTA(IN)+WTB(IN))
TMM=1.0D0/(RED*T)
KEQ(IN)=KEQ(IN)*(TMM**2/DSQRT(TMM))*EXP(DHO(IN)*WKA/RT)
CAPT(IN)=RATEMV(IN)*KEQ(IN)
STRED=6.8635D0*LOG10(RED)+11.4392D0*LOG10(T)-2.314D0
STRED=STRED*4.184D0
STREAC=STOTC(IN)-REXSRC(IN)-SQC+STRED
5350 CONTINUE
IF(.NOT.JAVX)GO TO 5420
C
C RE-GENERATE THERMODYNAMICS WITHOUT INACTIVE ROTORS
C
CALL THERM(QX,QROX,RT,NM,JM,NF,CX,SX,SROX,SX,STOX,WT1,NINTR,
1 BVECM,SIGVCM,IRTDMM,T,.TRUE.)
CORRAT = GI*DBLE(INC)/(QX*QROX)
DO 5410 IN=1,NCHAN
RATHP(IN)=RATE(IN)*CORRAT
5410 CONTINUE
PPLX1=PPL1*6.237D+4*TEM/(PRS*QROT*QV)
PPLX2=PPL2*6.237D+4*TEM/(PRS*QROT*QV)
PPLX1=PPLX1*STLP(2)*CORRAT/PPLX2
C
C CORRAT HERE CORRECTS THE HIGH PRESSURE RATE FOR NUMERICAL
C ERROR IN THE MOLECULAR PARTITION FUNCTION.
C
5420 WRITE(6,400) TITLE,T
WRITE(6,405) (I,I=1,NCHAN)
WRITE(6,410) (EOK(I),I=1,NCHAN)
IF(.NOT.JAVX)GO TO 5470
DO 5440 IN=1,NCHAN
RATE(IN)=RATEMV(IN)
RATHP(IN)=RATEMV(IN)*CORRAT
5440 CONTINUE
IF(JRECOM.NE.1)WRITE(6,430) (RATHP(IN),IN=1,NCHAN)
IF(JRECOM.NE.1)WRITE(6,470)CORRAT,(RATHP(IN),IN=1,NCHAN)
IF(JRECOM.NE.1)GO TO 5460
WRITE(6,615)KEQ(INCHAN)
WRITE(6,620)CAPT(INCHAN)*CORRAT
IF(NCHAN.EQ.1)GO TO 5460
WRITE(6,625)
DO 5450 IN=1,NCHAN
IF(IN.EQ.INCHAN)GO TO 5450
WRITE(6,630)IN,BRATIO(IN)
5450 CONTINUE
IF(RADST.GT.(1.0D-4))WRITE(6,635)RSTAB
5460 CONTINUE
IF(JRECOM.EQ.1)WRITE(6,640)PLIFE*1.0D+6
IF(JRECOM.EQ.1)WRITE(6,645)(UNIR(IN),IN=1,NCHAN)
IF(RADST.GT.(1.0D-4))WRITE(6,650)RADST
5470 WRITE(6,475) (IN,IN=1,NCHAN)
C WRITE(6,480) CP,(CPC(IN),IN=1,NCHAN)
WRITE(6,485) SV, (SVC(IN),IN=1,NCHAN)
WRITE(6,490) SROT, (SROTC(IN),IN=1,NCHAN)
WRITE(6,495)
WRITE(6,500) ST,(STC(IN),IN=1,NCHAN)
WRITE(6,505) STOT,(STOTC(IN),IN=1,NCHAN)
WRITE(6,510) (ALNAE(IN),IN=1,NCHAN)
C
C GENERATE ACTIVATION ENERGY FROM NUMERICAL RATE COEFFICIENT
C AND THERMODYNAMIC A FACTOR
C
IF(.NOT.JAVX)GO TO 5490
DO 5480 IN=1,NCHAN
EAC(IN)=-1.987D-3*T*(LOG(RATE(IN))
1 -2.303D0*ALNAE(IN))
5480 CONTINUE
5490 CONTINUE
WRITE(6,515) QV,(QVC(IN),IN=1,NCHAN)
WRITE(6,520) QROT,(QROTC(IN),IN=1,NCHAN)
WRITE(6,495)
WRITE(6,660)
C WRITE(6,480)REXCPM,(REXCPC(IN),IN=1,NCHAN)
WRITE(6,490)REXSRM,(REXSRC(IN),IN=1,NCHAN)
WRITE(6,520)REXQM,(REXQC(IN),IN=1,NCHAN)
IF(JRECOM.EQ.1)WRITE(6,665)STREAC
IF(.NOT.JAVX)GO TO 4999
WRITE(6,405)(IN,IN=1,NCHAN)
WRITE(6,415)(EAC(IN),IN=1,NCHAN)
C
C HAVING FINISHED COMPUTATION OF HIGH PRESSURE PARAMETERS, FIND FALL-OFF
C REACTION RATES, USING STRONG COLLISION FORMULAE.
C ADJUST CORRAT TO INCLUDE NUMERICAL ERROR INCURRED IN THE AVERAGE
C OVER K(E,J) AND ENSUING INTEGRATION OVER ENERGY E.(THE HIGH
C PRESSURE RATE IS CALCULATED DIRECTLY AND IS ESSENTIALLY EXACT
C EXCEPT FOR THE MOLECULAR PARTITION FUNCTION ERROR.)
C
C INTEGRATE USING TRAPEZOIDAL RULE TAKING HALF OF FIRST TERM. MJTJ 8.1.93
C
IF(WKJ1)GO TO 5600
WRITE(6,525)
IF(JRECOM.NE.1)WRITE(6,670)
IF(JRECOM.EQ.1)WRITE(6,675)
5500 CONTINUE
IF(WKJ1)GO TO 5600
DO 5510 IN=1,NCHAN
RATE(IN)=0.0D0
5510 CONTINUE
BOT=0.0D0
FEXP=DELT/RT
DO 5550 II=1,NI
G=0.0D0
IF(RHOMOL(II).GT.0.0D0) G=EXP(LOG(RHOMOL(II))-DBLE(II)*FEXP)
IF(II.LE.NREACT) GO TO 5540
AKTOT=0.0D0
DO 5520 IN=1,NCHAN
AKTOT=AKTOT+K(IN,II-NREACT)
5520 CONTINUE
G=G*OMEGA/(OMEGA+AKTOT)
DO 5530 IN=1,NCHAN
RATE(IN)=RATE(IN)+G*K(IN,II-NREACT)
IF(II.EQ.NREACT+1)RATE(IN)=0.5D0*RATE(IN)
5530 CONTINUE
5540 BOT=BOT+G
IF(II.EQ.1)BOT=0.5D0*BOT
5550 CONTINUE
DO 5560 IN=1,NCHAN
RATE(IN)=RATE(IN)*CORRAT
RATE(IN)=RATE(IN)/BOT
RATIO(IN)=RATE(IN)/RATHP(IN)
5560 CONTINUE
IF(JRECOM.NE.1)GO TO 5580
DO 5570 IN=1,NCHAN
IF(IN.EQ.INCHAN)RATE(IN)=RATE(IN)*KEQ(INCHAN)
IF(IN.NE.INCHAN)RATE(IN)=0.0D0
IF(IN.NE.INCHAN)RATIO(IN)=0.0D0
5570 CONTINUE
5580 CONTINUE
WRITE(6,530) PRS,OMEGA,(RATE(IN),RATIO(IN),IN=1,NCHAN)
5600 IF(J.LT.NP)GO TO 5610
IF(JRECOM.NE.1)WRITE(6,540)STLP(1)*CORRAT
IF(JRECOM.EQ.1)WRITE(6,680)PPLX1*KEQ(INCHAN)/6.023D+23
5610 CONTINUE
C
C PREPARE INPUT FILE FOR MASTER EQUATION PROGRAM
C
IF((J.EQ.1).AND.(I4999.EQ.1))GO TO 5700
GO TO 5710
5700 CONTINUE
IDELT=400
REWIND 10
WRITE(10,550) TITLE,IDELT,NCHAN,INC
ERR6=1.0D-3
ERR2=1.0D-3
WRITE(10,555) ERR6,ERR2,ERR6
WRITE(10,685) EOK(1)
WRITE(10,690)GAMMA(1)
WRITE(10,570)NP
WRITE(10,555) (PRESS(IP),IP=1,NP)
WRITE(10,565)
WRITE(10,575)NT,(TEMP(I),I=1,NT)
WRITE(10,695)-KCAPT
IOPTHT=0
IOPTPR=NCHAN
WRITE(10,585) IOPTHT,IOPTPR
WRITE(10,570) NI
WRITE(10,555) (RHOMOL(I),I=1,NI)
WRITE(10,570)NTHR
WRITE(10,595)(ROTHR(I),I=1,NTHR)
5710 CONTINUE
EOEFF=(DBLE(NJC*INC)+5.0D0)/WKA
WRITE(10,600)(RATHP(I),I=1,NCHAN),CORRAT,CORRAT
WRITE(10,605)STLP(1)*CORRAT
WRITE(10,590)EOEFF
WRITE(10,570)NMAX4
WRITE(10,555) (K(1,I),I=1,NMAX4)
DO 5720 I5720=2,NCHP
WRITE(10,555) (K(I5720,I),I=1,NMAX4)
5720 CONTINUE
4990 CONTINUE
4999 CONTINUE
CLOSE(5,STATUS='KEEP')
CLOSE(6,STATUS='KEEP')
CLOSE(10,STATUS='KEEP')
STOP
END
C
C
SUBROUTINE BSWINE(NM,JM,RHO,NI,NF,INC,BVEC,SIGVC,IROT,NINTR,
1 SUM)
C
C SUBROUTINE FOR DIRECT COUNT TO GIVE DENSITY OF STATES
C OR SUM OF STATES.
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
PARAMETER (NMAX8=8000, NMAX5=5000, NMAX80=80000)
DIMENSION NM(100),JM(100),BVEC(20),SIGVC(20),IROT(20),RHO(NMAX8)
C
LOGICAL REORD,SUM,JAVX,IHINDX1,IONX1,STERICX
C
COMMON /ALL/ R,RT,H,WKA
COMMON /BIG2/ PFT(NMAX8)
COMMON /BS/ FRE(200),TI(NMAX80),TII(NMAX80),DELTAE
COMMON /INTBS/ M,IR(200),IS
COMMON /ISTP/ ISTEP
COMMON /JPARTF/ PF(NMAX8,3),PFC(3),PFM(NMAX8)
COMMON /JSUM/ TI2(NMAX8,3,1)
COMMON /LOOP/ IN,NCHAN
COMMON /LOGVAR/ IONX1,IHINDX1,JAVX,STERICX
C
C DELTAE IS THE BEYER-SWINEHART ENERGY GRAIN
C
IS=0
REORD=.TRUE.
DO 6000 I6000=1,NMAX8
RHO(I6000)=0.00D+0
6000 CONTINUE
EMAX=DFLT(NI*INC)
ISTEP=INT(DFLT(INC)/DELTAE)
CN=1.0D0/(DELTAE*DFLT(ISTEP))
IF(.NOT.REORD) GO TO 6050
SPACE=DELTAE*0.5D0
C
C SEPARATE ANY DEGENERATE FREQUENCIES TO AVOID BEYER-SWINEHART
C BUNCHING PROBLEMS.
C
DO 6030 I=1,NF
NTEMP=JM(I)
J=INT(DFLT(NTEMP)*0.5D0)
JR=NTEMP-2*J
VAL=DFLT(NM(I))
IF(JR.EQ.0) GO TO 6010
IS=IS+1
FRE(IS)=VAL
6010 CONTINUE
IF(J.LT.1) GO TO 6030
DO 6020 IX=1,J
IS=IS+1
BIT=SPACE*DFLT(IX)
FRE(IS)=VAL+BIT
IS=IS+1
FRE(IS)=VAL-BIT
6020 CONTINUE
6030 CONTINUE
C
C HAVING UNSCRAMBLED FREQUENCIES, REPLACE THE OLD WITH THE NEW VALUES.
C
NF=IS
DO 6040 I6040=1,NF
NM(I6040)=INT(FRE(I6040))
6040 CONTINUE
GO TO 6070
6050 CONTINUE
DO 6060 I6060=1,NF
FRE(I6060)=DFLT(NM(I6060))
6060 CONTINUE
IS=NF
6070 CONTINUE
CALL GRAIN
M=1.0D0+EMAX/DELTAE
IF(M.GT.NMAX80) WRITE(6,70)
70 FORMAT(' MAXIMUM ENERGY TOO GREAT FOR DIRECT COUNT',
1 ' ARRAY, ABORT')
IF(M.GT.NMAX80) STOP
C
C THE TNONH SUBROUTINE USES THE LOGICAL VARIABLE 'SUM'
C TO CALCULATE THE ANALYTIC ROTATIONAL SUM OF STATES
C EXPRESSION.
C
C TO SPEED CALCULATION THERE ARE TWO VERSIONS OF THE TNONH
C SUBROUTINE; INCORPORATING STERIC INTERACTIONS EXPLICITLY
C SIGNIFICANTLY SLOWS DOWN THE CALCULATION.
C
IF(STERICX)CALL TNONHS(BVEC,SIGVC,IROT,NINTR,SUM)
IF(.NOT.STERICX)CALL TNONHF(BVEC,SIGVC,IROT,NINTR,SUM)
C
C THE COUNT SUBROUTINE INCLUDES THE VIBRATIONAL DENSITY OF STATES
C INTO THE TI MATRIX
C
CALL COUNT(SUM)
IF((IONX1).AND.(.NOT.SUM))GO TO 6080
IF((.NOT.IONX1).AND.(JAVX))GO TO 6080
IF((.NOT.IONX1).AND.(.NOT.SUM))GO TO 6190
IF(.NOT.IONX1)GO TO 6140
6080 CONTINUE
GCI=0.0D0
FEXP1=DELTAE/RT
DO 6090 I6090=1,NI
J6090=NI+1-I6090
GC=0.0D0
JT10=J6090*ISTEP
6085 IF(TI(JT10).LE.0.0D0)GO TO 6090
GC=EXP(LOG(TI(JT10))-DFLT(JT10)*FEXP1)
GCI= GCI+GC
PFT(J6090)=GCI*ISTEP
6090 CONTINUE
PFCT=GCI*ISTEP
IF(SUM)GO TO 6120
DO 6110 I6110=1,NI
PFM(I6110)=PFT(I6110)
6110 CONTINUE
GO TO 6190
6120 CONTINUE
DO 6130 I6130=1,NI
PF(I6130,IN)=PFT(I6130)
6130 CONTINUE
PFC(IN)=PFCT
IF(JAVX.AND.(NINTR.GT.0))GO TO 6160
6140 CONTINUE
SUMS=0.0D0
DO 6150 I6150=1,M
IF(NINTR.GT.0)TI(I6150)=TII(I6150)
IF(NINTR.GT.0)GO TO 6150
SUMS=SUMS+TI(I6150)
TII(I6150)=DELTAE*SUMS
IF(.NOT.JAVX)TI(I6150)=TII(I6150)
6150 CONTINUE
IF(.NOT.JAVX)GO TO 6190
6160 DO 6170 I6170=1,NI
TI2(I6170,IN,1)=TII(I6170*ISTEP)
6170 CONTINUE
IF(IN.LT.NCHAN)GO TO 6180
CALL JAVRGE(RHO,DELTAE,M,NI,ISTEP,INC)
6180 RETURN
6190 CONTINUE
ISUM=0
DO 6210 I6210=1,M,ISTEP
ISUM=ISUM+1
IP=I6210+ISTEP-1
AV=0.0D0
DO 6200 J6200=I6210,IP
AV=AV+TI(J6200)
6200 CONTINUE
RHO(ISUM)=AV*CN
6210 CONTINUE
IF(ISUM.GE.NI) RETURN
ISUM=ISUM+1
DO 6220 I6220=ISUM,NI
RHO(I6220)=0.0D0
6220 CONTINUE
RETURN
END
C
C
SUBROUTINE TNONHF(BVEC,SIGVC,IROT,NINTR,SUM)
C
C INITIALIZES ARRAY FOR BEYER-SWINEHART ALGORITHM, INCLUDING SEMI-CLASSICAL
C FORM FOR FREE INTERNAL ROTORS, AS SUGGESTED BY TROE.
C
C INCLUDES EXTENSION TO SEMI-CLASSICAL FREE ROTOR EXPRESSION TO
C INCLUDE SINUSOIDALLY HINDERED ROTOR, (JORDAN,SMITH,GILBERT).
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
PARAMETER (NMAX8=8000, NMAX5=5000, NMAX80=80000)
C
DIMENSION BVEC(20),SIGVC(20),IROT(20)
C
LOGICAL SUM,IHINDX1,IONX1,JAVX,STERICX
C
COMMON /BS/ FRE(200),TI(NMAX80),TII(NMAX80),DELTAE
COMMON /INTBS/ M,IR(200),IS
COMMON /ROTF/ GAMON2(40),ROTTM(3),V01,OMEGA1,OMEGA2
COMMON /ROTINT/ ISYMX1,ISYMX2,NH(3,1)
COMMON /LOOP/ IN,NCHAN
COMMON /LOGVAR/ IONX1,IHINDX1,JAVX,STERICX
C
IIR=0
Q=1.0D0
KZ1=0.0D0
KZ2=0.0D0
IF(NINTR.LE.0) GO TO 6600
DO 6500 I=1,NINTR
IRX=IROT(I)
IIR=IIR+IRX
HD=0.5D0*DFLT(IRX)
Q=Q*GAMON2(IRX)/(SIGVC(I)*(BVEC(I)**HD))
6500 CONTINUE
HR=0.5D0*DFLT(IIR)
EROT=0.0D0
Q=DELTAE*Q/GAMON2(IIR)
6600 CONTINUE
DO 6900 I6900=1,M
TI(I6900)=0.0D0
TII(I6900)=0.0D0
IF(NINTR.LE.0) GO TO 6900
EROT=EROT+DELTAE
TI(I6900)=Q*(EROT**HR)/EROT
IF(.NOT.SUM)GO TO 6900
C
C TI IS THE DENSITY OF STATES TERM THAT IS USED TO WORK OUT PARTITION
C FUNCTIONS IN THE BSWINE SUBROUTINE.
C
C TII IS THE ANALYTIC FORM OF THE ROTATIONAL SUM OF STATES. TII IS
C CALCULATED IF SUM=.TRUE.
C USING THE TII MATRIX SAVES DOING A NUMERICAL INTEGRATION IN THE
C BSWINE SUBROUTINE.
C
TII(I6900)=TI(I6900)*EROT*GAMON2(IIR)/GAMON2(IIR+2)
IF(.NOT.IHINDX1)GO TO 6900
IF(NH(IN,1).LE.0)GO TO 6900
IF(IONX1)NH(IN,1)=1
C
C THE FOLLOWING IS THE CORRECTION FACTOR FOR HINDERED ROTATIONS
C IN THE DENSITY AND SUM OF STATES TERMS FOR THE TRANSITION STATE.
C
C IF NH=1 : ONE 2-DIM SINUSOIDALLY HINDERED ROTOR
C IF NH=2 : TWO 2-DIM SINUSOIDALLY HINDERED ROTORS
C
IF(V01.LE.DELTAE)GO TO 6900
C
C IE THE SINUSOIDALLY HINDERED ROTOR WILL ONLY HAVE AN EFFECT IF ITS
C ENERGY, V01, IS LARGER THAN THE ENERGY GRAIN SIZE, DELTAE.
C
KZ1=INT((V01/DELTAE)+0.5D0)
IF(NH(IN,1).EQ.2)KZ2=INT(((V01+V01)/DELTAE)+0.5D0)
TI(I6900)=TI(I6900)*EROT*GAMON2(IIR)/(GAMON2(IIR+2)*V01)
TII(I6900)=TII(I6900)*EROT*GAMON2(IIR+2)/(GAMON2(IIR+4)*V01)
C
C DIFFERENT FORM FOR NH=2:
C
IF(NH(IN,1).EQ.2)TI(I6900)=TI(I6900)*EROT*GAMON2(IIR+2)/
1 (GAMON2(IIR+4)*V01)
IF(NH(IN,1).EQ.2)TII(I6900)=TII(I6900)*EROT*GAMON2(IIR+4)/
1 (GAMON2(IIR+6)*V01)
C
C INCLUDING THE STEP FUNCTIONS IN THE DENSITY AND SUM OF STATES TERMS:
C
IF(I6900.LE.KZ1)GO TO 6800
IF(NH(IN,1).EQ.1)TEMPQ=Q*((EROT-V01)**HR)*GAMON2(IIR)/
1 (GAMON2(IIR+2)*V01)
IF(NH(IN,1).EQ.1)TEMPQ1=Q*((EROT-V01)**(HR+1.0D0))*
1 GAMON2(IIR)/(GAMON2(IIR+4)*V01)
C
C DIFFERENT FORM FOR NH=2:
C
IF(NH(IN,1).EQ.2)TEMPQ=Q*((EROT-V01)**(HR+1.0D0))*
1 GAMON2(IIR)/(GAMON2(IIR+4)*V01*V01)
IF(NH(IN,1).EQ.2)TEMPQ1=Q*((EROT-V01)**(HR+2.0D0))*
1 GAMON2(IIR)/(GAMON2(IIR+6)*V01*V01)
C
TI(I6900)=TI(I6900)-NH(IN,1)*TEMPQ
TII(I6900)=TII(I6900)-NH(IN,1)*TEMPQ1
C
6800 IF(NH(IN,1).EQ.1)GO TO 6900
C
C IF WE HAVE 2 D-DIM SINUSOIDALLY HINDERED ROTORS:
C
IF(I6900.LE.KZ2)GO TO 6900
TEMPR=Q*((EROT-V01)**(HR+1.0D0))*GAMON2(IIR)/
1 (GAMON2(IIR+4)*V01*V01)
TEMPR1=Q*((EROT-V01)**(HR+2.0D0))*GAMON2(IIR)/
1 (GAMON2(IIR+6)*V01*V01)
C
TI(I6900)=TI(I6900)+TEMPR
TII(I6900)=TII(I6900)+TEMPR1
C
6900 CONTINUE
IF(NINTR.LE.0) TI(1)=1.0D0
RETURN
END
C
C
SUBROUTINE TNONHS(BVEC,SIGVC,IROT,NINTR,SUM)
C
C INITIALIZES ARRAY FOR BEYER-SWINEHART ALGORITHM, INCLUDING SEMI-CLASSICAL
C FORM FOR FREE INTERNAL ROTORS, AS SUGGESTED BY TROE.
C
C INCLUDES EXTENSION TO SEMI-CLASSICAL FREE ROTOR EXPRESSION TO
C INCLUDE SINUSOIDALLY HINDERED ROTOR, (JORDAN,SMITH,GILBERT).
C
C INCLUDES EXPLICITLY STERIC INTERACTIONS IN THE THETA EULER ANGLE FOR ONE OR
C BOTH FRAGMENTS: HINDRANCE ANGLES ARE LABELLED OMEGA1 AND OMEGA2, SYMMETRY
C NUMBERS FOR THE 2-DIM ROTATIONS LABELLED ISYMX1 AND ISYMX2,
C (JORDAN,SMITH,GILBERT).
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
PARAMETER (NMAX8=8000, NMAX5=5000, NMAX80=80000)
C
DIMENSION BVEC(20),SIGVC(20),IROT(20)
LOGICAL SUM,IHINDX1,IONX1,JAVX,STERICX
C
COMMON /BS/ FRE(200),TI(NMAX80),TII(NMAX80),DELTAE
COMMON /INTBS/ M,IR(200),IS
COMMON /ROTF/ GAMON2(40),ROTTM(3),V01,OMEGA1,OMEGA2
COMMON /ROTINT/ ISYMX1,ISYMX2,NH(3,1)
COMMON /LOOP/ IN,NCHAN
COMMON /LOGVAR/ IONX1,IHINDX1,JAVX,STERICX
C
IIR=0
Q=1.0D0
KZ1=0.0D0
KZ2=0.0D0
KZ3=0.0D0
VSCALE1=0.0D0
VSCALE2=0.0D0
IF(NINTR.LE.0) GO TO 6600
DO 6500 I=1,NINTR
IRX=IROT(I)
IIR=IIR+IRX
HD=0.5D0*DFLT(IRX)
Q=Q*GAMON2(IRX)/(SIGVC(I)*(BVEC(I)**HD))
6500 CONTINUE
HR=0.5D0*DFLT(IIR)
EROT=0.0D0
Q=DELTAE*Q/GAMON2(IIR)
IF(.NOT.SUM)GO TO 6600
IF(NH(IN,1).EQ.0)GO TO 6600
C
C VSCALE1 AND VSCALE2 ARE THE SCALED BARRIER HEIGHTS THAT ARE EQUIVALENT
C TO INCORPORATING THE STERIC HINDRANCE ANGLES.
C
VSCALE1=V01*ISYMX1*0.5D0*(1.0D0-COS(OMEGA1))
IF(NH(IN,1).EQ.1)GO TO 6600
VSCALE2=V01*ISYMX2*0.5D0*(1.0D0-COS(OMEGA2))
C
6600 CONTINUE
DO 6900 I6900=1,M
TI(I6900)=0.0D0
TII(I6900)=0.0D0
IF(NINTR.LE.0) GO TO 6900
EROT=EROT+DELTAE
TI(I6900)=Q*(EROT**HR)/EROT
IF(.NOT.SUM)GO TO 6900
C
C TI IS THE DENSITY OF STATES TERM THAT IS USED TO WORK OUT PARTITION
C FUNCTIONS IN THE BSWINE SUBROUTINE.
C
C TII IS THE ANALYTIC FORM OF THE ROTATIONAL SUM OF STATES. TII IS
C CALCULATED IF SUM=.TRUE.
C USING THE TII MATRIX SAVES DOING A NUMERICAL INTEGRATION IN THE
C BSWINE SUBROUTINE.
C
TII(I6900)=TI(I6900)*EROT*GAMON2(IIR)/GAMON2(IIR+2)
IF(.NOT.IHINDX1)GO TO 6900
IF(NH(IN,1).LE.0)GO TO 6900
IF(IONX1)NH(IN,1)=1
C
C THE FOLLOWING IS THE CORRECTION FACTOR FOR HINDERED ROTATIONS
C IN THE DENSITY AND SUM OF STATES TERMS FOR THE TRANSITION STATE.
C
C IF NH=1 : ONE 2-DIM SINUSOIDALLY HINDERED ROTOR
C IF NH=2 : TWO 2-DIM SINUSOIDALLY HINDERED ROTORS
C
IF(V01.LE.DELTAE)GO TO 6900
C
C IE THE SINUSOIDALLY HINDERED ROTOR WILL ONLY HAVE AN EFFECT IF ITS
C ENERGY, V01, IS LARGER THAN THE ENERGY GRAIN SIZE, DELTAE.
C
KZ1=INT((VSCALE1/DELTAE)+0.5D0)
IF(NH(IN,1).EQ.2)KZ2=INT((VSCALE2/DELTAE)+0.5D0)
IF(NH(IN,1).EQ.2)KZ3=INT(((VSCALE1+VSCALE2)/DELTAE)+0.5D0)
TI(I6900)=TI(I6900)*EROT*GAMON2(IIR)/(GAMON2(IIR+2)*V01)
TII(I6900)=TII(I6900)*EROT*GAMON2(IIR+2)/(GAMON2(IIR+4)*V01)
C
C DIFFERENT FORM FOR NH=2:
C
IF(NH(IN,1).EQ.2)TI(I6900)=TI(I6900)*EROT*GAMON2(IIR+2)/
1 (GAMON2(IIR+4)*V01)
IF(NH(IN,1).EQ.2)TII(I6900)=TII(I6900)*EROT*GAMON2(IIR+4)/
1 (GAMON2(IIR+6)*V01)
C
C INCORPORATING THE STEP FUNCTION IN THE DENSITY AND SUM OF STATES TERMS -
C THIS IS DONE USING PARAMETERS KZ1, KZ2, KZ3.
C
IF(I6900.LE.KZ1)GO TO 6800
IF(NH(IN,1).EQ.1)TEMPQ=Q*((EROT-VSCALE1)**HR)*GAMON2(IIR)/
1 (GAMON2(IIR+2)*V01)
IF(NH(IN,1).EQ.1)TEMPQ1=Q*((EROT-VSCALE1)**(HR+1.0D0))*
1 GAMON2(IIR)/(GAMON2(IIR+4)*V01)
C
C DIFFERENT FORM FOR NH=2:
C
IF(NH(IN,1).EQ.2)TEMPQ=Q*((EROT-VSCALE1)**(HR+1.0D0))*
1 GAMON2(IIR)/(GAMON2(IIR+4)*V01*V01)
IF(NH(IN,1).EQ.2)TEMPQ1=Q*((EROT-VSCALE1)**(HR+2.0D0))*
1 GAMON2(IIR)/(GAMON2(IIR+6)*V01*V01)
C
TI(I6900)=TI(I6900)-TEMPQ
TII(I6900)=TII(I6900)-TEMPQ1
6800 IF(NH(IN,1).EQ.1)GO TO 6900
C
C IF WE HAVE 2 D-DIM SINUSOIDALLY HINDERED ROTORS:
C
IF(I6900.LE.KZ2)GO TO 6900
TEMPR=Q*((EROT-VSCALE2)**(HR+1.0D0))*GAMON2(IIR)/
1 (GAMON2(IIR+4)*V01*V01)
TEMPR1=Q*((EROT-VSCALE2)**(HR+2.0D0))*GAMON2(IIR)/
1 (GAMON2(IIR+6)*V01*V01)
TI(I6900)=TI(I6900)-TEMPR
TII(I6900)=TII(I6900)-TEMPR1
IF(I6900.LE.KZ3)GO TO 6900
TEMPS=Q*((EROT-(VSCALE1+VSCALE2))**(HR+1.0D0))*GAMON2(IIR)/
1 (GAMON2(IIR+4)*V01*V01)
TEMPS1=Q*((EROT-(VSCALE1+VSCALE2))**(HR+2.0D0))*GAMON2(IIR)/
1 (GAMON2(IIR+6)*V01*V01)
TI(I6900)=TI(I6900)+TEMPS
TII(I6900)=TII(I6900)+TEMPS1
6900 CONTINUE
IF(NINTR.LE.0) TI(1)=1.0D0
RETURN
END
C
C
SUBROUTINE THERM(QV,QR,RT,NM,JM,NF,CP,SV,SROT,ST,STOT,WT1,NINTR,
1 BVEC,SVECM,IRDIM,T,MOLTH)
C
C CALCULATES ALL THERMODYNAMIC QUANTITIES.
C EXPANDED TO ALSO CALCULATE EXPRESSIONS USING THE VARIOUS SINUSOIDALLY
C HINDERED (AND STERICALLY HINDERED) ROTOR RESULTS (JORDAN,SMITH,GILBERT).
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
PARAMETER (NMAX8=8000, NMAX5=5000, NMAX80=80000)
C
COMMON /ROTF/ GAMON2(40),ROTTM(3),V01,OMEGA1,OMEGA2
COMMON /ROTINT/ ISYMX1,ISYMX2,NH(3,1)
COMMON /LOOP/ IN,NCHAN
COMMON /LOGVAR/ IONX1,IHINDX1,JAVX,STERICX
COMMON /THML/ HRCORR,SQC
C
LOGICAL MOLTH,IHINDX1,IONX1,JAVX,STERICX
DIMENSION NM(100),JM(100),BVEC(30),SVECM(20),IRDIM(20)
C
CV=0.0D0
QV=1.0D0
SV=0.0D0
VSCALE1=0.0D0
VSCALE2=0.0D0
VSCALE1=0.50D0*V01*ISYMX1*(1.0D0-COS(OMEGA1))
IF(NH(IN,1).EQ.2)VSCALE2=0.5D0*V01*ISYMX2*(1.0D0-COS(OMEGA2))
DO 7005 I7005=1,NF
X=DFLT(NM(I7005))/(RT)
EX=EXP(X)
REX=1.0D0/EX
EXL1=EX-1
C
C NOTE: AFTER THE BSWINE SUBROUTINE THE FREQUENCIES HAVE ALL
C BEEN SEPARATED AND THEIR DEGENERACIES SET TO 1
C
CV=CV+JM(I7005)*X*X*EX/(EXL1*EXL1)
SV=SV+JM(I7005)*(-LOG(1.0D0-REX)+(X/EXL1))
QV=QV/((1.0D0-REX)**JM(I7005))
7005 CONTINUE
SV=SV*8.314D0
ST=6.8635D0*LOG10(WT1)+11.4392D0*LOG10(T)-2.314D0
ST=ST*4.184D0
SROT=0.0D0
QR=1.0D0
IIR=0
IF(NINTR.LE.0) GO TO 7035
DO 7015 I7015=1,NINTR
IND=IRDIM(I7015)
IIR=IIR+IND
PWR=0.5D0*DFLT(IND)
QR=QR*((RT/BVEC(I7015))**PWR)*ROTTM(IND)/SVECM(I7015)
C
C ROTTM IS ACTUALLY JUST THE GAMMA FUNCTION EVALUATED AT ONE HALF
C OF THE ARGUMENT VALUE, IE FOR 1/2, 1 AND 3/2.
C
7015 CONTINUE
IF(MOLTH)GO TO 7025
C
C MOLTH IS A FLAG FOR GENERATING THE THERMODYNAMIC PARAMETERS OF
C JUST THE MOLECULE.
C MOLTH=.FALSE. IMPLIES THAT WE ARE LOOKING AT THE COMPLEX.
C
IF(.NOT.IHINDX1)GO TO 7025
C
C THE EXPRESSION FOR THE ROTATIONAL PARTITION FUNCTION IS
C DIFFERENT IF SINUSOIDALLY HINDERED ROTORS ARE USED.
C
IF(NH(IN,1).LE.0)GO TO 7025
VQ=V01/RT
VQ1=VSCALE1/RT
EXPV1=1.0D0/EXP(VQ1)
HRCORR=(1.0D0-EXPV1)/VQ
QR=QR*HRCORR
SQC=(1.0D0-EXPV1*(1.0D0+VQ1))/(1.0D0-EXPV1)
SQC=8.314D0*SQC
SROT=SQC
IF(NH(IN,1).NE.2)GO TO 7025
VQ2=VSCALE2/RT
EXPV2=1.0D0/EXP(VQ2)
HRCORR2=(1.0D0-EXPV2)/VQ
QR=QR*HRCORR2
SQC2=(1.0D0-EXPV2*(1.0D0+VQ2))/(1.0D0-EXPV2)
SQC2=8.314D0*SQC2
SROT=SROT+SQC2
7025 CONTINUE
SROT=SROT+8.314D0*(LOG(QR) + 0.5D0*DFLT(IIR))
7035 CONTINUE
STOT=SROT+ST+SV
CP=8.314D0*(4.0D0+CV+DFLT(IIR)*0.5D0)
RETURN
END
C
C
SUBROUTINE GRAIN
C
C GRAIN CALCULATES FREQUENCIES IN GRAIN UNITS
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
PARAMETER (NMAX8=8000, NMAX5=5000, NMAX80=80000)
C
COMMON /BS/ FRE(200),TI(NMAX80),TII(NMAX80),DELTAE
COMMON /INTBS/ M,IR(200),IS
C
DO 8000 I8000=1,IS
FFR=FRE(I8000)/DELTAE
IR(I8000)=NINT(FFR)
8000 CONTINUE
RETURN
END
C
C
SUBROUTINE COUNT(SUM)
C
C CALCULATES DENSITY OF STATES USING BEYER-SWINEHART ALGORITHM
C
C THE TI MATRIX IS THE ROTATIONAL DENSITY OF STATES MATRIX, WHEREAS
C THE TII MATRIX IS THE ROTATIONAL SUM OF STATES MATRIX.
C USING THE TII MATRIX ALLEVIATES THE NEED FOR A NUMERIAL INTEGRATION.
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
PARAMETER (NMAX8=8000, NMAX5=5000, NMAX80=80000)
C
COMMON /BS/ FRE(200),TI(NMAX80),TII(NMAX80),DELTAE
COMMON /INTBS/ M,IR(200),IS
C
LOGICAL SUM
DO 8001 I8001=1,IS
IR8001=1+IR(I8001)
DO 8002 I8002=IR8001,M
II8002=I8002-IR8001+1
TI(I8002)=TI(I8002)+TI(II8002)
IF(SUM)TII(I8002)=TII(I8002)+TII(II8002)
8002 CONTINUE
8001 CONTINUE
RETURN
END
C
C
SUBROUTINE MINIM(BETAX,ICOUNT,DBETA,ACB,CHECK)
C
C FINDS BEST VALUE OF BETA TO FIT INPUT POTENTIAL\
C TO MORSE CURVE.
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
PARAMETER (NMAX8=8000, NMAX5=5000, NMAX80=80000)
LOGICAL CHECK
CHECK=.FALSE.
FUN=DIFF(BETAX)
DO 8020 I8020=1,ICOUNT
FUNL=FUN
BETAX=BETAX+DBETA
FUN=DIFF(BETAX)
IF(FUN.LE.FUNL)GO TO 8010
DBETA=-DBETA/2.0D0
GO TO 8020
8010 IF(DABS(DBETA).LT.ACB)GO TO 8030
8020 CONTINUE
GO TO 8040
8030 CHECK=.TRUE.
8040 CONTINUE
RETURN
END
C
C SUBROUTINE WEAKJ CALCULTES THE CORRECTION TERM IF WEAK ROTATIONAL
C RELAXATION IS USED
C
SUBROUTINE WEAKJ(BSTAR,EILO,EISM,PHIL,EROT,
* INC,DELTAE,N)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DOUBLE PRECISION KEJT,BSTAR(3),EILO(3),EISM(3)
C
PARAMETER (NMAX8=8000, NMAX5=5000, NMAX80=80000)
COMMON /ALL/ R,RT,H,WKA
COMMON /JCOLL/ DELTA,GAMMA1,OMEGA
COMMON /JMICRO/ RHOMOL(NMAX8),ROTH(NMAX5),WE(NMAX8,3)
COMMON /JPARAM/ ROTINC1,ERR1
COMMON /JSUM/ TI2(NMAX8,3,1)
COMMON /LOGVAR/ IONX1,IHINDX1,JAVX,STERICX
COMMON /LOOP/ IN,NCHAN
C
LOGICAL IONX1,IHINDX1,JAVX,STERICX
C
T11=0.0D0
T11B=0.0D0
ERR3=1.0D-2
DO 8630 I8630=1,150
R1=EROT-I8630*ROTINC1
IF(R1.LT.-1.0D0)GO TO 8635
KEJT=0.0D0
SUMCH1=0.0D0
DO 8620 I8620=1,NCHAN
WD=0.0D0
ISUM1=0
THRSH=EILO(I8620)+BSTAR(I8620)*(EROT-R1)
IF(DBLE(N*INC).LE.THRSH)GO TO 8620
THRCP=EISM(I8620)+BSTAR(I8620)*(EROT-R1)
DEC11=DBLE(N*INC)-THRCP
IF(DEC11.LE.0.0D0)DEC11=1.0D0
ISUM1=INT(DEC11/INC)
REM1=MOD(DEC11,DBLE(INC))/DBLE(INC)
IF(ISUM1.LT.1)GO TO 8600
WD=((1.0D0-REM1)*TI2(ISUM1,I8620,1)
1 +REM1*TI2(ISUM1+1,I8620,1))/DELTAE
GO TO 8610
8600 WD=REM1*TI2(1,I8620,1)/DELTAE
8610 SUMCH1=SUMCH1+WD
8620 CONTINUE
KEJT=SUMCH1/(H*RHOMOL(N))
TT=EXP((R1-EROT)/GAMMA1)/(OMEGA+KEJT)
TTB=EXP((R1-EROT)/GAMMA1)
IF(R1.LT.5.0D0)TT=TT/2.0D0
IF(R1.LT.5.0D0)TTB=TTB/2.0D0
T11=T11+TT
T11B=T11B+TTB
IF((TTB.LT.ERR3*T11B).AND.(TT.LT.ERR3*T11))GO TO 8635
8630 CONTINUE
WRITE(6,863)
863 FORMAT(/,' PHI(E,R) INTEGRAL NOT CONVERGED ON DOWNWARD ',
1 'SIDE. ABORT.')
STOP
8635 CONTINUE
T22=0.0D0
T22B=0.0D0
DO 8670 I8670=1,150
SUMCH2=0.0D0
KEJT=0.0D0
R1=EROT+(I8670-1)*ROTINC1
DO 8660 I8660=1,NCHAN
THRSH=EILO(I8660)-BSTAR(I8660)*(R1-EROT)
IF(DBLE(N*INC).LE.THRSH)GO TO 8660
THRCP=EISM(I8660)-BSTAR(I8660)*(R1-EROT)
DEC11=DBLE(N*INC)-THRCP
IF(DEC11.LE.0.0D0)DEC11=1.0D0
ISUM1=INT(DEC11/INC)
REM1=MOD(DEC11,DBLE(INC))/DBLE(INC)
IF(ISUM1.LT.1)GO TO 8640
WD=((1.0D0-REM1)*TI2(ISUM1,I8660,1)
1 +REM1*TI2(ISUM1+1,I8660,1))/DELTAE
GO TO 8650
8640 WD=REM1*TI2(1,I8660,1)/DELTAE
8650 SUMCH2=SUMCH2+WD
8660 CONTINUE
KEJT=SUMCH2/(H*RHOMOL(N))
TT=EXP((EROT-R1)/DELTA)/(OMEGA+KEJT)
TTB=EXP((EROT-R1)/DELTA)
T22=T22+TT
T22B=T22B+TTB
IF((TTB.LT.ERR3*T22B).AND.(TT.LT.ERR3*T22))GO TO 8680
8670 CONTINUE
WRITE(6,867)EROT,N
867 FORMAT(/,' PHI(E,R) INTEGRAL NOT CONVERGED ON UPWARD SIDE.'
1 ,/,' ABORT AT ROTATIONAL ENERGY =',E10.3,', N =',I4)
STOP
8680 CONTINUE
PHIL=(T11+T22)/(T11B+T22B)
RETURN
END
C
C
SUBROUTINE CENT(IN,EROT,JEJ,XXMU,BMOLEC,EOK1,EILOW,ECB)
C
C FINDS CENTRIFUGAL BARRIER.
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DOUBLE PRECISION MORSE
PARAMETER (NMAX8=8000, NMAX5=5000, NMAX80=80000)
LOGICAL BINPUT
C
COMMON /ALL/ R,RT,H,WKA
COMMON /BFIT/ RR(3,30),BB(3,30),AA1(3),AA2(3)
COMMON /BINTEG/ NB(3)
COMMON /BLOGIC/ BINPUT
COMMON /POT/ D,BE,REQ1
C
DATA SINC/0.1D0/,NTOT/200/
C
ECB=0.0D0
IF(JEJ.EQ.0)GO TO 8080
RSTART=REQ1+LOG(0.2D+1)/BE
DO 8060 I8060=1,NTOT
RCOM=RSTART+DFLT(I8060)*SINC
IF(BINPUT)THEN
BT=AA1(IN)*RCOM+AA2(IN)
ELSE
BT=16.8477D0/(XXMU*(RCOM**2))
ENDIF
EFFV=WKA*MORSE(RCOM)+EROT*(BT/BMOLEC)
C
C TEST TO ENSURE A GENUINE MAXIMUM IN EFFV (ADDED CA. 20/3/89)
C
IF(ECB.LT.EFFV)GO TO 8055
RCOM1=RCOM+0.7D0
IF(BINPUT)THEN
BT1=AA1(IN)*RCOM1+AA2(IN)
ELSE
BT1=16.8477D0/(XXMU*(RCOM1**2))
ENDIF
EFFV1=WKA*MORSE(RCOM1)+EROT*(BT1/BMOLEC)
IF(EFFV1.GT.EFFV)GO TO 8055
GO TO 8070
8055 ECB=EFFV
8060 CONTINUE
WRITE(6,806)EROT
806 FORMAT(/,' WARNING : NO CENTRIFUGAL BARRIER FOUND FOR EROT',
1 ' =',F6.2,' CM-1',/,' ABORT')
STOP
C
C DETERMINE THE (J-DEPENDENT) THRESHOLD ENERGY
C FOR THIS VALUE OF THE ROTATIONAL ENERGY.THIS IS EVALUATED
C AT THE POSITION OF THE CENTRIFUGAL BARRIER.
C
8070 RCB=RCOM-SINC
EILOW=ECB-EROT
GO TO 8090
C
C EO(J)=EO WHEN J=0.
C
8080 EILOW=EOK1*WKA
ECB=EOK1*WKA
8090 CONTINUE
RETURN
END
C
SUBROUTINE CENDIP(EROT,XXMU,BMOLEC,DHO1,EILOW,ECB,
1 DIP1,POLZ1,UNSTJ,REQ1)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
PARAMETER (NMAX8=8000, NMAX5=5000, NMAX80=80000)
LOGICAL UNSTJ
C
COMMON /ALL/ R,RT,H,WKA
C
DATA NTOT/116/,SINC/0.25D0/
C
ECB=-1.0D+6
C
C EVALUATE EFFECTIVE (R**-2) CONSTANT.
C
C=16.8477D0*EROT/(BMOLEC*XXMU)-2.4185D+4*DIP1
IF(C.LE.0.0D0)GO TO 8130
RSTART=30.0D0
DO 8100 I8100=1,NTOT
II=I8100-1
RCOM=RSTART-DBLE(II)*SINC
IF(RCOM.LT.(REQ1+0.1D0))GO TO 8110
IF(RCOM.GE.(REQ1+3.0D0))BT=16.8477D0/(XXMU*(RCOM**2))
IF(RCOM.GE.(REQ1+3.0D0))BT1=BT
IF(RCOM.LT.(REQ1+3.0D0))BT=BMOLEC*(1.0D0-(RCOM-REQ1)/3.0D0)
1 +BT1*(RCOM-REQ1)/3.0D0
RONREQ=RCOM/REQ1
EFFV=WKA*POTDIP(RONREQ,DHO1,DIP1,POLZ1,REQ1)+EROT*(BT/BMOLEC)
IF(ECB.GT.EFFV)GO TO 8120
ECB=EFFV
8100 CONTINUE
8110 UNSTJ=.TRUE.
RETURN
C
C DETERMINE THE (J-DEPENDENT) THRESHOLD ENERGY
C FOR THIS VALUE OF THE ROTATIONAL ENERGY.THIS IS
C EVALUATED AT THE CENTRIFUGAL BARRIER.
C
8120 RCB=RCOM+SINC
EILOW=ECB-EROT
GO TO 8140
C
C E0(J)=E0 WHEN C<0.
C
8130 EILOW=DHO1*WKA-EROT
ECB=DHO1*WKA
RCB=30.0D0
8140 IF(RCB.LT.1.5D0*REQ1)UNSTJ=.TRUE.
CONTINUE
RETURN
END
C
C
SUBROUTINE BESTFIT(NB1,RR1,BB1,A1,A2)
C
C THIS SUBROUTINE FITS A STRAIGHT LINE OF BEST FIT TO THE SET OF
C POINTS {RR1,BB1}
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DOUBLE PRECISION NUM1,NUM2
PARAMETER (NMAX8=8000, NMAX5=5000, NMAX80=80000)
DIMENSION RR1(30),BB1(30)
DATA TEST/0.9D0/
C
C INITIALISE SUMS
C
SUMR=0.0D0
SUMB=0.0D0
SUMRR=0.0D0
SUMBB=0.0D0
SUMRB=0.0D0
C
DO I=1,NB1
SUMR=SUMR+RR1(I)
SUMB=SUMB+BB1(I)
SUMRR=SUMRR+RR1(I)*RR1(I)
SUMBB=SUMBB+BB1(I)*BB1(I)
SUMRB=SUMRB+RR1(I)*BB1(I)
ENDDO
C
C CALCULATE NUMERATORS AND DENOMINATORS
C
DEN=NB1*SUMRR-SUMR*SUMR
NUM1=NB1*SUMRB-SUMR*SUMB
NUM2=SUMB*SUMRR-SUMR*SUMRB
C
C CALCULATE SLOPE, A1 AND INTERCEPT, A2
C
A1=NUM1/DEN
A2=NUM2/DEN
C
C CALCULATE CORRELATION COEFFICIENT, R
C
DENR=SQRT(DEN*(NB1*SUMBB-SUMB*SUMB))
R=NUM1/DENR
C
C CHECK CORRELATION OF LINE
C
ABSR=ABS(R)
IF(ABSR.LT.TEST)THEN
WRITE(6,10)R,TEST
STOP
ELSE
WRITE(6,20)A1,A2,R
ENDIF
10 FORMAT(/,1X,'STRAIGHT LINE CORRELATION BETWEEN R AND B IS WEAK'/,
* 1X,'CORRELATION = ',F6.4,/,
* 1X,'OPTIONS: CHECK VALUES INPUTTED',/,
* 1X,' ALTER VALUE OF TEST IN SUBROUTINE BESTFIT',/,
* 1X,' (CURRENTLY ',F5.3,')',/,
* 1X,' USE A HIGHER ORDER POLYNOMIAL FIT',/,
* 1X,'ABORT')
20 FORMAT(/,1X,'STRAIGHT LINE FIT FOR B AS A FUNCTION OF R',/,
* 1X,' SLOPE = ',F10.6,/,
* 1X,' INTERCEPT = ',F10.6,/,
* 1X,'CORRELATION = ',F6.4,/)
RETURN
END
C
C
BLOCK DATA
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DOUBLE PRECISION K
PARAMETER (NMAX8=8000, NMAX5=5000, NMAX80=80000, NMAX24=24000)
C NMAX24 SHOULD BE 3 TIMES NMAX8
C
C CORRECTION BY GILBERT NOV 8 1992, BY JORDAN 23 NOV 92
C
COMMON /ALL/ R,RT,H,WKA
COMMON /BIG1/ K(3,NMAX8),AK(NMAX8),AK2(NMAX8),ROTHR(NMAX5)
COMMON /BS/ FRE(200),TI(NMAX80),TII(NMAX80),DELTAE
COMMON /JINT/ INCHAN,JRECOM,NJC
COMMON /JMICRO/ RHOMOL(NMAX8),ROTH(NMAX5),WE(NMAX8,3)
COMMON /JPARTF/ PF(NMAX8,3),PFC(3),PFM(NMAX8)
COMMON /JPOT/ BETA(3),DHO(3),DIP(3),EOK(3),POLZ(3)
COMMON /JROTC/ BCMPLX(3,1),RCPL(3,1),SRC(3,1)
COMMON /JROTM/ BMOLEC,BRATIO(3),REQ(3),SRM
COMMON /STOR/ BVEC(3,1,20),SIGVC(3,1,20)
COMMON /ROTF/ GAMON2(40),ROTTM(3),V01,OMEGA1,OMEGA2
C
DATA R/0.694944D+0/,H/0.333566D-10/,WKA/349.689D0/,
* K/NMAX24*0.0D0/,AK2/NMAX8*0.0D0/,ROTHR/NMAX5*0.0D+0/,
* TI/NMAX80*0.0D+0/,TII/NMAX80*0.0D+0/,DELTAE/10.0D+0/,
* INCHAN/1/,
* RHOMOL/NMAX8*0.0D+0/,
* PF/NMAX24*0.0D+0/,PFM/NMAX8*0.0D0/,
* DIP/3*0.0D0/,EOK/3*0.0D0/,POLZ/3*0.0D0/,
* SRC/3*1.0D0/,
* SRM/1.0D0/,
* BVEC/60*0.0D0/,SIGVC/60*0.0D0/,
* ROTTM/1.772454D+0,1.0D+0,1.772454D+0/
END
C
C
DOUBLE PRECISION FUNCTION DIFF(BETA1)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DOUBLE PRECISION MORSE
C
COMMON /VFIT/ VCH1(30),RVCH1(30)
COMMON /INTVFIT/ NV1
COMMON /POT/ D,BE,REQ1
C
BE=BETA1
DIFFX=0.0D0
DO 8050 I8050=1,NV1
VTEMP=MORSE(RVCH1(I8050))
DIFFX=DIFFX+((VTEMP-VCH1(I8050))**2)
8050 CONTINUE
DIFF=DIFFX
RETURN
END
C
C
DOUBLE PRECISION FUNCTION DFLT(I)
DFLT = DBLE(FLOAT(I))
RETURN
END
C
C
DOUBLE PRECISION FUNCTION MORSE(RVAL)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
COMMON /POT/ D,BE,REQ1
C
TR1=1.0D0/EXP((RVAL-REQ1)*BE)
MORSE=((1.0D0-TR1)**2)*D
RETURN
END
C
C
DOUBLE PRECISION FUNCTION POTDIP(RVAL,DHO1,DIP1,POLZ1,REQ1)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DD=69.1235D0*DIP1/(REQ1**2)
CC=1.6609D+26*POLZ1/(REQ1**4)
BB=2.0D0*DHO1-(4.0D0*CC+5.0D0*DD)/3.0D0
AA=DHO1-CC/3.0D0-2.0D0*DD/3.0D0
POTD=DD/(RVAL**2)
POTIND=CC/(RVAL**4)
POTCH=AA/(RVAL**12)-BB/(RVAL**6)
POTDIP=DHO1+(POTCH-POTD-POTIND)
RETURN
END
C
C
C
SUBROUTINE JAVRGE(RHO,DELTAE,M,NI,ISTEP,INC)
C
C SETS UP J-AVERAGED K(E).
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C
INTEGER NE0(3),NLO(3),ITP1(3)
DOUBLE PRECISION KEJ,MORSE
C
PARAMETER (NMAX8=8000, NMAX5=5000, NMAX80=80000)
DIMENSION BR1(3),BR(3),BSTAR(3),
* ECNB(3),EOCPLX(3),EILO(3),EISM(3),ITA(NMAX8),
* PHI(NMAX8),RHO(NMAX8),SCLX(3),SIM(3),XKH(3)
C
LOGICAL HPC(3),HPCONV,CONV,TRUNC,UNSTJ,
* JCBS,WKJ1,IHINDX1,IONX1,JAVX,STERICX
C
COMMON /ALL/ R,RT,H,WKA
COMMON /BLOK21/ WS1(NMAX8,3),BOTLN(NMAX8)
COMMON /JCOLL/ DELTA,GAMMA1,OMEGA
COMMON /JHIGHP/ XKHP(NMAX8,3)
COMMON /JINT/ INCHAN,JRECOM,NJC
COMMON /JLIFE/ PLIFE,PPL1,PPL2,RADST,RSTAB,UNIR(3)
COMMON /JLOGIC/ JCBS(3),WKJ1
COMMON /JMICRO/ RHOMOL(NMAX8),ROTH(NMAX5),WE(NMAX8,3)
COMMON /JPARAM/ ROTINC1,ERR1
COMMON /JPARTF/ PF(NMAX8,3),PFC(3),PFM(NMAX8)
COMMON /JPOT/ BETA(3),DHO(3),DIP(3),EOK(3),POLZ(3)
COMMON /JRATES/ SCLE,SCLOW(2),RTHI(3),RTHI2(3),FWR(3)
COMMON /JROTC/ BCMPLX(3,1),RCPL(3,1),SRC(3,1)
COMMON /JROTM/ BMOLEC,BRATIO(3),REQ(3),SRM
COMMON /JSUM/ TI2(NMAX8,3,1)
COMMON /LOOP/ IN,NCHAN
COMMON /LOGVAR/ IONX1,IHINDX1,JAVX,STERICX
COMMON /POT/ D,BE,REQ1
C
ZERO=0.0D0
C
IDC=0
ION=0
IF(IONX1)ION=5
ISUML=0
ITPLST=0
NE0MX=0
IF(IONX1)NJC=1
NMAXT=NMAX8
PLIFE=ZERO
PPL1=ZERO
PPL2=ZERO
TEMP1=ZERO
TEMP2=ZERO
TEMP3=ZERO
THPSUM=ZERO
TSCSUM=ZERO
UNSTJ=.FALSE.
DO N=1,NMAX8
BOTLN(N)=ZERO
ITA(N)=0
DO I=1,NCHAN
WE(N,I)=ZERO
XKHP(N,I)=ZERO
ENDDO
ENDDO
DO N=1,NMAX5
ROTH(N)=ZERO
ENDDO
DO I=1,NCHAN
BR(I)=ZERO
BSTAR(I)=1.0D0-BCMPLX(I,1)/BMOLEC
FWR(I)=ZERO
HPC(I)=.FALSE.
NE0(I)=NINT(EOK(I)*WKA/DFLT(INC))
NE0MX=MAX0(NE0MX,NE0(I))
RTHI(I)=ZERO
RTHI2(I)=ZERO
SCLX(I)=ZERO
ENDDO
C
C SET ROTATIONAL ENERGY GRAIN SIZE AND MAXIMUM.
C
JMAX=6000
C
C original default was rotinc = 50
C rotinc now 10 if inc=10 (gpk 7/94)
C
IF((.NOT.IHINDX1) .AND. (INC .EQ. 10)) ROTINC1=10.0D0
IF((.NOT.IHINDX1) .AND. (INC .NE. 10)) ROTINC1=50.0D0
RSPACE=ROTINC1
JSTEP=INT((ROTINC1+1.0D0)/10.0D0)
C
C SET CONVERGENCE PARAMETER FOR INTEGRAL OF W(E,J) OF COMPLEX
C OVER J & SET CONVERGENCE PARAMETERS TO FALSE.
C
ERR1=0.250D-02
HPCONV=.FALSE.
TRUNC=.FALSE.
JCOUNT=0
DO 7210 IEJ=1,JMAX,JSTEP
JCOUNT=JCOUNT+1
JEJ=IEJ-1
EROT=JEJ*DELTAE
C
C EROT IS THE ROTATIONAL ENERGY OF THE MOLECULE
C
C FOR EACH ROTATIONAL ENERGY CONSIDERED LOCATE THE POSITION OF
C THE CENTRIFUGAL BARRIER AND THE EFFECTIVE THRESHOLD ENERGY
C FOR THE REACTION.
C
BR2=ZERO
EILOM=1.0D+7
PLIF1=ZERO
PP1=ZERO
PP2=ZERO
DO N=1,NMAX8
PHI(N)=1.0D0
IF(WKJ1)PHI(N)=ZERO
DO I=1,NCHAN
WS1(N,I)=ZERO
ENDDO
ENDDO
DO 7020 I=1,NCHAN
ECNB(I)=ZERO
EILO(I)=ZERO
EISM(I)=ZERO
EOCPLX(I)=ZERO
IF(JCBS(I))GO TO 7010
BE=BETA(I)
D=EOK(I)
REQ1=REQ(I)
XXM=16.8477D0/(BCMPLX(I,1)*RCPL(I,1)**2)
IF(.NOT.IONX1)GO TO 7000
DHO1=DHO(I)
DIP1=DIP(I)
POLZ1=POLZ(I)
RCPL1=RCPL(I,1)/REQ1
C
C CENDIP & CENT CALCULATE THE CENTRIFUGAL BARRIER.
C EILOW = EFFECTIVE THRESHOLD ENERGY &
C ECB = ENERGY AT THE CENTRIFUGAL BARRIER
C
CALL CENDIP(EROT,XXM,BMOLEC,DHO1,EILOW,ECB,
* DIP1,POLZ1,UNSTJ,REQ1)
UNSTJ=EILOW.LE.0.0D0
EILO(I)=EILOW
EILOM=DMIN1(EILOM,EILO(I))
ECNB(I)=ECB
EOCPLX(I)=POTDIP(RCPL1,DHO1,DIP1,POLZ1,REQ1)
GO TO 7010
7000 EOK1=EOK(I)
CALL CENT(IN,EROT,JEJ,XXM,BMOLEC,EOK1,EILOW,ECB)
EILO(I)=EILOW
ECNB(I)=ECB
EOCPLX(I)=MORSE(RCPL(I,1))
7010 IF(UNSTJ)GO TO 7020
DEJ=BSTAR(I)*EROT
IF(JCBS(I))EOCPLX(I)=EOK(I)
EISM(I)=EOCPLX(I)*WKA-DEJ
C
C CHECK THAT BARRIER IS NOT NEGATIVE ADDED 5.1.93 MJTJ
C
IF(EILO(I).LT.0.0D0)EILO(I)=ZERO
IF(EISM(I).LT.0.0D0)EISM(I)=ZERO
IF(JCBS(I))EILO(I)=EISM(I)
NLO(I)=INT(EILO(I)/DFLT(INC))
C
C DEJ = ROT. ENERGY OF MOLECULE - ROT. ENERGY OF COMPLEX
C EOCPLX = J=0 ENERGY AT COMPLEX
C EISM = DIFFERENCE IN EFFECTIVE POTENTIAL BETWEEN COMPLEX & MOLECULE
C
C IF CHANNEL I HAS A CHEMICAL BARRIER THEN THE CENTRIFUGAL BARRIER
C EFFECT IS IGNORED AND EILOW IS EVALUATED AT THE TRANSITION STATE
C THUS EILO = EISM
C
7020 CONTINUE
C
C EVALUATE TERMS IN THE HIGH PRESSURE LIMIT EXPRESSION
C
C VEFCPL = THE EFFECTIVE POTENTIAL AT THE COMPLEX
C RTHI = NUMERATOR FOR HIGH PRESSURE RATE / RSPACE
C HPC & HPCONV = HIGH PRESSURE CONVERGENCE PARAMETERS
C
DO 7050 I=1,NCHAN
IF(UNSTJ)GO TO 7220
IF(JCBS(I))GO TO 7050
VEFCPL=EOCPLX(I)*WKA+EROT*(BCMPLX(I,1)/BMOLEC)
DLEM=EILO(I)-EISM(I)
IDLEM=INT(DLEM/INC)
REM=MOD(DLEM,DFLT(INC))/DFLT(INC)
C
C EXPRESSION FOR RAD2 CORRECTED 23.11.92 MJTJ
C LIMIT ON VALUE OF IDLEM CORRECTED 5.1.93 MJTJ
C
IF(IDLEM.GE.1)GO TO 7030
RAD1=REM*TI2(1,I,1)/(DELTAE*EXP(ECNB(I)/RT))
RAD2=((1.0D0-REM)*PFC(I)+REM*PF(1,I))/EXP(VEFCPL/RT)
GO TO 7040
7030 RAD1=((1.0D0-REM)*TI2(IDLEM,I,1)+REM*TI2(IDLEM+1,I,1))
* /(DELTAE*EXP(ECNB(I)/RT))
RAD2=((1.0D0-REM)*PF(IDLEM,I)+REM*PF(IDLEM+1,I))
* /EXP(VEFCPL/RT)
7040 HPC(I)=(RAD1+RAD2).LT.(ERR1*RTHI(I))
C
C USE TRAPEZOIDAL RULE TO INTEGRATE OVER J TAKING HALF OF THE FIRST
C TERM MJTJ 6.1.93
C
RTHI(I)=RTHI(I)+RAD1+RAD2
FWR(I)=FWR(I)+1.0D0/EXP((ECNB(I)-EOK(I)*WKA)/RT)
IF(JCOUNT.EQ.1)THEN
RTHI(I)=0.5D0*RTHI(I)
FWR(I)=0.5D0*FWR(I)
ENDIF
7050 CONTINUE
HPCONV=.TRUE.
DO I=1,NCHAN
IF(JCBS(I))HPC(I)=.TRUE.
HPCONV=HPC(I).AND.HPCONV
ENDDO
C
C EVALUATE TERM FOR STRONG COLLISION (E,J) LOW PRESSURE LIMITING
C RATE COEFFICIENT SCLX.
C
IIX=1
IF(IONX1.AND.(JRECOM.EQ.1))IIX=2
DO 7075 JJ=1,IIX
IF(JJ.EQ.1)ELO1=EILO(1)
IF(JJ.EQ.2)ELO1=EILO(INCHAN)
NLO1=INT(ELO1/DFLT(INC))
REM=MOD(ELO1,DFLT(INC))/DFLT(INC)
C
C EXPRESSION FOR SCLX FOR NLO1.LT.1 CORRECTED 5.1.93 MJTJ
C
IF(NLO1.GE.1)GO TO 7060
SCLX(JJ)=SCLX(JJ)+REM*PFM(1)/EXP(EROT/RT)
GO TO 7070
7060 SCLX(JJ)=SCLX(JJ)+((1.0D0-REM)*PFM(NLO1)
* +REM*PFM(NLO1+1))/EXP(EROT/RT)
7070 CONTINUE
C
C USE TRAPEZOIDAL RULE TO INTEGRAT SCLX. MJTJ 8.1.93
C
IF(JCOUNT.EQ.1)SCLX(JJ)=0.5D0*SCLX(JJ)
7075 CONTINUE
NRCTM=nmax80
DO I=1,NCHAN
NRCTM=MIN0(NRCTM,NINT(EISM(I)/INC))
ENDDO
C
C IE NEED TO ADD EISM TO EOCPLX TO REACH THE THRESHOLD ENERGY
C
NMAXT=MIN0(NMAXT,(NRCTM+NI))
IF(NMAXT.GT.(NE0MX+100))GO TO 7080
WRITE(6,7)
C
7 FORMAT(' NUMBER OF ENERGY INCREMENTS (NN) IS NOT LARGE',
1 ' ENOUGH : INCREASE BY 50.',//,' ABORT')
C
STOP
7080 CONTINUE
IF(.NOT.TRUNC)NMAX=NMAXT
C
C EVALUATE W(E,J) OF THE COMPLEX, WS1
C
ITPM=nmax80
DO I=1,NCHAN
ITP1(I)=NLO(I)+1
ITPM=MIN0(ITPM,ITP1(I))
DEL1=ITP1(I)*DFLT(INC)-EISM(I)
IF(DEL1.LE.0.0D0)DEL1=1.0D0
ISUM=INT(DEL1/INC)-1
REM=MOD(DEL1,DFLT(INC))/DFLT(INC)
DO 7110 N=ITP1(I),NMAX
ISUM=ISUM+1
IF(ISUM.GE.1)GO TO 7090
WSTEMP=REM*TI2(1,I,1)/DELTAE
WS1(N,I)=WSTEMP*SRM/SRC(I,1)
GO TO 7100
7090 WSTEMP=((1.0D0-REM)*TI2(ISUM,I,1)+
* REM*TI2(ISUM+1,I,1))/DELTAE
WS1(N,I)=WSTEMP*SRM/SRC(I,1)
7100 CONTINUE
IF(.NOT.IONX1)GO TO 7110
IF((WS1(N,I).LE.0.0D0).OR.
* (WSTEMP.GT.WS1(N,I)))WS1(N,I)=WSTEMP
7110 CONTINUE
ENDDO
IF(ITPM.GT.NE0(1))THEN
ROTH(ITPM)=ZERO
ELSEIF(ITPLST.EQ.ITPM)THEN
ROTH(ITPM)=(ROTH(ITPM)+EROT)*0.5D0
ELSE
ROTH(ITPM)=EROT
ENDIF
ITPLST=ITPM
C
C ITPM = ENERGY INCREMENT EQUIVALENT TO THE ENERGY THRESHOLD
C & NE0 = ENERGY INCREMENT EQUIVALENT TO THE J=O DISSOCIATION
C ENERGY
C
IF(.NOT.WKJ1)GO TO 7120
C
C CALCULATE THE CORRECTION TERM IF WEAK ROTATIONAL RELAXATION
C IS USED
C
DO N=1,NMAX
CALL WEAKJ(BSTAR,EILO,EISM,PHIL,EROT,INC,DELTAE,N)
PHI(N)=PHIL
ENDDO
7120 CONTINUE
CONV=.FALSE.
C
C PROGRESSIVELY EVALUATE (USING COMPOSITE SIMPSON RULE)
C THE INTEGRAL OF W(E,J) OF THE COMPLEX OVER A STRONG
C COLLISION DISTRIBUTION OF ROTATIONAL ENERGIES.
C
C IE: USE SIMPSON'S RULE TO INTEGRATE K(E) OVER J
C
C FOR NEUTRAL/NEUTRAL WE START WITH NJC=0.5*NE0
C FOR ION/DIPOLE WE START WITH NJC=1
C
IF(MOD(JCOUNT,2).EQ.0)THEN
NSP=4
ELSEIF(JCOUNT.EQ.1)THEN
NSP=1
ELSE
NSP=2
ENDIF
C
C XKHP IS FOR CALCULATING HIGH PRESSURE LIMIT OF K(E).
C
CONV=.FALSE.
FEXP=DFLT(INC)/RT
TERM1=DEXP(-EROT/RT)
DO 7150 N=NJC,NMAX
SUMCH=ZERO
DO I=1,NCHAN
KEJ=WS1(N,I)/(H*RHOMOL(N))
SUMCH=SUMCH+KEJ
ENDDO
TM2=1.0D0/(OMEGA+SUMCH)
TERM=TERM1*TM2
C
C SET UP "TEMPORARY" VARIABLES:
C * SIM IS USED FOR NUMERATOR OF INTEGRAL
C * XKH IS USED FOR INTEGRAL HIGH PRESSURE RATE
C * BL IS USED FOR DENOMINATOR OF INTEGRAL
C
DO I=1,NCHAN
SIM(I)=NSP*WS1(N,I)*TERM
XKH(I)=NSP*WS1(N,I)*TERM1
ENDDO
BL=NSP*TERM
C
C IF HIGH PRESSURE RATE AT ENERGY N HAS ALREADY CONVERGED SKIP TO
C THE NEXT ENERGY
C
IF(ITA(N).EQ.1)GO TO 7150
C
C CHECK CONVERGENCE ON HIGH PRESSURE RATE T ENERGY N.
C CHECK IS ONLY MADE ON CHANNEL 1 AS CRITICAL ENERGIES WERE INPUT
C IN DESCENDING ORDER
C
C NB WE ARE CHECKING EACH INDIVIDUAL TERM'S CONTRIBUTION TO THE
C OVERALL INTEGRAL RATHER THAN THE CONTRIBUTION OF EACH INTERVAL
C THE EXTRA FACTOR OF 4 TIGHTENS THE CONVERGENCE CRITERIA
C
IF((4.0D0*XKH(1)).GE.(ERR1*XKHP(N,1)))GO TO 7130
ITA(N)=1
C
C TEST CONVERGENCE FOR TOTAL HIGH PRESSURE RATE: THPSUM
C WRT HIGH PRESSURE RATE AT ENERGY N: THP
C AND TOTAL STRONG COLLISION RATE: TSCSUM
C WRT STRONG COLLISION RATE AT ENERGY N: TSCG
C AGAIN CHECK IS ONLY MADE FOR CHANNEL 1
C
THP=XKHP(N,1)*RSPACE/(3.0D0*RT)
THP=THP/(H*DEXP(DFLT(N)*FEXP))
TSC=WE(N,1)/BOTLN(N)
TSC=TSC/(H*RHOMOL(N))
G=EXP(LOG(RHOMOL(N))-DFLT(N)*FEXP)
G=G*OMEGA/(OMEGA+TSC)
TSCG=TSC*G
C
C CONVERGENCE CRITERIA:
C * TOTAL HIGH PRESSURE RATE IS CONVERGED AT ENERGY N
C * STRONG COLLISION RATE IS CONVERGED AT ENERGY N
C * ENERGY N IS LESS THAN NE0, IE THRESHOLD IS LESS THAN E0
C (E0+5 IN THE CASE OF ION-DIPOLE CALCULATIONS)
C
CONV=(THP.LT.ERR1*THPSUM)
CONV=CONV.AND.(TSCG.LT.ERR1*TSCSUM)
CONV=CONV.AND.(N.LE.NE0(1)+ION)
IF(CONV.AND.HPCONV)NJC=MIN0(NE0(1),N-1)
THPSUM=THPSUM+THP
TSCSUM=TSCSUM+TSCG
C
C CALCULATE J-AVERAGED QUANTITIES AT ENERGY N:
C * SUM OF STATES WE(N,I)
C * HIGH PRESSURE RATE XKHP(N,I)
C
7130 CONTINUE
DO I=1,NCHAN
WE(N,I)=WE(N,I)+SIM(I)
XKHP(N,I)=XKHP(N,I)+XKH(I)
C
C RESET TEMPORARY PARAMETERS USED IN THE INTEGRATION
C
SIM(I)=ZERO
XKH(I)=ZERO
ENDDO
C
C CALCULATE DENOMINATOR FOR
C
BOTLN(N)=BOTLN(N)+BL
BL=ZERO
C
IF(.NOT.IONX1)GO TO 7150
C
C CALCULATE BRANCHING RATIOS & LIFETIMES USING STRAIGHT
C SUM OVER E AND J
C
IF(N.LE.ITP1(INCHAN))GO TO 7150
BR2=BR2+WS1(N,INCHAN)*EXP(-DBLE(N)*FEXP)
PTEMP=WS1(N,INCHAN)*EXP(-DBLE(N)*FEXP)
* /((RADST+SUMCH)*H)
PP1=PP1+PTEMP
PLIF1=PLIF1+PTEMP/SUMCH
PP2=PP2+RHOMOL(N)*EXP(-DBLE(N)*FEXP)
DO 7140 I=1,NCHAN
IF(I.EQ.INCHAN)GO TO 7140
IF(N.LE.ITP1(I))GO TO 7140
FRR=WS1(N,I)/((RADST+SUMCH)*H*RHOMOL(N))
BR1(I)=BR1(I)+FRR*WS1(N,INCHAN)*EXP(-DBLE(N)*FEXP)
7140 CONTINUE
7150 CONTINUE
C
C CONDITIONS FOR A REACTION WITH A CHEMICAL BARRIER
C IF ALL CHANNELS HAVE A BARRIER CHECK IS ON CHANNEL 1
C ELSE CHECK IS ON CHANNEL WITHOUT BARRIER WITH HIGHEST E0
C
IF(.NOT.JCBS(1))GO TO 7170
NX=1
DO I=2,3
NX=NX+1
IF(.NOT.JCBS(I))GO TO 7160
ENDDO
IF(BSTAR(1).GE.0.2D0)NNM=35
IF(BSTAR(1).LT.0.2D0)NNM=5
IF(BSTAR(1).LE.0.0D0)NNM=NE0(1)-NLO(1)
IF(ITA(NE0(1)-NNM).EQ.1)CONV=.TRUE.
IF(CONV)NJC=NE0(1)-NNM-1
GO TO 7170
C
C CHANNEL NX DOES NOT HAVE A CHEMICAL BARRIER
C
7160 IF(EILO(NX)/WKA.LE.(EOK(NX)-12.0D0))CONV=.TRUE.
IF(CONV)NJC=INT((EOK(NX)-12.0D0)*WKA/dflt(inc))+1
7170 CONTINUE
IF(CONV.AND.HPCONV)GO TO 7220
C
C CALCULATE BRANCHING RATIOS AND LIFETIMES FOR THE ION-DIPOLE
C CALCULATIONS:
C
IF(.NOT.IONX1)GO TO 7190
DO 7180 I=1,NCHAN
IF(I.EQ.INCHAN)GO TO 7180
BR(I)=BR(I)+BR1(I)*EXP(-EROT/RT)
7180 CONTINUE
BR2T=BR2T+BR2*EXP(-EROT/RT)
PPL1=PPL1+PP1*DBLE(INC)*EXP(-EROT/RT)
PPL2=PPL2+PP2*DBLE(INC)*EXP(-EROT/RT)
PLIFE=PLIFE+PLIF1*DBLE(INC)*EXP(-EROT/RT)
7190 CONTINUE
IF(TRUNC)GO TO 7210
C
C IF ALL W(E,J) INTEGRALS ABOVE E0 ARE CONVERGED,THEN
C LEAVE OUT FOR HIGHER EROT.
C
TRUNC=.TRUE.
DO 7200 N=NE0(1),NMAX
IF(ITA(N).EQ.1)GO TO 7200
TRUNC=.FALSE.
7200 CONTINUE
IF(.NOT.TRUNC)GO TO 7210
NMAX=NE0(1)+1
7210 CONTINUE
C
C AVERAGE OVER J HAS NOT CONVERGED
C
WRITE(6,77)ERR1,JMAX
C
77 FORMAT(/,' K(E,J) AVERAGE OVER STRONG COLLISION '
* ,'J DISTRIBUTION NOT CONVERGED.',/,' EITHER: 1.INCREASE',
* ' ERROR TOLERANCE PARAMETER ERR1',/,
* ' (CURRENT VALUE ',F10.4,')',
* /,' OR : 2.INCREASE MAXIMUM ROTATIONAL ENERGY FOR',
* ' CALCULATION',/,' (CURRENT VALUE ',I5,').',/,' NOTE: EXTRA ',
* 'NUMERICAL ERROR INCURRED BY 1. WILL BE APPROXIMATELY ',
* /,' ACCOUNTED FOR BY CALIBRATION OF CALCULATED WITH EXACT',
* ' HIGH PRESSURE RATE.')
STOP
C
C AVERAGE HAS CONVERGED
C
7220 CONTINUE
IF(UNSTJ)NJC=MAX0(1,ITPM)
C
C EVALUATE W(E) [ = ] FOR THE COMPLEX.
C EVALUATE THE HIGH PRESSURE LIMIT, XKHP, OF W(E).
C
DO N=NJC,NMAXT
RFAC=1.0D0/BOTLN(N)
DO I=1,NCHAN
WE(N,I)=WE(N,I)*RFAC
XKHP(N,I)=XKHP(N,I)*RSPACE/(3.0D0*RT)
ENDDO
ENDDO
C
C EVALUATE NUMERATOR FOR HIGH PRESSURE RATE (RTHI) AND THE HIGH
C PRESSURE RATE USING THE EFFECTIVE ROTATIONAL PARTITION
C FUNCTION RATIO I+/I (RTHI2).
C
DO I=1,NCHAN
RTHI(I)=RTHI(I)*RSPACE
IF(JCBS(I))RTHI(I)=RT*BMOLEC*PFC(I)/
* (BCMPLX(I,1)*EXP(EOCPLX(I)*WKA/RT))
FWR(I)=FWR(I)*RSPACE/RT
IF(JCBS(I))FWR(I)=BMOLEC/BCMPLX(I,1)
RTHI2(I)=FWR(I)*RT*PFC(I)/EXP(EOK(I)*WKA/RT)
ENDDO
C
C BRATIO(I) IS FOR LOW DENSITY BRANCHING RATIO IN CHEMICAL ACTIVATION
C
IF(.NOT.IONX1)GO TO 7240
SUM=ZERO
DO 7230 I=1,NCHAN
IF(I.EQ.INCHAN)GO TO 7230
BRATIO(I)=BR(I)/BR2T
UNIR(I)=BR(I)*DBLE(INC)/(H*PPL1)
SUM=SUM+BR(I)
7230 CONTINUE
UNIR(INCHAN)=(BR2T-SUM)*DBLE(INC)/(H*PPL1)
RSTAB=RADST*PPL1*H/(BR2T*DBLE(INC))
7240 CONTINUE
C
C EVALUATE EXACT STRONG COLLISION (E,J) LOW PRESSURE RATE
C COEFFICIENT, SCLOW
C CORRECTED INDEX IIX 16.5.94 (MJTJ)
C
DO I=1,IIX
SCLOW(I)=SCLX(I)*RSPACE/RT
SCLOW(I)=SCLOW(I)*OMEGA
ENDDO
C
C EVALUATE STRONG COLLISION (E,J=0) LOW PRESSURE RATE COEFFICIENT.
C
IF(IONX1)GO TO 7250
NR2=INT(EOK(1)*WKA/DFLT(INC))
REM=MOD((EOK(1)*WKA),DFLT(INC))/DFLT(INC)
SCLE=(1.0D0-REM)*PFM(NR2)+REM*PFM(NR2+1)
SCLE=SCLE*OMEGA
NI=NMAXT-NJC
GO TO 7260
7250 CONTINUE
C
C EVALUATE LIFETIMES
C
NI=MIN0(NI,NMAXT-NJC)
PLIFE=PLIFE/PPL1
PPL1=OMEGA*PPL1*RSPACE/BMOLEC
PPL2=OMEGA*PPL2*RSPACE/BMOLEC
7260 CONTINUE
C
C SET W(E) TO ZERO BELOW THE LOWEST ENERGY, (NJC).
C
DO N=1,NJC
DO I=1,NCHAN
WE(N,I)=ZERO
ENDDO
ENDDO
DO N=1,NI
RHO(N)=ZERO
DO I=1,NCHAN
WE(N,I)=WE(N+NJC,I)
XKHP(N,I)=XKHP(N+NJC,I)
ENDDO
ENDDO
DO NN=NI,NMAXT
DO I=1,NCHAN
XKHP(NN,I)=ZERO
ENDDO
ENDDO
DO II=1,MIN0(NE0(1),ITPLST+1)
ROTH(II)=1.0D+6
ENDDO
RETURN
END
|
