C PROGRAM TO COMPUTE MATRIX OF RATE CONSTANTS C FOR SIO - H2 AT TEMPERATURES FROM 300 - 1500 K C WRITTEN BY SHELDON GREEN, NASA/GISS, FEB 1995 C INTERNAL DATA FOR GENERALIZED IOS CROSS SECTIONS, QL (ANG**2), C FROM SIO-H2 IOS CALCS USING THE RIGID ROTOR C VERSION OF THE BIENIEK & GREEN, CHEM PHYS LETT 84, 380 (1981) C POTENTIAL AS DESCRIBED IN TURNER ET AL, AP J 399, 114 (1992) C USES ADIABATIC-CORRECTED-SUDDEN METHODS OF DEPRISTO ET AL., C J. CHEM. PHYS. 71, 850 (1979) C REDUCES TO IOS APPROX. WHEN PARAMETER LC=0, XC=1. C LC=0 (DEFAULT) AND XC=1. (DEFAULT) ARE RECOMMENDED FROM C LIMITED COMPARISONS WITH COUPLED STATES VALUES AT 200-300 K C USES "AVERAGE" QL'S AS RECOMMENDED IN CHAPMAN & GREEN, C CHEM PHYS LETT 98, 467 (1983); 112, 436 (1984) C C ************************************************************** C * NAMELIST INPUT (&INPUT ): * C * ============== * C * TEMP (EQUIVALENCED TO T) -- IN DEGREES KELVIN * C * JMAX -- HIGHEST J-VALUE FOR RATES * C * LC, XC -- ECS PARAMETERS * C ************************************************************** C * LIMITATIONS: * C * =========== * C * OWING TO LIMITED QL VALUES INCLUDED I RECOMMEND * C * 300.LE.TEMP.LE.1500 AND JMAX.LE.30 * C * THESE LIMITS ARE CURRENTLY ENFORCED BY THE PROGRAM * C ************************************************************** C * REQUIRED SUBROUTINES: * C * ==================== * C * INTEG - DOES BOLTZMANN AVEARGE * C * EORDR - DOES ENERGY ORDERING (NOT ACTUALLY NEEDED HERE) * C * THREEJ(J1,J2,J3) - THREE-J SYMBOL WITH M1=M2=M3=0 * C ************************************************************** C C DIMENSIONS FOR UP TO MXQL L-VALUES AND NNRG ENERGIES PARAMETER (NNRG=13,MXQL=31,MXQLE=NNRG*MXQL) DOUBLE PRECISION RATES(MXQL),QL(NNRG,MXQL),ENERGY(NNRG),Q0S(NNRG) DOUBLE PRECISION INTEG,KBOLTZ,ERR,LC CHARACTER*8 UNITS(2) EQUIVALENCE(T,TEMP) C LOGICAL LAVG = .TRUE. GETS (GEOMETRIC) AVERAGE OF UP/DOWN IOS C LOGICAL LTHRSH = .TRUE. TO START BOLTZ INTEG AT 0.5*DELTAE LOGICAL LAVG, LTHRSH C DATA LAVG/.TRUE./, LTHRSH/.FALSE./ C DATA UNITS/' CM**3/S','EC '/ DATA RATES/MXQL*0.D0/ DATA URED/1.91/, BE/0.7242324/ C EPS CONTROLS ERROR MESSAGE ON INTEG DATA EPS/0.1D0/ C DATA ENERGY/ 1.0000D+02, * 2.0000D+02, 4.0000D+02, 6.0000D+02, 9.0000D+02, 1.2000D+03, * 1.5000D+03, 1.8000D+03, 2.2000D+03, 2.6000D+03, 3.0000D+03, * 3.5000D+03, 4.0000D+03/ DATA Q0S/ -4.0244D+01, * -3.5766D+01, -3.1865D+01, -3.0963D+01, -3.0039D+01, -2.9540D+01, * -2.9203D+01, -2.8933D+01, -2.8658D+01, -2.8429D+01, -2.8214D+01, * -2.8021D+01, -2.7825D+01/ DATA (QL(I, 1),I=1,NNRG)/ 3.4374D+00, * 1.8111D+00, 1.0359D+00, 8.6692D-01, 6.7818D-01, 5.6596D-01, * 5.0612D-01, 4.5888D-01, 4.0452D-01, 3.6082D-01, 3.3727D-01, * 3.0928D-01, 2.8984D-01/ DATA (QL(I, 2),I=1,NNRG)/ 2.6856D+01, * 1.8933D+01, 1.1750D+01, 9.7400D+00, 7.7939D+00, 6.8150D+00, * 6.2173D+00, 5.7752D+00, 5.4086D+00, 5.1683D+00, 4.9317D+00, * 4.7939D+00, 4.6459D+00/ DATA (QL(I, 3),I=1,NNRG)/ 1.5284D+00, * 1.3893D+00, 1.1391D+00, 9.6026D-01, 8.4336D-01, 7.5271D-01, * 6.6839D-01, 6.1930D-01, 5.7836D-01, 5.2574D-01, 5.0454D-01, * 4.5998D-01, 4.3923D-01/ DATA (QL(I, 4),I=1,NNRG)/ 6.3140D+00, * 8.5980D+00, 6.9056D+00, 5.3959D+00, 4.7656D+00, 4.1431D+00, * 3.7638D+00, 3.5669D+00, 3.2999D+00, 3.1081D+00, 2.9922D+00, * 2.8465D+00, 2.7534D+00/ DATA (QL(I, 5),I=1,NNRG)/ 3.5575D-01, * 5.1447D-01, 9.0782D-01, 8.4945D-01, 6.8411D-01, 6.4798D-01, * 5.9667D-01, 5.3922D-01, 4.8562D-01, 4.4266D-01, 4.2256D-01, * 3.8834D-01, 3.6386D-01/ DATA (QL(I, 6),I=1,NNRG)/ 1.2281D+00, * 3.1735D+00, 5.2068D+00, 4.5485D+00, 3.6427D+00, 3.5632D+00, * 3.2956D+00, 3.0167D+00, 2.8576D+00, 2.7089D+00, 2.5474D+00, * 2.4327D+00, 2.3239D+00/ DATA (QL(I, 7),I=1,NNRG)/ 1.2828D-01, * 4.8499D-01, 1.2665D+00, 1.5630D+00, 1.2582D+00, 1.0729D+00, * 9.9074D-01, 9.0057D-01, 7.8780D-01, 6.9571D-01, 6.4543D-01, * 5.8211D-01, 5.4144D-01/ DATA (QL(I, 8),I=1,NNRG)/ 8.4768D-02, * 4.9756D-01, 2.1189D+00, 3.2222D+00, 3.0205D+00, 2.6478D+00, * 2.7169D+00, 2.6722D+00, 2.4339D+00, 2.3168D+00, 2.2340D+00, * 2.1010D+00, 2.0052D+00/ DATA (QL(I, 9),I=1,NNRG)/ 4.9527D-02, * 9.5276D-02, 6.4777D-01, 1.2825D+00, 1.4685D+00, 1.2060D+00, * 1.0635D+00, 1.0006D+00, 8.9815D-01, 7.9711D-01, 7.1686D-01, * 6.4409D-01, 5.8777D-01/ DATA (QL(I,10),I=1,NNRG)/ 5.0491D-02, * 5.9125D-02, 4.1180D-01, 1.1879D+00, 2.1085D+00, 2.0273D+00, * 1.8368D+00, 1.9111D+00, 1.9749D+00, 1.8496D+00, 1.7578D+00, * 1.7188D+00, 1.6297D+00/ DATA (QL(I,11),I=1,NNRG)/ 4.9286D-02, * 3.1165D-02, 1.6538D-01, 6.1261D-01, 1.2778D+00, 1.3216D+00, * 1.1195D+00, 1.0024D+00, 9.4521D-01, 8.6317D-01, 7.7605D-01, * 6.8785D-01, 6.2115D-01/ DATA (QL(I,12),I=1,NNRG)/ 2.4550D-02, * 1.8846D-02, 3.7765D-02, 2.2112D-01, 9.5139D-01, 1.4782D+00, * 1.4758D+00, 1.3556D+00, 1.4324D+00, 1.5405D+00, 1.4657D+00, * 1.4047D+00, 1.3941D+00/ DATA (QL(I,13),I=1,NNRG)/ 6.5306D-02, * 5.1613D-02, 7.3183D-02, 1.9971D-01, 6.7877D-01, 1.0764D+00, * 1.1023D+00, 9.6423D-01, 8.5514D-01, 8.3500D-01, 7.7553D-01, * 6.8867D-01, 6.1944D-01/ DATA (QL(I,14),I=1,NNRG)/ 1.5235D-02, * 1.9739D-02, 2.5037D-02, 3.0474D-02, 2.4096D-01, 7.2783D-01, * 1.1019D+00, 1.1673D+00, 1.0763D+00, 1.1515D+00, 1.2647D+00, * 1.2490D+00, 1.2015D+00/ DATA (QL(I,15),I=1,NNRG)/ 3.1689D-02, * 4.0211D-02, 5.6628D-02, 8.1414D-02, 2.6257D-01, 6.2024D-01, * 8.9106D-01, 9.3799D-01, 8.1899D-01, 7.4080D-01, 7.2565D-01, * 6.7168D-01, 6.0058D-01/ DATA (QL(I,16),I=1,NNRG)/ 9.1991D-03, * 1.4312D-02, 2.8693D-02, 3.7373D-02, 4.5671D-02, 2.0742D-01, * 5.3143D-01, 8.2682D-01, 9.4924D-01, 8.9655D-01, 9.4362D-01, * 1.1010D+00, 1.1035D+00/ DATA (QL(I,17),I=1,NNRG)/ 9.1054D-03, * 1.7307D-02, 3.7053D-02, 4.9502D-02, 9.8996D-02, 2.6291D-01, * 5.1317D-01, 7.1695D-01, 7.8366D-01, 6.8954D-01, 6.4262D-01, * 6.1261D-01, 5.6847D-01/ DATA (QL(I,18),I=1,NNRG)/ 2.5050D-03, * 4.6920D-03, 1.1617D-02, 2.7570D-02, 3.9390D-02, 5.2568D-02, * 1.5554D-01, 3.6607D-01, 6.6960D-01, 8.0028D-01, 7.6743D-01, * 8.1821D-01, 9.6823D-01/ DATA (QL(I,19),I=1,NNRG)/ 2.2047D-03, * 5.3691D-03, 1.9238D-02, 3.3816D-02, 5.5094D-02, 1.1063D-01, * 2.3246D-01, 4.0176D-01, 5.9382D-01, 6.3109D-01, 5.9143D-01, * 5.2743D-01, 5.0887D-01/ DATA (QL(I,20),I=1,NNRG)/ 1.1290D-03, * 2.7749D-03, 5.5922D-03, 1.0623D-02, 3.1849D-02, 4.7439D-02, * 6.1977D-02, 1.2081D-01, 3.0927D-01, 5.5931D-01, 6.9246D-01, * 6.7875D-01, 7.2326D-01/ DATA (QL(I,21),I=1,NNRG)/ 2.4474D-04, * 8.3620D-04, 6.5375D-03, 1.8147D-02, 3.6475D-02, 6.2325D-02, * 1.1258D-01, 1.9752D-01, 3.4898D-01, 4.7160D-01, 5.1847D-01, * 4.6525D-01, 4.3227D-01/ DATA (QL(I,22),I=1,NNRG)/ 5.3474D-04, * 1.3583D-03, 3.6664D-03, 5.8487D-03, 1.2256D-02, 3.4419D-02, * 5.5655D-02, 7.1745D-02, 1.2348D-01, 2.6707D-01, 4.6842D-01, * 6.2761D-01, 6.1935D-01/ DATA (QL(I,23),I=1,NNRG)/ 1.2464D-04, * 1.0547D-04, 1.1713D-03, 6.0307D-03, 2.0017D-02, 3.9020D-02, * 6.5298D-02, 1.0491D-01, 1.8397D-01, 2.8326D-01, 3.7325D-01, * 3.9369D-01, 3.6227D-01/ DATA (QL(I,24),I=1,NNRG)/ 7.4786D-05, * 3.4566D-04, 1.8576D-03, 5.0629D-03, 4.9682D-03, 1.2919D-02, * 3.5068D-02, 6.0274D-02, 8.6227D-02, 1.2704D-01, 2.2795D-01, * 4.4091D-01, 5.7797D-01/ DATA (QL(I,25),I=1,NNRG)/ 1.6451D-04, * 1.5340D-04, 1.4820D-04, 1.0972D-03, 8.0317D-03, 2.1632D-02, * 3.9470D-02, 6.1914D-02, 1.0206D-01, 1.5614D-01, 2.2227D-01, * 2.8368D-01, 2.8997D-01/ DATA (QL(I,26),I=1,NNRG)/ 4.1954D-05, * 9.0377D-05, 8.0608D-04, 3.3404D-03, 4.0844D-03, 3.9336D-03, * 1.1555D-02, 3.1613D-02, 6.7966D-02, 9.8799D-02, 1.2705D-01, * 2.3224D-01, 4.0952D-01/ DATA (QL(I,27),I=1,NNRG)/ 9.8841D-05, * 1.2035D-04, 1.3137D-04, 1.3085D-04, 2.1858D-03, 9.5674D-03, * 2.1594D-02, 3.6524D-02, 6.0484D-02, 9.0509D-02, 1.2436D-01, * 1.7485D-01, 2.0484D-01/ DATA (QL(I,28),I=1,NNRG)/ 3.6341D-05, * 2.7770D-05, 2.8802D-04, 1.4998D-03, 2.6311D-03, 2.5604D-03, * 2.7441D-03, 8.9593D-03, 3.3960D-02, 7.3852D-02, 1.0567D-01, * 1.4260D-01, 2.2902D-01/ DATA (QL(I,29),I=1,NNRG)/ 4.3478D-05, * 5.2743D-05, 9.6529D-05, 5.6084D-05, 3.7175D-04, 3.1190D-03, * 9.6960D-03, 1.9185D-02, 3.4907D-02, 5.4859D-02, 7.3039D-02, * 1.0324D-01, 1.2915D-01/ DATA (QL(I,30),I=1,NNRG)/ 1.6678D-05, * 7.8750D-06, 7.5765D-05, 4.8197D-04, 9.4247D-04, 1.3646D-03, * 1.2995D-03, 1.7977D-03, 9.8320D-03, 3.5835D-02, 7.5268D-02, * 1.1808D-01, 1.5254D-01/ DATA (QL(I,31),I=1,NNRG)/ 1.1361D-05, * 1.3766D-05, 3.3101D-05, 4.2831D-05, 1.1888D-04, 8.4078D-04, * 3.6106D-03, 8.6813D-03, 1.8213D-02, 3.1919D-02, 4.4278D-02, * 6.3265D-02, 7.9650D-02/ C NAMELIST /INPUT/ LC,T,TEMP,JMAX,XC C C STATEMENT FUNCTION . . . . F(J)=2.*FLOAT(J)+1. C *** C BELOW SUPPRESSES UNDERFLOWS C Replace with an appropriate system-dependent routine, if necessary C CALL XUFLOW(0) C *** C C SET DEFAULT VALUES, READ NAMELIST INPUT, ECHO NQL=MXQL IF (.NOT.LAVG) LTHRSH=.FALSE. LC=0. XC=1. JMAX=0 TEMP=0. 9000 READ(5,INPUT,END=9999) WRITE(6,600) URED,BE,LC,XC,TEMP 600 FORMAT(' PROGRAM TO COMPUTE ADIABATICALLY-CORRECTED-SUDDEN', & ' MATRIX OF RATE COEFFICIENTS.'/ & /' SIO-H2 IOS (S. GREEN, FEB 95) '/ & /' URED =',F12.5,' A.M.U.'/ & /' BE (ROTATION CONSTANT) =',F12.5,' (1/CM)'/ & /' LC (CRITICAL DISTANCE) =',F6.2,' ANGSTROMS'/ & /' XC (LONG-RANGE PARAMETER) =',F6.2/ & /' TEMPERATURE REQUESTED =',F8.2) IF (LAVG) WRITE(6,631) 631 FORMAT(/' (GEOMETRIC) AVERAGE OF UP/DOWN IOS RATES USED.') C IF (TEMP.LT.300. .OR. TEMP.GT.1500.) THEN WRITE(6,699) TEMP 699 FORMAT(/' ***** WARNING. REQUESTED TEMPERATURE =',F8.2/ 1 ' ***** OUTSIDE RANGE RANGE CONSIDERED RELIABLE'/ 2 ' ***** YOU CAN FORCE PROGRAM TO RUN BY MODIFYING CODE') STOP ENDIF C C CONVERSION FACTORS TO GIVE RATES IN CM**3/SEC 9002 KBOLTZ=1.3805E-16 AVOGAD=6.024E23 PI=3.1415926 ERGPCM=KBOLTZ*300./208.4 ERGPCM=ERGPCM*ERGPCM FACT=(PI*URED/AVOGAD)*8.*KBOLTZ**3 TFACT=8.E-16*ERGPCM/SQRT(FACT*TEMP**3) SFACT=(208.4/300.)*TEMP BEKT2=BE*(300./208.4)/(2.*TEMP) C WRITE(6,604) NQL,NNRG,TEMP 604 FORMAT(//' BOLTZMANN AVERAGE OF',I4,' QL AND',I4,' ENERGIES' & ,' AT TEMP =',F8.2) C C ORDER QL ON INCREASING ENERGY CALL EORDR(NNRG,NNRG,ENERGY,NQL,MXQL,QL,Q0S) C C IF Q0S ARE SET, BOLTZMANN AVERAGE THESE FIRST -- NOT USED HERE C DO 2900 I=1,NNRG C IF (Q0S(I).GT.0.) GO TO 2902 C2900 CONTINUE C GO TO 2901 C2902 I0=-999 C WRITE(6,605) I0,(ENERGY(IE),Q0S(IE),IE=1,NNRG) C RAVG=INTEG(NNRG,ENERGY,Q0S,0.,TEMP,ERR,ERR1) C R0S = TFACT*RAVG C RAVG=RAVG/(SFACT*SFACT) C WRITE(6,606) RAVG,R0S C IF (ERR.GT.EPS.OR.ERR1.GT.EPS) WRITE(6,696) ERR,ERR1 C 2901 DO 2000 L=1,NQL WRITE(6,605) L,(ENERGY(IE),QL(IE,L),IE=1,NNRG) 605 FORMAT(/' FOR L=',I4,' ENERGY QL(ANGSTROMS) ARE'/ & (14X,F10.2,F10.4) ) THRESH=0. C OPTIONAL THRESHOLD ADDED 30 NOVEMBER 1984. IF (LTHRSH) THRESH=0.5*BE*FLOAT(L*(L+1)) IF (LTHRSH) WRITE(6,684) THRESH 684 FORMAT(/' THRESHOLD ENERGY IN INTEG =',F20.4) RAVG=INTEG(NNRG,ENERGY,QL(1,L),THRESH,TEMP,ERR,ERR1) RATES(L)=TFACT*RAVG RAVG=RAVG/(SFACT*SFACT) WRITE(6,606) RAVG,RATES(L) 606 FORMAT('0 AVERAGE IS',F12.4,' ANG**2 =',1PE12.4,' CM**3/SEC') IF (ERR.GT.EPS.OR.ERR1.GT.EPS) WRITE(6,696) ERR,ERR1 696 FORMAT( ' *** WARNING. FRACTION IN EXTRAP AT HIGH/LOW ENERGY =' & ,2F8.3) C GET DOWNWARD RATE, R(L->0) IS QL/(2*L+1), DO "LAVG" CORRECTION IF (LAVG) RATES(L)=RATES(L)*EXP(FLOAT(L*(L+1))*BEKT2) 2000 RATES(L)=RATES(L)/F(L) C C OUTPUT RATES(L->0) / / USE ECS TO GET MATRIX OF RATES WRITE(6,601) TEMP,(L,RATES(L),UNITS,L=1,NQL) 601 FORMAT(//' RATES FOR ENERGY CORRECTED ', & 'SUDDEN SCALING AT TEMP =',F8.2/ & (' RATE(',I3,' -> 0 ) =' ,1PE12.4,2A8)) C C JMAX.LT.0 INPUT TO SKIP PROCESSING OF STATE-TO-STATE RATES IF (JMAX.LE.0) THEN WRITE(6,*) WRITE(6,*) ' *** JMAX .LE. 0 ==> RATES NOT COMPUTED.' STOP ELSEIF (JMAX.GT.30) THEN WRITE(6,*) WRITE(6,*) ' *** VALUES MAY NOT BE ACCURATE FOR JMAX =',JMAX WRITE(6,*) ' *** SUGGESTED LIMIT IS JMAX=30' WRITE(6,*) ' *** YOU CAN FORCE A RUN BY MODIFYING CODE' STOP ELSE WRITE(6,602) JMAX 602 FORMAT(/' RATES REQUESTED FOR LEVELS TO JMAX =',I4) ENDIF C C SCALING FACTORS FOR ADIABATIC CORRECTION . . . C OMEGA = 0.13*BE*LC/SQRT(T/URED) * J C FOR UNITS BE IN (1/CM), LC IN ANG., T IN KELVIN, URED IN A.M.U. OFACT=0.13*BE*LC/SQRT(T/URED) OFACT=OFACT*XC OFACT=OFACT*OFACT DO 1000 J=1,JMAX WRITE(6,*) OD=6.+OFACT*FLOAT(J*J) DO 1000 J1=0,J-1 RDOWN=0. LMIN=MIN0(J-J1,NQL) LTOP=MIN0(NQL,J+J1) DO 1001 L=LMIN,LTOP,2 ON=6.+OFACT*FLOAT(L*L) A=ON/OD TJ=THREEJ(J,J1,L) 1001 RDOWN=RDOWN+TJ*TJ*A*A*RATES(L) * F(L) DE=(300./208.4)*BE*(J*(J+1)-J1*(J1+1)) C GET UP RATE BY DETAILED BALANCE RUP=RDOWN* EXP(-DE/TEMP)*F(J) RDOWN=RDOWN*F(J1) 1991 WRITE(6,610) J,J1,RDOWN,RUP 610 FORMAT(' ',I3,' TO',I3,' RATE =',1PE10.2,5X,'UP RATE =',E10.2) 1000 CONTINUE C END OF LOOP OVER J, J1 C 9999 STOP END SUBROUTINE EORDR(NNRG,MXNRG,E,NQL,MXQL,QL,Q0S) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C ROUTINE TO ORDER QL VALUES ON ENERGY DIMENSION E(MXNRG),QL(MXNRG,MXQL),Q0S(MXNRG) IF (NNRG.LE.1) RETURN NM1=NNRG-1 DO 1000 I=1,NM1 IST=I+1 DO 1000 J=IST,NNRG IF (E(I).LE.E(J)) GO TO 1000 T=E(I) E(I)=E(J) E(J)=T DO 1001 L=1,NQL T=QL(I,L) QL(I,L)=QL(J,L) 1001 QL(J,L)=T T=Q0S(I) Q0S(I)=Q0S(J) Q0S(J)=T 1000 CONTINUE RETURN END FUNCTION INTEG(N,E,SIG,TH,T,ERR,ERR1) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION INTEG,KBT DIMENSION E(N),SIG(N) C C BOLTZMAN INTEGRATOR. VERSION OF 26 MAY 93. C INPUT: E(I) ARE ENNERGIES IN CM**-1 C SIG(I) ARE CROSS SECTIONS IN ANGSTROMS**2 C N IS NUMBER OF INPUT E, SIG C T IS TEMPERATURE IN KELVINS C TH IS THRESHOLD ENERGY IN 1/CM C ERR RETURNS THE FRACTION FROM EXTRAPOLATION AT HIGH E. C ERR1 RETURNS THE FRACTION FROM EXTRAPOLATION AT LOW E C C NUMERICALLY CALCULATES INTEGRAL OF DE*E*EXP(-E/KT)*SIG(E) ----- C IN EACH INTERVAL FITS LINE TO SIG(E)=B+SLOPE*E C LOW ENERGY EXTRAP ASSUMES CONSTANT SIG FROM TH TO E(1) C HIGH ENERGY EXTRAP USES STRAIGHT LINE FROM LAST SEGMENT, UNLESS C THIS GIVES NEGATIVE CONTRIB, IN WHICH CASE IT FITS EXPONENTIAL C TO LAST TWO POINTS. C LOGICAL FIRST LOGICAL LPRT/.FALSE./ C ERR=0. ERR1=0. IF (N.LE.1) GO TO 9876 FIRST=.TRUE. KBT=208.4*T/300. XINT=0. EL=TH XL=-EL/KBT WL=EXP(XL) C SET VALUE AT THRESHOLD TO FIRST SIGMA, I.E. ASSUME CONST AT LOW E SL=SIG(1) DO 1000 I=1,N ER=E(I) DE=ER-EL C SKIP UNTIL ABOVE THRESHOLD. C IF SEVERAL CARDS HAVE SAME ENERFY, THIS AVERAGES FIRST AND LAST. IF (DE.LE.0.0) GO TO 1000 SR=SIG(I) XR=-ER/KBT WR=EXP(XR) C FIT STRAIGHT LINE TO (EL,SL), (ER,SR). . . SLOPE=(SR-SL)/(ER-EL) B=SL-EL*SLOPE C ADD INTEGRAL OVER THIS SEGMENT XADD=(KBT*KBT)*(B*(WR*(XR-1.)-WL*(XL-1.)) 1 -KBT*SLOPE*(WR*(XR*XR-2.*XR+2.)-WL*(XL*XL-2.*XL+2.))) XINT=XINT+XADD IF (.NOT.FIRST) GO TO 1001 ERR1=XINT FIRST=.FALSE. C SET UP LEFT VALUES FOR NEXT POINT. 1001 EL=ER XL=XR WL=WR SL=SR IF (LPRT.OR.XADD.LT.0.) WRITE(6,609) EL,SL,XADD 609 FORMAT(' INTERVAL W/ ER,SR =',2F8.2,' CONTRIB =',F10.2) 1000 CONTINUE C C ------------ HIGH-ENERGY EXTRAPOLATION -------------- C EXTEND LINE FROM FINAL SEGMENT TO INFINITY C NOTE: ON EXIT FROM LOOP, B AND SLOPE ARE SET FROM LAST INTEVAL C AND EL,XL,WL,SL ARE SET TO LAST POINT C SET WR TO EXP(-INFINITY)=ZERO WR=0.D0 DEL=(KBT*KBT)*(B*(WR*(XR-1.)-WL*(XL-1.)) 1 -KBT*SLOPE*(WR*(XR*XR-2.*XR+2.)-WL*(XL*XL-2.*XL+2.))) C IF THIS EXTRAPOLATION GIVES NEGATIVE CONTRIB, FIT EXPONENTIAL IF (DEL.LT.0.D0) THEN DEL=0. NM1=N-1 C BACK UP IF DUPLICATE LAST ENERGIES. 2300 IF (E(NM1).LT.E(N)) GO TO 2200 NM1=NM1-1 GO TO 2300 2200 RAT=SIG(NM1)/SIG(N) IF (RAT.GE.1.) GO TO 2000 C INCREASING SIG SHOULDN'T HAVE GIVEN NEGATIVE EXTRAP CONTRIB. WRITE(6,*) ' *** INTEG. IMPROPER LOGICAL PATH. TERMINATING.' STOP 2000 B=DLOG(RAT)/(E(N)-E(NM1)) XXXXX= B * E(N) C IF XXXXX IS LARGE, THE EXP FUNC. WILL OVERFLOW. THIS ARISES FOR C EITHER LARGE E(N) OR RAPIDLY DECREASING SIG. IN EITHER CASE, C ASSUME DEL = 0. IF (XXXXX.LE.1.5D2) THEN A = SIG(N) * EXP(XXXXX) B=B+1./KBT DEL=A*(1.+E(N)*B)*EXP(-E(N)*B)/(B*B) ENDIF ENDIF C ADD EXTRAPOLATION, CALC. HIGH-E AND LOW-E ERRS, AND RETURN TEMP =XINT + DEL ERR=ABS(DEL/TEMP ) ERR1=ABS(ERR1/TEMP) INTEG=TEMP RETURN C 9876 INTEG=0. C WRITE(6,653) N 653 FORMAT(' ***** INTEG CALLED WITH ONLY',I2,' POINTS -- SET TO 0') ERR=1.0 RETURN END FUNCTION THREEJ (J1,J2,J3) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C C COMPUTATION OF SPECIAL WIGNER 3J COEFFICIENT WITH C VANISHING PROJECTIONS. SEE EDMONDS, P. 50. C C THIS VERSION EVALUATES BINOM AND PARITY IN-LINE C SHOULD IMPROVE EFFICIENCY, ESPECIALLY ON CRAY; C ALSO GIVES IMPROVEMENT ON AMDAHL (SG: 20 DEC 92) C C STATEMENT FUNCTION FOR DELTA ASSOCIATED W/ RACAH AND SIXJ SYMBOLS C DELTA(I,J,K)= SQRT(1.D0/ ( BINOM(I+J+K+1,I+J-K) * C 1 BINOM(K+K+1,I-J+K) * DBLE(K+J-I+1) ) ) C I1=J1+J2+J3 IF (I1-2*(I1/2).NE.0) GO TO 8 1 I2=J1-J2+J3 IF (I2) 8,2,2 2 I3=J1+J2-J3 IF (I3) 8,3,3 3 I4=-J1+J2+J3 IF (I4) 8,4,4 4 I5=I1/2 I6=I2/2 SIGN=1.D0 IF (I5-2*(I5/2).NE.0) SIGN=-SIGN C 7 THREEJ=SIGN*DELTA(J1,J2,J3)*BINOM(I5,J1)*BINOM(J1,I6) C B1,B2 ARE BINOM ASSOCIATED W/ DELTA N=J1+J2+J3+1 M=J1+J2-J3 NM = N-M MNM = MIN(NM,M) IF(MNM.LE.0) THEN B1=1.D0 ELSE FN = N+1 F = 0.D0 B = 1.D0 DO 101 I = 1,MNM F = F+1.D0 C = (FN-F)*B 101 B = C/F B1 = B ENDIF N=J3+J3+1 M=J1-J2+J3 NM = N-M MNM = MIN(NM,M) IF(MNM.LE.0) THEN B2=1.D0 ELSE FN = N+1 F = 0.D0 B = 1.D0 DO 102 I = 1,MNM F = F+1.D0 C = (FN-F)*B 102 B = C/F B2 = B ENDIF DELTA=SQRT(1.D0/(B1*B2*(J3+J2-J1+1))) C B3=BINOM(I5,J1), B4=BINOM(J1,I6) N=I5 M=J1 NM = N-M MNM = MIN(NM,M) IF(MNM.LE.0) THEN B3=1.D0 ELSE FN = N+1 F = 0.D0 B = 1.D0 DO 103 I = 1,MNM F = F+1.D0 C = (FN-F)*B 103 B = C/F B3 = B ENDIF N=J1 M=I6 NM = N-M MNM = MIN(NM,M) IF(MNM.LE.0) THEN B4=1.D0 ELSE FN = N+1 F = 0.D0 B = 1.D0 DO 104 I = 1,MNM F = F+1.D0 C = (FN-F)*B 104 B = C/F B4 = B ENDIF THREEJ=SIGN*DELTA*B3*B4 RETURN 8 THREEJ=0.D0 RETURN END