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