BLOCK DATA H2HE C C ============================================================================ C Copyright (C) 1991 Aleksandra Borysow C ============================================================================ C C Copyright Notice: C You may use this program (H2HEOVER) for your scienftific applications, C but please do not distribute it yourself. C Direct all your inquires to the author: e-mail aborysow@nbi.dk. C Final request: If you publish your work which benefited from this C program, please acknowledge using this program and QUOTE C the original papers describing the models. C ============================================================================ C C Program written by Aleksandra Borysow; c Previously at: c Michigan Technological University - Physics Department c Houghton, MI 49931 c e-mail: aborysow@phy.mtu.edu c Currently at: c Niels Bohr Institute for Astronomy, Physics and Geophysics, c University Observatory, c Juliane Maries Vej 30, c DK-2100 Copenhagen, c Denmark c e-mail: aborysow@nbi.dk c ***********************************************************yy c c Program compatible with two papers published in c The Astrophysical Journal; c A. Borysow, L. Frommhold and M. Moraldi: Ap.J., 336, 495-503 (1989) c (fundamental band, 18-7000K ); c and c A. Borysow and L. Frommhold: Ap. J. 341, 549-555 (1989) c (overtones + hot bands) c c Covered range of temperatures: 20-7000K c C COMPUTES THE ROTO/VIBRATIONAL CIA SPECTRUM OF H2-HE FOR C V (=NVIB_INI) --> V'(=NVIB_FIN) TRANSITION, V < OR = V', V,V'={0,1,2,3} C I.E. FOR 0-1, 0-2, 0-3, 1-1, 1-2, 1-3, 2-2, 2-3, 3-3. c NOTE: c files with the parameters for X and Y functions ("parxy..") c for (NVIB_INI --> NVIB_FIN) MUST reside in the same directory C IMPLICIT REAL*8 (A-H,O-Z) COMMON/SPEKTRM/ NF,LIKE, NVIB_FIN, NVIB_INI, FREQLO, DFREQ, 1 NTYPE, FMAX common/temp/temp DATA LIKE / 0 / DATA NTYPE /0/ data FMAX / 5000.d0/ c FMAX is not used any more... data temp / 1000.d0/ C NTYPE=0 for equilibrium hydrogen; for "normal" H2 set NTYPE=1 ! c fmax - the limit od integration for moments and for alfa(omega) END PROGRAM H2HEOVER C Generates detailed listings of the roto-vibrational spectra C Try to run it and input parameters are asked interactively... IMPLICIT REAL*8 (A-H,O-Z) CHARACTER*10 MESSAGE, MARKER, AL1, AL2 CHARACTER*4 PAR1, PAR2, PAR3 PARAMETER (ITE=801) COMMON /PRTSUM/ Q1,JRANGE1,W1(2),Q2,JRANGE2,W2(2),MARKER COMMON /SPECIT/ LIST,FREQ(ITE),ABSCOEF(ITE),ALFATOT(ITE) COMMON/SPEKTRM/ NF,LIKE, NVIB_FIN,NVIB_INI, FREQLO, Deltastep, 1 NTYPE, FMAX COMMON/TEMP/TEMP COMMON/ALPHA/ALPHA1,GAMMA1 COMMON/YMAX/YMAXXR, YMAXXL C DATA FOR 0--> 1 TRANSITION GIVEN HERE, FOR OTHER TRANSITION READ C THE DATA FROM CORRESPONDING FILES c data "01" from fitting the Q.M. lineshapes data A01A, A01B, A01C, A01D, A01E /0.33040D+03, 1 -.21649D+01, 0.52749D+00, -.67195D-01, 0.30651D-02/ data B01A, B01B, B01C, B01D, B01E /0.54319D+01, 1 0.20547D+01,-.84063D+00, 0.10838D+00, -.48703D-02/ data C01A, C01B, C01C, C01D, C01E /0.83228D+02, 1 -.15051D+01, 0.37391D+00, -.17139D-01, -.93260D-04/ data D01A, D01B, D01C, D01D, D01E /0.40478D-02, 1 0.18367D+01, -.42343D+00, 0.41328D-01, -.16901D-02/ data E01A, E01B, E01C, E01D, E01E / 0.49231D-03, 1 0.56766D+01, -.11506D+01, 0.83468D-01, -.18655D-02/ data F01A, F01B, F01C, F01D, F01E /0.18547D+04, 1 -.12944D+01, -.20196D+00, 0.60123D-01, -.35471D-02/ c ================================================ data A23A, A23B, A23C, A23D, A23E /0.12235D+03, 1 -0.61810, 0.18320D-01, -0.61509D-02, 0.43208D-03/ data B23A, B23B, B23C, B23D, B23E /0.26987D+02, 1 -0.70068D-01, -0.14427D+00, 0.20516D-01, -0.10215D-02/ data C23A, C23B, C23C, C23D, C23E / 0.95270D+01, 1 -0.61553D+00, 0.11515D+00, 0.41458D-02, -0.50238D-03/ data D23A, D23B, D23C, D23D, D23E /0.43532D+01, 1 -0.20351D+01, 0.24114, -0.48096D-03 , -0.99351D-03/ data E23A, E23B, E23C, E23D, E23E / 0.53493D+04, 1 -0.34232D+01, 0.47824D+00 , -0.24050D-01, 0.25984D-04/ data F23A, F23B, F23C, F23D, F23E / 0.12477D+02, 1 0.26576D+00, -0.35979D-01, -0.11942D-01, 0.92788D-03/ data a21_k0,b21_k0,c21_k0,d21_k0,E21_k0,F21_k0/0.65849D-66 , 1 0.14311D+00,-0.20553D+00, 0.79330D-01,-0.79669D-02, 0.24910D-03/ data a21_k1,b21_k1,c21_k1,d21_k1,E21_k1,F21_k1/0.19681D-53 , 1 0.31620D+01,-0.14486D+01, 0.31465D+00,-0.27135D-01, 0.83022D-03/ data a21_k2,b21_k2,c21_k2,d21_k2,E21_k2,F21_k2/0.72730D-39 , 1 0.16919D+01,-0.96719D+00 ,0.24603D+00,-0.22554D-01,0.70624D-03/ data a211,b211,c211,d211,e211,f211/0.65884d-12, -1.10380d0, 1 0.12066d0, 0.00997d0, -0.00491d0, 0.00035d0/ c data for tau1 ("Y", "21") data a212,b212,c212,d212,e212,f212/0.54418d-12, -0.83587d0, 1 0.20806d0, -0.05608d0, 0.00713d0, -0.00034d0/ data a21_F0,b21_F0,c21_F0,d21_F0,E21_F0,F21_F0/0.12524D-64, 1 -0.20424D+00,-0.12240D+00, 0.76295D-01,-0.85239D-02,0.30576D-03/ data a21_F1,b21_F1,c21_F1,d21_F1,E21_F1,F21_F1/0.29219D-51, 1 -0.19336D+00,-0.11755D+00, 0.73404D-01,-0.81282D-02,0.28816D-03/ data a21_F2,b21_F2,c21_F2,d21_F2,E21_F2,F21_F2/0.90029D-39 , 1 0.28128D+01, -0.13524D+01,0.30663D+00,-0.27084D-01,0.86170D-03/ C ============================================================= C STATEMENT FUNCTIONS: C ============================================================= S(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 ) S7(T,A,B,C,D,E,F,G) = A*dexp(b*dlog(T) + C*(dlog(T))**2 + 1 D*(dlog(T))**3 + E*(dlog(T))**4 + F*(dlog(T))**5 + 1 G*(DLOG(T))**6 ) C ============================================================= open(unit=1, file='addspec.dat', status='unknown',form='formatted') print 9234 9234 format(' Program H2HEOVER',//, 1 ' Write the INITIAL vibrational level (I1):') accept 1235, nvib_ini 1235 format(i1) print 1236 1236 format(' Write the FINAL vibrational level (I1):') accept 1235, nvib_fin PRINT 1266 1266 FORMAT(' Frequency STEP in cm^{-1}') ACCEPT 1257, DELTASTEP if (nvib_ini .eq. nvib_fin) go to 7676 C CALCULATE OMEGA_0 C DETERMINING THE PURELY VIBRATIONAL TRANSITION FREQUENCY:(CM^-1) JJ=0 OMEGA_0 = H2ELEV(NVIB_FIN,JJ) - H2ELEV(NVIB_INI,JJ) PRINT 6998, NVIB_INI,NVIB_FIN, OMEGA_0 6998 FORMAT(' VIBRATIONAL TRANSITION FREQUENCY FOR ',I1,' --> ', 1 I1,' IS EQUAL TO',/,F12.5,' CM-1') PRINT 1256 1256 FORMAT(' HOW MANY CM^{-1} ON BOTH SIDES OF THIS FREQUENCY?') ACCEPT 1257, DELTAFREQ 1257 FORMAT (F10.4) FREQLO = OMEGA_0 - DELTAFREQ PRINT 7799, FREQLO 7799 FORMAT(' THE LOWEST FREQ = ', F12.5, ' CM-1',/, 1 ' ADJUST IT TO THE CLOSEST ROUND NUMBER:') ACCEPT 7665, FREQLO 7665 FORMAT(F12.1) FREQMAX = OMEGA_0 + DELTAFREQ NF = INT (FREQMAX - FREQLO)/DELTASTEP IF(NF.GT.ITE) STOP 963 C ==> CHANGE "ITE>801" OR ADJUST DELTASTEP FOR A LARGER ONE! PRINT 1233, FREQLO, FREQMAX, DELTASTEP, NF 1233 FORMAT (' THE SPECTRUM WILL BE COMPUTED FOR FREQUENCIES',/, 1 ' FROM', F12.3,' TO ', F12.3, ' EVERY', F10.1, ' IN ',I4, 1 ' STEPS ') go to 9994 7676 PRINT 1344 1344 FORMAT(' HOW MANY CM^{-1} from 0.0 to:?') accept 7665, freqmax nf = int (freqmax)/deltastep freqlo = 0.0D0 PRINT 1233, FREQLO, FREQMAX, DELTASTEP, NF 9994 LIST = NF DO 10 I=1,NF FREQ(I)=FREQLO+DFLOAT(I-1)*DELTASTEP 10 ALFATOT(I)=0.D0 c NOTE: c files with the parameters for X and Y functions c for (NVIB_INI --> NVIB_FIN MUST reside in your current directory !!! c THESE FILES WILL CONTAIN INFO ABOUT ALL INDUCTION TERMS IN A SPECIFIED C ORDER: "01", "21", "23", "45". C LABEL WILL BE INSERTED AT THE TOP OF EACH AND CHECKED EACH TIME; if(nvib_ini.eq.0 .and. nvib_fin.eq.2) open (unit=2, 1 file='parxy.02', status='old') if(nvib_ini.eq.0 .and. nvib_fin.eq.3) open (unit=2, 1 file='parxy.03', status='old' ) if(nvib_ini.eq.1 .and. nvib_fin.eq.1) open (unit=2, 1 file='parxy.11',status='old') if(nvib_ini.eq.1 .and. nvib_fin.eq.2) open (unit=2, 1 file='parxy.12',status='old') if(nvib_ini.eq.1 .and. nvib_fin.eq.3) open (unit=2, 1 file='parxy.13',status='old') if(nvib_ini.eq.2 .and. nvib_fin.eq.2) open (unit=2, 1 file='parxy.22',status='old') if(nvib_ini.eq.2 .and. nvib_fin.eq.3) open (unit=2, 1 file='parxy.23',status='old') if(nvib_ini.eq.3 .and. nvib_fin.eq.3) open (unit=2, 1 file='parxy.33',status='old') C ============================================================= open(unit=6, file='addspec.out', status='unknown',form='formatted') write(6,2) NVIB_INI, NVIB_FIN 2 FORMAT(' OUTPUT of H2HEOVER.for ',/, 1 ' CIA ROTO/VIBRATIONAL SPECTRUM CORRESPONDING TO A TRANSITION', 1 I3,' ---> ',I3,/) dfreq = freq(2)-freq(1) write(1,333) TEMP, NF, LIKE, NVIB_FIN, NVIB_INI, NTYPE, 1 FREQLO, DFREQ 333 format(' Temp(K),NF,LIKE,NVIB_FIN,NVIB_INI,NTYPE,FREQLO,DFREQ' 1 ,/, ' FREQ(I)=FREQLO+DFLOAT(I-1)*DFREQ',/, 1 ' NTYPE=0 for equil.H2; NTYPE=1 for "normal" H2',/, 1 f10.3, 5i10 ,/, 2F15.4) C ============================================================= CALL PARTFCT (TEMP, Q1, NTYPE, JRANGE1, MARKER, W1) C PARTITION FUNCTION IN THE VIBRATIONAL STATE =0 (BY DEFAULT) ALFA=1./(0.69519*TEMP) C FOR TESTING PURPOSE C ROTO/VIB/TRANSLATIONAL MOMENTS, SUM RULES BUILT IN: ALPHA1 = 0.0 GAMMA1 = 0.0 C******************************************************************** C "01" TERM COMPUTED HERE C******************************************************************** IBGAMA=0 c for (01) I will use K0 + Y c Moments for K0 (==X) first, for "01" term IF( NVIB_INI.EQ.0 .AND. NVIB_FIN .EQ. 1) GO TO 1991 GO TO 1992 C READING THE PARAMETERS A-F FOR FUNCTIONS X AND Y C READ FIRST, WHETHER LABEL IS '' 1992 READ (2, 1456) AL1,AL2,NNN 1456 FORMAT(A10,/,A10,/,I5) 1245 FORMAT(A4, 7E12.5) read (2, 1245) PAR1, X0A,X0B,X0C,X0D,X0E,X0F read (2, 1245) PAR2, X1A,X1B,X1C,X1D,X1E,X1F read (2, 1245) PAR3, X2A,X2B,X2C,X2D,X2E,X2F GF0=S(TEMP,X0A,X0B,X0C,X0D,X0E,X0F)*1.D-65 IF(PAR2.EQ.'G1 ') GF1=S(TEMP,X1A,X1B,X1C,X1D,X1E,X1F)*1.D-52 IF(PAR2.EQ.'TAU1') TAUF1=S(TEMP,X1A,X1B,X1C,X1D,X1E,X1F) IF(PAR3.EQ.'G2 ') GF2=S(TEMP,X2A,X2B,X2C,X2D,X2E,X2F)*1.D-38 IF(PAR3.EQ.'TAU2') TAUF2=S(TEMP,X2A,X2B,X2C,X2D,X2E,X2F) IF( PAR2.EQ. 'TAU1' .AND. PAR3.EQ.'TAU2') GO TO 5266 c WRITE(6, 6577) GF0, GF1, GF2 6577 FORMAT(/,' TEST OF MOMENTS OF K0 "01":',/,3E15.5,/) c WRITE(6, 8854) 8854 FORMAT(' Call to PARAMK0 :',/) CALL PARAMK0(GF0,GF1,GF2,TEMP,tauF1,tauF2) 5266 IF(NVIB_INI.EQ.NVIB_FIN) GO TO 1996 READ (2, 1456) AL1,AL2,NNN read (2, 1245) PAR1, Y0A,Y0B,Y0C,Y0D,Y0E,Y0F read (2, 1245) PAR2, Y1A,Y1B,Y1C,Y1D,Y1E,Y1F read (2, 1245) PAR3, Y2A,Y2B,Y2C,Y2D,Y2E,Y2F C GY0=S(TEMP,Y0A,Y0B,Y0C,Y0D,Y0E,Y0F)*1.D-65 IF(PAR2.EQ.'G1 ') GY1=S(TEMP,Y1A,Y1B,Y1C,Y1D,Y1E,Y1F)*1.D-52 IF(PAR2.EQ.'TAU1') TAUK1=S(TEMP,X1A,X1B,X1C,X1D,X1E,X1F) IF(PAR3.EQ.'G2 ') GY2=S(TEMP,Y2A,Y2B,Y2C,Y2D,Y2E,Y2F)*1.D-38 IF(PAR3.EQ.'TAU2') TAUK2=S(TEMP,X1A,X1B,X1C,X1D,X1E,X1F) IF(PAR2.EQ.'TAU1' .AND. PAR3.EQ.'TAU2') GO TO 1993 c WRITE(6, 1577) GY0, GY1, GY2 1577 FORMAT(/,' TEST OF MOMENTS OF Y "01":',/,3E15.5,/) c write (6,1551) 1551 format(/' Call to PARAM_Y ("01"):') CALL PARAM_Y(GY0,GY1,GY2,TEMP,tauK1,tauK2) tauk1=dabs(tauk1) GO TO 1993 1996 TAUK1=1.D-14 TAUK2=1.D-13 GY0 = 0.D0 GO TO 1993 c use fitted parameters (A,B,C,...) 1991 TAUF1=S(TEMP, A01A, A01B, A01C, A01D, A01E,0.D0)*1.D-14 TAUF2=S(TEMP, B01A, B01B, B01C, B01D, B01E,0.D0)*1.D-14 TAUK1=S(TEMP, E01A, E01B, E01C, E01D, E01E,0.D0)*1.D-14 TAUK2=S(TEMP, F01A, F01B, F01C, F01D, F01E,0.D0)*1.D-14 EPSK =S(TEMP, D01A, D01B, D01C, D01D, D01E,0.D0) GTOTAL=S(TEMP, C01A, C01B, C01C, C01D, C01E,0.D0)*1.D-65 GF0=GTOTAL/(1.D0+EPSK) GY0=GTOTAL* EPSK/(1.D0+EPSK) c WRITE(6, 1999) 1999 FORMAT(' COMPUTING THE "01" TERM',/) 1993 continue c choosing YMAXX here: STEP = 200.D0 XS = 1000.D0 GMAX = GF0*BGAMA0(0.D0, TAUF1, TAUF2, TEMP) 4441 GLOWR = GF0* BGAMA0(XS, TAUF1, TAUF2, TEMP) + 1 GY0*BGAMA_Y(XS, TAUK1, TAUK2, TEMP) XS = XS + STEP IF(GLOWR. GT. GMAX/3000.D0) GO TO 4441 YMAXXR = XS-STEP XS = -1000.D0 3441 GLOWL = GF0*BGAMA0(XS, TAUF1, TAUF2, TEMP) + 1 GY0*BGAMA_Y(XS, TAUK1, TAUK2, TEMP) XS = XS - STEP IF(GLOWL. GT. GMAX/1000.D0) GO TO 3441 YMAXXL = XS + STEP C ibgama will tell in addspec whever to call bgama1 (BC) or bgama0 (K0) CALL ADDSPEC(ibgama,TEMP,LIKE,GF0,tauF1,tauF2,GY0,tauK1,tauK2, 1 -1,-1,0,1, ymaxx , NVIB_FIN, NVIB_INI) DO 22 I=1,NF 22 ALFATOT(I) = ABSCOEF(I) c write(1,334) (abscoef(i), i=1,nf) 334 format( 81(10e12.4,/)) c WRITE (6, 5588) c ***************************************************************** C "21" TERM COMPUTED HERE: c ***************************************************************** C FOR "21" USE (K0 + Y) IBGAMA=0 IF( (NVIB_INI.EQ.0) .AND. (NVIB_FIN.EQ.1)) GO TO 991 C READING THE PARAMETERS A-F FOR FUNCTIONS X AND Y 992 READ (2, 1456) AL1,AL2,NNN read (2, 1245) PAR1, X0A,X0B,X0C,X0D,X0E,X0F read (2, 1245) PAR2, X1A,X1B,X1C,X1D,X1E,X1F read (2, 1245) PAR3, X2A,X2B,X2C,X2D,X2E,X2F GF0=S(TEMP,X0A,X0B,X0C,X0D,X0E,X0F)*1.D-65 IF(PAR2.EQ.'G1 ') GF1=S(TEMP,X1A,X1B,X1C,X1D,X1E,X1F)*1.D-52 IF(PAR2.EQ.'TAU1') TAUF1 = S(TEMP,X1A,X1B,X1C,X1D,X1E,X1F) IF(PAR3.EQ.'G2 ') GF2=S(TEMP,X2A,X2B,X2C,X2D,X2E,X2F)*1.D-38 IF(PAR3.EQ.'TAU2') TAUF2 =S(TEMP,X2A,X2B,X2C,X2D,X2E,X2F) IF(PAR2.EQ.'TAU1' .AND. PAR3.EQ.'TAU2') GO TO 4466 CALL PARAMK0(GF0,GF1,GF2,TEMP,tauF1,tauF2) 4466 IF(NVIB_INI.EQ.NVIB_FIN) GO TO 996 READ (2, 1456) AL1,AL2,NNN read (2, 1245) PAR1, Y0A,Y0B,Y0C,Y0D,Y0E,Y0F read (2, 1245) PAR2, Y1A,Y1B,Y1C,Y1D,Y1E,Y1F read (2, 1245) PAR3, Y2A,Y2B,Y2C,Y2D,Y2E,Y2F GY0=S(TEMP,Y0A,Y0B,Y0C,Y0D,Y0E,Y0F)*1.D-65 IF(PAR2.EQ.'G1 ') GY1=S(TEMP,Y1A,Y1B,Y1C,Y1D,Y1E,Y1F)*1.D-52 IF(PAR2.EQ.'TAU1') TAUK1=S(TEMP,Y1A,Y1B,Y1C,Y1D,Y1E,Y1F) IF(PAR3.EQ.'G2 ') GY2=S(TEMP,Y2A,Y2B,Y2C,Y2D,Y2E,Y2F)*1.D-38 IF(PAR3.EQ.'TAU2') TAUK2=S(TEMP,Y2A,Y2B,Y2C,Y2D,Y2E,Y2F) IF(PAR2.EQ.'TAU1' .AND. PAR3.EQ.'TAU2') GO TO 993 CALL PARAM_Y(GY0,GY1,GY2,TEMP,tauK1,tauK2) tauk1=dabs(tauk1) GO TO 993 996 TAUK1=1.D-14 TAUK2 = 1.D-13 GY0 = 0.0D0 GO TO 993 C HERE: c Moments for K0 (==X) first, for "21" term (0-->1) TRANSITION ONLY 991 GF0=S(temp,A21_F0,B21_F0,C21_F0,D21_F0,E21_F0,F21_F0) GF1=S(temp,A21_F1,B21_F1,C21_F1,D21_F1,E21_F1,F21_F1) GF2=S(temp,A21_F2,B21_F2,C21_F2,D21_F2,E21_F2,F21_F2) CALL PARAMK0(GF0,GF1,GF2,temp,tauF1,tauF2) GY0=S(temp,A21_K0,B21_K0,C21_K0,D21_K0,E21_K0,F21_K0) tauk1=S(temp,a211,b211,c211,d211,e211,f211) tauk2=S(temp,a212,b212,c212,d212,e212,f212) 993 CONTINUE c choosing YMAXX here: STEP = 200.D0 XS = 1000.D0 GMAX = GF0*BGAMA0(0.D0, TAUF1, TAUF2, TEMP) 2441 GLOWR = GF0* BGAMA0(XS, TAUF1, TAUF2, TEMP) + 1 GY0*BGAMA_Y(XS, TAUK1, TAUK2, TEMP) XS = XS + STEP IF(GLOWR. GT. GMAX/1000.D0) GO TO 2441 YMAXXR = XS-STEP XS = -1000.D0 1441 GLOWL = GF0*BGAMA0(XS, TAUF1, TAUF2, TEMP) + 1 GY0*BGAMA_Y(XS,TAUK1, TAUK2, TEMP) XS = XS - STEP IF(GLOWL. GT. GMAX/1000.D0) GO TO 1441 YMAXXL = XS + STEP c write (6,7555) 7555 format(' COMPUTING THE "21" TERM',/ ) C IBGAMA will tell in ADDSPEC whever to call BGAMA1 (BC) or BGAMA0 (K0) CALL ADDSPEC(ibgama,TEMP,LIKE,GF0,tauF1,tauF2,GY0,tauK1,tauK2, 1 -1,-1,2,1, FMAX, NVIB_FIN, NVIB_INI) DO 20 I=1,NF 20 ALFATOT(I)=ABSCOEF(I) + ALFATOT(I) c write(1,334) (abscoef(i), i=1,nf) c WRITE (6, 5588) 5588 FORMAT(1X,/) c ***************************************************************** C "23" TERM COMPUTED HERE: c ***************************************************************** IBGAMA=1 c For (23) use BC + Y IF( NVIB_INI.EQ.0 .AND. NVIB_FIN .EQ. 1) GO TO 2991 GO TO 2992 C READING THE PARAMETERS A-F FOR FUNCTIONS X AND Y 2992 READ (2, 1456) AL1,AL2,NNN read (2, 1245) PAR1, X0A,X0B,X0C,X0D,X0E,X0F read (2, 1245) PAR2, X1A,X1B,X1C,X1D,X1E,X1F read (2, 1245) PAR3, X2A,X2B,X2C,X2D,X2E,X2F GF0=S(TEMP,X0A,X0B,X0C,X0D,X0E,X0F)*1.D-65 IF(PAR2.EQ.'G1 ') GF1=S(TEMP,X1A,X1B,X1C,X1D,X1E,X1F)*1.D-52 IF(PAR2.EQ.'TAU1') TAUF1 = S(TEMP,X1A,X1B,X1C,X1D,X1E,X1F) IF(PAR3.EQ.'G2 ') GF2=S(TEMP,X2A,X2B,X2C,X2D,X2E,X2F)*1.D-38 IF(PAR3.EQ.'TAU2') TAUF2=S(TEMP,X2A,X2B,X2C,X2D,X2E,X2F) if(par2.eq.'tau1' .and. par3.eq.'tau2') go to 5577 CALL PARAMBC(GF0,GF1,GF2,TEMP,tauF1,tauF2) 5577 IF(NVIB_INI.EQ.NVIB_FIN) GO TO 2996 READ (2, 1456) AL1,AL2,NNN read (2, 1245) PAR1, Y0A,Y0B,Y0C,Y0D,Y0E,Y0F read (2, 1245) PAR2, Y1A,Y1B,Y1C,Y1D,Y1E,Y1F read (2, 1245) PAR3, Y2A,Y2B,Y2C,Y2D,Y2E,Y2F GY0=S(TEMP,Y0A,Y0B,Y0C,Y0D,Y0E,Y0F)*1.D-65 IF(PAR2.EQ.'G1 ') GY1=S(TEMP,Y1A,Y1B,Y1C,Y1D,Y1E,Y1F)*1.D-52 IF(PAR2.EQ.'TAU1') TAUK1=S(TEMP,Y1A,Y1B,Y1C,Y1D,Y1E,Y1F) IF(PAR3.EQ.'G2 ') GY2=S(TEMP,Y2A,Y2B,Y2C,Y2D,Y2E,Y2F)*1.D-38 IF(PAR3.EQ.'TAU2 ') TAUK2=S(TEMP,Y2A,Y2B,Y2C,Y2D,Y2E,Y2F) IF(PAR2.EQ.'TAU1' .AND. PAR3.EQ.'TAU2') GO TO 2993 CALL PARAM_Y(GY0,GY1,GY2,TEMP,tauK1,tauK2) tauk1=dabs(tauk1) GO TO 2993 2996 TAUK1 = 1.D-14 TAUK2 = 1.D-13 GY0 = 0.D0 GO TO 2993 C FITTED PARAMETERS USED ("BC+Y") (FOR (0-->1) TRANSITION ONLY) 2991 TAUF1=S(TEMP, A23A, A23B, A23C, A23D, A23E,0.D0)*1.D-14 TAUF2=S(TEMP, B23A, B23B, B23C, B23D, B23E,0.D0)*1.D-14 TAUK1=S(TEMP, E23A, E23B, E23C, E23D, E23E,0.D0)*1.D-14 TAUK2=S(TEMP, F23A, F23B, F23C, F23D, F23E,0.D0)*1.D-14 EPSK=S(TEMP, D23A, D23B, D23C, D23D, D23E,0.D0) GTOTAL=S(TEMP, C23A, C23B, C23C, C23D, C23E,0.D0)*1.D-65 GF0=GTOTAL/(1.D0+EPSK) GY0=GTOTAL* EPSK/(1.D0+EPSK) 2993 CONTINUE c choosing YMAXX here: STEP = 200.D0 XS = 1000.D0 GMAX = GF0*BGAMA1(0.D0, TAUF1, TAUF2, TEMP) 6441 GLOWR = GF0* BGAMA1(XS, TAUF1, TAUF2, TEMP) + 1 GY0*BGAMA_Y(XS, TAUK1, TAUK2, TEMP) XS = XS + STEP IF(GLOWR. GT. GMAX/1000.D0) GO TO 6441 YMAXXR = XS-STEP XS = -1000.D0 7441 GLOWL = GF0*BGAMA1(XS, TAUF1, TAUF2, TEMP) 1 + GY0*BGAMA_Y(XS, TAUK1, TAUK2, TEMP) XS = XS - STEP IF(GLOWL. GT. GMAX/1000.D0) GO TO 7441 YMAXXL = XS + STEP c WRITE(6, 7255) 7255 FORMAT(' Computing the "23" TERM',/ ) CALL ADDSPEC(IBGAMA,TEMP,LIKE,GF0,tauF1,tauF2,GY0,tauK1,tauK2, 1 -1,-1,2,3, FMAX, NVIB_FIN, NVIB_INI) DO 24 I=1,NF 24 ALFATOT(I)=ALFATOT(I)+ABSCOEF(I) c write(1,334) (abscoef(i), i=1,nf) c WRITE (6,5588) C *********************************************************************** C "45" TERM COMPUTED HERE: c ***************************************************************** C NEGLECTED, TOO SMALL. c ***************************************************************** 3997 CONTINUE MESSAGE = 'GRANDTOTAL' CALL ALFAMOM (ALFATOT,TEMP,MESSAGE,6,LIKE,-1,-1,-1,-1) c write(1,334) (alfatot(i), i=1,nf) write(1,9999) (freq(i), alfatot(i), i=1,nf) 9999 format(f10.2, e12.4) WRITE (6, 8867) GAMMA1, ALPHA1 8867 FORMAT(/,' THE ROTO/VIB/TRANSL. MOMENTS UNDER THE TOTAL', 1 ' SPECTRUM (SUM RULES):',/,' GAMMA1=',E12.5,/, 1 ' ALPHA1=',E12.5,/ ) 11 continue STOP END SUBROUTINE ADDSPEC(ibgama,temp,like,gF0,tauF1,tauF2, 1 GY0,tauK1,tauK2,LAMBDA1,LAMBDA2,LAMBDA,LVALUE, 1 FREQCUT, NVIB_FIN, NVIB_INI) C C This program generates detailed listing of CIA ALFA(OMEGA). C If both LAMBDA1 and LAMBDA2 are negative: single transitions; C Double transitions are assumed otherwise. C FREQCUT serves THE purpose to skip interpolation of G(OMEGA) for C OMEGA > FREQCUT. ACTUALLY, YMAXXR AND YMAXXL DO IT. C LIKE=1 for like systems (as H2-H2) (zero else) c IBGAMA=0 (X:==K0); IBGAMA=1 (X:==BC) C IMPLICIT REAL*8 (A-H,O-Z) CHARACTER*10 MESSAGE, MARKER CHARACTER*20 NHOL PARAMETER (ITE=801) COMMON /PRTSUM/ Q1,JRANGE1,W1(2),Q2,JRANGE2,W2(2),MARKER COMMON /SPECIT/ LIST,FREQ(ITE),ABSCOEF(ITE),ALFATOT(ITE) COMMON/ALPHA/ ALPHA1, GAMMA1 COMMON/CHECK_X/ XM1,ICHKX COMMON/CHECK_Y/ YM1,ICHKY COMMON/YMAX/YMAXXR, YMAXXL DIMENSION BETA(3) DATA SCALEF/1.D80/, BOLTZWN/.6950304D0/ DATA HBAR,PI,CLIGHT/1.054588757D-27,3.141592653589797D0, 1 2.9979250D10/ DATA BOLTZK /1.380658D-16/ C C DETERMINING THE PURELY VIBRATIONAL TRANSITION FREQUENCY:(CM^-1) JJ=0 OMEGA_0 = H2ELEV(NVIB_FIN,JJ) - H2ELEV(NVIB_INI,JJ) PRINT 7766, NVIB_INI, NVIB_FIN, OMEGA_0 7766 FORMAT(/,' VIBRATIONAL TRANSITION FREQUENCY FOR ',I1,' -->', 1 I1,' IS EQUAL TO ',/, 1X, F12.6,' CM-1',/) if ( (IBGAMA.eq.0).or.(ibgama.eq.1)) go to 1113 1112 write(6,1124) 1124 format(/, ' Ibgama undefined, STOP 123',/) stop 123 1113 continue XM0 = GF0 YM0 = GY0 ICHX = 0 ICHY = 0 DO 10 I=1,LIST 10 ABSCOEF(I)=0.D0 TWOPIC=2.D0*PI*CLIGHT MESSAGE='NON-COMUL.' IF (LIKE.NE.1) LIKE=0 CALIB=TWOPIC*((4.*PI**2)/(3.*HBAR*CLIGHT))*(2.68676D19)**2 CALIB=CALIB/FLOAT(1+LIKE) BETA(1)=1.D0/(BOLTZWN*TEMP) WRITE (6,270) LAMBDA1, LAMBDA2, LAMBDA, LVALUE, JRANGE1, 1 nvib_ini, nvib_fin, MARKER, W1(1), W1(2) 270 FORMAT (//,' LAMBDA1,LAMBDA2, LAMBDA,LVALUE=',2I3,2X,2I3,/, 1 ' JRANGE1=', I2, 5X, 'NVIB_INITIAL, NVIB_FINAL', 2I5, 7X, A10, 1 ' W(1), W(2)=', 2F10.2, /, 1X,46(1H=) ,/ ) C ROTATIONAL SPECTRUM FOR THE DETAILED LISTING ************ C c IF ((LAMBDA1.LT.0).AND.(LAMBDA2.LT.0)) GO TO 170 C 170 JPLUSL=JRANGE1+LAMBDA C C Molecule interacting with an atom (but not another molecule): C DO 210 I=1,JRANGE1 J=I-1 DO 210 IP=1,JPLUSL JP=IP-1 CGS=CLEBSQR(J,LAMBDA,JP) IF (CGS) 210,210,180 180 P=H2POPL(NVIB_INI,J,TEMP,W1)/Q1 OMEGA1=H2ELEV(NVIB_FIN,JP) - H2ELEV(NVIB_INI,J) c PRINT 444, NVIB_FIN, JP, NVIB_INI, J, OMEGA1 444 FORMAT(' (v=',I2,' J=',I2,'), (v=',I2,' J=',I2,'); OMEGA=', 1 F12.5) FAC=CALIB*P*CGS DO 200 IQ=1,LIST FRQ=FREQ(IQ)-OMEGA1 C IF (DABS(FRQ)-FREQCUT) 190,200,200 IF(FRQ.GT.YMAXXR) GO TO 200 IF(FRQ.LT.YMAXXL) GO TO 200 190 WKF=FREQ(IQ)*(1.-DEXP(-BETA(1)*FREQ(IQ)))*FAC if(Ibgama.eq.0) XX=GF0*BGAMA0(frq,tauF1,tauF2,temp) if(Ibgama.eq.1) XX=GF0*BGAMA1(frq,tauF1,tauF2,temp) IF(GY0.NE.0.0) XX = XX + GY0 * bgama_y(FRQ,tauK1,tauK2,temp) ABSCOEF(IQ)=ABSCOEF(IQ) + XX*WKF 200 CONTINUE 210 CONTINUE C C ALL SINGLE AND DOUBLE TRANSITIONS ACCOUNTED FOR ********* C 220 CALL ALFAMOM (ABSCOEF,TEMP,MESSAGE,6,LIKE, 1 LAMBDA1,LAMBDA2,LAMBDA,LVALUE) C NOW, CHECKING THE SUM RULES: ICHKX = 1 ICHKY = 1 if(Ibgama.eq.0) X=BGAMA0(0.,tauF1,tauF2,temp) if(Ibgama.eq.1) X=BGAMA1(0.,tauF1,tauF2,temp) C I KNOW NOW XM1 (REAL M1 FOR "X" IS XM0 * XM1) X = BGAMA_Y(0.,TAUK1,TAUK2,TEMP) C I NOW KNOW YM1 (REAL M1 FOR "Y" IS YM0 * YM1) c PRINT 993, XM0, XM1*XM0, YM0, YM1*YM0 993 FORMAT(/,' X0, X1, Y0, Y1: ', 4E12.4,/) C CHANGE THE UNITS TO [1/S] OMEGA_0 = OMEGA_0 * 2.*PI*CLIGHT C DEFINE B_ROT IN VIBRATIONAL STATE NVIB_INI: (IN CM^-1) C E(J+1) - E(J) = 2 B (J+1); I CHOOSE J=0, J+1 = 1: B_ROT = (H2ELEV(NVIB_INI,1) - H2ELEV(NVIB_INI,0))/2.D0 IF(LIKE.EQ.0) XLIKE = 2. IF(LIKE.EQ.1) XLIKE = 1. GM0 = (XM0+YM0) GM1 = ( XM0*XM1 + YM0*YM1) XW1 = GM1 + 1 2.*PI*CLIGHT*B_ROT*DFLOAT(LAMBDA*(LAMBDA+1))*GM0 1 + OMEGA_0 * GM0 ALPHA = 2.*PI**2/(3.*HBAR*CLIGHT) * XW1*XLIKE GAMMA = PI * PI/(3.*CLIGHT*BOLTZK*TEMP) * GM0*XLIKE WRITE(6, 1446) LAMBDA1, LAMBDA2, LAMBDA, LVALUE, ALPHA, GAMMA 1446 FORMAT(//,' ROTO/VIB/TRANSL MOMENTS FROM THE SUM RULES',/, 1 ' TERM: ',4I2,/, 1 ' ALPHA1 =', E12.5,/,' GAMMA1 =',E12.5,/) ALPHA1 = ALPHA1 + ALPHA GAMMA1 = GAMMA1 + GAMMA RETURN C 290 FORMAT (1X,10E13.4) 300 FORMAT (1X,10F13.1) END SUBROUTINE PARTFCT (TEMP, Q, NTYPE, JRANGE, MARKER, W) C C NTYPE = 0 for equilibrium hydrogen (or other gas in therm. eq.) C NTYPE = 1 redefines the weights W for "normal" hydrogen with a C ratio of ortho:para H2 of 3:1. C IMPLICIT REAL*8 (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/900.) GO TO 30 J=-1 S=0. IF (NTYPE.EQ.1) GO TO 60 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.996D0) GO TO 40 JRANGE=J+3 c print 8555, jrange 8555 format(' Jrange=', i5) MARKER='EQUILIBRM.' GO TO 90 C C "NORMAL" HYDROGEN REQUIRES REDEFINITION OF W(2): C 60 SEV=0. 70 J=J+1 DS=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+DS SEV=SEV+DS J=J+1 S=S+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) IF (S.LE.0.996*Q) GO TO 70 JRANGE=J+3 SODD=S-SEV W(2)=W(2)*(3.*SEV)/SODD C C DEFINITION OF "NORMAL" HYDROGEN: 3*S(EVEN) = S(ODD) C Q=4.*SEV MARKER='NORMAL (!)' C 90 WRITE (6,310) TEMP,Q,W(1),W(2) 310 FORMAT (/, ' PARTITION SUM OF H2 AT Temperature=',F8.2, 1 8H K IS Q=,E12.4,20X,'W(1)=',F8.3,' W(2)=',f8.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/, INIT /1/ C IF (INIT) 10,20,10 10 INIT=0 W(1) = 1. W(2) = 3. C set W(2) = 0. as A TEMPORARY MODIFICATION TO SIMULATE para-H2 C 20 DDD = H2ELEV(NVIB,J) H2POPL = DFLOAT(2*J+1)*W(1+MOD(J,2))*DEXP(-ddd/(BOLTZWN*T)) RETURN END REAL*8 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 REAL*8 (A-H,O-Z) C 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 ====================================================== C ACCURACY OF FIT: 0.04% !!! 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 C ACCURACY OF FIT: 0.03% !!! 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 C ACCURACY OF FIT: 0.034% !!! 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(ALFATOT,TEMP,MESSAGE,K,LIKE, 1 LAMBDA1,LAMBDA2,LAMBDA,LVALUE) C C Writes out the ALFA (absorption) array and computes moments C IMPLICIT REAL*8 (A-H,O-Z) CHARACTER*10 MESSAGE PARAMETER (ITE=801) COMMON /SPECIT/ LIST,FREQ(ITE),A1(ITE),A2(ITE) DIMENSION ABS2(ITE), XI(1), FX(1), ABSP(ITE), ALFATOT(ITE) 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 TWOPIC= 2.D0*PI*CLIGHT CALIB = TWOPIC*((4.D0*PI**2)/(3.D0*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',F10.1,' CM-1.', 1 ' every',F8.1,' CM-1' 10x, ' Temperature:', F8.2, 'K',/, 1 2X,4I2,5X,' IN UNITS OF CM-1 AMAGAT-2.',4X,A10/) WRITE (K,340) (ALFATOT(I),I=1,LIST) 340 FORMAT ( 90(10E13.5,/)) C C MOMENTS OF THE TRANSL/ROTATIONAL CIA SPECTRUM C CALL SPLINE (LIST,1,2,EPS,FREQ,ALFATOT,XI,FX,ALFA1,NR,ABS2) ALFA1=ALFA1*TWOPIC/CLOSCHM**2 DO 250 I=2,LIST E2=DEXP(HBAR*TWOPIC*FREQ(I)/(BK*TEMP)) DIV=FREQ(I)*(E2-1.)/(E2+1.) ABSP(I)=ALFATOT(I)/DIV 250 CONTINUE FR=(FREQ(2)/FREQ(3))**2 ABSP(1)=(ABSP(2)-FR*ABSP(3))/(1.-FR) CALL SPLINE (LIST,1,2,EPS,FREQ,ABSP,XI,FX,GAMA1,NR,ABS2) GAMA0=GAMA1*TWOPIC/CALIB GAMA1=GAMA1*HBAR/(2.*BK*TEMP*CLOSCHM**2) DO 260 I=1,LIST 260 ABSP(I)=ABSP(I)*FREQ(I)**2 CALL SPLINE (LIST,1,2,EPS,FREQ,ABSP,XI,FX,DELTA1,NR,ABS2) DELTA1=DELTA1*TWOPIC**3/CALIB CORR=-1. WRITE (K,350) TEMP,LVALUE,LAMBDA,GAMA1,CORR,ALFA1,CORR, 1 DELTA1,CORR 350 FORMAT (' MOMENTS OF TRANSL/ROTATIONAL SPECTRUM AT',F6.1, 1 ' K FOR L=',I3,' AND LAMBDA=',I3, 7H EQUAL:,/, 2 22H ZEROTH MOMENT GAMA1 =,E12.4,8H S.CM**5,20X,F10.3,/, 3 22H FIRST MOMENT ALFA1 =,E12.4,12H S**-1.CM**5,16X,F10.3,/, 4/ 22H SECOND MOMENT DELTA =,E12.4,15H ERG*CM**6/S**2,13X, 5 F10.3) WRITE (K,360) GAMA0 360 FORMAT (10H GAMMA 0 =,E12.4,10H ERG*CM**6) RETURN 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 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.333333333D-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 SUBROUTINE PARAMK0(G0,G1,G2,TEMP,TAU1,TAU2) c ***************************************************************** c K0 c ****************************************************************** c computes the parameters for K0 function from 3 first moments implicit double precision (a-h,o-z) DATA HBAR, BOLTZK / 1.05458875D-27, 1.380662D-16 / T=temp c WRITE(6, 17) g0, g1, g2 17 FORMAT(' MOMENTS FOR K0 FUNCTION ', 3e15.5) TAU0 = HBAR/(2.D0 * BOLTZK * T) DELT=(TAU0*G1)**2 -4.*(G1*G1/G0+G1/TAU0-G2)*TAU0*TAU0*G0 if(delt.le.0.d0) go to 999 tau1=(-DSQRT(DELT)-TAU0*G1)/(2.*(G1*G1/G0+G1/TAU0-G2)) if(tau1.le.0.d0) go to 999 tau1=DSQRT(tau1) tau2=TAU0*tau1*G0/(G1*tau1*tau1-TAU0*G0) c WRITE(6, 21) tau0, tau1,tau2 21 FORMAT(' tau0=',e13.7,/, 1 ' TAU1= ' ,E13.7,/,' TAU2= ' ,E13.7) go to 889 999 write (6,1177) delt, tau1 1177 format(' One of the following is negative for K0, ', 1 'K0 skipped ',/,' delt, tau1:', 2e13.4,/,' STOP 111') stop 111 889 return end SUBROUTINE PARAM_Y(G0,G1,G2,TEMP,TAU1,TAU2) c ***************************************************************** c K_function (bgama_y) c ****************************************************************** implicit double precision (a-h,o-z) DATA HBAR, BOLTZK / 1.05458875D-27, 1.380662D-16 / c g0, g1 g2 moments for Y function AA=g1/g0 BB=g2/g0 c WRITE(6, 17) g0, g1, g2, aa, bb 17 FORMAT(' MOMENTS FOR Y FUNCTION ', 3e15.5,/, 1 ' Reduced moments: ', 2e15.5) T = Temp TAU0 = HBAR/(2.D0 * BOLTZK * T) t0 = tau0 X = AA - 1.d0/t0 cc = 4.d0* BB - 3.d0*AA**2 - 6.d0*AA/t0 + 9.d0/t0**2 IF(CC.LT.0.0) GO TO 999 GO TO 998 999 G0=0.0 C SET PARAMETERS FOR ZERO-VALUED Y_FUNCTION (CAN NOT BE DETERMINED) TAU1=1.E-14 TAU2=1.E-14 WRITE(6,222), X, CC, 4.d0* BB, 3.d0*AA**2, 6.d0*AA/t0, 1 9.d0/t0**2 222 FORMAT(/, ' FUNCTION "Y" CAN NOT BE DETERMINED, ', 1 ' TAKEN AS ZERO ',/, ' X= ', E12.5,/,' CC= ', E12.5,/, 1 ' 4.d0* BB, 3.d0*AA**2, 6.d0*AA/t0, 9.d0/t0**2:',4E15.5,/) RETURN 998 Y = Dsqrt(CC) c write(6,9555) TAU0, G0,G1,G2,aa,bb,x,y,CC 9555 format(' Test:',/,' Tau_0 (in PARAM_Y) = ',e12.6,/, 1 ' G0,G1,G2 (Y):', 3e15.5,/, 1 ' AA,BB:', 2e15.5,/,' x,y: ', 2e15.5 ,/, 1 ' CC(y=sqrt(CC))=', e13.6 /) iset=0 c Set #1: if((-X-Y).lt.0.d0) go to 333 iset=1 tau2=dsqrt(2.d0*t0/(-X-Y)) tau1=2.d0*t0/(3.d0*X+Y)/dabs(tau2) c write(6, 226) tau1, tau2 226 formaT(' Solutions for Y: tau1,tau2:', 2e15.5) c Set #2: 333 if( (-X+Y) .lt.0.d0 ) go to 335 iset=2 tau2=dsqrt(2.d0*t0/(-X+Y)) tau1= 2.d0*t0/(3.d0*X-Y)/dabs(tau2) c write(6, 225) tau1, tau2 225 format(' Solutions for Y: tau1,tau2:', 8X,2e15.5) 335 write(6,1999) iset, tau1,tau2 1999 format(' Parameters for Y function (set #',i2,'): ',2e12.5) if (iset.eq.1 .or. iset.eq.2) Return if(iset .eq.0) write(6,994) 994 format(/,' Parameters for Y function not found, STOP 155',/) stop 155 end SUBROUTINE PARAMBC ( G0,G1,G2,TEMP,TAU1,TAU2) c ************************************************************************ c BC c ************************************************************************ implicit double precision (a-h,o-z) DATA HBAR, BOLTZK / 1.05458875D-27, 1.380662D-16 / c WRITE(6, 17) g0, g1, g2 17 FORMAT(' MOMENTS FOR BC FUNCTION ', 3e15.5) T=temp TAU0 = HBAR/(2.D0 * BOLTZK * T) 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) stop 777 122 format (/, ' xxx<0 in PARAMBC, STOP 777',/) 98 TAU1=DSQRT(xxx) TAU2=TTA/TAU1 c write(6,1999) tau0,tau1,tau2 1999 format(' Parameters for BC function (tau0,tau1,tau2):',/, 1 3e13.7) 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) c b1 is a value of a B-C lineshape at fnu omega = twopic * fnu tau0 = hbok /(2.*Temp) bgama_y = dexp( -(tau2-dabs(tau2))/tau1) * tau1 * abs(tau2)/tau0 1 * omega * b1 c has units of bgama1 (BC) 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) 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 END