BLOCK DATA H2H2 c c Written by Aleksandra Borysow, c Previously at: Michigan Technological University, Physics Department c Houghton, Mi 49931 - 1295 c e-mail: aborysow@phy.mtu.edu. c Currently at: Niels Bofr Institute for Astronomy, Physics and Geophysics, c University Observatory, c Juliane Maries vej 30, c DK-2100 Copenhagen, Denmark c e-mail: aborysow@nbi.dk C ============================================================================ C Copyright (C) 1991 Aleksandra Borysow C ============================================================================ C Copyright Notice: C You may use this program for your scientific applications, C but please do not distribute it yourself. C Direct all your inquires to the author: e-mail aborysow@phy.mtu.edu C Final request: If you publish your work which benefited from this C program, please acknowledge using this program and QUOTE C the original paper describing the procedure: c A. Borysow & L. Frommhold; Ap.J.Letters, 348, L41-L43, (1990) C ============================================================================ c For other details refer to the paper: c W. Meyer, A. Borysow, L. Frommhold, c on H2-H2 fundamental band, Phys.Rev.A, 40, 6931, 1989. c ***************************************************************** c ================================================ c Compile with a G_fl qualifier on VMS machine which handles large c exponentials, double precision assumed c ================================================================ c Program asks interactively for: c TEMPERATURE [K], choose from 600-5000K; c min. frequency (FREQLO) (in cm^[-1]) c max. frequency (FREQMAX); c frequency interval (in cm^{-1}) c c Program computes absorption spectrum at temperatures c between 600 and 5000K, H2-H2 fundamental band, c frequency range roughly 1000 - 8000 cm^{-1}. c Note: it handles 0-->1 vibrational transition only!!! c ***************************************************************** c At this time it will handle thermal hydrogen only C IMPLICIT double precision (A-H,O-Z) COMMON/SPEKTRM/ LIKE, nvib_F,nvib_I, FREQLO, DFREQ, NTYPE, 1 FMAX common /wagi/ w(2) DATA LIKE, nvib_F, nvib_I / 1, 1, 0 / c nvib_F - upper state, nvib_I - lower state DATA NTYPE / 0 / DATA W(1), W(2) /1., 3./ ! rotational weights DATA FMAX / 5000.d0/ ! to be kept this way C NTYPE=0 for equilibrium hydrogen; c NF - sets the maximum number of points where the spectrum can be c requested, it can be changed as long as a statement c "NF=" is changed in ALL common blocks! END PROGRAM ADDEM IMPLICIT double precision (A-H,O-Z) CHARACTER*10 MESSAGE, MARKER CHARACTER*11 labeta c PARAMETER NF = 600 c Adjust NF in all subroutines the same; take dimension of the c arrays equal to number of needed frequency points; c common/ integ / rlo, rup, frac, Ninteg common/temp/temp COMMON/YMAX/YMAXXR, YMAXXL COMMON/BC/ G0_BC, TAU1_BC, TAU2_BC COMMON/K0/ G0_K0, TAU10, TAU20 COMMON/Y_BC/ G0_YBC, TAU1_YBC, TAU2_YBC COMMON/Y_EXP/ G0_y, TAU_YEXP COMMON /PRTSUM/ Q1,JRANGE1,W1(2),MARKER common /wagi/ w(2) common/ bb/ labeta, j1, j1p, j2, j2p, xxm0 COMMON /SPECIT/ LIST, FREQ(600), ABSCOEF(600), ALFATOT(600) COMMON/SPEKTRM/ LIKE, nvib_F, nvib_I, FREQLO, DFREQ, NTYPE, 1 FMAX common/ read/ go common/om/omega_0 Z(T,a,b,c,d,e) = A * DEXP( B * dlog(T) + C * (dlog(T))**2 + 1 D * (dlog(T))**3 + E * (dlog(T))**4 ) Z6(T,a,b,c,d,e,f) = A * DEXP( B * dlog(T) + C * (dlog(T))**2 + 1 D * (dlog(T))**3 + E * (dlog(T))**4 + F * (dlog(T))**5 ) NF = 600 write(*,*) ' Which temperature? (choose between 600 and 5000K) ' read (*,*) temp if ( (temp. lt. 600.d0) .or. (temp.gt.5000.d0)) stop 1299 c Temperatures should be within the range 600 - 5000K C CALCULATE OMEGA_0 C DETERMINING THE PURELY VIBRATIONAL TRANSITION FREQUENCY:(CM^-1) JJ=0 OMEGA_0 = H2ELEV(NVIB_F,JJ) - H2ELEV(NVIB_I,JJ) PRINT 6998, NVIB_I,NVIB_F, OMEGA_0 6998 FORMAT(' VIBRATIONAL TRANSITION FREQUENCY FOR ',I1,' --> ', 1 I1,' IS EQUAL TO',/,F12.5,' CM-1') PRINT 7799 7799 FORMAT(' The spectrum will be computed from FREQLO', 1 ' to FREQMAX (cm^{-1})',/,' FREQLO=?') read(*,*) FREQLO 7665 FORMAT(F12.1) write (*,*)' FREQMAX=?' read(*,*) freqmax write(*,*) ' Currently, you are allowed to request',nf,' points' write(*,*) ' Frequency STEP in cm^{-1} = ? ' read(*,*) DELTASTEP list = INT (FREQMAX - FREQLO)/DELTASTEP IF(list.GT.NF) STOP 963 C ==> CHANGE ' NF = ' OR ADJUST DELTASTEP FOR A LARGER ONE! PRINT 1233, FREQLO, FREQMAX, DELTASTEP, list 1233 FORMAT (' THE SPECTRUM WILL BE COMPUTED FOR FREQUENCIES',/, 1 ' FROM', F12.3,' TO ', F12.3, ' EVERY', F10.1, ' IN ',I4, 1 ' STEPS ') write(*,*) write(*,*) ' ==> Final results in file: output.h2h2' DO 10 I=1,list FREQ(I)=FREQLO+DFLOAT(I-1)*DELTASTEP 10 ALFATOT(I)=0. open(unit=6, file='output.h2h2', 1 status='unknown',form='formatted') write(6,*) ' Program computing H2-H2 RV CIA spectrum ' write(6,*) ' in the fundamental band of hydrogen' write(6,*) write(6,*) ' Written by Aleksandra Borysow, refer to published work: ' write(6,777) 777 format (' A.Borysow & L.Frommhold; Ap.J.Letters, ', 1 '348, L41-L43, (1990)' ) write(6,*) write(6, 1212) list, freqlo, deltastep 1212 format(' Absorption will be given in NF=',i4,' points',/, 1 ' From ', f10.2, ' cm-1, every ', f10.2, ' cm-1',/) write(6,*) write(6,*) ' For final results, look at the GRANDTOTAL listing' write(6,*) ' at the end of this output ' c =============================================================== w1(1) = w(1) w1(2) = w(2) c Compute the partition function Q1: CALL PARTFCT (TEMP, Q1, NTYPE, JRANGE1, MARKER, W ) ALFA =1./(0.69519*TEMP) c ***************************************************************** c compute radial distribution function g(R) at temperature TEMP ldist=1 labeta =' call g(R)' rlo = 0.5 rup = 8. Ninteg = 400 frac = 100. CALL MOMTS (ldist, labeta, RLO,RUP,FRAC,Ninteg,TEMP,XXM0) c ***************************************************************** labeta = '<0|B0001|1>' XM0= Z6(temp,0.63902D-67,0.50455D01,-0.14016D01,0.20789D00, 1 -0.13105D-01,0.26577D-03) XM1 = Z6(temp,0.43091D-54,0.94752D01,-0.35582D01,0.62993D00, 1 -0.50446D-01,0.15025D-02) XM2 = Z6(temp,0.18491D-41,0.12692D02,-0.49708D01,0.90444D00, 1 -0.73471D-01,0.22248D-02) c Moments of function Y fitted 00|0001|01 12-JUN-1989 YM0 = Z6(temp,0.31099d-69,0.44458D01,-0.55863D00,-0.82697d-02, 1 0.86826d-02,-0.54053d-03) YM1 = Z6(temp,0.15904d-56,0.67978D01,-0.14943D01,0.16360D00, 1 -0.44924d-02,-0.17136d-03) YM2 = Z6(temp,0.84551d-44,-0.33042D01,0.44784D01,-0.10762D01, 1 0.10815D00,-0.39964D-02) c write (*,*) ' XM2:', xm2 call ReadMom( 2, 4, temp, xm0, xm1, xm2, ym0, ym1, ym2) gf0 = g0_k0 tauf1 = tau10 tauf2 = tau20 GY0 = G0_Y tauk1 = tau_Yexp tauk2 = 0. ! dummy parameter for Y=bgama_exp call range (2, 4) ibgama = 0 ! means, X==K0 iy = 2 ! Y based on BC(1) or EXP(2) freqcut = 5000.d0 CALL ADDSPEC(TEMP,LIKE, 0, 0, 0, 1, freqcut,nvib_F,nvib_I, 2 ibgama, iy, GF0, tauf1,tauf2, GY0, tauk1, tauk2 ) DO 920 I=1,list 920 ALFATOT(I) = 2.* ABSCOEF(I) c * * * * ** * * * * * * * * * * ** * * * * * * * * * * ** * * * * * c 0221 Labeta='<0|B0221|1>' c Moments of function X fitted 00|0221|01 12-JUN-1989 12:29:03 XM0 = Z6(temp,0.43194D-66,0.33857D01,-0.12462D01,0.23933D00, 1 -0.19857D-01,0.61580D-03) xm1=Z6(temp,0.58831D-51,0.21533D01,-0.13399D01,0.32085D00, 1 -0.31226D-01,0.11251D-02) xm2=z6(temp,0.21773D-38,0.13955D01,-0.54992D00,0.14333D00, 1 -0.13238D-01,0.45764D-03) C Moments of function Y fitted 00|0221|01 12-JUN-1989 (M0 & M1 only) YM0 = Z6(temp,0.16252D-68,-0.11550D02,0.76317D01,-0.16525D01, 1 0.15655D00,-0.55169D-02) YM1 = Z6(temp,0.43540D-56,-0.86291D01,0.65347D01,-0.14617D01, 1 0.14257D00,-0.51469D-02) YM2 = Z(temp, 0.3596D-47, 0.2443D+2, -0.7896D+1, 0.1021D+1, 1 -0.4517D-01) call ReadMom( 2, 4, temp, xm0, xm1, xm2, ym0, ym1, ym2) gf0 = g0_k0 tauf1 = tau10 tauf2 = tau20 GY0 = G0_Y tauk1 = tau_Yexp tauk2 = 0. ! dummy parameter for Y=bgama_exp call range (2, 4) ibgama = 0 ! means, X==K0 iy = 2 ! Y based on EXP=2, iy=1 (BC for Y) c write(6, 78) YMAXXR, YMAXXL 78 format(' Found, YMAXXR, YMAXXL:', 2f10.2 ) CALL ADDSPEC(TEMP,LIKE, 0, 2, 2, 1, FMAX, nvib_F,nvib_I, 2 ibgama, iy, GF0, tauf1,tauf2, GY0, tauk1, tauk2 ) DO 24 I=1,list 24 ALFATOT(I)= ALFATOT(I)+ 2.* ABSCOEF(I) c * * * * * * * * * * * * * ** * * * * * * * * * * * labeta='<0|B0223|1>' C Moments of function X fitted 00|0223|01 XM0 = Z(temp, 0.5446D-63, 0.2518 d0, -0.1763d0, 0.3394D-01, 1 -0.1456D-2 ) XM1 = Z(TEMP, 0.2274D-50, 0.5840d0, -0.2414d0, 0.4345D-01, 1 -0.1869D-2 ) XM2 = Z(temp, 0.4188D-38, 0.1598D+1, -0.4015d0, 0.6851D-01, 1 -0.2990D-2) c Moments of function Y fitted 00|0223|01 12-JUN-1989 YM0 = Z6(temp,0.28455d-66,-0.88758D01,0.56035D01,-0.12136D01, 1 0.11552D00, -0.40873d-02) YM1 = Z6(temp,0.34097d-54,-0.56887D01,0.43820D01,-0.98537D00, 1 0.96919d-01,-0.35211d-02) YM2 = Z6(temp,0.11633d-42,-0.24079D02,0.15595D02,-0.33886D01, 1 0.32260D00, -0.11418d-01) call ReadMom( 1, 4, temp, xm0, xm1, xm2, ym0, ym1, ym2) gf0 = g0_bc tauf1 = tau1_bc tauf2 = tau2_bc GY0 = G0_Y tauk1 = tau_Yexp tauk2 = 0. ! dummy parameter for Y=bgama_exp call range ( 1, 4) c write(6, 78) YMAXXR, YMAXXL ibgama = 1 ! means, X==BC iy = 2 ! Y based on BC(1) or EXP(2) CALL ADDSPEC(TEMP,LIKE,0, 2, 2, 3, FMAX, nvib_F,nvib_I, 2 ibgama, iy, GF0, tauf1,tauf2, GY0, tauk1, tauk2 ) DO 124 I=1,list 124 ALFATOT(I)= ALFATOT(I)+ 2.* ABSCOEF(I) c * * * * * * * * * * * * * ** * * * * * * * * * * * labeta = '<0|B2233|1>' C Moments of function X fitted 00|2233|01 XM0 = Z(temp, 0.7788D-64, -0.1366D+1, 0.3287d0, -0.3172D-01, 1 0.1357D-2 ) XM1 = Z(temp, 0.8431D-50, -0.3149D+1, 0.7554d0, -0.7748D-01, 1 0.3380D-2 ) XM2 = Z(temp, 0.8145D-38, -0.2132D+1, 0.6398d0, -0.5909D-01, 1 0.2545D-2) C Moments of function Y fitted 00|2233|01 12-JUN-1989 M0&M1 new YM0 = Z6(temp,0.30931D-66,-0.11140D02,0.60326D01,-0.12367D01, 1 0.11312D00, -0.38684D-02) YM1 = Z6(temp,0.26828D-54,-0.10187D02,0.60208D01,-0.12485D01, 1 0.11563D00,-0.39985D-02) YM2 = Z(temp, 0.2166D-44, 0.1417D+2, -0.4704D+1, 0.6150d0, 1 -0.2706D-01) call ReadMom( 1, 4, temp, xm0, xm1, xm2, ym0, ym1, ym2) gf0 = g0_bc tauf1 = tau1_bc tauf2 = tau2_bc GY0 = G0_Y tauk1 = tau_Yexp tauk2 = 0. ! dummy parameter for Y=bgama_exp call range ( 1, 4) c write(6, 78) YMAXXR, YMAXXL ibgama = 1 ! means, X==BC iy = 2 ! Y based on BC(1) or EXP(2) CALL ADDSPEC(TEMP,LIKE,2, 2, 3, 3, FMAX, nvib_F,nvib_I, 2 ibgama, iy, GF0, tauf1,tauf2, GY0, tauk1, tauk2 ) DO 324 I=1,list 324 ALFATOT(I)= ALFATOT(I)+ 2.* ABSCOEF(I) c * * * * * * * * * * * ** * * * * * * * * * * * labeta='<0|B2023|1>' C Moments of function X fitted 00|2023|01 new 12-JUN-1989 XM0 = Z6(temp,0.82694D-64,0.13580D01,-0.61129D00,0.11585D00, 1 -0.92305D-02, 0.27314D-03) XM1 = Z6(temp,0.10307D-50,0.15066D01,-0.85921D00,0.18654D00, 1 -0.16878D-01, 0.57104D-03) XM2 = Z6(temp,0.94523D-39,0.32469D01,-0.13433D01,0.27780D00, 1 -0.24235D-01, 0.79216D-03) C Moments of function Y fitted 00|2023|01 new 12-JUN-1989 YM0 = Z6(temp,0.23319D-66,-0.34530D01,0.23290D01,-0.50258D00, 1 0.47254D-01, -0.16557D-02) YM1 = Z6(temp,0.19958D-54,-0.13442D01,0.17088D01,-0.39385D00, 1 0.39121D-01,-0.14325D-02) YM2 = Z6(temp,0.15699D-41,-0.21459D02,0.13081D02,-0.27548D01, 1 0.25578D00,-0.88617D-02) call ReadMom( 1, 4, temp, xm0, xm1, xm2, ym0, ym1, ym2) gf0 = g0_bc tauf1 = tau1_bc tauf2 = tau2_bc GY0 = G0_Y tauk1 = tau_Yexp tauk2 = 0. ! dummy parameter for Y=bgama_exp call range ( 1, 4) c write(6, 78) YMAXXR, YMAXXL ibgama = 1 ! means, X==BC iy = 2 ! Y based on BC(1) or EXP(2) CALL ADDSPEC( TEMP,LIKE, 2, 0, 2, 3, FMAX, nvib_F,nvib_I, 2 ibgama, iy, GF0, tauf1,tauf2, GY0, tauk1, tauk2 ) DO 724 I=1,list 724 ALFATOT(I)= ALFATOT(I)+ 2.* ABSCOEF(I) c * * * * * * * * * * * ** * * * * * * * * * * * c For high temperatures, I neglect weak terms like 2211, 0445,4045,2021 c * * * * * * * * * * * ** * * * * * * * * * * ** * * * * * * * * * * * MESSAGE = 'GRANDTOTAL' CALL ALFAMOM (TEMP,MESSAGE,6,LIKE,-1,-1,-1,-1) 11 continue write(6,*) write(6, 1111) (freq(i), alfatot(i), i=1, list) 1111 format( f10.3, e12.4) stop END SUBROUTINE ADDSPEC( TEMP, $ LIKE,LAMBDA1,LAMBDA2,LAMBDA,LVALUE, CUT,nvib_F,nvib_I, 2 ibgama, iy, GF0, tauf1,tauf2, GY0, tauk1, tauk2 ) C C This subroutine generates a detailed listing of CIA ALFA(OMEGA). C FREQCUT (or, CUT) serves two purposes: cut-off for moment integrals, and C to skip interpolation of G(OMEGA) for OMEGA > FREQCUT. c Actually, Ymaxxr and Ymaxxl do it. c IBGAMA=0 means X==K0, IBGAMA=1 (X==BC) c IY will tell which "Y" to take (3 or 2 parameters) c IY = 1 (Y based on BC, BGAMA_Y) c IY=2 ( Y based on EXP, BGAMA_exp ) --> give ONLY tauF1 C LIKE=1 for like systems (as H2-H2) C IMPLICIT double precision (A-H,O-Z) CHARACTER*10 MESSAGE, MARKER CHARACTER*11 labeta character*6 note CHARACTER*20 NHOL c PARAMETER NF = 600 COMMON /PRTSUM/ Q1,JRANGE1,W1(2),MARKER COMMON /SPECIT/ LIST, FREQ(600),ABSCOEF(600), ALFATOT(600) COMMON/YMAX/YMAXXR, YMAXXL common /BB/ labeta, J1, Jp1, J2, Jp2, XM00 common/gam/ gamma(3) common/om/omeg_0 common/ check_x / xm1, ichkx common/ check_y / ym1, ichky common/ integ / rlo, rup, frac, ninteg DATA SCALEF/1.D80/, BOLTZWN/ 0.6950304D0/ DATA HBAR,PI,CLIGHT/1.054588757D-27, 3.141592653589797D0, 1 2.9979250D10/ DATA Boltzk / 1.380658d-16 / DO 10 I=1,LIST 10 ABSCOEF(I)=0.D0 TWOPIC=2.*PI*CLIGHT MESSAGE='NON-COMUL.' CALIB=TWOPIC*((4.*PI**2)/(3.*HBAR*CLIGHT))*(2.68676D19)**2 CALIB=CALIB/DFLOAT(1+LIKE) FREQCUT = CUT ! do not alter the outside parameter CUT BETA1 = 1.D0/(BOLTZWN*TEMP) WRITE (6,270) LAMBDA1,LAMBDA2,LAMBDA,LVALUE,MARKER, 1 W1(1),W1(2) 270 FORMAT (//,' LAMBDA1,LAMBDA2, LAMBDA,LVALUE=',2I3,2X,2I3,8X, 1 /, ' OUTPUT ADDSPEC',10X,A10,'HYDROGEN',2F10.2,/, 2 1X,45(1H=) ) C Molecule interacting with molecule (as opposed to atom): XM00 = GF0 + GY0 ! pure quentum result c Checking, M0 semiclassical, how different from quantum ??? ldist = 0 j1=0 jp1=0 j2=0 jp2=0 CALL MOMTS (ldist, labeta, RLO,RUP,FRAC,Ninteg,TEMP, XM00p) c write(6, 7551) LABETA, XM00, XM00p 7551 format(/,' For ', a11, ' and all J=0, quantum M0=',e12.5,/, 1 12x, ' ... and semiclassical M0 =', e12.5 ,/) JPLUSL=JRANGE1+MAX0(LAMBDA1,LAMBDA2) DO 160 I1=1,JRANGE1 J1=I1-1 DO 160 IP1=1,JPLUSL JP1=IP1-1 CG1S=CLEBSQR(J1,LAMBDA1,JP1) IF (CG1S) 160,160,110 110 P1=H2POPL(nvib_I,J1,TEMP,W1)/Q1 OMEGA1=H2ELEV(nvib_I,JP1)-H2ELEV(nvib_I,J1) C "FIRST" H2 does not undergo vibrational transition DO 150 I2=1,JRANGE1 J2=I2-1 DO 150 IP2=1,JPLUSL JP2=IP2-1 CG2S=CLEBSQR(J2,LAMBDA2,JP2) IF (CG2S) 150,150,120 120 P2=H2POPL(nvib_I,J2,TEMP,W1)/Q1 OMEGA2=H2ELEV(nvib_F,JP2)-H2ELEV(NVIB_i,J2) c "second" H2 does undergo vibrational transition ldist = 0 CALL MOMTS (ldist, labeta, RLO,RUP,FRAC,Ninteg,TEMP, XM0j ) if(j1.eq.jp1) note = 'Single' if(j1.ne.jp1) note = 'Double' c Below a correction made for different M0(J1,J1',J2,J2'): FAC=CALIB*P1*P2*CG1S*CG2S * ( Xm0j/Xm00) c test:OK write(*,*) xm0j, xm00 , j2,jp2,j1,jp1 DO 140 I=1,LIST FRQ=FREQ(I)-OMEGA1-OMEGA2 c IF (DABS(FRQ)-FREQCUT) 130,140,140 if (frq.gt.ymaxxr) go to 140 if (frq.lt.ymaxxl) go to 140 130 WKF=FREQ(I) * (1.-DEXP(- BETA1 * FREQ(I))) * FAC if (ibgama.eq.0) xx = gf0 * bgama0(frq,tauf1,tauf2,temp) if (ibgama.eq.1) xx = gf0 * bgama1(frq,tauf1,tauf2,temp) yy = 0.d0 if( gy0.eq.0.0) go to 444 if(iy.eq.1) yy = gy0 * bgama_y(frq,tauk1,tauk2,temp) if(iy.eq.2) yy = gy0 * bgama_exp(frq,tauk1,temp) 444 gg = xx + yy ABSCOEF(I) = ABSCOEF(I) + gg * WKF 140 CONTINUE 150 CONTINUE 160 CONTINUE CALL ALFAMOM (TEMP,MESSAGE,6,LIKE,LAMBDA1,LAMBDA2,LAMBDA,LVALUE) C 290 FORMAT (1X,10E13.4) 300 FORMAT (1X,10F13.1) RETURN END SUBROUTINE PARTFCT (TEMP, Q, NTYPE, JRANGE, MARKER, W) C C NTYPE = 0 for equilibrium hydrogen C IMPLICIT double precision (A-H,O-Z) CHARACTER*10 MARKER DIMENSION W(2) DATA BOLTZWN /.6950304D0/ C Q=0. J=0 30 DQ=H2POPL(0,J,TEMP,W) + h2popl(1,j,temp,w) + 1 h2popl(2,j,temp,w) + h2popl(3,j,temp,w) + 1 h2popl(4,j,temp,w) + h2popl(5,j,temp,w) Q = Q + DQ J = J + 1 IF (DQ.GT.Q/800.) GO TO 30 J=-1 S=0. 40 J=J+1 DDD = H2POPL(0, J, TEMP, W) + h2popl(1,j,temp,w) + 1 h2popl(2,j,temp,w) + h2popl(3,j,temp,w) + 1 h2popl(4,j,temp,w) + h2popl(5,j,temp,w) S= S + DDD/Q IF (S.LE.0.99D0) GO TO 40 JRANGE=J MARKER='EQUILIBRM ' WRITE (6,310) marker, TEMP,Q,W(1),W(2) 310 FORMAT (1x,a10, ' HYDROGEN',/, 1 ' THE ROTATION PARTITION SUM OF H2 ',/, 1 ' AT TEMPERATURE =', F8.2, 8H K IS Q=, F12.4,/, 1 ' Weights J_even/odd:', 2F10.3,/ ) RETURN END FUNCTION H2POPL (NVIB, J, T, W) C UNNORMALIZED (!) Boltzmann factor of the roto- vibrational levels C of the H2 molecule. (To be normalized by Q = partition sum) C IMPLICIT double precision (A-H,O-Z) DIMENSION W(2) DATA BOLTZWN /.6950304D0/ C 20 DDD = H2ELEV(nvib,J) H2POPL = DFLOAT(2*J+1)*W(1+MOD(J,2))*DEXP(-DDD/(BOLTZWN*T)) RETURN END FUNCTION H2ELEV (NVIB, JROT) C Roto- vibrational energy levels of the H2 molecule C Needed for the computation of the partition sum Q of H2 ****** C IMPLICIT double precision (A-H,O-Z) C c The following results are from LEVELS (levels.lev file) of c Kolos et all. (HSINGLX.FOR function), fitted as function of J(J+1) c Fitting program FITLEVELS.FOR (.for), below only for v=0, 1, 2 C NOTE: ROTATIONAL CONSTANTS B AND D HAVE BEEN FIXED WHILE FITTING THE C ROTATIONAL LEVELS, AFTER S.L.BRAGG, J.W.BRAULT,W.H.SMITH, AP.J.,VOL.263, C P.999-1004, (1982). EH2(I, A, B, D, H, XL, G, P, R)= A + B*DFLOAT(I) - 1 D*1.D-2*DFLOAT(I)**2 + H * 1.D-5 * DFLOAT(I)**3 - 1 XL*1.D-8*DFLOAT(I)**4 + G * 1.D-11*DFLOAT(I)**5 - 1 P * 1.D-14 * DFLOAT(I)**6 + R*1.D-17*DFLOAT(I)**7 C A - THE ROTATIONAL ENERGY (NVIB, J=0) c |E(v=0,J=0) - E(v=1,J=0)| = 4162.14 cm-1 if(nvib.eq.0) h2elev=eh2(jrot*(jrot+1), -0.361132D5, 1 59.33451D0, 4.56510D0, 4.568777D0, 4.366951D0, 2.863721D0, 1 1.051977D0, 0.156903D0) + 36113.d0 c add 36113.d0 to each energy to prevent h2elev/(kT) being too c big at small temperatures c ====================================================== if(nvib.eq.1) h2elev=eh2(jrot*(jrot+1), -0.319511D5, 1 56.3742D0, 4.4050D0, 4.515341D0, 4.614924D0, 3.301549D0, 1 1.32896D0, 0.212922D0) + 36113.d0 if(nvib.eq.2) h2elev=eh2(jrot*(jrot+1),-0.280238D5, 53.482D0, 1 4.28D0, 4.552623D0, 4.957454D0, 3.695921D0, 1.469646D0, 1 0.199078D0) + 36113.d0 if(nvib.eq.3) h2elev=eh2(jrot*(jrot+1), -.243268d+05, 1 0.5050d+02, 0.3818d+01, 0.2393d+01, -.8313d+00, -.4652d+01, 1 -.4660d+01, -.1628d+01) + 36113.d0 if(nvib.eq.4) h2elev=eh2(jrot*(jrot+1), -.208566d+05, 1 0.4762d+02, 0.3605d+01, 0.1511d+01, -.3982d+01, -.1106d+02, 1 -.1108d+02, -.4167d+01) + 36113.d0 if(nvib.eq.5) h2elev=eh2(jrot*(jrot+1), -.176133d+05, 1 0.4481d+02, 0.3511d+01, 0.1099d+01, -.6643d+01, -.1913d+02, 1 -.2174d+02, -.9433d+01) + 36113.d0 if(nvib.gt.5) stop 888 C RETURN END SUBROUTINE ALFAMOM(TEMP,MESSAGE,K,LIKE, 1 LAMBDA1,LAMBDA2,LAMBDA,LVALUE) C C Writes out the ALFA (absorption) array C IMPLICIT double precision (A-H,O-Z) CHARACTER*10 MESSAGE c PARAMETER NF = 600 COMMON /SPECIT/ LIST,FREQ(600),Abscoef(600),Alfatot(600) DIMENSION ABS2(600), XI(1), FX(1), ABSP(600), aa(600), a2(600) DATA EPS, XI(1) /1.D-8, 1.D0/, PI /3.141592653589797D0/ DATA CLIGHT,CLOSCHM /2.9979250D10, 2.68675484D19/ DATA HBAR,BK / 1.054588757D-27, 1.38066244D-16/ C IF(MESSAGE.EQ. 'GRANDTOTAL' ) go to 11 if(message.eq. 'NON-COMUL.') go to 12 stop 987 11 do 22 i=1, list 22 aa(i)=alfatot(i) go to 13 12 do 23 i=1, list 23 aa(i)=abscoef(i) 13 continue TWOPIC= 2.*PI*CLIGHT CALIB = TWOPIC*((4.*PI**2)/(3.*HBAR*CLIGHT))*CLOSCHM**2 1 /DFLOAT(1+LIKE) WRITE (K,330) FREQ(1),FREQ(LIST),FREQ(2)-FREQ(1),TEMP, 1 LAMBDA1,LAMBDA2,LAMBDA,LVALUE,MESSAGE 330 FORMAT (/,' ABSORPTION FROM',F10.1,' CM-1 TO',F7.1,' CM-1.', 1 ' every',F8.2, 10x, /, ' Temperature:', f8.2, /, 1 2X,4I2,5X,' IN UNITS OF CM-1 AMAGAT-2.',4X,A10/) WRITE (K,340) ( AA(I),I=1,LIST) 340 FORMAT ( 200(6E12.5,/)) RETURN END FUNCTION XK0(X) C MODIFIED BESSEL FUNCTION K0(X) C ABRAMOWITZ AND STEGUN P.379 implicit double precision (a-h,o-z) IF(X-2.) 10,10,20 10 T=(X/3.75d0)**2 FI0=(((((.0045813*T+.0360768)*T+.2659732)*T 1 +1.2067492)*T+3.0899424)*T+3.5156229)*T+1. T=(X/2.)**2 P=(((((.00000740*T+.00010750)*T+.00262698)*T 1 +.03488590)*T+.23069756)*T+.42278420)*T+ 2 (-.57721566) X=ABS(X) XK0=-dLOG(X/2.)*FI0+P RETURN 20 T=(2./X) P=(((((.00053208*T-.00251540)*T+.00587872)*T 1 -.01062446)*T+.02189568)*T-.07832358)*T+ 2 1.25331414 X=dMIN1(X,330.d0) XK0=dEXP(-X)*P/dSQRT(X) RETURN END FUNCTION XK1(X) C MODIFIED BESSEL FUNCTION K1(X) TIMES X C PRECISION IS BETTER THAN 2.2E-7 EVERYWHERE. C ABRAMOWITZ AND S,TEGUN, P.379; TABLES P.417. implicit double precision (a-h,o-z) IF(X-2.) 10,10,20 10 T=(X/3.75)**2 FI1=X*((((((.00032411*T+.00301532)*T+.02658733)*T+.15084934) 1 *T+.51498869)*T+.87890594)*T+.5) T=(X/2.)**2 P=(((((-.00004686*T-.00110404)*T-.01919402)*T-.18156897)*T- 1 .67278579)*T+.15443144)*T+1. XK1=X*dLOG(X/2)*FI1+P RETURN 20 T=2./X P=(((((-.00068245*T+.00325614)*T-.00780353)*T+.01504268)*T- 1 .03655620)*T+.23498619)*T+1.25331414 X=dMIN1(X,330.d0) XK1=dSQRT(X)*dEXP(-X)*P RETURN END FUNCTION bgama_y(FNU,TAU1,TAU2,TEMP) C NORMALIZATION SUCH THAT 0TH MOMENT EQUALS UNITY, OK C FNU IS THE FREQUENCY IN CM-1; TEMP IN KELVIN. IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON/CHECK_Y/ G1,ICHY DATA TWOPIC,BKW,HBOK,PI/1.883651568D11,.6950304256D0, 1 7.638280918D-12, 3.141592654D0/ B1 = BGAMA1(FNU,TAU1,TAU2,TEMP) omega = twopic * fnu tau0 = hbok /(2.*Temp) bgama_y = dexp( -(tau2-dabs(tau2))/tau1) * tau1 * abs(tau2)/tau0 1 * omega * b1 IF(ICHY.EQ.1) G1 = TAU0 * DABS(TAU2)/(TAU1 * TAU2**2) + 1 TAU0/TAU2**2 + 1.D0/TAU0 RETURN END FUNCTION BGAMA1(FNU,TAU1,TAU2,TEMP) C BIRNBAUM'S CIA LINE SHAPE MODEL (K1) C NORMALIZATION SUCH THAT 0TH MOMENT EQUALS UNITY C FNU IS THE FREQUENCY IN CM-1; TEMP IN KELVIN. implicit double precision (a-h,o-z) COMMON/CHECK_X/G1,ICHX DATA TWOPIC,BKW,HBOK,PI/1.883651568d11,.6950304256d0, 1 7.638280918d-12, 3.141592654d0/ TAU0 = HBOK /(2.*TEMP) TAU3=DSQRT(TAU2*TAU2+(TAU0)**2) OMEGA=TWOPIC*FNU DENOM=1.d0+(OMEGA*TAU1)**2 X=(TAU3/TAU1)*dSQRT(DENOM) AAA=TAU2/TAU1 AAA=dMIN1(AAA,430.d0) BGAMA1=(TAU1/PI)*dEXP(AAA+FNU/(2.d0*BKW*TEMP))* 1 XK1(X)/DENOM IF(ICHX.EQ.1) G1 = TAU0/(TAU1*TAU2) RETURN END FUNCTION BGAMA0(FNU,TAU1,TAU2,TEMP) C K0 LINE SHAPE MODEL C NORMALIZATION SUCH THAT 0TH MOMENT EQUALS UNITY C FNU IS THE FREQUENCY IN CM-1; TEMP IN KELVIN. implicit double precision (a-h,o-z) COMMON/CHECK_X/G1,ICHX DATA TWOPIC,BKW,HBOK,PI/1.883651568d11,.6950304256d0, 1 7.638280918d-12, 3.141592654d0/ TAU0 = HBOK /(2.*TEMP) TAU4=DSQRT(TAU1*TAU1+(TAU0)**2) OMEGA=TWOPIC*FNU XNOM=1.d0/(TAU2*TAU2)+OMEGA*OMEGA X=TAU4*dSQRT(XNOM) TAU12=TAU1/TAU2 TAU12=DMIN1(TAU12,430.d0) BGAMA0=(TAU1/PI)*DEXP(TAU12+FNU/(2.d0*BKW*TEMP))* 1 XK0(X) IF(ICHX.EQ.1) G1 = TAU0/TAU1 * (1./TAU1 + 1./TAU2) end Subroutine ReadMom( IFUNx, ifuny, temp, xm0, xm1, xm2, 1 ym0, ym1, ym2) C AS FOR NOW, CHOOSE IFUN 1,2,3,4 c c The subroutine is computing the PARAMETERS OF X(omega)+Y(omega) from c two sets of moments, M_(n) (v-->v'), n=0,1,2; c and M_(n) (v'-->v), n=0,1,2 c c returns parameters of X or Y depending on the value of IFUN c implicit double precision (a-h,o-z) CHARACTER*10 LABEL,note COMMON/BC/ G0_BC, TAU1_BC, TAU2_BC COMMON/K0/ G0_K0, TAU10, TAU20 COMMON/Y_BC/ G0_YBC, TAU1_YBC, TAU2_YBC COMMON/Y_EXP/ G0_y, TAU_YEXP DATA HBAR, BOLTZK / 1.05458875D-27, 1.380662D-16 / DATA TWOPIC /1.883651568D11/ IBC=0 IK0=0 IYBC=0 IYEXP = 0 C FOR "X" FUNCTION, 3 PARAMETER LINESHAPE, BC OR K0 if (ifunx.eq.1) IBC=1 if (ifunx.eq.2) IK0=1 C FOR "Y" FUNCTION, 3-PARAMETER: Y(BC) OR 2-PARAMETER: Y_EXP if (ifuny.eq.3) IYBC=1 IF( IFUNy.EQ.4) IYEXP=1 c ============================== T = TEMP TAU0 = HBAR/(2.D0 * BOLTZK * T) IF(IBC.EQ.1) GO TO 111 IF(IK0.EQ.1) GO TO 222 c ************************************************************************ c BC ... for X function (BC+Y) c ************************************************************************ 111 g0 = xm0 g1 = xm1 g2 = xm2 TTA = G0 * TAU0/ G1 XXX = (G2*TTA-G0*(1.+TAU0**2/TTA))/(G0*(TAU0/TTA)**2) if(xxx.le.0.d0) go to 99 go to 98 99 write (6,122) 122 format (' xxx<0 in K1, skipped BC for X') stop 1555 98 TAU1_BC=DSQRT(xxx) TAU2_BC=TTA/TAU1_BC G0_BC = G0 GO TO 5566 c ***************************************************************** c K0... for X function (K0+Y) c ****************************************************************** 222 g0 = xm0 g1 = xm1 g2 = xm2 DELT=(TAU0*G1)**2 -4.*(G1*G1/G0+G1/TAU0-G2)*TAU0*TAU0*G0 if(delt.le.0.d0) go to 999 TAU10=(-DSQRT(DELT)-TAU0*G1)/(2.*(G1*G1/G0+G1/TAU0-G2)) if(tau10.le.0.d0) go to 999 TAU10=DSQRT(TAU10) TAU20=TAU0*TAU10*G0/(G1*TAU10*TAU10-TAU0*G0) G0_K0 = G0 go to 889 999 write (6,1177) delt, tau10 1177 format(' One of the following is negative for K0 (as X)', 1 ' K0 skipped ',/,' delt, tau10:', 2e13.4) stop 9667 889 CONTINUE 5566 if(ifuny.eq.0) go to 444 IF(IYBC.EQ.1) GO TO 333 IF(IYEXP.EQ.1) GO TO 555 c ***************************************************************** c "Y" function (BGAMA_Y) c ****************************************************************** 333 g0 = ym0 g1 = ym1 g2 = ym2 G0_YBC = G0 AA=g1/g0 BB=g2/g0 t0 = tau0 x = AA - 1.d0/t0 DD = 4.d0* BB - 3.d0*AA**2 - 6.d0*AA/t0 +9.d0/t0**2 IF(DD.LE.0.D0) GO TO 991 GO TO 1774 991 WRITE(6,445) 445 FORMAT( ' NO SOLUTIONS FOR Y FUNCTION, DD set to zero!!!') if(dd.le.0.d0) dd = 0.d0 1774 y = Dsqrt(DD) c Set #1: if((-X-Y).lt.0.d0) go to 1333 ISET1=1 tau2=dsqrt(2.d0*t0/(-X-Y)) tau1=2.d0*t0/(3.d0*X+Y)/dabs(tau2) c Set #2: 1333 if( (-X+Y) .lt.0.d0 ) go to 335 ISET2=1 tau2=dsqrt(2.d0*t0/(-X+Y)) tau1= 2.d0*t0/(3.d0*X-Y)/dabs(tau2) 335 continue TAU1_YBC = TAU1 TAU2_YBC = TAU2 go to 444 C *********************************************************** c Y as Y_exp here IYFUN = 4; 2 parameter function c =========================================================== 555 g0 = ym0 g1 = ym1 g2 = ym2 G0_Y = G0 G1_RED=G1/G0 if ( ( G1_RED * TAU0 - 1.).lt.0.d0) go to 866 TAU_Yexp = 2. * TAU0 / DSQRT( G1_RED * TAU0 - 1.) GO TO 966 866 CONTINUE 966 continue c ================================================ 444 RETURN END FUNCTION BGAMA_EXP(FNU,TAU,TEMP) IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON/CHECK_Y/ G1,ICHY C DISPERSION, "Y" TYPE DATA TWOPIC,BKW,HBOK,PI/1.883651568D11,.6950304256, 1 7.638280918d-12, 3.141592654/ data hbar / 1.0545887d-27 / data boltzk / 1.380662d-16 / OMEGA=TWOPIC*FNU TAU0 = HBAR /(2.d0 * BOLTZK *TEMP) SC = 2.* TAU0/TAU**2 BGAMA_EXP = BGAMAE(FNU,TAU,TEMP) * OMEGA/SC if(ichy.eq.1) G1 = 4.*tau0/tau**2 + 1./tau0 RETURN END FUNCTION BGAMAE(FNU,tau ,TEMP) c "EXPONENTIAL" one parameter model lineshape 14-JUL-1988 c with Egelstaff time "X" (BELL) TYPE C NORMALIZATION SUCH THAT 0TH MOMENT EQUALS UNITY C FNU IS THE FREQUENCY IN CM-1; TEMP IN KELVIN. IMPLICIT DOUBLE PRECISION (A-H,O-Z) DATA TWOPIC,BKW,HBOK,PI/1.883651568E11,.6950304256, 1 7.638280918D-12, 3.141592654/ TAU3=DSQRT( TAU**2 + (HBOK/(TEMP*2.))**2 ) OMEGA=TWOPIC*FNU ARG = DABS(OMEGA)*TAU3 BGAMAE = DEXP(FNU/(2.*BKW*TEMP))* DEXP(-ARG) 1 * (TAU)**2 / 2. / TAU3 RETURN END SUBROUTINE RANGE(IFUNX, IFUNY ) implicit double precision (a-h,o-z) COMMON/YMAX/YMAXXR, YMAXXL COMMON/BC/ G0_BC, TAU1_BC, TAU2_BC COMMON/K0/ G0_K0, TAU10, TAU20 COMMON/Y_BC/ G0_YBC, TAU1_YBC, TAU2_YBC COMMON/Y_EXP/ G0_y, TAU_YEXP common/temp/temp C FOR "X" FUNCTION, 3 PARAMETER LINESHAPE, BC OR K0 C if (ifunx.eq.1) IBC=1 C if (ifunx.eq.2) IK0=1 C FOR "Y" FUNCTION, 3-PARAMETER: Y(BC) OR 2-PARAMETER: Y_EXP C if (ifuny.eq.3) IYBC=1 C IF( IFUNy.EQ.4) IYEXP=1 c choosing YMAXX here: if(ifunx.eq.1) gmax =g0_bc* bgama1(0.,tau1_bc,tau2_bc,temp) if(ifunx.eq.2) gmax =g0_k0* bgama0(0.,tau10,tau20,temp) STEP = 200.D0 XS = 1000.D0 4441 if(ifunx.eq.1) gx =g0_bc* bgama1(xs,tau1_bc,tau2_bc,temp) if(ifunx.eq.2) gx =g0_K0* bgama0(xs,tau10,tau20,temp) gy=0.d0 if( ifuny.eq.3 .and. g0_ybc.ne.0.d0 ) gy = 1 g0_ybc*bgama_y(xs,tau1_ybc,tau2_ybc,temp) if( ifuny.eq.4 .and. g0_y.ne.0.d0) gy = 1 g0_y*bgama_exp(xs,tau_yexp,temp) 975 GLOWR = GX + GY XS = XS + STEP if(xs.gt.5000. ) go to 6565 IF(GLOWR. GT. GMAX/5000.D0) GO TO 4441 6565 YMAXXR = XS-STEP if(ymaxxr.lt.2500.) ymaxxr = 3000. XS = -1000.D0 if(ifunx.eq.1) gx =g0_bc* bgama1(0.,tau1_bc,tau2_bc,temp) if(ifunx.eq.2) gx =g0_K0* bgama0(0.,tau10,tau20,temp) 6441 XS = XS - STEP if(xs.le.-5000.) go to 7766 if(ifunx.eq.1) gx =g0_bc* bgama1(xs,tau1_bc,tau2_bc,temp) if(ifunx.eq.2) gx =g0_K0* bgama0(xs,tau10,tau20,temp) gy=0.d0 if(ifuny.eq.3 .and. g0_ybc.ne.0.) gy = 1 g0_ybc*bgama_y(xs,tau1_ybc,tau2_ybc,temp) if(ifuny.eq.4 .and. g0_y.ne.0.) gy = 1 g0_y*bgama_exp(xs,tau_yexp,temp) 978 GLOWL = GX + GY IF(GLOWL. GT. GMAX/2500.D0) GO TO 6441 7766 YMAXXL = XS + STEP return end SUBROUTINE MOMTS (IFIRST, labeta, RLO,RUP,FRAC,N,TEMP, G0) c c If IFIRST = 1 then do only preliminary set-up; set ifirst=0 c when calling moments, c G0 is the zeroth semiclassical moment c Recommended: RLO=0.5, RUP=7.5 (Angstroms) c FRAC = 100., N = 400 C COMPUTES 0-TH, SEMICLASSICAL TRANSLATIONAL SPECTRAL MOMENT C WORKS WITH "DOUBLE" POTENTIAL FUNCTIONS IMPLICIT DOUBLE PRECISION (A-H,O-Z) CHARACTER*10 LGAS, LABEL, NOTE(2) COMMON /COMUNIC/ NOTE,LGAS,LABEL,FFMP,EPSP(2),RMN(2), 1 LVALUE COMMON/EPOT/EPS(2),RM(2),FM,RMIN,RMAX,HHH DATA ANGAU / 0.529917706 / C RLO,RUP - LOWER AND UPPER LIMIT OF INTEGRATION IN ANG C N - NUMBER OF STEPS, FRAC REQUIRED BY DERIVAS, USUALLY SET=100. C FMP - REDUCED MASS IN HYDROGEN MASS, LVALUE = L C LGAS - 10 CHARACTER DESCRIPTION OF SYSTEM if (ifirst.eq.1) go to 222 go to 333 222 BBB=BETACOM(0.) ! call BETA_JJ.for VVV=VACOMUN(0.) ! call VA(R,k) RM(1)=RMN(1)/ANGAU EPS(1)=EPSP(1) RM(2)=RMN(2)/ANGAU EPS(2)=EPSP(2) 10 FORMAT(3F10.5,I5) 11 FORMAT(F10.5) RLO1=RLO/ANGAU RUP1=RUP/ANGAU H=(RUP1-RLO1)/DFLOAT(N) T= TEMP * 1.380662D-23 FMP=FFMP/2.d0 FMPI=2.*FMP*1.67265D-27 return 333 CALL MOMENTS(RLO1,H,N,T,LVALUE,FRAC,FMPI,Gamma, gammap) G0=GAMMA RETURN END SUBROUTINE MOMENTS(RLO,H,N,T,LVALUE,FRAC,FMP,G0, g0p) C COMPUTES ZEROTH MOMENT (G0) FOR ANY L VALUE c G0p is the confirmation of the accuracy; lower order diff. IMPLICIT DOUBLE PRECISION(A-H,O-Z) COMMON/ORDER/Y1P,Y2P,Y3P,Y4P EXTERNAL VTILDE,BETA DATA ANGAU, FOURPI, HBAR/ 0.52917706, 12.566370614359, 1 1.05458876D-34/ SP=0. SQ = 0. R=RLO-H F2=HBAR*HBAR/(12.*FMP*T*T*ANGAU**2*1.D-20) DO 10 I=1,N R=R+H DR=R/FRAC B = beta(R) CALL DERIVAS(VTILDE,R,DR,V,DV,DDV,Z3,Z4,JFLAG) DVP=Y1P DDVP=Y2P G = DEXP(-V/T)*R*R G2=G*(1.+F2*((DV/T-4./R)*DV-2.*DDV)) G2P=G*(1.+F2*((DVP/T-4./R)*DVP-2.*DDVP)) SP=SP+B*B*G2 SQ=SQ+B*B*G2P 10 CONTINUE G0 = FOURPI*SP*H*(ANGAU*1.D-8)**3 G0p = FOURPI*SQ*H*(ANGAU*1.D-8)**3 RETURN END SUBROUTINE DERIVAS(FCT,X0,DX1, Y,Y1,Y2,Y3,Y4,IFLAG) C COMPUTES OTH..4TH DERIVATIVES OF Y=FCT(X) AT X=X0 USING FIVE C ABSCISSAE AND CENTRAL DIFFERENCES. FCT MUST BE DECLARED EXTERNAL C IN CALLING ROUTINE. NOTE THAT FCT(X) MUST BE DEFINED FOR C X0-2DX<=X<=X0+2DX C A TEST OF THE DEFINITION OF THE FIRST DERIVATIVE: IFLAG=1 OUTPUT C SEE ABRAMOWITZ AND STEGUN, P. 914 IMPLICIT DOUBLE PRECISION(A-H,O-Z) COMMON/ORDER/Y1P,Y2P,Y3P,Y4P DIMENSION F(7) 9 X=X0 H=DX1 HSQ=H*H IFLAG=0 DO 10 I=1,7 10 F(I)=FCT(X-DFLOAT(4-I)*H) Y=F(4) A=F(1)+F(7) A1=F(1)-F(7) B=F(2)+F(6) B1=F(2)-F(6) C=F(3)+F(5) C1=F(3)-F(5) C FOLLOWING ARE LOWER ORDER DERIVATIVES FOR TESTING Y1P=(B1-8.*C1)/(12.*H) Y2P=(16.*C-B-30.*Y)/(12.*HSQ) Y3P=(2.*C1-B1)/(2.*H*HSQ) Y4P=(B-4.*C+6.*Y)/(HSQ*HSQ) C THESE ARE HIGHEST ORDER DERIVATIVES Y1 = (-A1+9.*B1-45.*C1)/(60.*H) Y2 = (A-13.5*B+135.*C-245.*Y)/(90.*HSQ) Y3 = (A1-8.*B1+13.*C1)/(8.*HSQ*H) Y4 = (-A+12.*B-39.*C+56.*Y)/(6.*HSQ*HSQ) RETURN END FUNCTION VTILDE(R) IMPLICIT DOUBLE PRECISION(A-H,O-Z) K=0 VTILDE = VA(R,K) RETURN END FUNCTION BETA(R) IMPLICIT REAL*8 (A-H,O-Z) character*11 labeta CHARACTER*10 NOTE,LGAS,LABEL COMMON /BB/ LABETA, J1, J1p, J2, J2p , XM0 COMMON /COMUNIC/ NOTE(2),LGAS,LABEL,FMP,EPSP(2),RMP(2),LVALUE DATA DEBYE / 2.541769713d-18 / E(X) = B0 * DEXP(A*X + B*X**2) E1(X) = B01 * DEXP(A1*X + B1*X**2) E2(X) = B02 * DEXP(A2*X + B2*X**2) if (labeta .eq. '<0|B0001|1>') go to 10 if (labeta .eq. '<0|B2023|1>') go to 20 if (labeta .eq. '<0|B0223|1>') go to 30 if (labeta .eq. '<0|B2233|1>') go to 40 if (labeta .eq. '<0|B2021|1>') go to 50 if (labeta .eq. '<0|B0221|1>') go to 60 if (labeta .eq. '<0|B2211|1>') go to 70 if (labeta .eq. '<0|B4045|1>') go to 80 if (labeta .eq. '<0|B0445|1>') go to 90 print 555, labeta 555 format(' Unpredicted LABEL in BETA:', a11) C stop 999 ! then, what term is that??? 10 continue ! 0001 here LAMDA = 0 LVALUE = 1 LABEL = '00|0001|01' N = 7 BN = -38.20D0 B0 = 0.652d-3 A = -1.521D0 B = -0.033D0 N1=1 ! really 0 but does not matter BN1=0.0d0 B01 = 0.883d-6 A1=-1.517D0 B1 = -0.063D0 N2=7 BN2 =0.161d-1 B02 =-0.318d-7 A2 =-2.65D0 B2 = -0.232D0 go to 100 20 continue ! 2023 here LAMDA = 2 LVALUE =3 LABEL = '00|2023|01' N =4 BN =-0.622D0 B0 =0.811d-5 A =-2.268D0 B =-0.074D0 N1= 4 BN1= -0.152d-1 B01= -0.645d-6 A1= -1.307D0 B1= -0.044D0 N2= 4 BN2= 0.148d-1 B02= 0.667d-6 A2= -1.339D0 B2= -0.036D0 go to 100 30 continue ! 0223 LAMDA = 2 LVALUE =3 LABEL = '00|0223|01' N = 4 BN = 0.853D0 B0 = 0.118d-3 A = -1.479D0 B = -0.026D0 N1 = 4 BN1 = 0.154d-1 B01 = 0.883d-6 A1 = -1.402D0 B1= -0.035D0 N2 = 4 BN2 = -0.149d-1 B02 = -0.616d-6 A2 = -1.311D0 B2 = -0.036D0 go to 100 40 continue ! 2233 here LAMDA = 2 LVALUE =3 LABEL = '00|2233|01' N =4 BN =0.163D0 B0 =-0.107d-4 A =-1.507D0 B =-0.028D0 N1 =4 BN1 =0.244d-2 B01 =-0.122d-6 A1 =-1.557D0 B1 =-0.019D0 N2 =4 BN2 =-0.228d-2 B02 =0.105d-6 A2 =-1.576D0 B2 =-0.016D0 go to 100 50 continue ! 2021 here LAMDA =2 LVALUE = 1 LABEL = '00|2021|01' N = 7 BN =15.5D0 B0 =-0.148d-4 A =-1.219D0 N1 = 1 ! really 0 BN1 = 0.0D0 B01 =0.152d-5 A1 =-1.597D0 B1 = -0.014D0 N2 = 1 BN2 =0.d0 B02 =-0.150d-5 A2 =-1.617D0 B2 = -0.017D0 go to 100 60 continue ! 0221 here LAMDA = 2 LVALUE = 1 LABEL = '00|0221|01' N = 7 BN = -0.161d2 B0 = -0.153d-3 A = -1.806D0 B =-0.130D0 N1 = 1 ! really 0 BN1 = 0.d0 B01 = -.19d-5 A1 = -1.598D0 B1 =-0.018D0 N2 = 1 BN2 = 0.d0 B02 = 0.141d-5 A2 = -1.629D0 B2 = -0.017D0 go to 100 70 continue ! 2211 LAMDA = 2 LVALUE =1 LABEL = '00|2211|01' N = 7 BN = -0.105d1 B0 = 0.862d-5 A = -1.431D0 B = -0.002D0 N1 = 1 BN1 = 0.0d0 B01 = 0.572d-7 A1 = -1.665D0 B1 =-0.070D0 N2 = 7 BN2 = 0.534d-2 B02 = -0.679d-7 A2 = -1.569D0 B2 = -0.022D0 go to 100 80 continue ! 4045 here LAMDA =4 LVALUE = 5 LABEL ='00|4045|01' N = 6 BN = -0.779D0 B0 = 0.316d-6 A = -3.755D0 B = -0.341D0 N1 =6 BN1 = -0.234d-1 B01 = 0. A1 = 0. B1= 0. N2 = 6 BN2 = 0.176d-1 B02 = 0.99d-7 A2 = -1.91D0 B2 = -0.332D0 go to 100 90 continue ! 0445 here LAMDA = 4 LVALUE = 5 LABEL = '00|0445|01' N = 6 BN = 0.220d1 B0 = 0.295d-4 A = -1.430D0 B =-0.131D0 N1 = 6 BN1 = 0.311d-1 B01 = 0. A1 = 0. B1 = 0. N2 = 6 BN2 = -0.17d-1 B02 = -0.602d-7 A2 = -2.227D0 B2 =-.497D0 100 BETAjj = (E1(R-6.D0)+ BN1/R**n1) * dfloat(J2*(J2+1)) 1 + (E2(R-6.D0) + BN2/R**N2) * dfloat(J2p*(J2p+1)) BETA00 = E(R-6.D0) + BN/R**N beta = (betajj + beta00) * debye RETURN ENTRY BETACOM ! the same for all terms BETA = 0.D0 RETURN END FUNCTION VA(R,K) C IMPLICIT DOUBLE PRECISION (A-H,O-Z) CHARACTER*10 NOTE(2), LGAS, LABEL COMMON /EPOT/ EPS(2),RV(2),FM,RMIN,RMAX,HQ COMMON /COMUNIC/ NOTE,LGAS,LABEL,FMP,EPSP(2),RM(2),LVALUE DIMENSION A(2), GAMA(2), ALPHA(2), D(2), RMIN1(2), E(2), 1 C6(2), C8(2), C10(2), C6S(2), C8S(2), C10S(2) K1=K+1 X = R/RV(K1) V_REP = A(K1) * X**GAMA(K1) * DEXP(- ALPHA(K1)*X) F = 1.D0 IF(X .LE. D(K1) ) F = DEXP(- ( D(K1)/X-1.)**2) V_DISP = F * (C6S(K1)/X**6 + C8S(K1)/X**8 + C10S(K1)/X**10) FF = V_REP + V_DISP VA = FF * EPS(K1) RETURN ENTRY VACOMUN FACTOR = (4.803250D-10)**2/(0.52916607D-8) *(1.D-7) FMP=2.D0 LGAS = 'H2-H2 GAS ' C POTENTIAL H2(v=0) - H2(v=0) NOTE(1)='00 - 11 H2' RM(1) = 6.5*0.52917706D0 RMIN1(1)= 6.5D0 E(1)= 0.11163D-3 EPSP(1)= 0.11163D-3 * FACTOR A(1)= 0.37996D7 GAMA(1)= 2.39588D0 ALPHA(1)=14.84169D0 D(1) = 1.02354D0 C6(1)= -12.134D0 C8(1)= -221.26D0 C10(1)= -0.50358D4 C POTENTIAL H2(v=1) - H2(v=0) NOTE(2)='00 - 11 H2' RM(2) = 6.5*0.52917706D0 EPSP(2)= 0.11658D-3 * FACTOR RMIN1(2)=6.5D0 E(2)= 0.11658D-3 A(2) = 0.38000D7 GAMA(2)=2.52498 ALPHA(2)=14.76805 D(2)= 1.03796D0 C6(2)= -12.396D0 C8(2)= -290.878D0 C10(2)= -4370.587D0 DO 10 I=1,2 C6S(I) = C6(I) /( E(I) * RMIN1(I)**6 ) C8S(I) = C8(I) /( E(I) * RMIN1(I)**8 ) 10 C10S(I) = C10(I) /( E(I) * RMIN1(I)**10 ) VA=0.D0 RETURN END SUBROUTINE SPLINE (L,M,K,EPS,X,Y,T,SS,SI,NR,S2) C C Spline interpolation and quadrature, third order after Greville. C Input arguments L...Y, output SS...NR. (Except NR=-1 allows to C read in S2(0) and S2(L) if desired) C L data points X(1), Y(1) ... X(L),Y(L) C EPS=error criterion, typically EPS=1.E-5 for 5 deci. places accuracy C M arguments T(1)..T(M) for which function values SS(1)..SS(M), for C K=0; or first or second derivative for K=1 or -1, respectively. C Note that M has to be at least equal to 1. C SI=integral (over whole interval) for K=2 only. C 'NATURAL' spline functions (if NR.NE.-1) with S2(0)=S2(L)=0 UNLESS C NR was set initially to -1; in that case, read in S2(0),S2(L). C NR indicates the number of out-of-range calls. X(1)@T(I)@X(L) C Extrapolate with caution. (ASSUMPTION D2Y/DX2 = 0.) C S2(I) is the 2nd derivative at X=X(I) and is computed internally. C IMPLICIT REAL*8 (A-H,O-Z) DIMENSION X(L), Y(L), T(M), SS(M), S2(L) DELY(I)=(Y(I+1)-Y(I))/(X(I+1)-X(I)) B(I)=(X(I)-X(I-1))*0.5/(X(I+1)-X(I-1)) C(I)=3.*(DELY(I)-DELY(I-1))/(X(I+1)-X(I-1)) N=L N1=N-1 DO 10 I=2,N1 10 S2(I)=C(I)/1.5 OMEGA=1.0717968 IC=0 C C 'NATURAL' spline functions of third order unless NR=-1: C IF (NR.EQ.-1) GOTO 20 S2(N)=0. S2(1)=0. 20 ETA=0. IC=IC+1 SM=DABS(S2(1)) DO 30 I=2,N IF (DABS(S2(I)).GT.SM) SM=DABS(S2(I)) 30 CONTINUE EPSI=EPS*SM DO 50 I=2,N1 W=(C(I)-B(I)*S2(I-1)-(0.5-B(I))*S2(I+1)-S2(I))*OMEGA IF (DABS(W)-ETA) 50,50,40 40 ETA=DABS(W) 50 S2(I)=S2(I)+W IF (ETA-EPSI) 60,20,20 ENTRY IXPOLAT C C THIS ENTRY USEFUL WHEN ITERATION PREVIOUSLY COMPLETED C N=L N1=N-1 IC=-1 60 IF (K.EQ.2) GO TO 260 NR=0 DO 250 J=1,M I=1 IF (T(J)-X(1)) 110,210,80 80 IF (T(J)-X(N)) 100,190,150 90 IF (T(J)-X(I)) 200,210,100 100 I=I+1 GO TO 90 110 NR=NR+1 HT1=T(J)-X(1) HT2=T(J)-X(2) YP1=DELY(1)+(X(1)-X(2))*(2.*S2(1)+S2(2))/6. IF (K) 140,130,120 120 SS(J)=YP1+HT1*S2(1) GO TO 250 130 SS(J)=Y(1)+YP1*HT1+S2(1)*HT1*HT1/2. GO TO 250 140 SS(J)=S2(I) GO TO 250 150 HT2=T(J)-X(N) HT1=T(J)-X(N1) NR=NR+1 YPN=DELY(N1)+(X(N)-X(N1))*(S2(N1)+2.*S2(N))/6. IF (K) 180,170,160 160 SS(J)=YPN+HT2*S2(N) GO TO 250 170 SS(J)=Y(N)+YPN*HT2+S2(N)*HT2*HT2/2. GO TO 250 180 SS(J)=S2(N) GO TO 250 190 I=N 200 I=I-1 210 HT1=T(J)-X(I) HT2=T(J)-X(I+1) PROD=HT1*HT2 S3=(S2(I+1)-S2(I))/(X(I+1)-X(I)) SS2=S2(I)+HT1*S3 DELSQS=(S2(I)+S2(I+1)+SS2)/6. IF (K) 240,220,230 220 SS(J)=Y(I)+HT1*DELY(I)+PROD*DELSQS GO TO 250 230 SS(J)=DELY(I)+(HT1+HT2)*DELSQS+PROD*S3/6. GO TO 250 240 SS(J)=SS2 250 CONTINUE 260 SI=0. DO 270 I=1,N1 H=X(I+1)-X(I) 270 SI=SI+0.5*H*(Y(I)+Y(I+1))-H**3*(S2(I)+S2(I+1))/24. IF (K.EQ.2) NR=IC RETURN C END FUNCTION CLEBSQR (L,LAMBDA,LP) C SQUARE OF CLEBSCH-GORDAN COEFFICIENT (L,LAMBDA,0,0;LP,0) C FOR INTEGER ARGUMENTS ONLY C NOTE THAT LAMBDA SHOULD BE SMALL, MAYBE @10 OR SO. C IMPLICIT REAL*8 (A-H,O-Z) FC=DFLOAT(2*LP+1) GO TO 10 C ENTRY THREEJ2 C C THIS ENTRY RETURNS THE SQUARED 3-J SYMBOL L LAMBDA LP C 0 0 0 C INSTEAD OF THE CLEBSCH-GORDAN COEFFICIENT C (LIMITATION TO INTEGER ARGUMENTS ONLY) C C NOTE THAT THE THREE-J SYMBOLS ARE COMPLETELY SYMMETRIC IN THE C ARGUMENTS. IT WOULD BE ADVANTAGEOUS TO REORDER THE INPUT ARGUMENT C LIST SO THAT LAMBDA BECOMES THE SMALLEST OF THE 3 ARGUMENTS. C FC=1. 10 CLEBSQR=0. IF (((L+LAMBDA).LT.LP).OR.((LAMBDA+LP).LT.L).OR.((L+LP).LT.LAMBDA) 1) RETURN IF (MOD(L+LP+LAMBDA,2).NE.0) RETURN IF ((L.LT.0).OR.(LP.LT.0).OR.(LAMBDA.LT.0)) RETURN F=1./DFLOAT(L+LP+1-LAMBDA) IF (LAMBDA.EQ.0) GO TO 30 I1=(L+LP+LAMBDA)/2 I0=(L+LP-LAMBDA)/2+1 DO 20 I=I0,I1 20 F=F*DFLOAT(I)/DFLOAT(2*(2*I+1)) 30 P=FC*F*FCTL(LAMBDA+L-LP)*FCTL(LAMBDA+LP-L) CLEBSQR=P/(FCTL((LAMBDA+L-LP)/2)*FCTL((LAMBDA+LP-L)/2))**2 RETURN END FUNCTION FCTL (N) C Simple FACTORIALS program IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER*4 (I-N) P(Z)=((((-2.294720936D-4)/Z-(2.681327160D-3))/Z+(3.472222222D-3))/ 1Z+(8.3333333333D-2))/Z+1. FCTL=1 IF (N.LE.1) RETURN IF (N.GT.15) GO TO 20 J=1 DO 10 I=2,N 10 J=J*I FCTL=DFLOAT(J) RETURN 20 Z=DFLOAT(N+1) FCTL=(DEXP(-Z)*(Z**(Z-0.5))*P(Z)*2.5066282746310) RETURN END