BLOCK DATA INPUT C c ============================================= c Copyright, Aleksandra Borysow, 1988 c ============================================= c If you use this program, please cite original work: c A. Borysow and L. Frommhold, c "Theoretical collision induced rototranslational absorption c spectra for the outer planets: H2-CH4 pairs", c Astrophysical Journal, vol. 304, p.849-865, (1986) c ============================================= C THIS PROGRAM MUST BE ACCOMPANIED BY FILE corrch4 C THESE DATA SERVE AS THE INPUT TO THE MAIN PROGRAM ADDEM. DAT00003 C TEMP = TEMPERATURE IN KELVIN, SHOULD BE BETWEEN 40. AND 300. DAT00004 C FNUMIN = LOWEST FREQUENCY IN CM-1, FOR LISTING OF ALPHA(FNU) DAT00005 C FNUMAX = HIGHEST FREQUENCY IN CM-1, FOR LISTING OF ALPHA(FNU) DAT00006 C DNU = FREQUENCY INCREMENT IN CM-1. DNU SHOULD BE CHOSEN SO DAT00007 C THAT NOT MORE THAN 600 STEPS ARE NEEDED TO GO FROM DAT00008 C FNUMIN TO FNUMAX (ELSE ARRAY DIMENSIONS OF FREQ,ABSCOEF DAT00009 C MUST BE ADJUSTED IN ADDEM). DAT00010 C DAT00011 IMPLICIT DOUBLE PRECISION (A-H,O-Z) character*10 LGAS DAT00012 COMMON /BLOCKIN/ TEMP,FNUMIN,FNUMAX,DNU DAT00013 COMMON /LIKE/ LIKE,LGAS DAT00014 COMMON /APP3/ SLIT,DX,NSRI,WNRMAX3,NS,NSRIUP,JSLIT DAT00015 DATA TEMP/70.D0/ DATA FNUMIN/10.D0/ DATA FNUMAX/400.d0/ DAT00019 DATA SLIT/4.3d0/ DAT00018 DATA DNU/5.0D0/ DAT00020 DATA LIKE/0/ DAT00021 DATA LGAS/' H2 - CH4' / DAT00022 END DAT00024 PROGRAM ADDEM C ADD00026 C ****************************************************************ADD00027 C PROGRAM PREPARED BY ALEKSANDRA BORYSOW, UNIVERSITY OF TEXAS @AUSADD00028 C AND JOINT INSTITUTE FOR LABORATORY ASTROPHYSICS, UNIV. COLORADO ADD00029 C LAST CORRECTION DATE: 14 NOVEMBER 1988 ADD00030 C THIS PROGRAM IS COMPATIBLE WITH PAPER: A. BORYSOW & L. FROMMHOLDADD00031 C ASTROPHYSICAL JOURNAL; VOL. 304, PP.849-865; 1986. ADD00032 C ****************************************************************ADD00033 C PROGRAM GENERATES H2-CH4 FREE-FREE, BOUND-FREE & BOUND-BOUND C COLLISION INDUCED ABSORPTION SPECTRA C ADD00035 C OUTPUT: file output ADD00036 C File corrch4 contains STATISTICAL (NUCLEAR) CORRECTIONS ADD00037 C FOR METHANE INDUCTION ADD00038 C ADD00039 IMPLICIT DOUBLE PRECISION (A-H,O-Z) Character*10 LGAS COMMON /BLOCKIN/ TEMP,FNUMIN,FNUMAX,DNU ADD00040 COMMON /APP3/ SLIT,DX,NSRI,WNRMAX3,NS,NSRIUP,JSLIT ADD00041 COMMON /RSILO/ RSILO(201) ADD00042 COMMON /BOU43/ INITB ADD00043 COMMON /BB/ OMEG,RSI,RSIGG,NSOL,ALFA,SCAL ADD00045 COMMON /BF/ G0BF,DELBF,OM0 ADD00046 COMMON /STAT/ QW3(51,3),QW4(51,4) ADD00047 COMMON /LIKE/ LIKE,LGAS ADD00048 DIMENSION FREQ(601), ABSCOEF(601), ALFATOT(601) ADD00049 DIMENSION RSI(201), RSIGG(201), TT(2), SS(1), OMEG(201), AIG(201) ADD00050 Y(X,A,B,C)=A*DEXP((C*X+B)*X) ADD00051 C ADD00052 OPEN( UNIT=6, FILE='output', STATUS='unknown') OPEN(UNIT=33, FILE='corrch4', STATUS='OLD') NF=INT((FNUMAX-FNUMIN)/DNU+0.5)+1 ADD00053 IF (NF.GT.601) NF=601 ADD00054 FNUMAX=FNUMIN+DFLOAT(NF-1)*DNU ADD00055 WRITE (6,140) LGAS,TEMP,FNUMIN,FNUMAX,DNU,NF-1 ADD00056 CALL PARTSUM (TEMP) ADD00057 C ADD00058 C THE H2 - CH4 SPECTRA FOR 50-300K ADD00059 C ================= ADD00060 C ADD00061 C ONLY B-F + F-F TERMS INCLUDED IN MODELLING ADD00062 C TEMPERATURE INTERPOLATION GOOD FOR 45-300K ADD00063 C ADD00064 X=DLOG(TEMP) ADD00065 DO 10 I=1,NF ADD00066 FREQ(I)=FNUMIN+DFLOAT(I-1)*DNU ADD00067 ALFATOT(I)=0.0 ADD00068 10 ABSCOEF(I)=0. ADD00069 C ADD00070 C READ FROM FILE corrch4 /HERE TAPE 33/ THE QUANTUM CORRECTIONS ADD00071 C (12 Y(J,JP,LAM)-1) ADD00072 C ADD00073 REWIND 33 ADD00074 DO 30 N=1,3 ADD00076 DO 20 I=1,5 ADD00077 K=10*(I-1)+1 ADD00078 K1=10*I ADD00079 READ (33,160) (QW3(J,N),J=K,K1) ADD00080 20 CONTINUE ADD00081 READ (33,150) QW3(51,N) ADD00082 30 CONTINUE ADD00083 DO 50 N=1,4 ADD00084 DO 40 I=1,5 ADD00085 K=(I-1)*10+1 ADD00086 K1=10*I ADD00087 READ (33,160) (QW4(J,N),J=K,K1) ADD00088 40 CONTINUE ADD00089 READ (33,160) QW4(51,N) ADD00090 50 CONTINUE ADD00091 CLOSE(33) C ADD00092 C THESE STATISTICAL FACTORS ARE VALID ONLY FOR DELTA J>0 ADD00093 C FOR DELTA J<0 THESE WILL BE EQUAL TO 1. ADD00094 C HYDROGEN'S QUADRUPOLE ADD00095 C ADD00096 EPS=1.D-5 ADD00097 TT(1)=10.D0 ADD00098 CALL BOUND32 (TEMP,RSI,NSOL) ADD00099 DO 60 I=1,NSOL ADD00100 RSILO(I)=DLOG(RSI(I)*1.D80) ADD00101 60 OMEG(I)=DFLOAT(I-1)*DX ADD00102 CALL SPLINE (NSOL,1,0,EPS,OMEG,RSILO,TT,SS,SI,NR,RSIGG) ADD00103 SCAL=1.D0 ADD00104 S=Y(X,.22916D-59,-1.27134D0, .12423D0) E=Y(X,-.2019D1,-.02259d0,-.05891D0) ADD00106 T1=Y(X,.48254D-13,.6421d0,-.10109D0) ADD00107 T2=Y(X,.97826D-12,-.48654d0,-.0146D0) ADD00108 T3=Y(X,.35247D-12,.10164d0,-.07879D0) ADD00109 T4=Y(X,.13961D-13,1.11146d0,-.09561D0) ADD00110 C ADD00111 C THIS PART FOR MODELING A LOW FREQUENCY PART OF BOUND-FREE ADD00112 C TRANSLATIONAL SPECTRAL FUNCTION, BY A DESYMMETRIZED GAUSSIAN ADD00113 C PROFILE ADD00114 C ADD00116 G0BF=Y(X,0.73536D-72,-.79815D0,-.0585D0) ADD00117 DELBF=Y(X,3.123D0,-.00178D0,0.00021D0) ADD00118 OM0=Y(X,8.6922D0,0.00785D0,-0.00054D0) ADD00119 CALL ADDSPEC (S,E,T1,T2,T3,T4,TEMP,NF,FREQ,ABSCOEF,0,LIKE,2,0,2,3)ADD00121 DO 70 I=1,NF ADD00122 70 ALFATOT(I)=ABSCOEF(I)+ALFATOT(I) ADD00123 C ADD00124 C PARAMETERS FOR 4045 (HYDROGEN HEXADECAPOLE) COMPONENTS ADD00125 C ADD00126 CALL BOUN54C (TEMP,RSI,NSOL) ADD00127 DO 80 I=1,NSOL ADD00128 RSILO(I)=DLOG(RSI(I)*1.D80) ADD00129 80 OMEG(I)=DFLOAT(I-1)*DX ADD00130 CALL SPLINE (NSOL,1,0,EPS,OMEG,RSILO,TT,SS,SI,NR,RSIGG) ADD00131 S=Y(X,(2.7321D-60/9.)*0.038,-2.32012D0,.23082D0) ADD00132 E=Y(X,-1.8198D0,-.00665D0,-.05626D0) ADD00133 T1=Y(X,1.3492D-13,.14472D0,-.06506D0) ADD00134 T2=Y(X,1.6596D-12,-.77679D0,.01401D0) ADD00135 T3=Y(X,5.9914D-13,-.16208D0,-.05895D0) ADD00136 T4=Y(X,1.9405D-14,.95965D0,-.10246D0) ADD00137 SCAL=0.038D0 ADD00138 CALL ADDSPEC (S,E,T1,T2,T3,T4,TEMP,NF,FREQ,ABSCOEF,0,LIKE,4,0,4,5)ADD00140 DO 90 I=1,NF ADD00141 90 ALFATOT(I)=ALFATOT(I)+ABSCOEF(I) ADD00142 C ADD00143 C METHANE INDUCED COMPONENTS ADD00144 C =============================== ADD00145 C OCTOPOLE- INDUCED TERM (43) ADD00146 C ADD00147 EPS=1.D-5 ADD00148 TT(1)=10.D0 ADD00149 CALL BOUND43 (TEMP,RSI,NSOL) ADD00150 DO 100 I=1,NSOL ADD00151 RSILO(I)=DLOG(RSI(I)*1.D80) ADD00152 100 OMEG(I)=FLOAT(I+INITB-2)*DX ADD00153 CALL SPLINE (NSOL,1,0,EPS,OMEG,RSILO,TT,SS,SI,NR,RSIGG) ADD00154 SCAL=1.D0 ADD00155 S=Y(X,.91196D-60/7.,-1.56529D0,.15284D0) ADD00156 E=Y(X,-.5D0,0.D0,0.D0) ADD00157 T1=Y(X,.51456D-13,.60523D0,-.10856D0) ADD00158 T2=Y(X,.62514D-12,-.51384D0,.00754D0) ADD00159 T3=Y(X,.55346D-12,-.40381D0,.00208D0) ADD00160 T4=Y(X,.13804D-13,1.9307D0,-.27921D0) ADD00161 CALL ADSPEC1 (S,E,T1,T2,T3,T4,TEMP,NF,FREQ,ABSCOEF,0,LIKE,0,3,3,4)ADD00163 DO 110 I=1,NF ADD00164 110 ALFATOT(I)=ABSCOEF(I)+ALFATOT(I) ADD00165 C ADD00166 C HEXADECAPOLE-INDUCED TERM ADD00167 C ============================= ADD00168 C ADD00169 EPS=1.D-5 ADD00170 TT(1)=10.D0 ADD00171 CALL BOUN54C (TEMP,RSI,NSOL) ADD00172 DO 120 I=1,NSOL ADD00173 RSILO(I)=DLOG(RSI(I)*1.D80) ADD00174 120 OMEG(I)=DFLOAT(I-1)*DX ADD00175 CALL SPLINE (NSOL,1,0,EPS,OMEG,RSILO,TT,SS,SI,NR,RSIGG) ADD00176 SCAL=1.D0 ADD00177 S=Y(X,2.7321D-60/9.,-2.32012D0,.23082D0) ADD00178 E=Y(X,-1.8198D0,-.00665D0,-.05626D0) ADD00179 T1=Y(X,1.3492D-13,.14472D0,-.06506D0) ADD00180 T2=Y(X,1.6596D-12,-.77679D0,.01401D0) ADD00181 T3=Y(X,5.9914D-13,-.16208D0,-.05895D0) ADD00182 T4=Y(X,1.9405D-14,.95965D0,-.10246D0) ADD00183 CALL ADSPEC1 (S,E,T1,T2,T3,T4,TEMP,NF,FREQ,ABSCOEF,0,LIKE,0,4,4,5)ADD00185 DO 130 I=1,NF ADD00186 130 ALFATOT(I)=ABSCOEF(I)+ALFATOT(I) ADD00187 WRITE(6, 170) FREQ(1),FREQ(NF),FREQ(2)-FREQ(1),TEMP ADD00190 WRITE(6,180) (ALFATOT(I),I=1,NF) ADD00191 C ADD00194 140 FORMAT (' ABSORPTION SPECTRA OF ', A10,' AT', F8.1, 1 'K',/1X, 1 ' MIN.FREQ.=', F8.1, ' CM-1', 10X, ' MAX.FREQ.=', F8.1, 2 ' CM-1',/, ' FREQ.INCREMENT=', F8.2, ' CM-1', 5X, 'IN', I5, 1 ' STEPS',/) ADD00198 150 FORMAT (F8.5) ADD00199 160 FORMAT (10F8.5) ADD00200 170 FORMAT (/, ' ABSORPTION COEFFICIENT ALPHA(FNU), FROM ',F5.1, ADD00201 1' CM-1 TO', F7.1, 9H CM-1, AT, F6.2, /, 1 23H CM-1 INCREMENTS, AT T=,F7.2 ADD00202 2, 29H K, IN UNITS OF CM-1 AMAGAT-2,/) ADD00203 180 FORMAT ( 1H ,10E13.5) ADD00204 190 FORMAT (4F10.5,I5) ADD00205 C ADD00206 END ADD00207 SUBROUTINE ADDSPEC (G0,EP,TAU1,TAU2,TAU5,TAU6,TEMP,NF,FREQ,ABSCOEFADD00208 1,MP,LIKE,LAMBDA1,LAMBDA2,LAMBDA,LVALUE) ADD00209 C ADD00210 C THIS ENTRY FOR L I N E A R M O L E C U L E!!!!! ADD00211 C SET LAMBDA1 NONZERO ADD00212 C THIS PROGRAM GENERATES LISTING OF CIA TR ALFA(OMEGA) ADD00213 C IF EITHER LAMBDA1 OR LAMBDA2 EQUAL TO ZERO - SINGLE TRANSITIONS; ADD00214 C LAMBDA1 CORRESPONDS TO H2, LAMBDA2 TO CH4 ADD00215 C LIKE=1 FOR LIKE SYSTEMS (AS H2-H2), SET LIKE=0 ELSE. ADD00216 C ADD00217 IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /BB/ OMEG(201),RSI(201),RSIGG(201),NSOL,BETA,SCAL ADD00219 COMMON /RSILO/ RSILO(201) ADD00220 COMMON /APP3/ SLIT,DX,NSRI,WNRMAX3,NS,NSRIUP,JSLIT ADD00221 COMMON /H2PART/ Q,WH2(2),B0,D0,JRANGE1 ADD00222 COMMON /BF/ G0BF,DELBF,OM0 ADD00223 COMMON /CHPART/ Q1,WCH(2),B01,D01,JRANGE2 ADD00224 DIMENSION ABSCOEF(NF), FREQ(NF) ADD00225 DATA CLOSCHM,BOLTZWN/2.68675484D19,.6950304/ ADD00226 DATA HBAR,PI,CLIGHT/1.054588757D-27,3.1415926535898,2.9979250D10/ ADD00227 EH2(I)=(B0-DFLOAT(I)*D0)*DFLOAT(I) ADD00228 PH2(J,T)=DFLOAT(2*J+1)*WH2(1+MOD(J,2))*DEXP(-1.4387859/T* 1 EH2(J*(J+1))) TWOPIC=2.*PI*CLIGHT ADD00231 CALIB=TWOPIC*((4.*PI**2)/(3.*HBAR*CLIGHT))*CLOSCHM**2 ADD00232 CALIB=CALIB/DFLOAT(1+LIKE) ADD00233 BETA=1./(BOLTZWN*TEMP) ADD00234 LIST=NF ADD00235 DO 10 I=1,LIST ADD00236 10 ABSCOEF(I)=0.0 ADD00237 C ADD00238 C ROTATIONAL SPECTRUM FOR THE DETAILED LISTING *******************ADD00239 C ADD00240 WRITE (6,50) LAMBDA1,LAMBDA2,LAMBDA,LVALUE,G0,EP,TAU1,TAU2, 1 TAU5,TAU6,G0*BGAMA(0.,TAU1,TAU2,EP,TAU5,TAU6,TEMP) IF (LAMBDA1.EQ.0) STOP 6666 ADD00246 C ADD00247 C FOR THIS ENTRY LAMBDA1 HAS TO BE NONZERO ADD00248 C SINGLE TRANSITIONS ON HYDROGEN'S ROTATIONAL FREQUENCIES. ADD00249 C LAMBDA IS EQUAL TO LAMBDA1 FOR SINGLE TRANSITIONS ADD00250 C ADD00251 JPLUSL=JRANGE1+LAMBDA ADD00252 DO 40 I=1,JRANGE1 ADD00253 J=I-1 ADD00254 DO 40 IP=1,JPLUSL ADD00255 JP=IP-1 ADD00256 CGS=CLEBSQR(J,LAMBDA,JP) ADD00257 IF (CGS) 40,40,20 ADD00258 20 P=PH2(J,TEMP)/Q ADD00259 OMEGA1=EH2(JP*IP)-EH2(J*I) ADD00260 FAC=CALIB*P*CGS ADD00261 DO 30 IQ=1,LIST ADD00262 FRQ=FREQ(IQ)-OMEGA1 ADD00263 WKI=FREQ(IQ)*(1.-DEXP(-BETA*FREQ(IQ))) ADD00264 WKF=WKI*FAC ADD00265 XBG=G0*BGAMA(FRQ,TAU1,TAU2,EP,TAU5,TAU6,TEMP) ADD00266 IF ( DABS(FRQ).LE.WNRMAX3 ) 1 XBG=XBG+SCAL*SPECFCT(FRQ,OMEG,RSILO,RSIGG,NSOL,BETA) IF (LVALUE.EQ.3 .AND. G0BF.NE.0.D0) 1 XBG = XBG + BGAUS(FRQ, G0BF, DELBF, OM0, TEMP) C ADD00272 C THIS MODELES THE PART OF THE BOUND-FREE SPECTRAL FUNCTION BY ADD00273 C MEANS OF THE DESYMMETRIZED GAUSSIAN PROFILE ADD00274 C ADD00275 ABSCOEF(IQ)=ABSCOEF(IQ)+XBG*WKF ADD00276 30 CONTINUE ADD00277 40 CONTINUE ADD00278 WRITE(6,60) (ABSCOEF(I),I=1,LIST) ADD00279 RETURN ADD00280 C ADD00281 50 FORMAT (/, 32H LAMBDA1,LAMBDA2, LAMBDA,LVALUE=,2I3,2X,2I3, 12H COMADD00282 1PONENT .,/15X, 22HLINE SHAPE PARAMETERS:,6E12.3,5X, 5HG(0)=,E12.3ADD00283 2/) ADD00284 60 FORMAT ((1X,10E12.4,/)) ADD00285 C ADD00286 END ADD00287 SUBROUTINE ADSPEC1 (G0,EP,TAU1,TAU2,TAU5,TAU6,TEMP,NF,FREQ,ABSCOEFADS00288 1,MP,LIKE,LAMBDA1,LAMBDA2,LAMBDA,LVALUE) ADS00289 C ADS00290 C FOR INDUCTION BY A T E T R A H E D R A L M O L E C U L E!! ADS00291 C THIS PROGRAM GENERATES LISTING OF CIA TR ALFA(OMEGA) ADS00292 C LAMBDA1 CORRESPONDS TO H2, LAMBDA2 TO CH4 ADS00293 C FOR THIS ENTRY SET LAMBDA2 NONZERO ADS00294 C LIKE=1 FOR LIKE SYSTEMS (AS H2-H2), SET LIKE=0 ELSE. ADS00295 C ADS00296 IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /BB/ OMEG(201),RSI(201),RSIGG(201),NSOL,BETA,SCAL ADS00298 COMMON /RSILO/ RSILO(201) ADS00299 COMMON /APP3/ SLIT,DX,NSRI,WNRMAX3,NS,NSRIUP,JSLIT ADS00300 COMMON /H2PART/ Q,WH2(2),B0,D0,JRANGE1 ADS00301 COMMON /CHPART/ Q1,WCH(2),B01,D01,JRANGE2 ADS00302 COMMON /STAT/ QQW3(51,3),QQW4(51,4) ADS00303 DIMENSION ABSCOEF(NF), FREQ(NF) ADS00304 DATA CLOSCHM,BOLTZWN/2.68675484D19,.6950304/ ADS00305 DATA HBAR,PI,CLIGHT/1.054588757D-27,3.1415926535898,2.9979250D10/ ADS00306 ECH4(I)=B01*DFLOAT(I) ADS00307 PCH4(J,T)=DEXP(-1.4387859/T*ECH4(J*(J+1))) ADS00308 TWOPIC=2.*PI*CLIGHT ADS00309 CALIB=TWOPIC*((4.*PI**2)/(3.*HBAR*CLIGHT))*CLOSCHM**2 ADS00310 CALIB=CALIB/DFLOAT(1+LIKE) ADS00311 BETA=1./(BOLTZWN*TEMP) ADS00312 LIST=NF ADS00313 DO 10 I=1,LIST ADS00314 10 ABSCOEF(I)=0.0 ADS00315 C ADS00316 C ROTATIONAL SPECTRUM: DETAILED LISTING ******************* ADS00317 C ADS00318 WRITE (6,70) LAMBDA1,LAMBDA2,LAMBDA,LVALUE,G0,EP,TAU1,TAU2, 1 TAU5,TAU6, G0*BGAMA(0.,TAU1,TAU2,EP,TAU5,TAU6,TEMP) C ADS00323 DO 60 I=1,JRANGE2 ADS00324 J=I-1 ADS00325 P=DFLOAT(2*J+1)*PCH4(J,TEMP)/Q1 ADS00326 DO 40 IP=1,LAMBDA ADS00327 C ADS00328 C POSITIVE DELTA J ADS00329 C ADS00330 JP=J+IP ADS00331 OMEGA1=ECH4(JP*(JP+1))-ECH4(J*I) ADS00332 IF (LAMBDA.EQ.3) CC=(1.+QQW3(I,IP)/4.) ADS00333 IF (LAMBDA.EQ.4) CC=(1.+QQW4(I,IP)/4.) ADS00334 FAC=CALIB*P*CC ADS00335 DO 20 IQ=1,LIST ADS00336 FRQ=FREQ(IQ)-OMEGA1 ADS00337 WKF=FREQ(IQ)*(1.-DEXP(-BETA*FREQ(IQ)))*FAC* 1 (2.*DFLOAT(JP)+1.) XBG=G0*BGAMA(FRQ,TAU1,TAU2,EP,TAU5,TAU6,TEMP) ADS00340 IF ( DABS(FRQ).LE.WNRMAX3 ) 1 XBG= XBG + SPECFCT(FRQ,OMEG,RSILO,RSIGG,NSOL,BETA) ABSCOEF(IQ)=ABSCOEF(IQ)+XBG*WKF ADS00343 20 CONTINUE ADS00344 C ADS00345 C NEGATIVE DELTA J ADS00346 C ADS00347 JP=J-IP ADS00348 IF (JP.LT.0) GO TO 40 ADS00349 OMEGA1=ECH4(JP*(JP+1))-ECH4(J*I) ADS00350 FAC=CALIB*P ADS00351 DO 30 IQ=1,LIST ADS00352 FRQ=FREQ(IQ)-OMEGA1 ADS00353 WKF=FREQ(IQ)*(1.-DEXP(-BETA*FREQ(IQ)))*FAC* 1 (2.*DFLOAT(JP)+1.) XBG=G0*BGAMA(FRQ,TAU1,TAU2,EP,TAU5,TAU6,TEMP) ADS00356 IF ( DABS(FRQ).LE.WNRMAX3 ) 1 XBG=XBG+SPECFCT(FRQ,OMEG,RSILO,RSIGG,NSOL,BETA) ABSCOEF(IQ)=ABSCOEF(IQ)+XBG*WKF ADS00359 30 CONTINUE ADS00360 40 CONTINUE ADS00361 C ADS00362 C DELTA J=0 ADS00363 C ADS00364 JP=J ADS00365 FAC=CALIB*P ADS00366 DO 50 IQ=1,LIST ADS00367 FRQ=FREQ(IQ) ADS00368 WKF=FREQ(IQ)*(1.-DEXP(-BETA*FREQ(IQ)))*FAC* 1 (2.*DFLOAT(JP)+1.) XBG=G0*BGAMA(FRQ,TAU1,TAU2,EP,TAU5,TAU6,TEMP) ADS00370 IF (DABS(FRQ).LE.WNRMAX3) 1 XBG=XBG+SPECFCT(FRQ,OMEG,RSILO,RSIGG,NSOL,BETA) ABSCOEF(IQ)=ABSCOEF(IQ)+XBG*WKF ADS00373 50 CONTINUE ADS00374 60 CONTINUE ADS00375 WRITE(6,80) (ABSCOEF(I),I=1,LIST) ADS00376 RETURN ADS00377 C ADS00378 70 FORMAT (/, 32H LAMBDA1,LAMBDA2, LAMBDA,LVALUE=,2I3,2X,2I3, 12H COMADS00379 1PONENT .,/15X, 22HLINE SHAPE PARAMETERS:,6E12.3,5X, 5HG(0)=,E12.3ADS00380 2/) ADS00381 80 FORMAT ((1X,10E12.4,/)) ADS00382 C ADS00383 END ADS00384 FUNCTION CLEBSQR (L,LAMBDA,LP) CLE00385 C CLE00386 C SQUARE OF CLEBSCH-GORDAN COEFFICIENT (L,LAMBDA,0,0;LP,0) CLE00387 C FOR INTEGER ARGUMENTS ONLY CLE00388 C NOTE THAT LAMBDA SHOULD BE SMALL, MAYBE @10 OR SO. CLE00389 C CLE00390 IMPLICIT DOUBLE PRECISION (A-H,O-Z) FC=DFLOAT(2*LP+1) CLE00391 GO TO 10 CLE00392 C CLE00393 ENTRY THREEJ2 CLE00394 C CLE00395 C THIS ENTRY RETURNS THE SQUARED 3-J SYMBOL L LAMBDA LP CLE00396 C 0 0 0 CLE00397 C INSTEAD OF THE CLEBSCH-GORDAN COEFFICIENT CLE00398 C (LIMITATION TO INTEGER ARGUMENTS ONLY) CLE00399 C CLE00400 C NOTE THAT THE THREE-J SYMBOLS ARE COMPLETELY SYMMETRIC IN THE CLE00401 C ARGUMENTS. IT WOULD BE ADVANTAGEOUS TO REORDER THE INPUT ARGUMENT CLE00402 C LIST SO THAT LAMBDA BECOMES THE SMALLEST OF THE 3 ARGUMENTS. CLE00403 C CLE00404 FC=1. CLE00405 10 CLEBSQR=0. CLE00406 IF (((L+LAMBDA).LT.LP).OR.((LAMBDA+LP).LT.L).OR.((L+LP).LT.LAMBDA)CLE00407 1) RETURN CLE00408 IF (MOD(L+LP+LAMBDA,2).NE.0) RETURN CLE00409 IF ((L.LT.0).OR.(LP.LT.0).OR.(LAMBDA.LT.0)) RETURN CLE00410 F=1./DFLOAT(L+LP+1-LAMBDA) CLE00411 IF (LAMBDA.EQ.0) GO TO 30 CLE00412 I1=(L+LP+LAMBDA)/2 CLE00413 I0=(L+LP-LAMBDA)/2+1 CLE00414 DO 20 I=I0,I1 CLE00415 20 F=F*DFLOAT(I)/FLOAT(2*(2*I+1)) CLE00416 30 P=FC*F*FCTL(LAMBDA+L-LP)*FCTL(LAMBDA+LP-L) CLE00417 CLEBSQR=P/(FCTL((LAMBDA+L-LP)/2)*FCTL((LAMBDA+LP-L)/2))**2 CLE00418 RETURN CLE00419 C CLE00420 END CLE00421 FUNCTION FCTL (N) FCT00422 IMPLICIT DOUBLE PRECISION (A-H,O-Z) P(Z)=((((-2.294720936D-4)/Z-(2.681327160D-3))/Z+(3.472222222D-3))/FCT00423 1Z+(8.333333333D-2))/Z+1. FCT00424 FCTL=1. FCT00425 IF (N.LE.1) RETURN FCT00426 IF (N.GT.15) GO TO 20 FCT00427 J=1 FCT00428 DO 10 I=2,N FCT00429 10 J=J*I FCT00430 FCTL=DFLOAT(J) FCT00431 RETURN FCT00432 20 Z=DFLOAT(N+1) FCT00433 FCTL=(DEXP(-Z)*(Z**(Z-0.5))*P(Z)*2.5066282746310) FCT00434 RETURN FCT00435 C FCT00436 END FCT00437 FUNCTION BGAMA (FNU,T1,T2,EPS,T3,T4,TEMP) BGA00438 C BGA00439 C EQUATION 13, SO-CALLED EBC MODEL, OF BORYSOW,TRAFTON,FROMMHOLD, BGA00440 C AND BIRNBAUM, ASTROPHYS.J., (1985) BGA00441 C BGA00442 IMPLICIT DOUBLE PRECISION (A-H,O-Z) REAL*8 K0 BGA00443 DATA PI,CLIGHT/3.1415926535898,2.99792458D10/ BGA00444 DATA HBAR,BOLTZ/1.0545887D-27,1.380662D-16/ BGA00445 P1(X)=((((((.0045813*X+.0360768)*X+.2659732)*X+1.2067492)*X+3.0899BGA00446 1424)*X+3.5156229)*X+1.) BGA00447 P2(X)=((((((.00000740*X+.00010750)*X+.00262698)*X+.03488590)*X+.23BGA00448 1069756)*X+.42278420)*X-.57721566) BGA00449 P3(X)=((((((.00032411*X+.00301532)*X+.02658733)*X+.15084934)*X+.51BGA00450 1498869)*X+.87890594)*X+.5) BGA00451 P4(X)=((((((-.00004686*X-.00110404)*X-.01919402)*X-.18156897)*X-.6BGA00452 17278579)*X+.15443144)*X+1.) BGA00453 P5(X)=((((((.00053208*X-.00251540)*X+.00587872)*X-.01062446)*X+.02BGA00454 1189568)*X-.07832358)*X+1.25331414) BGA00455 P6(X)=((((((-.00068245*X+.00325614)*X-.00780353)*X+.01504268)*X-.0BGA00456 13655620)*X+.23498619)*X+1.25331414) BGA00457 C BGA00458 OMEGA=2.*PI*CLIGHT*FNU BGA00459 T0=HBAR/(2.*BOLTZ*TEMP) BGA00460 Z=DSQRT((1.+(OMEGA*T1)**2)*(T2*T2+T0*T0))/T1 BGA00461 IF (Z-2.) 10,10,20 BGA00462 10 XK1=Z*Z*DLOG(Z/2.)*P3((Z/3.75)**2)+P4((Z/2.)**2) BGA00463 GO TO 30 BGA00464 20 XK1=DSQRT(Z)*DEXP(-Z)*P6(2./Z) BGA00465 30 ZP=DSQRT((1.+(OMEGA*T4)**2)*(T3*T3+T0*T0))/T4 BGA00466 IF (ZP-2.) 40,40,50 BGA00467 40 K0=-DLOG(ZP/2.)*P1((ZP/3.75)**2)+P2((ZP/2.)**2) BGA00468 GO TO 60 BGA00469 50 K0=DEXP(-ZP)*P5(2./ZP)/DSQRT(ZP) BGA00470 60 BGAMA=((T1/PI)*DEXP(T2/T1+T0*OMEGA)*XK1/ 1 (1.+(T1*OMEGA)**2)+EPS*(T3/PI)*DEXP(T3/T4+T0*OMEGA)*K0)/ 1 (1.+EPS) RETURN BGA00473 END BGA00475 SUBROUTINE PARTSUM (TEMP) PAR00476 C PAR00477 C H2 ROTATIONAL PARTITION SUM Q = Q(T). PAR00478 C PAR00479 IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /H2PART/ Q,WH2(2),B0,D0,JRANGE1 PAR00480 COMMON /CHPART/ Q1,WCH(2),B01,D01,JRANGE2 PAR00481 DATA B0,D0,WH2(1),WH2(2)/59.3392,0.04599,1.,3./ PAR00482 DATA B01,D01,WCH(1),WCH(2)/5.24,0.,1.,1./ PAR00483 EH2(I)=(B0-DFLOAT(I)*D0)*DFLOAT(I) PAR00484 PH2(J,T)=DFLOAT(2*J+1)*WH2(1+MOD(J,2))*DEXP(-1.4387859* 1 EH2(J*(J+1))/T) ECH4(I)=B01*DFLOAT(I) PAR00487 PCH4(J,T)=DEXP(-1.4387859/T*ECH4(J*(J+1))) PAR00488 C PAR00489 C Q,B01,D01,WCH - PARTITION FCT., ROT.CONSTANTS, WEIGHTS FOR CH4 PAR00490 C Q,B0,D0,WH2 - PARTITION FCT., ROT.CONSTANTS, WEIGHTS FOR H2 PAR00491 C PAR00492 Q=0. PAR00493 J=0 PAR00494 10 DQ=PH2(J,TEMP) PAR00495 Q=Q+DQ PAR00496 J=J+1 PAR00497 IF (DQ.GT.Q/900.) GO TO 10 PAR00498 JRANGE1=J PAR00499 C PAR00501 C *** PARTITION FUNCTION FOR CH4 ******************* PAR00502 C PAR00503 Q1=0. PAR00504 J=0 PAR00505 20 DQ=PCH4(J,TEMP)*DFLOAT((2*J+1)**2) PAR00506 Q1=Q1+DQ PAR00507 J=J+1 PAR00508 IF (DQ.GT.Q1/1000.) GO TO 20 PAR00509 JRANGE2=J+3 PAR00510 RETURN PAR00511 C PAR00512 30 FORMAT (/, 40H ROTATIONAL PARTITION FUNCTION OF H2: Q=,F8.2,10X, PAR00513 17HJ MAX =,I3/) PAR00514 C PAR00515 END PAR00516 SUBROUTINE PROFILE (X,Y) PRO00517 IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /BL3/ RSI(401) PRO00518 C PRO00519 C ATRIANGULAR SLIT FUNCTION IS USED. PRO00520 C PRO00521 COMMON /APP3/ SLIT,DX,NSRI,WNRMAX3,NS,NSRIUP,JSLIT PRO00522 IF (Y) 50,60,10 PRO00523 10 CONTINUE PRO00524 X0=(NSRI+1.)+X/DX PRO00525 NC=X0 PRO00526 N1=NC+1 PRO00527 SLOPE=Y/SLIT PRO00528 NU=X0-NS PRO00529 IF (NU.LT.1) NU=1 PRO00530 IF (NU.GT.NSRIUP) RETURN PRO00531 NO=X0+NS PRO00532 IF (NO.GT.NSRIUP) NO=NSRIUP PRO00533 IF (NO.LT.1) RETURN PRO00534 IF (NC.GT.NSRIUP) NC=NSRIUP PRO00535 IF (NC.LE.1) GO TO 30 PRO00536 DO 20 I=NU,NC PRO00537 XI=(I-1.)*DX-WNRMAX3 PRO00538 DR=SLOPE*(XI-(X-SLIT)) PRO00539 IF (DR.LE.0.) GO TO 20 PRO00540 RSI(I)=RSI(I)+DR PRO00541 20 CONTINUE PRO00542 30 IF (NC.GE.NSRIUP) RETURN PRO00543 IF (N1.LT.1) N1=1 PRO00544 DO 40 I=N1,NO PRO00545 XI=(I-1.)*DX-WNRMAX3 PRO00546 DR=Y-SLOPE*(XI-X) PRO00547 IF (DR.LE.0.) GO TO 40 PRO00548 RSI(I)=RSI(I)+DR PRO00549 40 CONTINUE PRO00550 RETURN PRO00551 50 WRITE(6,70) SLIT PRO00552 60 CONTINUE RETURN PRO00554 C PRO00555 70 FORMAT (/, 30H A TRIANGULAR SLIT FUNCTION OF,F6.3, 23H CM-1 HALFWIPRO00556 1DTH IS USED,/) PRO00557 C PRO00558 END PRO00559 FUNCTION SPECFCT (FREQ,OMEGA,PHI,PHI2,N,RTEMP) SPE00471 C SPE00472 C THIS INTERPOLATES THE SPECTRAL FUNCTION PHI(FREQ) DEFINED AT SPE00473 C OMEGA(N) AS PHI(N). PHI2 IS THE SECOND DERIVATIVE AT OMEGA SPE00474 C WHICH MUST BE OBTAINED FIRST (USE SPLINE FOR THAT PURPOSE). SPE00475 C RTEMP IS THE RECIPROCAL TEMPERATURE IN CM-1 UNITS. SPE00476 C NOTE THAT WE INTERPOLATE 1.E80 TIMES THE LOGARITHM OF PHI(OMEGA) SPE00477 C SPE00478 IMPLICIT DOUBLE PRECISION (A-h,o-z) DIMENSION PHI(N), PHI2(N), OMEGA(N) SPE00479 dimension f(1), gp(1) c correction above 1-MAR-1989 18:10:02 in all subroutine TFAC=0. SPE00480 F(1)=FREQ IF (F(1) ) 10,20,20 10 F(1)=dABS(F(1)) TFAC=(-RTEMP*F(1)) 20 IF (F(1).LE.OMEGA(N)) GO TO 30 SPECFCT=DEXP(-(PHI(N-1)-PHI(N))*(F(1)-OMEGA(N))/ 1 (OMEGA(N)-OMEGA(N-1))+ SPE00486 1PHI(N)+TFAC)*(1.D-80) SPE00487 RETURN SPE00488 c exchanged ixpolat to spline 1-MAR-1989 18:08:27 30 CALL spline (N,1,0,1.D-6,OMEGA,PHI,F,GP,SI,NR,PHI2) SPE00489 c f(1), gp(1) SPECFCT=DEXP(TFAC+GP(1))*(1.D-80) RETURN SPE00491 C SPE00492 END SPE00493 SUBROUTINE SPLINE (L,M,K,EPS,X,Y,T,SS,SI,NR,S2) SPL00583 C SPL00584 C SPLINE INTERPOLATION AND QUADRATURE, THIRD ORDER AFTER GREVILLE. SPL00585 C INPUT ARGUMENTS L...Y, OUTPUT SS...NR. SPL00586 C L DATA POINTS X(1), Y(1) ... X(L),Y(L) SPL00587 C EPS=ERROR CRITERION, TYPICALLY EPS=1.E-5 FOR 5 DECI. PLACES ACCURASPL00588 C M ARGUMENTS T(1)..T(M) FOR WHICH FUNCTION VALUES SS(1)..SS(M), FORSPL00589 C K=0; OR FIRST OR SECOND DERIVATIVE FOR K=1 OR -1, RESPECTIVELY. SPL00590 C NOTE THAT M HAS TO BE AT LEAST EQUAL TO 1. SPL00591 C SI=INTEGRAL (OVER WHOLE INTERVAL) FOR K=2 ONLY. SPL00592 C FOR 'NATURAL' SPLINE FUNCTIONS, S2(1)=S2(L)=0. MUST BE INPUT*NOTE*SPL00593 C N0 INDICATES THE NUMBER OF OUT-OF-RANGE CALLS. X(1)@T(I)@X(L) SPL00594 C EXTRAPOLATE WITH CAUTION. (ASSUMPTION D2Y/DX2 = 0.) SPL00595 C S2(I) IS THE 2ND DERIVATIVE AT X=X(I) AND IS COMPUTED INTERNALLY. SPL00596 C SPL00597 IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION X(L), Y(L), T(M), SS(M), S2(L) SPL00598 DELY(I)=(Y(I+1)-Y(I))/(X(I+1)-X(I)) SPL00599 B(I)=(X(I)-X(I-1))*0.5/(X(I+1)-X(I-1)) SPL00600 C(I)=3.*(DELY(I)-DELY(I-1))/(X(I+1)-X(I-1)) SPL00601 N=L SPL00602 N1=N-1 SPL00603 NR=0 SPL00604 DO 10 I=2,N1 SPL00605 10 S2(I)=C(I)/1.5 SPL00606 OMEGA=1.0717968 SPL00607 IC=0 SPL00608 C SPL00609 C 'NATURAL' SPLINE FUNCTIONS OF THIRD ORDER. SPL00610 C SPL00611 S2(N)=0. SPL00612 S2(1)=S2(N) SPL00613 20 ETA=0. SPL00614 IC=IC+1 SPL00615 SM=DABS(S2(1)) SPL00616 DO 30 I=2,N SPL00617 IF (DABS(S2(I)).GT.SM) SM=DABS(S2(I)) SPL00618 30 CONTINUE SPL00619 EPSI=EPS*SM SPL00620 DO 50 I=2,N1 SPL00621 W=(C(I)-B(I)*S2(I-1)-(0.5-B(I))*S2(I+1)-S2(I))*OMEGA SPL00622 IF (DABS(W)-ETA) 50,50,40 SPL00623 40 ETA=DABS(W) SPL00624 50 S2(I)=S2(I)+W SPL00625 IF (ETA-EPSI) 60,20,20 SPL00626 ENTRY IXPOLAT SPL00627 C SPL00628 C THIS ENTRY USEFUL WHEN ITERATION PREVIOUSLY COMPLETED SPL00629 C SPL00630 N=L SPL00631 N1=N-1 SPL00632 NR=0 SPL00633 IC=-1 SPL00634 60 IF (K.EQ.2) GO TO 260 SPL00635 GO TO 70 SPL00636 70 DO 250 J=1,M SPL00637 I=1 SPL00638 IF (T(J)-X(1)) 110,210,80 SPL00639 80 IF (T(J)-X(N)) 100,190,150 SPL00640 90 IF (T(J)-X(I)) 200,210,100 SPL00641 100 I=I+1 SPL00642 GO TO 90 SPL00643 110 NR=NR+1 SPL00644 HT1=T(J)-X(1) SPL00645 HT2=T(J)-X(2) SPL00646 YP1=DELY(1)+(X(1)-X(2))*(2.*S2(1)+S2(2))/6. SPL00647 IF (K) 140,130,120 SPL00648 120 SS(J)=YP1+HT1*S2(1) SPL00649 GO TO 250 SPL00650 130 SS(J)=Y(1)+YP1*HT1+S2(1)*HT1*HT1/2. SPL00651 GO TO 250 SPL00652 140 SS(J)=S2(I) SPL00653 GO TO 250 SPL00654 150 HT2=T(J)-X(N) SPL00655 HT1=T(J)-X(N1) SPL00656 NR=NR+1 SPL00657 YPN=DELY(N1)+(X(N)-X(N1))*(S2(N1)+2.*S2(N))/6. SPL00658 IF (K) 180,170,160 SPL00659 160 SS(J)=YPN+HT2*S2(N) SPL00660 GO TO 250 SPL00661 170 SS(J)=Y(N)+YPN*HT2+S2(N)*HT2*HT2/2. SPL00662 GO TO 250 SPL00663 180 SS(J)=S2(N) SPL00664 GO TO 250 SPL00665 190 I=N SPL00666 200 I=I-1 SPL00667 210 HT1=T(J)-X(I) SPL00668 HT2=T(J)-X(I+1) SPL00669 PROD=HT1*HT2 SPL00670 S3=(S2(I+1)-S2(I))/(X(I+1)-X(I)) SPL00671 SS2=S2(I)+HT1*S3 SPL00672 DELSQS=(S2(I)+S2(I+1)+SS2)/6. SPL00673 IF (K) 240,220,230 SPL00674 220 SS(J)=Y(I)+HT1*DELY(I)+PROD*DELSQS SPL00675 GO TO 250 SPL00676 230 SS(J)=DELY(I)+(HT1+HT2)*DELSQS+PROD*S3/6. SPL00677 GO TO 250 SPL00678 240 SS(J)=SS2 SPL00679 250 CONTINUE SPL00680 260 SI=0. SPL00681 DO 270 I=1,N1 SPL00682 H=X(I+1)-X(I) SPL00683 270 SI=SI+0.5*H*(Y(I)+Y(I+1))-H**3*(S2(I)+S2(I+1))/24. SPL00684 IF (K.EQ.2) NR=IC SPL00685 RETURN SPL00686 END SPL00688 SUBROUTINE BOUND32 (TEMP,RSI,NSOL) BOU00689 IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION EB(2,9), A(9,9), RSI(201) BOU00690 C BOU00691 C EB(I,K) - BOUND ENERGIES FOR VIB.NR.V=I-1, ROT.NR J=K-1 BOU00692 C A(I,J)=( C(I,L,J)**2) * (MTX.EL. (I,BETA,J)**2 ) BOU00693 C I,J-DENOTE ROTATIONAL STATES FOR VIB. LEVEL V=0 BOU00694 C BOU00695 COMMON /APP3/ SLIT,DX,NSRI,WNRMAX3,NS,NSRIUP,JSLIT BOU00696 COMMON /BL3/ RSIBB(401) BOU00697 DATA (EB(1,J),J=1,9)/-26.1308,-24.9681,-22.6532,-19.2080, 1 -14.669,-9.0915,-2.564,4.7286,12.4947/ DATA (EB(2,J),J=1,9)/-0.516,-0.1246,0.,0.,0.,0.,0.,0.,0./ BOU00700 DATA TWOPIC/1.88365183D11/ BOU00701 NSRI=190 BOU00702 WNRMAX3=25.D0 BOU00703 C BOU00704 C NSRI - HOW MANY POINTS FOR B-B SPECTRAL FUNCTION TO BE GIVEN BOU00705 C WNRMAX3 - THE FREQUENCY RANGE OF B-B CONTRIBUTION BOU00706 C SLIT- THE HALFWIDTH OF THE SPECTRAL PROFILE CONVOLUTED WITH BOU00707 C B-B SPECTRUM, IN [CM-1]. BOU00708 C BOU00709 NSRIUP=2*NSRI+1 BOU00710 DX=WNRMAX3/DFLOAT(NSRI) BOU00711 DO 10 I=1,401 BOU00712 10 RSIBB(I)=0.0 BOU00713 NS=INT(SLIT/DX) BOU00714 DO 20 I=1,9 BOU00715 DO 20 J=1,9 BOU00716 20 A(I,J)=0.0 BOU00717 C BOU00718 C A(I,J) FOR L=32 CONTRIBUTION H2-CH4 SYSTEM BOU00719 C BOU00720 A(1,4)=.14127D-39 BOU00721 A(2,3)=.18414D-39 BOU00722 A(2,5)=.23414D-39 BOU00723 A(3,4)=.18534D-39 BOU00724 A(3,6)=.30894D-39 BOU00725 A(4,5)=.21745D-39 BOU00726 A(4,7)=.36564D-39 BOU00727 A(5,6)=.24659D-39 BOU00728 A(5,8)=.40032D-39 BOU00729 A(6,7)=.26656D-39 BOU00730 A(7,8)=.27376D-39 BOU00731 A(8,9)=.26805D-39 BOU00732 A(6,9)=.41233D-39 BOU00733 ALFA=1./(0.69519*TEMP) BOU00734 RM=1.7768*1.672649D-24 BOU00735 PI=3.141592654 BOU00736 PF=((RM*(1.380662D-16)*TEMP*2.*PI/(6.626176D-27)**2)**1.5) BOU00737 C BOU00738 C DELTA V = 0 BELOW BOU00739 C BOU00740 DO 40 L=1,8 BOU00741 STOKE1=EB(1,L+1)-EB(1,L) BOU00742 C BOU00743 C FREQUENCY SHIFT FOR DELTA J=+1,DELTA V=0 BOU00744 C BOU00745 IF (L.GT.6) GO TO 30 BOU00746 STOKE3=EB(1,L+3)-EB(1,L) BOU00747 C BOU00748 C STOKE3- FREQUENCY SHIFT FOR DELTA J=+3, DELTA V=0 BOU00749 C BOU00750 STOKI=A(L,L+3)*DEXP(-ALFA*EB(1,L))/PF BOU00751 CALL PROFILE (STOKE3,STOKI) BOU00752 STOKIP=A(L,L+3)*DEXP(-ALFA*EB(1,L+3))/PF BOU00753 CALL PROFILE (-STOKE3,STOKIP) BOU00754 30 STOKI=A(L,L+1)*DEXP(-ALFA*EB(1,L))/PF BOU00755 CALL PROFILE (STOKE1,STOKI) BOU00756 STOKIP=A(L,L+1)*DEXP(-ALFA*EB(1,L+1))/PF BOU00757 CALL PROFILE (-STOKE1,STOKIP) BOU00758 40 CONTINUE BOU00759 C BOU00760 C DELTA V=1 BELOW BOU00761 C BOU00762 AV1=.39933D-41 BOU00763 STOKE=22.5287 BOU00764 STOKI=AV1*DEXP(-ALFA*EB(1,3))/PF BOU00765 CALL PROFILE (STOKE,STOKI) BOU00766 STOKIP=AV1*DEXP(-ALFA*EB(2,2))/PF BOU00767 CALL PROFILE (-STOKE,STOKIP) BOU00768 AV1=0.3079D-41 BOU00769 STOKE=18.6921 BOU00770 STOKI=AV1*DEXP(-ALFA*EB(1,4))/PF BOU00771 CALL PROFILE (STOKE,STOKI) BOU00772 STOKIP=AV1*DEXP(-ALFA*EB(2,1))/PF BOU00773 CALL PROFILE (-STOKE,STOKIP) BOU00774 AV1=0.4281D-41 BOU00775 STOKE=14.544 BOU00776 STOKI=AV1*DEXP(-ALFA*EB(1,5))/PF BOU00777 CALL PROFILE (STOKE,STOKI) BOU00778 STOKIP=AV1*DEXP(-ALFA*EB(2,2))/PF BOU00779 CALL PROFILE (-STOKE,STOKIP) BOU00780 DO 50 N=1,NSRIUP BOU00781 50 RSIBB(N)=RSIBB(N)/TWOPIC/SLIT BOU00782 NSOL=NSRI+1 BOU00783 DO 60 I=1,NSOL BOU00784 K=(NSRI+1)+(I-1) BOU00785 60 RSI(I)=RSIBB(K) BOU00786 C BOU00787 C RSI - BOUND-BOUND CONTRIBUTION FOR POSITIVE FREQUENCY SHIFTS BOU00788 C BOU00789 RETURN BOU00790 C BOU00791 END BOU00792 SUBROUTINE BOUND43 (TEMP,RSI,NSOL) BOU00793 IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION EB(2,9), A(9,9), RSI(201) BOU00794 C BOU00795 C EB(I,K) - BOUND ENERGIES FOR VIB.NR.V=I-1, ROT.NR J=K-1 BOU00796 C A(I,J)=( C(I,L,J)**2) * (MTX.EL. (I,BETA,J)**2 ) BOU00797 C I,J-DENOTE ROTATIONAL STATES FOR VIB. LEVEL V=0 BOU00798 C BOU00799 COMMON /APP3/ SLIT,DX,NSRI,WNRMAX3,NS,NSRIUP,JSLIT BOU00800 COMMON /BL3/ RSIBB(401) BOU00801 COMMON /BOU43/ INITB BOU00802 DATA (EB(1,J),J=1,9)/-26.1308,-24.9681,-22.6532,-19.2080,-14.669,-BOU00803 19.0915,-2.564,4.7286,12.4947/ BOU00804 DATA (EB(2,J),J=1,9)/-0.516,-0.1246,0.,0.,0.,0.,0.,0.,0./ BOU00805 DATA TWOPIC/1.88365183D11/ BOU00806 NSRI=190 BOU00807 WNRMAX3=30. BOU00808 INITB=1 BOU00809 C BOU00810 C NSRI - HOW MANY POINTS FOR B-B SPECTRAL FUNCTION TO BE GIVEN BOU00811 C WNRMAX3 - THE FREQUENCY RANGE OF B-B CONTRIBUTION BOU00812 C SLIT- THE HALFWIDTH OF THE SPECTRAL PROFILE CONVOLUTED WITH BOU00813 C B-B SPECTRUM, IN [CM-1]. BOU00814 C BOU00815 NSRIUP=2*NSRI+1 BOU00816 DX=WNRMAX3/DFLOAT(NSRI) BOU00817 DO 10 I=1,401 BOU00818 10 RSIBB(I)=0.0 BOU00819 NS=INT(SLIT/DX) BOU00820 DO 20 I=1,9 BOU00821 DO 20 J=1,9 BOU00822 20 A(I,J)=0.0 BOU00823 C BOU00824 C A(I,J) FOR L,Lambda={3,4} CONTRIBUTION H2-CH4 SYSTEM BOU00825 C MATRIX ELEMENTS AS IF M(ll')/7 BOU00826 C BOU00827 A(1,5)=.33038D-41 BOU00828 A(2,4)=.45246D-41 BOU00829 A(2,6)=.52213D-41 BOU00830 A(3,5)=.42125D-41 BOU00831 A(3,7)=.65899D-41 BOU00832 A(4,6)=.46727D-41 BOU00833 A(4,8)=.74203D-41 BOU00834 A(5,7)=.50413D-41 BOU00835 A(5,9)=.77286D-41 BOU00836 A(6,8)=.51657D-41 BOU00837 A(7,9)=.50297D-41 BOU00838 ALFA=1./(0.69519*TEMP) BOU00839 RM=1.7768*1.672649D-24 BOU00840 PI=3.141592654 BOU00841 PF=((RM*(1.380662D-16)*TEMP*2.*PI/(6.626176D-27)**2)**1.5) BOU00842 C BOU00843 C DELTA V = 0 BELOW BOU00844 C BOU00845 DO 40 L=1,7 BOU00846 STOKE1=EB(1,L+2)-EB(1,L) BOU00847 C BOU00848 C FREQUENCY SHIFT FOR DELTA J=+2,DELTA V=0 BOU00849 C BOU00850 IF (L.GT.5) GO TO 30 BOU00851 STOKE3=EB(1,L+4)-EB(1,L) BOU00852 C BOU00853 C STOKE3- FREQUENCY SHIFT FOR DELTA J=+4, DELTA V=0 BOU00854 C BOU00855 STOKI=A(L,L+4)*DEXP(-ALFA*EB(1,L))/PF BOU00856 CALL PROFILE (STOKE3,STOKI) BOU00857 STOKIP=A(L,L+4)*DEXP(-ALFA*EB(1,L+4))/PF BOU00858 CALL PROFILE (-STOKE3,STOKIP) BOU00859 30 STOKI=A(L,L+2)*DEXP(-ALFA*EB(1,L))/PF BOU00860 CALL PROFILE (STOKE1,STOKI) BOU00861 STOKIP=A(L,L+2)*DEXP(-ALFA*EB(1,L+2))/PF BOU00862 CALL PROFILE (-STOKE1,STOKIP) BOU00863 40 CONTINUE BOU00864 C BOU00865 C DELTA V=1 BELOW BOU00866 C BOU00867 AV1=.12153D-42 BOU00868 STOKE=19.0835 BOU00869 STOKI=AV1*DEXP(-ALFA*EB(1,4))/PF BOU00870 CALL PROFILE (STOKE,STOKI) BOU00871 STOKIP=AV1*DEXP(-ALFA*EB(2,2))/PF BOU00872 CALL PROFILE (-STOKE,STOKIP) BOU00873 AV1=0.89592D-43 BOU00874 STOKE=14.1528 BOU00875 STOKI=AV1*DEXP(-ALFA*EB(1,5))/PF BOU00876 CALL PROFILE (STOKE,STOKI) BOU00877 STOKIP=AV1*DEXP(-ALFA*EB(2,1))/PF BOU00878 CALL PROFILE (-STOKE,STOKIP) BOU00879 AV1=0.1182D-42 BOU00880 STOKE=8.9669 BOU00881 STOKI=AV1*DEXP(-ALFA*EB(1,6))/PF BOU00882 CALL PROFILE (STOKE,STOKI) BOU00883 STOKIP=AV1*DEXP(-ALFA*EB(2,2))/PF BOU00884 CALL PROFILE (-STOKE,STOKIP) BOU00885 DO 50 N=1,NSRIUP BOU00886 50 RSIBB(N)=RSIBB(N)/TWOPIC/SLIT BOU00887 NSOL=NSRI+1 BOU00888 DO 60 I=1,NSOL BOU00889 K=(NSRI+1)+(I-1) BOU00890 60 RSI(I)=RSIBB(K) BOU00891 C BOU00892 C RSI - BOUND-BOUND CONTRIBUTION FOR POSITIVE FREQUENCY SHIFTS BOU00893 C BOU00894 I=0 BOU00895 70 I=I+1 BOU00896 IF (RSI(I).EQ.0.) GO TO 70 BOU00897 IF (I.EQ.1) GO TO 90 BOU00898 NSOL=NSOL-I+1 BOU00899 DO 80 J=1,NSOL BOU00900 80 RSI(J)=RSI(J+I-1) BOU00901 INITB=I BOU00902 90 CONTINUE BOU00903 RETURN BOU00904 END BOU00906 SUBROUTINE BOUN54C (TEMP,RSI,NSOL) BOU00907 IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION EB(2,9), A(9,9), RSI(201) BOU00908 C BOU00909 C EB(I,K) - BOUND ENERGIES FOR VIB.NR.V=I-1, ROT.NR J=K-1 BOU00910 C A(I,J)=( C(I,L,J)**2) * (MTX.EL. (I,BETA,J)**2 ) BOU00911 C I,J-DENOTE ROTATIONAL STATES FOR VIB. LEVEL V=0 BOU00912 C BOU00913 COMMON /APP3/ SLIT,DX,NSRI,WNRMAX3,NS,NSRIUP,JSLIT BOU00914 COMMON /BL3/ RSIBB(401) BOU00915 DATA (EB(1,J),J=1,9)/-26.1308,-24.9681,-22.6532,-19.2080,-14.669,-BOU00916 19.0915,-2.564,4.7286,12.4947/ BOU00917 DATA (EB(2,J),J=1,9)/-0.516,-0.1246,0.,0.,0.,0.,0.,0.,0./ BOU00918 DATA TWOPIC/1.88365183D11/ BOU00919 NSRI=190 BOU00920 WNRMAX3=34. BOU00921 C BOU00922 C NSRI - HOW MANY POINTS FOR B-B SPECTRAL FUNCTION TO BE GIVEN BOU00923 C WNRMAX3 - THE FREQUENCY RANGE OF B-B CONTRIBUTION BOU00924 C SLIT- THE HALFWIDTH OF THE SPECTRAL PROFILE CONVOLUTED WITH BOU00925 C B-B SPECTRUM, IN [CM-1]. BOU00926 C BOU00927 NSRIUP=2*NSRI+1 BOU00928 DX=WNRMAX3/DFLOAT(NSRI) BOU00929 DO 10 I=1,401 BOU00930 10 RSIBB(I)=0.0 BOU00931 NS=INT(SLIT/DX) BOU00932 DO 20 I=1,9 BOU00933 DO 20 J=1,9 BOU00934 20 A(I,J)=0.0 BOU00935 C BOU00936 C A(I,J) FOR L=54(CH4) CONTRIBUTION H2-CH4 SYSTEM BOU00937 C MATRIX ELEMENTS AS IF M(LL')/9 BOU00938 C BOU00939 A(1,6)=.99184D-42 BOU00940 A(2,5)=.14168D-41 BOU00941 A(2,7)=.14969D-41 BOU00942 A(3,4)=.16093D-41 BOU00943 A(3,6)=.12455D-41 BOU00944 A(4,5)=.12587D-41 BOU00945 A(4,7)=.13120D-41 BOU00946 A(5,6)=.12714D-41 BOU00947 A(5,8)=.13379D-41 BOU00948 A(6,7)=.12818D-41 BOU00949 A(7,8)=.12378D-41 BOU00950 A(6,9)=.12963D-41 BOU00951 A(8,9)=.11352D-41 BOU00952 A(3,8)=.17927D-41 BOU00953 A(4,9)=.19144D-41 BOU00954 ALFA=1./(0.69519*TEMP) BOU00955 RM=1.7768*1.672649D-24 BOU00956 PI=3.141592654 BOU00957 PF=((RM*(1.380662D-16)*TEMP*2.*PI/(6.626176D-27)**2)**1.5) BOU00958 C BOU00959 C DELTA V = 0 BELOW BOU00960 C BOU00961 DO 50 L=1,8 BOU00962 STOKE1=EB(1,L+1)-EB(1,L) BOU00963 C BOU00964 C FREQUENCY SHIFT FOR DELTA J=+1,DELTA V=0 BOU00965 C BOU00966 IF (L.GT.6) GO TO 40 BOU00967 STOKE3=EB(1,L+3)-EB(1,L) BOU00968 C BOU00969 C STOKE3- FREQUENCY SHIFT FOR DELTA J=+3, DELTA V=0 BOU00970 C BOU00971 IF (L.GT.4) GO TO 30 BOU00972 STOKE5=EB(1,L+5)-EB(1,L) BOU00973 STOKI=A(L,L+5)*DEXP(-ALFA*EB(1,L))/PF BOU00974 CALL PROFILE (STOKE5,STOKI) BOU00975 STOKIP=A(L,L+5)*DEXP(-ALFA*EB(1,L+5))/PF BOU00976 CALL PROFILE (-STOKE5,STOKIP) BOU00977 30 STOKI=A(L,L+3)*DEXP(-ALFA*EB(1,L))/PF BOU00978 CALL PROFILE (STOKE3,STOKI) BOU00979 STOKIP=A(L,L+3)*DEXP(-ALFA*EB(1,L+3))/PF BOU00980 CALL PROFILE (-STOKE3,STOKIP) BOU00981 40 STOKI=A(L,L+1)*DEXP(-ALFA*EB(1,L))/PF BOU00982 CALL PROFILE (STOKE1,STOKI) BOU00983 STOKIP=A(L,L+1)*DEXP(-ALFA*EB(1,L+1))/PF BOU00984 CALL PROFILE (-STOKE1,STOKIP) BOU00985 50 CONTINUE BOU00986 C BOU00987 C DELTA V=1 BELOW BOU00988 C BOU00989 AV1=0.45814D-43 BOU00990 STOKE=14.5442 BOU00991 STOKI=AV1*DEXP(-ALFA*EB(1,5))/PF BOU00992 CALL PROFILE (STOKE,STOKI) BOU00993 STOKIP=AV1*DEXP(-ALFA*EB(2,2))/PF BOU00994 CALL PROFILE (-STOKE,STOKIP) BOU00995 AV1=.32581D-43 BOU00996 STOKE=8.5755 BOU00997 STOKI=AV1*DEXP(-ALFA*EB(1,6))/PF BOU00998 CALL PROFILE (STOKE,STOKI) BOU00999 STOKIP=AV1*DEXP(-ALFA*EB(2,1))/PF BOU01000 CALL PROFILE (-STOKE,STOKIP) BOU01001 AV1=.41005D-43 BOU01002 STOKE=2.4395 BOU01003 STOKI=AV1*DEXP(-ALFA*EB(1,7))/PF BOU01004 CALL PROFILE (STOKE,STOKI) BOU01005 STOKIP=AV1*DEXP(-ALFA*EB(2,2))/PF BOU01006 CALL PROFILE (-STOKE,STOKIP) BOU01007 DO 60 N=1,NSRIUP BOU01008 60 RSIBB(N)=RSIBB(N)/TWOPIC/SLIT BOU01009 NSOL=NSRI+1 BOU01010 DO 70 I=1,NSOL BOU01011 K=(NSRI+1)+(I-1) BOU01012 70 RSI(I)=RSIBB(K) BOU01013 C BOU01014 C RSI - BOUND-BOUND CONTRIBUTION FOR POSITIVE FREQUENCY SHIFTS BOU01015 C BOU01016 RETURN BOU01017 END BOU01019 FUNCTION BGAUS (FNU,G0,DELTA,OMEGA0,TEMP) BGA01020 C BGA01021 C THIS IS DESYMMETRIZED GAUSSIAN PROFILE, G0 IS A ZEROTH MOMENT C FNU,DELTA AND OMEGA0 ARE IN CM-1 BGA01024 C G0 IN CGS BGA01025 C BGA01026 IMPLICIT DOUBLE PRECISION (A-H,O-Z) DATA BKW,TWOPI/0.6950304256D0,6.2831853D0/ BGA01027 D=G0/(DELTA*DSQRT(TWOPI)) BGA01028 DESYM=2./(1. + DEXP(-FNU/(BKW*TEMP))) BGA01029 FEXP=(FNU-OMEGA0)**2/(2.*DELTA**2) BGA01030 FEXP=DMIN1(FEXP,300.D0) BGA01031 BGAUS=D*DESYM*DEXP(-FEXP) BGA01032 RETURN BGA01033 END BGA01035