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