C PROGRAM MASS2 C C (INPUT,OUTPUT,TAPE5=INPUT,TAPE6=OUTPUT) C C COMPUTE EXACT (NON-EQUILIBRIUM) RATES FOR INTERMEDIATE C AND (IF REQUIRED) LOW-PRESSURE REGIME BY C SOLUTION OF THE MASTER EQUATION (LARGEST EIGENVALUE ONLY). C THE GAS/GAS COLLISION PROBABILITY IS ASSUMED TO BE A FUNCTION C OF THE ENERGY DIFFERENCE OF THE INITIAL AND FINAL STATES,WITH A C PARAMETER WHICH MAY BE A FUNCTION OF THE INITIAL ENERGY. C C ANGULAR MOMENTUM CONSERVATION IN THE FALL-OFF AND LOW-PRESSURE C REGIMES IS TREATED USING THE METHOD OF SMITH AND GILBERT. C C THIS PROGRAM WILL HANDLE UP TO THREE REACTION CHANNELS. C C OPTION PARAMETER IOPTHT ENABLES ONE TO DO A FULL FALL-OFF C CALCULATION AND, IF REQUIRED, THE LOW-PRESSURE LIMIT. C DETERMINE EIGENVALUE BY NESBET METHOD. C C THE COMPLETE INPUT FILE FOR THIS PROGRAM IS PREPARED BY RRKM C C C************************************************************************* C COSMETIC IMPROVEMENTS TO PROGRAM: COMMENT LINES, IF BLOCKS AND DO LOOPS C PUTTING LOGICAL, INTEGER AND DOUBLE PRECISION VARIABLES INTO SEPARATE C COMMON BLOCKS C REMOVE SOME OF THE ARGUMENTS IN CALLS INTO COMMON BLOCKS C PUTTING EXTRA CHECKS TO AVOID TAKING LOGARITHM OF ZERO C PRINT OUT K0*[M] RATHER THAN INCORRECT K0*OMEGA C DON'T PRINT ANYTHING INTO COLUMN 1 C************************************************************************** C C ALTERATIONS BY Gary P Knight 9th September 1994. C C SINCE RRKM NEEDS AN INC VALUE OF LESS THAN 100 cm-1 FOR SOME ACCURATE C CALCULATIONS AT LOW TEMPERATUES. THIS PROGRAM WILL NOW WORK WITH C INC=100 cm-1 OR 10 cm-1. C THE VARIABLE TESTINC TAKES THE VALUES 1 AND 10 RESPECTIVELY TO C ACHIEVE THIS. C THIS NEW VERSION OF MASNEW (MAS55) IS DESIGNED TO WORK WITH ARRAYS C SET TO 10x BIGGER. (NOW PARAMETERISED) C NOTE THE RESULTING .EXE FILE WILL BE VERY LARGE. C C THE PROGRAM LOOPS HAVE BEEN TRUNCATED FOR INC=100 cm-1, USING THE C VARIABLE TESTINC TO SPEED UP THE CALCULATION. C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION KT PARAMETER (NMAX100=10000, NMAX50=5000, NMAX25=2500) COMMON /AMIN/ ID,NCUT,NILOC COMMON /BLOCK1/ NRATES,NI,NE0,NEFF COMMON /BLOCK2/ DELT,CORRF COMMON /BLOCK3/ NREACT,NBAND,NRPL1 COMMON /BLOCK6/ Q(NMAX100),H(NMAX100),RTRHB(NMAX100) COMMON /BLOCK7/ D(NMAX100) COMMON /BLK10/ AVRATE(NMAX50) COMMON /BLK12/ AVR2(NMAX50,2) COMMON /BLK13/ AROW(NMAX100),OMEGA,ANORM(NMAX100) COMMON /BLK14/ NRCTH COMMON /BLK15/ NCH2 COMMON /B2/ PROB(NMAX25) COMMON /DEOC/ EI,WE0 COMMON /DENS/ RHO(NMAX100) COMMON /EXP/ ERR1,ERR2,ERR3 COMMON /F1/ IFIRST,NSPEC,IMIN COMMON /F2/ ALMT,ALMT1,DET COMMON /JLP1/ JAVX COMMON /JLP2/ CRF2(NMAX100),EXL2(NMAX100) COMMON /JLP3/ CRF(NMAX100),EXL(NMAX100) COMMON /LOW/ LOWP COMMON /OPT/ IOPTPR COMMON /ROTTX/ ROTHR(NMAX100),ROTKT(NMAX100) COMMON /TRANS1/ RATTOT,RATE1,RATE2,RATEA(2) COMMON /TRANS2/ ITR COMMON /VERSCH/ IXV,JXV COMMON /NEWKIDS/ TESTINC C DIMENSION PR(20),R1(NMAX100),R2(NMAX100,2),R1S(NMAX100), * ALPHAV(10),TITLE(20),TEMPV(10),CORRAT(10),RATHP(3) C LOGICAL LOWP,JAVX,BRW C DATA R1/NMAX100*0./,R2/20000*0./,R1S/NMAX100*0./, * WKA/349.689/,RATHP/3*0./ C C READ IN AND PRINT OUT TITLE. C READ(5,1) TITLE C 1 FORMAT(1X,20A4) C WRITE(6,2) TITLE C 2 FORMAT(1H1,20X,20A4,/, * 10X,'MASTER EQUATION SOLUTION FOR FALL-OFF AND LOW PRESSURE', * ' REGIMES',/, * 10X,'LATEST UPDATE: Sept 1994') C C READ IN INITIAL GRAIN SIZE, NUMBER OF TIMES (+1) THAT C ONE WANTS TO HALVE THE GRAIN SIZE, STARTING AT A GRAIN SIZE (CM-1) WHICH C MUST BE 100,200,400,800,1600,... CM-1 C OR 10,20,40,80,160,... CM-1 (GPK 9/9/94) C C B IS A PARAMETER, NORMALLY BETWEEN 0.01 AND 0.9, WHICH HELPS C THE EIGENVALUE ROUTINE FROM GOING WILD; B GT 0.1 MEANS SLOWER C CONVERGENCE BUT LESS LIKELIHOOD OF NUMERICAL INSTABILITY. C A VALUE OF 0.1 GIVES BEST CONVERGENCE PROPERTIES, BUT THIS C MAY BE TESTED FOR EACH CASE BY VARYING B BETWEEN 0.01 AND 1. C READ(5,*) INC,NCHAN,INCCHK IDM=2 NCH2=MAX0(NCHAN,2)-1 C C CHANGED TO DELTS=DFLT(INCCHK) INSTEAD OF INC WHICH IS ALWAYS 400 C DELTS=DFLT(INCCHK) TESTINC=1.0 IF(INCCHK .EQ. 10) TESTINC=10. IF(INCCHK .EQ. 100) GO TO 100 IF(INCCHK .EQ. 10) GO TO 100 WRITE(6,3) C 3 FORMAT(' RRKM PROGRAM WAS NOT CALCULATED * WITH INC=10 OR 100 - ABORT') C STOP 100 CONTINUE C C THE FOLLOWING IS THE FIXED VALUE OF THE CONVERGENCE-GOVERNING PARAMETER. C B=0.999 C C READ IN ERRORS C READ(5,*) ERR1,ERR2,ERR3 WRITE(6,4) ERR1,ERR2,ERR3 C 4 FORMAT(1P,/,10X,'ERROR FOR TRUNCATION OF MATRIX',T50,E10.1, * /,10X,'EIGENVALUE TOLERANCE',T50,E10.1,/, * 10X,'ERROR FOR TRUNCATION OF P(E,E',1H',')',T50,E10.1) C WRITE(6,5) NCHAN C 5 FORMAT(/,10X,'NUMBER OF REACTION CHANNELS',T58,I2) C ITMAX=800 C C ITMAX IS THE MAXIMUM NUMBER OF NESBET ITERATIONS PERMITTED. C READ(5,*) E0 E0T=E0 WE0=E0T*WKA C C E0 IS THE CRITICAL ENERGY; THE MAXIMUM ENERGY CONSIDERED IS COMPUTED C BELOW AS THE ENERGY OF THE HIGHEST INPUT MICROSCOPIC RATE. C NALPHA IS THE NUMBER OF ALPHA VALUES,AND ALPHAV THE VECTOR OF THEIR C VALUES. C READ(5,*) NALPHA READ(5,*) (ALPHAV(I),I=1,NALPHA) C C READ IN PARAMETERS GIVING THE C FUNCTIONAL FORM OF P(E,E'). C READ(5,*) IXV,JXV C C IF THE FIRST VALUE OF ALPHA IS GT 0, AN EXPONENTIAL MODEL C FOR P(E,E') IS ASSUMED. IF FIRST ALPHA VALUE LT 0, THEN ADDITIONAL C PARAMETERS MUST BE READ IN TO SPECIFIY FUNCTIONAL FORM OF P. C READ IN PRESSURES (TORR) AT WHICH CALCS TO BE PERFORMED C READ(5,*) NP READ(5,*) (PR(I),I=1,NP) C C IF JAV=0, SINGLE FILE OF K(E)'S WITH AVERAGE (I+/I) J CORRECTION C FACTOR ASSUMED. C IF JAV.NE.0, K(E)= ASSUMED.(SEPARATE K(E) FILE FOR EACH C PRESSURE AND TEMPERATURE). C READ(5,*) PALMT C C PALMT IS THE PARAMETER SUCH THAT THE POPULATION VECTOR FOR C ALL ENERGIES LT E0*PALMT ARE GIVEN THEIR EQUILIBRIUM VALUES. C READ(5,*)JAV JAVX = JAV.NE.0 IF(JAVX) WRITE(6,7) C 7 FORMAT(10X,'FULL ANGULAR MOMENTUM CONSERVATION USED (SMITH-', * 'GILBERT METHOD)') C IF(.NOT.JAVX) WRITE(6,8) C 8 FORMAT(10X,'THIS CALCULATION DOES NOT USE FULL ANGULAR MOMENTUM', * ' CONSERVATION;',/,10X,'ANGULAR MOMENTUM CONSERVATION ONLY', * ' IN HIGH PRESSURE LIMIT') C IF((NALPHA.EQ.1).OR.(NP.EQ.1).OR.(.NOT.JAVX))GO TO 110 WRITE(6,9) C 9 FORMAT(/,' DUE TO THE STRUCTURE OF THE INPUT FILE,',/, * ' YOU MAY NOT LOOP SEVERAL ALPHAS AND SEVERAL PRESSURES ', * ' IN THE SAME RUN.',/,' MAKE EITHER NALPHA OR NP =1',/,' ABORT.') C STOP 110 CONTINUE READ(5,*)NTEMP,(TEMPV(I),I=1,NTEMP) IF(.NOT.JAVX)GO TO 120 IF(PR(1).GE.0)GO TO 120 WRITE(6,10) C 10 FORMAT(/,' OPTION TO GENERATE FALL-OFF AUTOMATICALLY IS NOT', * ' VALID WHEN USING J-AVERAGED K(E)S. ABORT') C STOP 120 CONTINUE C C READ IN COLLISION DIAMETER (ANGSTROM),COLLISION MASSES C READ(5,*) SGMA,WT1,WT2,EPS WRITE(6,11) SGMA,WT1,WT2 C 11 FORMAT(/,10X,'COLLISION DIAMETER =',T50,F10.2,' ANGSTROMS', * /,10X,'MASS OF REACTANT=',T50,F10.2,' A.M.U.' *,/,10X,'MASS OF BATH GAS =' * ,T50,F10.2,' A.M.U.',/) C IF(EPS.LE.0.) WRITE(6,12) C 12 FORMAT(10X,'HARD SPHERE COLLISION FREQUENCY USED',/) C IF(EPS.GT.0.) WRITE(6,13) EPS C 13 FORMAT(10X,'LENNARD JONES COLLISION FREQUENCY USED',/, * 10X,'WELL DEPTH =',T50,F10.1,' K',/) C AMASS=1./((1./WT1)+(1./WT2)) C C READ IN OPTIONS FOR DOING CALCULATION C THE OPTION PARAMETERS ARE AS FOLLOWS: C IOPTHT = 0, DOES LOW PRESSURE CALCULATION IN ADDITION TO C FALL-OFF CALCULATION, IOPT.NE.0, DOES FALL-OFF CALC. ONLY. C N.B., IF LOW PRESSURE CALCULATION DONE, LOWEST GRAIN SIZE IS MADE 100 CM-1. C C IOPTPR = 1 OR -1 WEIGHTS CONVERGENCE TESTS BY OVERALL RATE C (SUITABLE FOR 1-CHANNEL REACTION) C = 2 OR -2 WEIGHTS CONVEGENCE TESTS FOR BOTH CHANNELS EQUALLY C (SUITABLE FOR 2-CHANNEL SYSTEM). C READ(5,*) IOPTHT,IOPTPR IOPTPR=IABS(IOPTPR) C C READ IN NUMBER OF DENSITIES OF STATES THEN THE LIST OF DENSITIES OF STATES; C IT IS ASSUMED THAT THESE HAVE BEEN COMPUTED WITH A SPACING OF 100 CM-1. C WRITE(6,14) NP C 14 FORMAT(10X,'NO. OF PRESSURES =',T50,I10,/) C C FOR LOW PRESSURE CALCULATION IDM=2 DELTS = 200 OR 20 C IDM=2 DELTS=200./TESTINC INC=DELTS WRITE(6,15) INC C 15 FORMAT(10X,'INITIAL GRAIN SIZE =',T50,I10,' CM-1') C WRITE(6,16) NALPHA,IDM C 16 FORMAT(/,10X,'NO. OF VALUES OF ALPHA =',T50,I10,/, * 10X,'NO. OF TIMES,+1,GRAIN SIZE HALVED =',T50,I10) C BRW=IXV.LE.0.AND.JXV.LE.0 C C BRW IS LOGICAL VARIABLE SPECIFIYING THAT BIASED RANDOM WALK MODEL C USED FOR P(E,E'). C WRITE(6,17) PALMT C 17 FORMAT(10X,'CUTOFF PARAMETER FOR H(E)=1',T50,F10.2) C I=IABS(IOPTPR) WRITE(6,18) I C 18 FORMAT(/,10X,'NUMBER OF REACTION CHANNELS =',T50,I10) C READ(5,*) NDEGS IF(NDEGS.LE.NMAX100) GO TO 140 WRITE(6,20) NMAX100 C 20 FORMAT(' MORE THAN ',I6,' ELEMENTS IN RHO INPUT,ABORT') C STOP 140 CONTINUE READ(5,*) (RHO(I),I=1,NDEGS) C C RENORMALIZE RHO C DO I=1,NDEGS IF(RHO(I).GT.0.) GO TO 150 ENDDO 150 DNM=RHO(I) DO I=1,NDEGS RHO(I)=RHO(I)/DNM ENDDO IF(JAVX)GO TO 160 DO I=1,NTEMP READ(5,*)T,CORRAT(I) ENDDO NE0=INT(WE0*TESTINC/100.) GO TO 170 160 CONTINUE C C READ IN ARRAY OF ROTATIONAL THRESHOLD ENERGIES. C READ(5,*)NTHR READ(5,*)(ROTHR(I),I=1,NTHR) NE0=NTHR NEON2=(NE0/2)-1 170 CONTINUE C C INCREASE NUMBER OF PRESSURES CONSIDERED BY 1 SO THAT THE LOW PRESSURE C CALCULATION CAN BE PERFORMED THE FIRST TIME THE PRESSURE LOOP IS C ACCESSED FOR EACH TEMP AND EACH ALPHA C NP=NP+1 DO 400 ITEMPR=1,NTEMP T=TEMPV(ITEMPR) WRITE(6,23) T C 23 FORMAT(/////,10X,'TEMPERATURE =',T50,F10.1,' KELVINS',//) C KT=T/1.43879 IF(.NOT.JAVX)GO TO 180 DO I=NEON2,NE0 ROTKT(I) = ROTHR(I)/KT CRF(I)=0. IF(ROTKT(I).LE.30.) CRF(I)=DEXP(-ROTKT(I)) ENDDO DO I=NE0+1,NDEGS CRF(I)=1. ROTKT(I) = 0. ENDDO 180 CONTINUE ICHP=0 P1=PR(1) C C LOOP THROUGH ALL THE COLLISIONAL ENERGY TRANSFER PARAMETERS CONSIDERED C DO 399 IALPHA=1,NALPHA L=-1-(NP/2) BETA=1. ALPHA=ALPHAV(IALPHA) IF(.NOT.BRW) WRITE(6,24) ALPHA,IXV,JXV C 24 FORMAT(///,T30,' ALPHA =',F8.0,' CM-1', * //,10X,'PROBABILITY FUNCTION IS R(E-E',1H',')/N(E',1H', * '), WHERE',/,8X,12HR(E-E')=(X**,I2, * 10H) EXP(-X**,I2,21H), WHERE X=E-E'/ALPHA) C IPRT=0 DO 398 N=1,NP C C DO LOW PRESSURE CALCULATION FIRST TIME THROUGH PRESSURE LOOP C LOWP=N.EQ.1 IF(JAVX)GO TO 190 IF((N.NE.1).OR.(IALPHA.NE.1).OR.(ITEMPR.NE.1))GO TO 220 GO TO 200 190 CONTINUE IF((IALPHA.NE.1).OR.(N.EQ.2))GO TO 220 READ(5,*)(RATHP(I),I=1,NCHAN),CORAV,CORPF READ(5,*)STLPJ READ(5,*) E0EFF E0=E0EFF 200 CONTINUE C C READ IN NUMBER OF MICROSCOPIC RATES TO BE READ IN, AND THEN THE ACTUAL C RATES; IT IS ASSUMED THAT THESE HAVE BEEN COMPUTED WITH A SPACING OF C EITHER 100 OR 10 CM-1. C C MICROSCOPIC RATES AND DENSITIES OF STATES,FOR EITHER 10 CM-1 OR C 100 CM-1 SPACING,ARE READ IN R1 (K1), R2 (K2) AND RHO (RHO). C READ(5,*) NRATES READ(5,*) (R1(I),I=1,NRATES) DO I=1,NCH2 READ(5,*) (R2(II,I),II=1,NRATES) ENDDO WRITE(6,25) NRATES,NDEGS C 25 FORMAT(/,10X,'NUMBER OF INPUT MICROSCOPIC RATES =',T50,I10,/, * 10X,'NUMBER OF INPUT DENSITIES OF STATES =',T50,I10) C IF(NRATES.LE.NMAX50) GO TO 210 WRITE(6,26) NMAX50,NMAX50 C 26 FORMAT(' MAXIMUM OF ',I10,' INPUT RATES ONLY ALLOWED, ABORT.', * /,' INCREASE DIMENSION FROM',I10,' TO RUN.') C STOP 210 CONTINUE C C MAXIMUM ENERGY CONSIDERED CORRESPONDS WITH THE HIGHEST INPUT RATE C NRTSTM=NRATES EMAX=E0+(0.2859*DFLT(NRATES)/TESTINC) NIORIG=INT((EMAX*200./DELTS)/0.5719)-1 C C NI IS NUMBER OF LEVELS ALTOGETHER C NREACT IS NUMBER NOT REACTING C DELT IS ENERGY DIFFERENCE BETWEEN LEVELS IN CM-1 C NBAND IS THE BANDWIDTH OF THE MATRIX C EJM=EMAX*4.184 ECM=E0*4.184 WRITE(6,27) EMAX,EJM,E0,ECM C 27 FORMAT(//,10X,'MAXIMUM ENERGY CONSIDERED',T50,F10.1,' KCAL MOL-1' * ,/,T45,'=',T50,F10.1,' KJ MOL-1' * ,//,10X,'CRITICAL ENERGY',T50,F10.2,' KCAL MOL-1',/, * 45X,'=',T50,F10.2,' KJ MOL-1') C 220 CONTINUE NI=NIORIG PN=PR(MAX0(1,N-1)) NEFF=INT(E0EFF*WKA*TESTINC/100.) IF(P1.LE.0.) PN=1. L=L+1 IF(.NOT.LOWP.AND.P1.LE.0.) PN=CSET*(10.**L) OMEGA=(4.412E+7*SGMA*SGMA/DSQRT(AMASS))*PN/(DSQRT(T)) C C CALCULATE LENNARD-JONES COLLISION FREQUENCY IF REQUIRED. C IF(EPS.GT.0.) OMEGA=OMEGA/(0.636+0.567*DLOG10(T/EPS)) C C THE FOLLOWING SET OF STATEMENTS IS FOR THE SPECIAL CASE OF THE C LOW PRESSURE LIMIT CALCULATION, IE IF IOPTHT.NE.0 THEN THE LOW C PRESSURE CALCULATION IS NOT PERFORMED. C IF(IOPTHT.NE.0.AND.LOWP) GO TO 398 C C CALL SUBROUTINE TO CARRY OUT PREPARATION FOR LOW-PRESSURE CALCULATION. C IF(JAVX.AND.(.NOT.LOWP))GO TO 230 E0=E0T GO TO 240 230 E0=E0EFF 240 CONTINUE IF (LOWP) CALL SETUPL(E0,T,R1S,DELTS,NDEGS,ALPHA,IPRT, * IOPTEM) DELT=DELTS*2. C C END OF SEPARATE SECTION FOR LOW PRESSURE LIMIT. C ALL INITIALIZATION HAS NOW BEEN COMPLETED. C -------------------------------------------------------------------- C START LOOP OVER GRAIN SIZES DELTA E C NCUT=2 IF(LOWP) NCUT=1 C C LOOP OVER THE NUMBER OF TIMES THAT THE GRAIN SIZE, DELT, IS HALVED C DO 397 IDT=1,IDM ID=IDT DELT=DELT/2. IF((.NOT.LOWP).OR.(ID.EQ.1))GO TO 250 NE01=NE0*INT(101./(DELT*TESTINC)) NI=NE01 250 CONTINUE NREACT=INT(E0/(DELT*2.859E-3)) IF(LOWP)NREACT=NE0*100/(TESTINC*2*INT(DELT+0.1)) C C PALMT IS THE PARAMETER SUCH THAT THE POPULATION VECTOR FOR C ALL ENERGIES LT E0*PALMT ARE GIVEN THEIR EQUILIBRIUM VALUES. C ALMT1=PALMT*E0/(DELT*2.859E-3) EXPON=1.43879*DELT/T IF(NI.LE.(1000*TESTINC)) GO TO 260 WRITE(6,28) C 28 FORMAT(' NI TRUNCATED ') C NI=1000*TESTINC 260 CONTINUE SUMSCT=0. SUMSC2=0. C C GENDEG AND RATES SORT THE DENSITIES OF STATES (RHO) AND THE MICROSCOPIC C RATE COEFFICIENTS (R1,R2) - FOR REACTING LEVELS - INTO NEW ARRAYS D AND C AVRATE WHERE ENERGY SPACING BETWEEN ELEMENTS IS 100 CM-1. C BOTH SUBROUTINES INTERPOLATE OR TABLE SELECT C DO I=1,(500*TESTINC) AVRATE(I)=0.0 ENDDO IF(.NOT.LOWP) GO TO 270 CALL RATES(R1S,R2) C C GENCF SORTS ROTATIONAL THRESHOLDS, DEXP(-R0/KT), INTO EITHER C 10 OR 100 CM-1 SPACINGS DEPENDING ON VALUE OF TESTINC C IF(JAVX)CALL GENCF GO TO 280 270 CALL RATES(R1,R2) 280 CONTINUE CALL GENDEG(NDEGS) C C COMPUTE STRONG COLLISION RATES (SUMSC1 AND SUMSC2) AND (FOR THE FIRST C VALUE OF DELTA E) THE STRONG COLLISION VECTOR AS INITIAL GUESS TO C EIGENVECTOR Q(I). C ALMT=0. DENST=0. IF(LOWP) BETA=1. DO I=1,NI CI2=0. IF(D(I).GT.0.) CI2=DEXP(DLOG(D(I))-DFLT(I)*EXPON) DENST=DENST+CI2 RTRHB(I)=DSQRT(CI2) C C RTRHB IS BOLTZMANN VECTOR (SYMMETRIZED),Q IS EIGENVECTOR. C FOR ENERGIES BELOW NREACT Q(I) IS ASSIGNED THE VALUE RTRHB(I) C IF(ID.NE.1)GO TO 290 Q(I)=RTRHB(I) C C EXTRA SYMMETRIZATION IS REQUIRED FOR LOW PRESSURE J-AVERAGED C CALCULATION. C IF(JAVX.AND.LOWP)Q(I)=Q(I)/DSQRT(1.-CRF2(I)) 290 CONTINUE C C STRONG COLLISION VECTOR IS FIRST GUESS C AA=1. TT=1. IF(I.LE.NREACT) GO TO 310 AK=AVRATE(I-NREACT) IF(JAVX.AND.LOWP)AK=AK+EXL2(I) AA=OMEGA/(OMEGA+AK) IF(.NOT.LOWP) TT=AA C C IN THE FOLLOWING, THE FIRST GUESS AT THE EIGENVECTOR, Q(I), FOR ENERGIES C ABOVE NREACT IS THE STRONG COLLISION FORM. C C BETA*OMEGA C QSC(E)=F(E)------------------- C BETA*OMEGA + K(E) C C IF(ID.NE.1)GO TO 300 IF(JAVX.AND.LOWP)Q(I)=Q(I)*(1.-CRF2(I)) Q(I)=Q(I)*OMEGA*BETA/(OMEGA*BETA+AK) 300 CONTINUE C THE STRONG COLLISION LOW PRESSURE RATE IS COMPUTED IN SUBROUTINE STRONG. C N.B., STRONG ALWAYS USES A GRAIN SIZE OF 100 CM-1. C C __ OMEGA C \ F(E) -------------- C /_ OMEGA + K(E) C KSC = -------------------------- C __ C \ F(E) C /_ C C STERM=CI2*AK*AA SUMSCT=SUMSCT+STERM SUMCH=0. DO ISUM=1,NCH2 SUMCH=SUMCH+AVR2(I-NREACT,ISUM) ENDDO SUMSC2=SUMSC2+CI2*SUMCH*AA 310 CONTINUE C C ALMT IS USED AS A CORRETED NUMERATOR IN SUBROUTINE NESBET C WHEN KUNI IS CALCULATED, IE: C C C NREACT NI C __ __ OMEGA C ALMT = \ F(E) + \ F(E) -------------- C /_ /_ OMEGA + K(E) C E=0 NREACT C C ALMT=ALMT+TT*CI2 ENDDO SUMSCT=SUMSCT/DENST SUMSC2=SUMSC2/DENST SUMSC1=SUMSCT-SUMSC2 IF(ID.EQ.1) GO TO 320 C C FOR SUBSEQUENT GUESS TO Q(I), USE LINEARLY-INTERPOLATED RESULT FROM C PREVIOUS GRAIN SIZE. C AROW(1)=0.5*Q(1) DO I=2,NILST*2,2 ION2=I/2 AROW(I)=Q(ION2) AROW(I+1)=0.5*(Q(ION2)+Q(ION2+1)) ENDDO DO I=NILST*2+1,NI AROW(I)=AROW(I-1)*0.5 ENDDO DO I=1,NI Q(I)=AROW(I) IF(D(I).LE.0.) Q(I)=0. IF(I.LT.INT(ALMT1)) Q(I)=RTRHB(I) ENDDO C C INTERPOLATED EIGENVECTOR GUESS HAS NOW BEEN FOUND. C 320 CONTINUE C C STORE D(I)*B(I) IN AROW. C NILOC=NI IF(LOWP) NILOC=INT(DFLT(NDEGS)*100./DELT) if(lowp) niloc=niloc/testinc DO 330 J=1,NILOC IF(N.NE.1) AROW(J)=RTRHB(J)**2 IF(.NOT.LOWP) GO TO 330 AROW(J)=0. C C AROW IS NOW IN (100/TESTINC) CM-1 INTERVALS C IF(D(J).GT.0.) AROW(J)= * DEXP(DLOG(D(J))-DFLT(J)*EXPON) 330 CONTINUE NRPL1=NREACT+1 IPRT=IPRT+1 CALL ENEXPT(DELT,ALPHA,IPRT) C C THIS SUBROUTINE GENERATES PROBABILITY MATRIX C THE RESULTING MATRIX IS STORED IN ARRAY PROB. C IF(NBAND.LE.0) GO TO 397 C C NESBET IS SUBROUTINE WHICH CARRIES OUT COMPUTATION OF EIGENVALUE. C SET STARTING VALUE AT LOWEST NON-ZERO POPULATION. C DO I=1,(200*TESTINC) IF(AROW(I).GT.0.) GO TO 340 ENDDO 340 IMIN=I C C THIS LAST SET OF INSTRUCTIONS IS TO ALLOW FOR SOME RHO(E) ( D(I) ) BEING 0. C CALL NESBET(B,ITMAX) C C FIRST TIME THROUGH BOTH DELTAE HALVING LOOP AND PRESSURE LOOP PRINT OUT C LOW PRESSURE RESULTS: C IF(ID.EQ.1)THEN IF(LOWP)THEN WRITE(6,29) C 29 FORMAT(//,10X,'THE FOLLOWING IS FOR THE LOW PRESSURE', *' LIMIT') C WRITE(6,30) C 30 FORMAT(/,' NO. NOT',T12,'DELTA',T20,'BAND-',T30,'MATRIX', * T44,'TOTAL',T54,'NO.OF',/, ' REACTING',T13,'E', * T20,'WIDTH',T31,'SIZE',T40,' RATE (K0)',T52,'ITERATIONS') C ELSE PPAS=PN*133.3 WRITE(6,31) PN,PPAS,OMEGA C 31 FORMAT(1P,6X,'PRESSURE =',T20,E12.2,' TORR (=',E10.2,' PASCAL)' * ,/,10X,'OMEGA =',T20,E12.2,' SEC-1',///, * ' EXACT (EIGENVALUE) RATE STRONG COLLISION ', * 3X,'TOTAL',' DELTA ',' BAND-',' NO.NOT',4X,'NUMBER',/, * 5X,'K1',4X,'K2,3',4X,'TOTAL',6X,'K1',4X,'K2,3', * 5X,'DIMENSION',3X,'E',3X,'WIDTH',' REACTING',1X,'ITERATIONS') C ENDIF ENDIF IF(JAVX)THEN IF(LOWP)THEN CORRF=CORPF ELSE CORRF=CORAV ENDIF ELSE CORRF=CORRAT(ITEMPR) ENDIF RATE1=RATE1*CORRF RATE2=RATE2*CORRF RATTOT=RATTOT*CORRF SUMSC1=SUMSC1*CORRF SUMSC2=SUMSC2*CORRF SUMSCT=SUMSCT*CORRF IF(.NOT.LOWP) WRITE(6,32) RATE1,RATE2, * RATTOT,SUMSC1,SUMSC2,NI,DELT,NBAND, * NREACT,ITR C 32 FORMAT(1P,3E9.2,1X,E7.1,1X,E7.1,0P,I7,F10.0,I5, * I7,3X,I6) C C FOR THE LOW PRESSURE CASE NESBET CALCULATES [M]*K0 IN S-1 C NEED TO CONVERT TO K0 WITH UNITS CM3 MOL S-1 C C USE IDEAL GAS RELATION PV=NRT TO CALCULATE V/N: C IF(LOWP) RSEC=RATTOT*T/(1.603E-5*PN) IF(LOWP) WRITE(6,33) NREACT,DELT,NBAND,NI,RSEC,ITR C 33 FORMAT(1X,I4,5X,F5.0,2X,I5,5X,I5,6X,1P,E10.2,0P,I9) C NILST=NI NI=NI*2 397 CONTINUE IF(N.NE.1.OR.IOPTHT.NE.0) WRITE(6,34) PN,OMEGA,RATE1 C 34 FORMAT(1P,//,10X,25(' *'),/, * 10X,'FINAL RESULT: AT PRESSURE =',E10.2, * ' TORR, OMEGA =',E10.2,' S-1',/, * 20X,' Kuni (S-1) =',T40,E10.2,' (CHANNEL 1)') C IF(.NOT.LOWP.AND.NCHAN.GE.2) * WRITE(6,35)(RATEA(I23)*CORRF,I23=1,NCH2) C 35 FORMAT(1P,T40,E10.2, ' (CHANNEL 2)',/,T40,E10.2,' (CHANNEL 3)') C IF((.NOT.LOWP).AND.RATTOT.GT.1.05*SUMSCT) WRITE(6,36) C 36 FORMAT(5X,'WARNING,WEAK COLL. RESULT GREATER THAN STRONG' * ,', INDICATES NUMERICAL INSTABILITY,' * /,5X, ' REDUCE GRAIN SIZE AND/OR TOLERANCES') C C CARRY OUT HIGH PRESSURE AND STRONG COLLISION LOW PRESSURE CALCULATIONS. C NRCTH=INT(E0*WKA*TESTINC/100.) IF(ICHP.EQ.0) CALL STRONG(HP,ALPT,NDEGS,T,R1,R2) IF(.NOT.JAVX)THEN IF(ICHP.EQ.0)THEN HP=HP*CORRF ALPT=ALPT*CORRF ENDIF ELSE HP=RATHP(1)+RATHP(2)+RATHP(3) ALPT=STLPJ*PN/(T*6.237E+4) ALPT=ALPT/OMEGA ENDIF ICHP=ICHP+1 C C CALL SUBROUTINE TO CALCULATE TROE QUANTITIES FLH, FSC, FSC. C IF(.NOT.LOWP) CALL TROEF(HP,ALPT,SUMSCT,IOPTHT, * RATOTP,T,BETA,NDEGS,R1,R2) WRITE(6,37) C 37 FORMAT(/,25X,'*******************************',///) C NI=NI/2 IF(.NOT.LOWP) GO TO 398 C C CALCULATE AND PRINT OUT LOW PRESSURE RATE COEFFICIENT, RSEC C AND COLLISION EFFICIENCY, BETA C RATOTP=RATTOT/OMEGA CSET=HP/RATTOT NRATES=NRTSTM WRITE(6,38) RSEC C 38 FORMAT(1P,10X,'LOW PRESSURE K0 =',E10.2,' CM3 MOL-1 S-1') C BETA=RATOTP/ALPT WRITE(6,39) BETA C 39 FORMAT(1P,/,20X,'BETA =',T56,E10.2) C NEFF=INT(E0EFF*WKA/DELT) WRITE(6,*) IF(NCHAN.GT.1) CALL MULTCH(R1,R2,E0,NDEGS, * RSEC,T,PN) WRITE(6,40) C 40 FORMAT(////,19X,'-------FALL-OFF CALCULATION-------',///) C C LOW PRESSURE CALC. NOW DONE, SO START FALL-OFF CALCS. C IOPTPR=IOPTEM 398 CONTINUE 399 CONTINUE 400 CONTINUE STOP END C C**************************************************************************** C SUBROUTINE SETUPL(E0,T,R1S,DELTS,NDEGS,ALPHA,IPRT,IOPTEM) C C PREPARATION FOR LOW-PRESSURE CALCULATION C IMPLICIT DOUBLE PRECISION (A-H,O-Z) C PARAMETER (NMAX100=10000, NMAX50=5000, NMAX25=2500) COMMON /AMIN/ ID,NCUT,NILOC COMMON /BLOCK1/ NRATES,NI,NE0,NEFF COMMON /BLOCK2/ DELT,CORRF COMMON /BLOCK3/ NREACT,NBAND,NRPL1 COMMON /BLK10/ AVRATE(NMAX50) COMMON /BLK13/ AROW(NMAX100),OMEGA,ANORM(NMAX100) COMMON /B2/ PROB(NMAX25) COMMON /DENS/ RHO(NMAX100) COMMON /EXP/ ERR1,ERR2,ERR3 COMMON /F1/ IFIRST,NSPEC,IMIN COMMON /F2/ ALMT,ALMT1,DET COMMON /JLP1/ JAVX COMMON /OPT/ IOPTPR COMMON /TRANS1/ RATTOT,RATE1,RATE2,RATEA(2) COMMON /VERSCH/ IXV,JXV COMMON /NEWKIDS/ TESTINC C DIMENSION R1S(NMAX100) C LOGICAL JAVX C C SET UP RATE VECTOR AS COLLISION RATES C DELT=(100./TESTINC) NI=NE0 EXPON=1.43879*DELT/T ITEST=TESTINC*1000 NI=MIN0(ITEST,NI) NREACT=NI/2 NRPL1=NREACT+1 ID=5 IOPTEM=IOPTPR IOPTPR=-1 C C AROW IS THE BOLTZMANN POPULATION VECTOR AT ENERGY I*DELT C C IE AROW(E) = F(E) = RHO(E).EXP(-E/KT) C C NB: THE NORMALIZING FACTOR IS NOT INCLUDED IN AROW C DO I=1,NDEGS AROW(I)=0. IF(RHO(I).GT.0.) AROW(I)=DEXP(DLOG(RHO(I))-DFLT(I)*EXPON) ENDDO DO I=1,(20*TESTINC) IF(RHO(I).GT.0.) GO TO 1000 ENDDO 1000 IMIN=I NCUT=1 NILOC=NDEGS C C THE MAXIMUM NUMBER OF LEVELS, FOR THE PURPOSE OF SETTING UP C THE DIAGONAL COLLISION TERMS IN THE LOW PRESSURE CALCULATION, C MUST BE THE MAXIMUM NUMBER CONSIDERED (NDEGS) FOR CALCULATING C THE NORMALIZING FACTOR ANORM. C C ENEXPT CALCULATES THE PROBABILITY MATRIX FOR COLLISIONAL ENERGY C TRANSFER, IE MATRIX J0, WHERE ALL MICROSCOPIC REACTION RATES ARE C ZERO C C IN A LOW PRESSURE CALCULATION WE CONSIDER THE MATRIX K OF RATE C COEFFICIENTS AS BEING PERTURBED BY THE MATRIX J0 - BUT ALL K'S C ARE ZERO BELOW THE CRITICAL ENERGY E0 SO CONSIDERING THE LARGEST C EIGENVALUE ONLY WE CAN TRUNCATE THE MATRIX AT E0 C CALL ENEXPT(DELT,ALPHA,IPRT) C NIP1=NI+1 INDI=0 DO 1020 I=NRPL1,NI INDI=INDI+1 R1S(INDI)=0. IF(RHO(I).LE.0.) GO TO 1020 SUM=0. C C ONLY COUNT TERMS WITHIN THE BANDWIDTH OF PROB MATRIX C IF(NIP1-I.GT.NBAND) GO TO 1020 JMAX=MIN0(NDEGS,NIP1+NBAND) IF(NIP1.GE.JMAX) GO TO 1020 C C IN LOW PRESSURE LIMIT THE EIGENVALUE RELATIONSHIP IS: C C EO EMAX C __ __ C -K0 G(I) = \ PROB(J-I)G(J)-PROB(I-J)G(I) - \ PROB(J-I)G(I) C /_ /_ C J=0 J=E0 C C WHERE G(I) IS THE EIGENVECTOR OF THE MATRIX J0, = F(E) C AND THE NORMALISATION IS ASSUMED C DO 1010 J=NIP1,JMAX II=J-I+1 IF(II.LE.0.OR.II.GT.NBAND) GO TO 1010 C C TERM IS THE SECOND TERM IN THE ABOVE SUM DIVIDED THROUGH BY G(I) C THIS TERM REPRESENTS THE 'LEAKAGE' TO ENERGY LEVELS ABOVE THE C CRITICAL ENERGY AND IS EQUIVALENT TO THE RATE TERM IN THE NORMAL C FALLOFF PROCEDURE C TERM=(RHO(J)/RHO(I))*DEXP(DFLT(I-J)*EXPON)* * PROB(II)*ANORM(J) C C FOR THE LOW PRESSURE LIMIT, ONE DOES NOT USE THE CORRECT C TRAPEZOIDAL RULE (0.5*FIRST AND LAST TERMS) AS PER THE C FOLLOWING COMMENT CARD, SINCE IT CAN BE SHOWN BY APPROPRIATE C MATRIX PERTURBATION THEORY THAT, FOR THE LOW PRESSURE LIMIT C ONLY, CONVERGENCE IN GRAIN SIZE IS ACCELERATED BY OMITTING IT. C C IF(J.EQ.NIP1.OR.J.EQ.JMAX) TERM=TERM*0.5 C SUM=SUM+TERM 1010 CONTINUE C C IE: E0 EMAX C __ __ C \ \ DELTAE.PROB(J-I).ANORM(J).G(I) C /_ /_ C I=0 J=E0 C [M].K0 = ---------------------------------------- C E0 C __ C \ G(I) C /_ C I=0 C R1S(INDI)=SUM*OMEGA*DELT 1020 CONTINUE C C COLLF GENERATES EXPRESSION FOR EXL: C IF(JAVX)CALL COLLF(DELT) ITEST=TESTINC*1000 NRATES=MIN0(ITEST,INDI) NI=INT(E0/(DELTS*2.859E-3)) RETURN END C C**************************************************************************** C SUBROUTINE RATES(R1,R2) C IMPLICIT DOUBLE PRECISION (A-H,O-Z) C PARAMETER (NMAX100=10000, NMAX50=5000, NMAX25=2500) COMMON /LOW/ LOWP COMMON /BLOCK1/ NRATES,NI,NE0,NEFF COMMON /BLOCK2/ DELT,CORRF COMMON /BLK10/ AVRATE(NMAX50) COMMON /BLK12/ AVR2(NMAX50,2) COMMON /BLK15/ NCH2 COMMON /DEOC/ EI,WE0 COMMON /NEWKIDS/ TESTINC C DIMENSION R1(1),R2(NMAX100,2),ALAST2(2),ANEXT2(2),SP2(2) C LOGICAL LOWP C C FORM VECTORS OF MICROSCOPIC RATES REQUIRED FOR GIVEN ENERGY SPACING C FROM THE INPUT MICRSCOPIC RATES FOR 100 CM-1 SPACING (R1,R2).IF C SPACING GE 100 CM-1, DO BY TABLE SELECTION, OTHERWISE BY LINEAR C INTERPOLATION. C BRANCH DEPENDING ON WHETHER DELT LT OR GE (100. OR 10.) C IF(DELT.LT.(99.99/TESTINC)) GO TO 2030 C C FIND INTEGER FOR TABLE SELECTION, EXTRA 1. IS TO ACCOUNT FOR ANY ROUNDOFF C N=INT((DELT+1.)*(0.01*TESTINC)) J=0 IF(DELT.GT.(101./TESTINC))GO TO 2000 N3=1 GO TO 2010 C C DELT GT (100/TESTINC) CM-1: C 2000 CONTINUE INTE0=INT(WE0*TESTINC/100.) IF(LOWP)INTE0=INT(WE0*TESTINC/200.) IF(MOD(INTE0,2).EQ.1)THEN N3=1 ELSE N3=2 ENDIF C C DELT GE 100/TESTINC CM-1 C WANT TO END UP WITH 100/TESTINC CM-1 SPACING REGARDLESS C 2010 CONTINUE DO 2020 II=N3,NRATES,N I=MIN0(II,NRATES) J=J+1 DO ISUM=1,NCH2 AVR2(J,ISUM)=0. ENDDO AVRATE(J)=R1(I) IF(LOWP) GO TO 2020 SUMCH=0. DO ISUM=1,NCH2 AVR2(J,ISUM)=R2(I,ISUM) SUMCH=SUMCH+AVR2(J,ISUM) ENDDO AVRATE(J)=R1(I)+SUMCH 2020 CONTINUE GO TO 2055 C C DELT LT (100/TESTINC) CM-1 C FORM RATE VECTOR USING LINEAR INTERPOLATION C 2030 N=INT(101./DELT) JJ=0 ALAST1=(R1(1)) DO ISUM=1,NCH2 ALAST2(ISUM)=0. ENDDO DO 2050 I=1,NRATES IF(I.LT.NRATES) ANEXT1=R1(I+1) DO ISUM=1,NCH2 ANEXT2(ISUM)=0. ENDDO IF(I.EQ.NRATES) GO TO 2040 IF(LOWP.OR.R2(I+1,1).LE.0.) GO TO 2040 DO ISUM=1,NCH2 ANEXT2(ISUM)=R2(I+1,ISUM) ENDDO 2040 CONTINUE SP1=(ANEXT1-ALAST1)/DFLT(N) DO ISUM=1,NCH2 SP2(ISUM)=(ANEXT2(ISUM)-ALAST2(ISUM))/DFLT(N) ENDDO DO J=1,N JJ=JJ+1 AVRATE(JJ)=(ALAST1+SP1*DFLT(J-1)) SUMCH=0. DO ISUM=1,NCH2 AVR2(JJ,ISUM)=ALAST2(ISUM)+SP2(ISUM)*DFLT(J-1) SUMCH=SUMCH+AVR2(JJ,ISUM) ENDDO AVRATE(JJ)=AVRATE(JJ)+SUMCH ENDDO ALAST1=ANEXT1 DO ISUM=1,NCH2 ALAST2(ISUM)=ANEXT2(ISUM) ENDDO 2050 CONTINUE 2055 RETURN END C C*********************************************************************** C SUBROUTINE GENCF C C SELECTS DEXP(-R0/KT) BY GRAIN SIZE, ONLY UTILISED IF JAVX. C IMPLICIT DOUBLE PRECISION (A-H,O-Z) C PARAMETER (NMAX100=10000, NMAX50=5000, NMAX25=2500) COMMON /BLOCK1/ NRATES,NI,NE0,NEFF COMMON /BLOCK2/ DELT,CORRF COMMON /JLP2/ CRF2(NMAX100),EXL2(NMAX100) COMMON /JLP3/ CRF(NMAX100),EXL(NMAX100) COMMON /NEWKIDS/ TESTINC C IF(DELT.LT.(99.99/TESTINC))GO TO 3000 N=INT((DELT+0.01)*(0.01*TESTINC)) J=0 C C NT IS THE MAXIMUM ENERGY FOR WHICH RATES WERE ENTERED C CRF IS EXPRESSION FOR EXP(-ROTHR/KT) C EXL IS CALCULATED IN SUBROUTINE COLLF CALLED FROM SETUPL C NT=NE0+NRATES DO II=N,NT,N I=II IF(I.GT.NT)I=NT J=J+1 CRF2(J)=CRF(I) EXL2(J)=EXL(I) ENDDO GO TO 3001 C C SPLIT FOR DELTA LT (100/TESTINC) CM-1, C FORM VECTORS CRF2 AND EXL2 BY LINEAR INTERPOLATION. C 3000 CONTINUE N=INT(101./(DELT*TESTINC)) JJ=0 ALAST2=CRF(1) ALAST3=EXL(1) DO I=1,NE0 IF(I.LT.NE0)ANEXT2=CRF(I+1) IF(I.LT.NE0)ANEXT3=EXL(I+1) SP2=(ANEXT2-ALAST2)/DFLT(N) SP3=(ANEXT3-ALAST3)/DFLT(N) DO J=1,N JJ=JJ+1 CRF2(JJ)=ALAST2+SP2*DFLT(J-1) EXL2(JJ)=ALAST3+SP3*DFLT(J-1) ENDDO ALAST2=ANEXT2 ALAST3=ANEXT3 ENDDO 3001 RETURN END C C************************************************************************** C SUBROUTINE GENDEG(NDEGS) C C FINDS VECTOR OF DEGENERACIES FOR ANY GRAIN SIZE BY INTERPOLATION C IN THIS CASE VECTOR IS THE DENSITY OF STATES, RHO C IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (NMAX100=10000, NMAX50=5000, NMAX25=2500) COMMON /BLOCK1/ NRATES,NI,NE0,NEFF COMMON /BLOCK2/ DELT,CORRF COMMON /BLOCK7/ D(NMAX100) COMMON /DENS/ RHO(NMAX100) COMMON /NEWKIDS/ TESTINC C C BRANCH DEPENDING ON WHETHER DELT IS GT OR LT 100. C IF(DELT.LT.(99.99/TESTINC)) GO TO 4000 C C FIND INTEGER FOR TABLE SELECTION C N=INT((DELT+0.01)*(0.01*TESTINC)) J=0 DO II=N,NDEGS,N I=II IF(I.GT.NDEGS) I=NDEGS J=J+1 D(J)=RHO(I) ENDDO GO TO 4001 4000 N=INT(101./(DELT*TESTINC)) JJ=0 ALAST=DLOG(RHO(1)) DO I=1,NDEGS IF(I.LT.NDEGS) ANEXT=DLOG(RHO(I+1)) SP=(ANEXT-ALAST)/DFLT(N) DO J=1,N JJ=JJ+1 D(JJ)=DEXP(ALAST+SP*DFLT(J-1)) ENDDO ALAST=ANEXT ENDDO 4001 RETURN END C C*********************************************************************** C SUBROUTINE ENEXPT(DELT,ALPHA,IPRT) C C COMPUTES COLLISION PROBABILITY ARRAY PROB. C TO CHANGE PROBABILITY FUNCTION, CHANGE FUNCTIONAL FORM C IMPLICIT DOUBLE PRECISION (A-H,O-Z) C PARAMETER (NMAX100=10000, NMAX50=5000, NMAX25=2500) COMMON /BLOCK1/ NRATES,NI,NE0,NEFF COMMON /VERSCH/ IXV,JXV COMMON /LOW/ LOWP COMMON /BLOCK3/ NREACT,NBAND,NRPL1 COMMON /AMIN/ ID,NCUT,NILOC COMMON /BLK13/ AROW(NMAX100),OMEGA,ANORM(NMAX100) COMMON /B2/ PROB(NMAX25) COMMON /EXP/ ERR1,ERR2,ERR3 COMMON /NEWKIDS/ TESTINC C LOGICAL LOWP,BRW C NBAND=0 BRW=IXV.LE.0.AND.JXV.LE.0 NE00=NREACT IF(LOWP) NE00=NI Z=DLOG(AROW(NE00+1)/AROW(NE00))*ALPHA*ALPHA/DELT FS2=4.*ALPHA*ALPHA NTEST=INT(249.9*TESTINC) NTOPX=MIN0(NTEST,NILOC) DO I=1,NTOPX PROB(I)=0. ENDDO C C GENERATE UN-NORMALIZED PROBABILITY MATRIX, J0, WITH ELEMENTS: C C J0(I,J)=R(I,J) C __ C J0(I,I)=- \ R(J,I) C /_ C I.NE.J C C WHERE R(I,J), THE RATE OF ENERGY TRANSFER FROM LEVEL J TO I, C DEPENDS ONLY ON THE ENERGY DIFFERENCE J-I AND IS GIVEN BY C PROB(J-I) - WHICH NEEDS TO BE NORMALISED BY ANORM(J) C BOT=0. DEDOWN=0. DO I=1,NTOPX DE=DFLT(I-1)*DELT X=DE/ALPHA IF(BRW) GO TO 5000 C C USE EXPONENTIAL DOWN MODEL FOR PROB(J-I) C XX=0. IF(X.GT.0.) XX=X**JXV PWR=1. IF((X.GT.0.).AND.IXV.NE.0) PWR=X**IXV IF((X.LE.0.).AND.IXV.NE.0) PWR=0. IF(XX.LT.20.) PROB(I)=PWR*DEXP(-XX) GO TO 5010 C C USE BIASED RANDOM WALK MODEL FOR PROB(J-I) C 5000 TOP=-(Z+DE) PROB(I)=DEXP(-TOP*TOP/FS2) C C FOLLOWING ADDED BY RGG JUN 29 1989 C IF(PROB(I).LT.ERR3) GO TO 5020 5010 CONTINUE IF(PROB(I).LE.ERR3.AND.X.GT.0.5) GO TO 5020 C C INTEGRATE USING THE TRAPEZOIDAL RULE, OMITTING THE FACTOR OF 0.5 C AT THE HIGH ENERGY END C IF(I.GT.1)THEN DEDOWN=DEDOWN+DE*PROB(I) BOT=BOT+PROB(I) ELSE DEDOWN=0.5*DE*PROB(1) BOT=0.5*PROB(1) ENDIF C C NBAND IS THE BANDWIDTH OF THE MATRIX PROB, IE HOW FAR YOU MUST GO C FROM THE DIAGONAL BEFORE THE OFF-DIAGONAL ELEMENTS ARE LESS THAN C ERR3 C NBAND=MAX0(NBAND,I) ENDDO C C E' C __ C \ (E'-E)PROB(E'-E) C /_ C E=0 C SET PROB(250*TESTINC) = = ------------------------ C E' C __ C \ PROB(E'-E) C /_ C E=0 C 5020 PROB(I)=0. PROB(TESTINC*250)=DEDOWN/BOT NBAND=NBAND-1 IF(NBAND.EQ.1) NBAND=0 IF(NBAND.EQ.0) GO TO 5125 IF(.NOT.LOWP) NILOC=NI C C COMMENCE FINITE DIFFERENCE SOLUTION OF INTEGRAL EQUATION FOR C. C JJ=NILOC+1 C C GENERATE ANORM, THE VECTOR OF NORMALIZERS OF THE ENERGY TRANSFER C RATE, R(I,J), INITIALLY AS C, ITS RECIPROCAL. C DO 5060 J=1,NILOC JJ=JJ-1 ANORM(JJ)=0. IF(AROW(JJ).LE.0.) GO TO 5060 JP=JJ+1 JM=JJ-1 IF(JJ.GE.(NREACT/NCUT)) GO TO 5030 C C TO AVOID OCCASIONAL PROBLEMS WITH NORMALIZATION ALGORITHM FOR STATES C FOR STATES OF LOW ENERGY (BELOW NREACT/NCUT), ALL ARE GIVEN THE C SAME NORMALIZATION. THIS HAS NO PHYSICAL EFFECT, SINCE THESE LEVELS C ALWAYS HAVE THEIR EQUILIBIRIUM POPULATIONS. C ANORM(JJ)=ANORM(JJ+1) GO TO 5060 5030 CONTINUE DE=PROB(1) IMIN=MAX0(1,JJ-NBAND) IF(JM.LT.IMIN) GO TO 5040 C C I.E. FOR THE FIRST LOOP THROUGH FIND NORMALISATION AT MAX ENERGY: C C NMAX C __ C ANORM(NMAX) = DELTAE * \ PROB(NMAX-I) C /_ C I=1 C DO I=IMIN,JM DE=DE+PROB(JP-I) ENDDO 5040 ANORM(JJ)=DE*DELT IF(JJ.EQ.NILOC) GO TO 5060 IUP=MIN0(NILOC,JJ+NBAND) DE=0. IF(IUP.LT.JP) GO TO 5060 C C SUBSEQUENTLY AT ENERGY JJ: C C C JJ C __ C DELTAE \ PROB(JJ-I) C /_ C I=1 C ANORM(JJ) = ----------------------------------------------- C NMAX C DELTAE __ C 1 - -------- \ AROW(I)*PROB(I-JJ)/ANORM(I) C AROW(JJ) /_ C I=JP C C DO 5050 I=JP,IUP IF(ANORM(I).LE.0.) GO TO 5050 DE=DE+(PROB(I-JM)*AROW(I)/ANORM(I)) 5050 CONTINUE ANORM(JJ)=ANORM(JJ)/(1.-(DELT*DE/AROW(JJ))) 5060 CONTINUE C C INVERT ANORM TO OBTAIN NORMALISATION VECTOR C JS=NILOC+1 DO J=1,NILOC JS=JS-1 IF(ANORM(JS).GT.0.) ANORM(JS)=1./ANORM(JS) IF(ANORM(JS).LE.0..AND.JS.LT.NILOC) ANORM(JS)=ANORM(JS+1) IF(AROW(J).LE.0.) ANORM(J)=0. ENDDO IF(LOWP.AND.ID.EQ.5) GO TO 5070 IF(IPRT.NE.1) GO TO 5125 IF(BRW) WRITE(6,5001) DABS(ALPHA) C 5001 FORMAT(//,10X,'BIASED RANDOM WALK MODEL FOR P(E,E',1H',')',/, * 20X,'S =',F10.1,' CM-1') C IF (TESTINC .EQ. 1)WRITE(6,5002) (PROB(I),I=1,NBAND) C C Print out only every 10 R(X) values in inc=10 C IF (TESTINC .EQ. 10)WRITE(6,5005) (PROB(I),I=1,NBAND,10) C 5002 FORMAT(/,' VALUES OF R(X) (FOR LARGEST GRAIN SIZE):', * /,1P,10(3X,6E12.2,/)) 5005 FORMAT(/,' VALUES OF R(X) (EVERY TEN ) * (FOR LARGEST GRAIN SIZE):', * /,1P,10(3X,6E12.2,/)) C GO TO 5125 C C COMPUTE AVERAGE ENERGY TRANSFER QUANTITIES 5070 WRITE(6,5003) C 5003 FORMAT(//,' ENERGY (KJ/MOL) DELTA E TOTAL (CM-1)', * ' SQRT(DE**2)',T56,'EQUILIBRIUM POPULATION') C NZ=NI/2 SUMP=0. DO J=1,NILOC SUMP=SUMP+AROW(J) ENDDO ATEST=1.0 DO 5120 II=NZ,NILOC,10 E=DFLT(II)*DELT*1.195E-2 ANM=ANORM(II) IF(ANM.LE.0.) GO TO 5120 JP=II+1 JM=II-1 DE=0. DE2=0. IF(JM.LT.1) GO TO 5090 DO 5080 I=1,JM IF(II-I.GE.NBAND) GO TO 5080 FAC=ANM*PROB(JP-I) DIF=DFLT(I-II) DE=DE+DIF*FAC DE2=DE2+DIF*DIF*FAC 5080 CONTINUE 5090 IF(JP.GT.NILOC.OR.AROW(II).LE.0.) GO TO 5110 RBJ=1./AROW(II) DO 5100 J=JP,NILOC IF((J-II.GE.NBAND).OR.(ANORM(J).LE.0.)) GO TO 5100 FAC=PROB(J-JM)*ANORM(J)*AROW(J)*RBJ DIF=DFLT(J-II) DE=DE+DIF*FAC DE2=DE2+DIF*DIF*FAC 5100 CONTINUE 5110 DE=DE*DELT*DELT DE2=DSQRT(DE2*DELT)*DELT AROWN=AROW(II)/(DELT*SUMP) C C If INC = 10cm-1 then print out only very ten sets of calculations. C i.e. if testinc=10 C IF (TESTINC .EQ. 1) WRITE(6,5004) E,DE,DE2,AROWN IF (TESTINC .EQ. 10 .AND. (ATEST/10) .EQ. * INT(ATEST/10)) WRITE(6,5004) E,DE,DE2,AROWN C 5004 FORMAT(1X,F9.1,10X,F10.1,6X,F10.1,T62,1P,E10.2) C ATEST=ATEST+1 5120 CONTINUE 5125 RETURN END C C************************************************************************ C SUBROUTINE NESBET(B,ITMAX) C C SUBROUTINE TO SOLVE EIGENVALUE PROBLEM. C IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARAMETER (NMAX100=10000, NMAX50=5000, NMAX25=2500) COMMON /AMIN/ ID,NCUT,NILOC COMMON /BLOCK1/ NRATES,NI,NE0,NEFF COMMON /BLOCK2/ DELT,CORRF COMMON /BLOCK3/ NREACT,NBAND,NRPL1 COMMON /BLOCK6/ Q(NMAX100),H(NMAX100),RTRHB(NMAX100) COMMON /BLK10/ AVRATE(NMAX50) COMMON /BLK12/ AVR2(NMAX50,2) COMMON /BLK13/ AROW(NMAX100),OMEGA,ANORM(NMAX100) COMMON /BLK15/ NCH2 COMMON /BLOK17/ AVEC(NMAX100) COMMON /B2/ PROB(NMAX25) COMMON /EXP/ ERR1,ERR2,ERR3 COMMON /F1/ IFIRST,NSPEC,IMIN COMMON /F2/ ALMT,ALMT1,DET COMMON /JLP1/ JAVX COMMON /JLP2/ CRF2(NMAX100),EXL2(NMAX100) COMMON /LOW/ LOWP COMMON /OPT/ IOPTPR COMMON /TRANS1/ RATTOT,RATE1,RATE2,RATEA(2) COMMON /TRANS2/ ITR COMMON /NEWKIDS/ TESTINC C DIMENSION RATK(NMAX50,2) C LOGICAL DONE,LOWP,JAVX C C INTEGRATION IS DONE USING COMPLETE TRAPEZOIDAL RULE, FOR INNER LOOPS C (J LE I OR J GE I). FOR OUTER LOOPS, THE FIRST AND LAST VALUES OF THE C INTEGRAND CONTRIBUTE NEGLIGIBLY TO THE INTEGRAL, SO THE 0.5*INTEGRAND C CORRECTION USED IN PROPER TRAPEZOIDAL INTEGRATION IS NOT INCLUDED. C AROW(I) IS RHO(I)*DEXP(-E/KT) = F(I) C GAMMA IS Q(I)/RTRHB(I) = Q(I)/DSQRT(F(I)) C Q(I) IS THE EIGENVECTOR AND INITIALLY WE START WITH: C C Q(I) = SQRT(F(I)) FOR I < ALMT1 C = F(I) FOR ALMT1 < I < NI C RATEL=0. IFIRST=INT(ALMT1) IF(IFIRST.LT.IMIN) IFIRST=IMIN DO I=1,2 DO II=1,(500*TESTINC) RATK(II,I)=0. ENDDO ENDDO RATEL2=0. DEOMGA=OMEGA*DELT DET=0. DET2=0. SS=0. C C THE FOLLOWING LOOP ALSO GENERATES FIRST GUESS AT RATE COEFFICIENTS C LOOP FROM THE FIRST NON-ZERO ELEMENT OF AROW C DO 6000 I=IMIN,NI H(I)=1. AVEC(I)=0. IF(AROW(I).LE.0.) GO TO 6000 AVEC(I)=AROW(I)*ANORM(I) CI2=Q(I)*RTRHB(I) IF(JAVX.AND.LOWP)CI2=CI2*DSQRT(1.-CRF2(I)) SS=SS+CI2 H(I)=Q(I)/RTRHB(I) IF(I.LT.IFIRST) H(I)=1. IF(I.LE.NREACT) GO TO 6000 EFR=AVRATE(I-NREACT) IF(JAVX.AND.LOWP)EFR=EFR+EXL2(I) DET=DET+EFR*CI2 SUMCH=0. DO ISUM=1,NCH2 SUMCH=SUMCH+AVR2(I-NREACT,ISUM) ENDDO DET2=DET2+SUMCH*CI2 6000 CONTINUE IF(IOPTPR.LT.0) SS=ALMT C C FIRST GUESS AT KUNI IS: C C __ C \ Q(I)*K(I)*SQRT(F(I)) C /_ C KUNI = ------------------------- C __ C \ Q(I)*SQRT(F(I)) C /_ C C C WHERE Q(I) IS THE FIRST GUESS AT THE EIGENVECTOR C RATE=DET/SS IF(ID.NE.1) RATE=RATTOT RATE2=DET2/SS C C INITIALIZATION PROCESS NOW COMPLETE. C C COMMENCE NESBET ITERATIVE PROCESS. C C LOOP THROUGH THE NUMBER OF NESBET ITERATIONS: C DO 6060 ITRX=1,ITMAX ITR=ITRX C C THE FOLLOWING CONVERGENCE TEST GIVES EQUAL WEIGHTING TO THE ACCURACY C OF THE OVERALL RATE AND THE RATE FROM THE SECOND CHANNEL. C ALTERNATIVELY,THE WEIGHTING CAN BE FOR OVERALL RATE,DEPENDING ON INPUT C OPTION. C DONE=ABS(1.-(RATEL/RATE)).LE.ERR2 IF(IABS(IOPTPR).EQ.2) DONE=DONE.AND.ABS(1.-(RATEL2/RATE2)).LE. * ERR2 IF(DONE) GO TO 6070 RATEL2=RATE2 RATEL=RATE C C LOOP THROUGH THE ENERGY: C DO 6050 I=IFIRST,NI IF(AROW(I).LE.0.) GO TO 6050 QI=Q(I) CRFI=(1.-CRF2(I)) RTI=RTRHB(I) IP=I+1 HI=H(I) IL1=I-1 JMIN=MAX0(IMIN,I-NBAND+1) S=0. IF(I.EQ.NI) GO TO 6010 JMAX=MIN0(NI,I+NBAND-1) C C COMMENCE LOOP OVER J > I. C C NMAX C __ C 1 \ F(J) * PROB(J-I) * ANORM(J) * (H(I)-H(J)) C ---------- /_ C SQRT(F(I)) J=I C C IF(JMAX-IL1.LE.0) GO TO 6010 HT=H(JMAX) IF(JAVX.AND.LOWP)HT=HT*DSQRT(CRFI*(1.-CRF2(JMAX))) S=0.5*(HI-HT)*AVEC(JMAX)*PROB(JMAX-IL1) JM=JMAX-1 IF(JM.LT.IP) GO TO 6010 DO J=IP,JM HT=H(J) IF(JAVX.AND.LOWP)HT=HT*DSQRT(CRFI*(1.-CRF2(J))) S=S+(HI-HT)*AVEC(J)*PROB(J-IL1) ENDDO 6010 CONTINUE S=S/RTI S1=0. IF(I.LE.JMIN) GO TO 6020 C C COMMENCE LOOP OVER J < I. C C I C __ C SQRT(F(I)) * ANORM(I) * \ PROB(I-J) * (H(I)-H(J)) C /_ C J=1 C C IF(IP-JMIN.LE.0) GO TO 6020 HL=H(JMIN) IF(JAVX.AND.LOWP)HL=HL*DSQRT(CRFI*(1.-CRF2(JMIN))) S1=0.5*PROB(IP-JMIN)*(HI-HL) JM=JMIN+1 IM=I-1 IF(IM.LT.JM) GO TO 6020 DO J=JM,IM HL=H(J) IF(JAVX.AND.LOWP)HL=HL*DSQRT(CRFI*(1.-CRF2(J))) S1=S1+PROB(IP-J)*(HI-HL) ENDDO 6020 AA=-RATE S=S+S1*RTI*ANORM(I) IF(JAVX.AND.LOWP)S=S+PROB(1)*QI*ANORM(I)*CRF2(I) IF(I.LE.NREACT) GO TO 6030 ILNR=I-NREACT AKE=AVRATE(ILNR) AA=AA+AKE 6030 CONTINUE C C ( __ __ ) C ( \ + \ ) C S(I) = (K(I) - KUNI)*Q(I) + OMEGA*DELTAE * ( /_ /_ ) C ( J>I J (BY NUMERICAL INTEGRATION) =',F10.1, * ' CM-1') C RATIOT=RATTOT/HP APOM=ALPT*OMEGA WRITE(6,8002) HP,RATIOT,APOM C 8002 FORMAT(1P,20X,'TOTAL K(INF) =',T56,E10.2,' S-1' * ,/,20X,'K(TOT)/K(INF) =' * ,T56,E10.2,/,20X,'STRONG COLLISION K0*[M] =',T56 * ,E10.2,' S-1') C IF(APOM.LT.SUMSCT) WRITE(6,8003) C 8003 FORMAT(' WARNING, STRONG COLLISION RATE GT STRONG COLLISION', * /,' VALUE FOR K0*[M] USE SMALLER GRAIN SIZE') C IF(IOPTHT.NE.0) RETURN C C COMPUTE LINDEMANN-HINSHELWOOD AND STRONG COLLISION BROADENING FACTORS. C WCK0=RATOTP*OMEGA IF(BETA.GT.1.) WRITE(6,8004) C 8004 FORMAT(' WARNING, BETA GREATER THAN 1, INDICATES' * ,' TOLERANCES SHOULD BE ALTERED') C WCRAT=WCK0/HP WRITE(6,8005) WCK0,WCRAT C 8005 FORMAT(1P,20X,'WEAK COLLISION K0*[M] =',T56,E10.2,' S-1', * /,20X,'WEAK COLLISION K0*[M]/K(INF)=',T56,E10.2) C WRITE(6,8006) BETA C 8006 FORMAT(1P,/,20X,'BETA =',T56,E10.2) C FLH=WCRAT/(1.+WCRAT) OMSC=RATOTP*OMEGA/ALPT AKSC=0. EXPON=143.879/(T*TESTINC) BOT=0. ICSC=0. IF(JAVX)NN11=NEFF IF(.NOT.JAVX)NN11=NE0 DO 8000 I=1,NDEGS TERM=0. IF(RHO(I).GT.0.) TERM=DEXP(DLOG(RHO(I))-DFLT(I)*EXPON) BOT=BOT+TERM IF(I.LE.(NN11)) GO TO 8000 ICSC=ICSC+1 AKE=R1(ICSC)+(R2(ICSC,1)+R2(ICSC,2)) AKSC=AKSC+(TERM*OMSC*AKE/(AKE+OMSC)) 8000 CONTINUE AKSC=AKSC*CORRF/BOT FSC=AKSC/(FLH*HP) FWC=RATTOT/(HP*FLH*FSC) WRITE(6,8007) FLH,FSC,FWC C 8007 FORMAT(1P,//,20X,'FLH =',E10.2,',FSC =',E10.2,',FWC =',E10.2) C IF(FWC.GT.1.5) WRITE(6,8008) C 8008 FORMAT(' WARNING,FWC SIGNIFICANTLY GREATER THAN 1 :' * ,/,' TOLERANCES SHOULD BE CHANGED', * ' OR THAT GRAIN SIZE SHOULD BE' * ,/,' REDUCED TO 100 CM-1 (OR 10 CM-1 AT LOW TEMPERATURES)') C RETURN END C C**************************************************************************** C SUBROUTINE MULTCH(R1,R2,E0,NDEGS,RSEC,T,PN) C C CALCULATES INDIVIDUAL LOW-PRESSURE RATES FOR MULTIPLE CHANNEL REACTIONS. C IMPLICIT DOUBLE PRECISION (A-H,O-Z) C PARAMETER (NMAX100=10000, NMAX50=5000, NMAX25=2500) COMMON /AMIN/ ID,NCUT,NILOC COMMON /BLOCK1/ NRATES,NI,NE0,NEFF COMMON /BLOCK2/ DELT,CORRF COMMON /BLOCK3/ NREACT,NBAND,NRPL1 COMMON /BLOCK6/ Q(NMAX100),H(NMAX100),RTRHB(NMAX100) COMMON /BLK13/ AROW(NMAX100),OMEGA,ANORM(NMAX100) COMMON /BLK15/ NCH2 COMMON /B2/ PROB(NMAX25) COMMON /DENS/ RHO(NMAX100) COMMON /F1/ IFIRST,NSPEC,IMIN COMMON /F2/ ALMT,ALMT1,DET COMMON /JLP1/ JAVX COMMON /JLP2/ CRF2(NMAX100),EXL2(NMAX100) COMMON /OPT/ IOPTPR COMMON /NEWKIDS/ TESTINC C DIMENSION R1(NMAX100),R2(NMAX100,2),RTE2(2),SUM(2), * SUMR(NMAX100,2),CHMRT(2) C LOGICAL JAVX C DO ISUM=1,NCH2 RTE2(ISUM)=0. CHMRT(ISUM)=0. ENDDO RTOT=0. DELT=100./TESTINC EXPON = 1.43879*DELT/T NIT=INT(NE0*100/DELT) NIP1=NIT+1 S=0. DO I=1,NRATES SUMX=R1(I) DO ISUM=1,NCH2 SUMX=SUMX+R2(I,ISUM) ENDDO DO ISUM=1,NCH2 SUMR(I,ISUM)=R2(I,ISUM)/SUMX ENDDO ENDDO IF(.NOT.JAVX)NEFF=NIT DO 9040 I=IMIN,NIT BBI=Q(I)*RTRHB(I) IF(JAVX )BBI=BBI*DSQRT(1.D0-CRF2(I)) S=S+BBI DO ISUM=1,NCH2 SUM(ISUM)=0. ENDDO SUMX=0. IF(I.LE.NREACT) GO TO 9040 IF(NIP1-I.GT.NBAND) GO TO 9010 JMAX=MIN0(NDEGS,NIP1+NBAND) IF(NIP1.GE.JMAX) GO TO 9040 DO 9000 J=NIP1,JMAX II=J-I+1 IF(II.LE.0) GO TO 9000 TERM=(RHO(J)/RHO(I))*DEXP(DFLT(I-J)*EXPON)* * PROB(II)*ANORM(J) SUMX=SUMX+TERM IF(R2(J-NEFF,1).LE.0.) GO TO 9000 DO ISUM=1,NCH2 SUM(ISUM)=SUM(ISUM)+TERM*SUMR(J-NEFF,ISUM) ENDDO 9000 CONTINUE 9010 CONTINUE IF(.NOT.JAVX)GO TO 9030 IB=MAX0(NEFF+1,I-NBAND+1) IT=MIN0(NIT,I+NBAND-1) IF(I.GE.IT)GO TO 9020 IF(I.LE.NEFF)GO TO 9030 DO II=I+1,IT IK=II-I+1 TERM=(RHO(II)/RHO(I))*DEXP(DFLT(I-II)*EXPON)*PROB(IK) * *ANORM(II)*CRF2(II) DO ISUM=1,NCH2 SUM(ISUM)=SUM(ISUM)+TERM*SUMR(II-NEFF,ISUM) ENDDO SUMX=SUMX+TERM ENDDO 9020 CONTINUE DO II=IB,I IK=I+1-II TERM=PROB(IK)*ANORM(I)*CRF2(II) DO ISUM=1,NCH2 SUM(ISUM)=SUM(ISUM)+TERM*SUMR(II-NEFF,ISUM) ENDDO SUMX=SUMX+TERM ENDDO 9030 CONTINUE DO ISUM=1,NCH2 RTE2(ISUM)=RTE2(ISUM)+BBI*SUM(ISUM) ENDDO RTOT=RTOT+BBI*SUMX 9040 CONTINUE C C RSEC IS TOTAL SECOND-ORDER LOW PRESSURE RATE COEFFT; HERE WE HAVE C CALCULATED RATE 2, RATE 1 IS FOUND BY DIFFERENCE. C SUMCH=0. DO ISUM=1,NCH2 CHMRT(ISUM)=RTE2(ISUM)/RTOT SUMCH=SUMCH+CHMRT(ISUM) RTE2(ISUM)=CHMRT(ISUM)*RSEC ENDDO RATE1=(1.-SUMCH)*RSEC WRITE(6,9001) RATE1,(RTE2(I7),I7=1,NCH2) C 9001 FORMAT(//,' MULTIPLE CHANNEL LOW-PRESSURE RATE COEFFICIENTS', * ' (CM3 MOL-1 S-1):',/,15X,1P,3E10.2) C IF(JAVX)WRITE(6,9002) C 9002 FORMAT(/,' NOTE: THIS RATIO WILL ONLY BE VALID IF THE ', * 'FIRST FILE OF MULTICHANNEL K(E) ',/,' INPUT CORRESPONDS ', * 'TO A SUFFICIENTLY LOW PRESSURE.(SEE PROGRAM NOTES).',/) C RETURN END C C*********************************************************************** C SUBROUTINE COLLF(DELT) C C FINDS TOTAL COLLISION RATE OUTSIDE ROTATIONAL THRESHOLD. C IMPLICIT DOUBLE PRECISION (A-H,O-Z) C PARAMETER (NMAX100=10000, NMAX50=5000, NMAX25=2500) COMMON /BLOCK1/ NRATES,NI,NE0,NEFF COMMON /BLOCK3/ NREACT,NBAND,NRPL1 COMMON /BLK13/ AROW(NMAX100),OMEGA,ANORM(NMAX100) COMMON /B2/ PROB(NMAX25) COMMON /JLP3/ CRF(NMAX100),EXL(NMAX100) COMMON /NEWKIDS/ TESTINC C DO I=1,(TESTINC*1000) EXL(I)=0. ENDDO DO N=NREACT,NI IMAX=MIN0(NI,N+NBAND-1) IMIN=MAX0(NREACT,N-NBAND+1) IP=N+1 IM=N-1 S33=0. IF(N.EQ.NI)GO TO 1021 DO I=IP,IMAX IK=I-IM S33=S33+PROB(IK)*AROW(I)*CRF(I)*ANORM(I) ENDDO S33=S33*DELT*OMEGA/AROW(N) 1021 CONTINUE S44=0. IF(N.EQ.NREACT)GO TO 1022 DO I=IMIN,IM IK=IP-I S44=S44+PROB(IK)*CRF(I) ENDDO S44=S44+PROB(1)*CRF(N) S44=S44*OMEGA*DELT*ANORM(N) 1022 CONTINUE EXL(N)=S44+S33 ENDDO RETURN END C C************************************************************************* C BLOCK DATA C IMPLICIT DOUBLE PRECISION (A-H,O-Z) C PARAMETER (NMAX100=10000, NMAX50=5000, NMAX25=2500) COMMON /BLOCK6/ Q(NMAX100),H(NMAX100),RTRHB(NMAX100) COMMON /BLOCK7/ D(NMAX100) COMMON /BLK10/ AVRATE(NMAX50) COMMON /BLK12/ AVR2(NMAX50,2) COMMON /BLK13/ AROW(NMAX100),OMEGA,ANORM(NMAX100) COMMON /B2/ PROB(NMAX25) COMMON /DENS/ RHO(NMAX100) COMMON /JLP2/ CRF2(NMAX100),EXL2(NMAX100) COMMON /JLP3/ CRF(NMAX100),EXL(NMAX100) COMMON /ROTTX/ ROTHR(NMAX100),ROTKT(NMAX100) C LOGICAL JAVX C DATA Q/NMAX100*0./,RTRHB/NMAX100*0./,D/NMAX100*0./ DATA AVRATE/NMAX50*0./,AVR2/NMAX100*0./,AROW/NMAX100*0./ DATA ANORM/NMAX100*0./,PROB/NMAX25*0./,ROTHR/NMAX100*0./ DATA CRF/NMAX100*0./,RHO/NMAX100*0./ END C C************************************************************************** C DOUBLE PRECISION FUNCTION DFLT(I) C DFLT = DBLE(FLOAT(I)) RETURN END