BLOCK DATA INPUT DAT00001 c ==================================== c Copyright, Aleksandra Borysow, 1989 c ==================================== c Clean F77 version of an old CYBER version c Changed to DOUBLE PRECISION TO HANDLE LARGE EXPONENSES c Program tested on VAX, results identical with CYBER, including b-b c 1-MAR-1989 : WORKS OK C DAT00002 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 modeling Titan's atmosphere: H2-N2 pairs", c Astrophysical Journal, vol. 303, p.495-510, (1986). C THESE DATA SERVE AS 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 c You'll find the results in file: output c IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /BLOCKIN/ TEMP,FNUMIN,FNUMAX,DNU DAT00012 COMMON /APP3/ SLIT,DX,NSRI,WNRMAX3,NS,NSRIUP,JSLIT DAT00013 DATA TEMP/70./ DATA SLIT/4.3/ c The parameters taken in subroutines BOUNDxx are specially chosen c for SLIT=4.3 cm-1; if a much smaller slit is chosen, parameters c WNRMAX, DX and NRSI have to be adjusted c NOTE: Take care that all b-b intensities are larger than zero c in case of separated lines, (extremelty small slit-width) put c RSI(i)=1.d-99 or so, to prevent a failure in LOG(RSI). DATA FNUMIN/10./ DATA FNUMAX/400./ DATA DNU/10./ END DAT00020 PROGRAM ADDEM C ADD00022 C ****************************************************************ADD00023 C PROGRAM PREPARED BY ALEKSANDRA BORYSOW; UNIVERSITY OF TEXAS AT AADD00024 C AND JOINT INSTITUTE FOR LABORATORY ASTROPHYSICS, UNIVERSITY OF CADD00025 C LAST revision DATE: 14 NOVEMBER 1988 ADD00026 C PROGRAM COMPATIBLE WITH PAPER: A. BORYSOW & L. FROMMHOLD; ADD00027 C ASTROPHYSICAL JOURNAL; VOL. 303, PP. 495-510, (1986). ADD00028 C ****************************************************************ADD00029 C PROGRAM GENERATES THE H2-N2, FREE-FREE, BOUND-FREE & B-B CONTRIBUTADD00030 C TO THE CIA SPECTRA ADD00031 C TAPE3 IS OUTPUT: HEADER PLUS ABSORPTION COEFF. ALPHA(NU) ADD00032 C ADD00033 IMPLICIT DOUBLE PRECISION (A-h,o-z) COMMON /BLOCKIN/ TEMP,FNUMIN,FNUMAX,DNU ADD00034 COMMON /APP3/ SLIT,DX,NSRI,WNRMAX3,NS,NSRIUP,JSLIT ADD00035 COMMON /RSILO/ RSILO(201) ADD00036 COMMON /BB/ OMEG,RSI,RSIGG,NSOL,ALFA,SCAL ADD00037 DIMENSION FREQ(601), ABSCOEF(601), ALFATOT(601) ADD00038 DIMENSION RSI(201), RSIGG(201), TT(2), SS(1), OMEG(201), AIG(201) ADD00039 Y(X,A,B,C)=A*DEXP((C*X+B)*X) ADD00040 C ADD00041 OPEN(UNIT=6,FILE='output', STATUS='unknown') NF=INT((FNUMAX-FNUMIN)/DNU+0.5)+1 ADD00042 IF (NF.GT.601) NF=601 ADD00043 FNUMAX=FNUMIN+DFLOAT(NF-1)*DNU ADD00044 WRITE (6, 180) TEMP,FNUMIN,FNUMAX,DNU,NF-1 ADD00046 CALL PARTSUM (TEMP) ADD00047 C ADD00048 C THE H2-N2 SPECTRA FOR 50-300K ADD00049 C ================= ADD00050 C ADD00051 C ONLY B-F + F-F TERMS INCLUDED IN MODELLING ADD00052 C TEMPERATURE INTERPOLATION GOOD FOR 45-300K ADD00053 C ADD00054 X=DLOG(TEMP) ADD00055 DO 10 I=1,NF ADD00056 FREQ(I)=FNUMIN+DFLOAT(I-1)*DNU ADD00057 ALFATOT(I)=0.0 ADD00058 10 ABSCOEF(I)=0. ADD00059 EPS=1.D-5 ADD00060 TT(1)=10. ADD00061 CALL BOUND32 (TEMP,RSI,NSOL) ADD00062 c print 666, (rsi(i), i=1, nsol) OK! as in Cyber 1-MAR-1989 18:12:10 666 format(' rsi(i), for L=32 H2',/,60(5e12.5,/)) DO 20 I=1,NSOL ADD00063 RSILO(I)=DLOG(RSI(I)*1.D80) ADD00064 20 OMEG(I)=DFLOAT(I-1)*DX ADD00065 CALL SPLINE (NSOL,1,0,EPS,OMEG,RSILO,TT,SS,SI,NR,RSIGG) ADD00066 C ADD00067 C THE LAMBDA1,LAMBDA2,LAMBDA,L = 2023 AND 0223 COMPONENTS: ADD00068 C 2023 - CORRESPONDS TO HYDROGEN'S SINGLE TRANSITIONS ADD00069 C 0223 - CORRESPONDS TO NITROGEN'S TRANSITIONS ADD00070 C HYDROGEN'S QUADRUPOLE ADD00071 C ADD00072 SCAL=33.74d0 ADD00073 IF (TEMP-145.d0) 30,30,40 ADD00074 30 S=Y(X,33.74*5.9340D-62,-1.5429D0,.15677D0) ADD00075 E=-Y(X,4.239838D-3,2.76328D0,-.38222D0) ADD00076 T1=Y(X,8.66209D-14,.59129D0,-.11062D0) ADD00077 T2=Y(X,4.576639D-13,-.20236D0,-.03646D0) ADD00078 T3=Y(X,1.564106D-13,.42334D0,-.10565D0) ADD00079 T4=Y(X,1.37791D-11,-1.61496D0,.19558D0) ADD00080 GO TO 50 ADD00081 40 S=Y(X,33.74*5.5440D-63,-0.58016D0,.05844D0) ADD00082 T1=Y(X,4.4024D-12,-.92834D0,.03489D0) ADD00083 T2=Y(X,2.0536D-13,-.0549D0,-.03709D0) ADD00084 E=-Y(X,3.79637D-6,4.41619D0,-.43833D0) ADD00085 T3=Y(X,1.1773D-13,.22391D0,-.06086D0) ADD00086 T4=Y(X,1.28155D-4,-6.30489D0,.47603D0) ADD00087 50 CONTINUE ADD00088 CALL ADDSPEC (S,E,T1,T2,T3,T4,TEMP,NF,FREQ,ABSCOEF,0,0,2,0,2,3) ADD00089 DO 60 I=1,NF ADD00090 60 ALFATOT(I)=ABSCOEF(I) ADD00091 C ADD00092 C NITROGEN'S QUADRUPOLE (0,4,4,5) ADD00093 C ADD00094 SCAL=34.9472D0 ADD00095 IF (TEMP-145.) 70,70,80 ADD00096 70 S=Y(X,34.9472*5.9340D-62,-1.5429D0,.15677D0) ADD00097 GO TO 90 ADD00098 80 S=Y(X,34.9472*5.5440D-63,-0.58016D0,.05844D0) ADD00099 90 CONTINUE ADD00100 CALL ADDSPEC (S,E,T1,T2,T3,T4,TEMP,NF,FREQ,ABSCOEF,0,0,0,2,2,3) ADD00101 DO 100 I=1,NF ADD00102 100 ALFATOT(I)=ALFATOT(I)+ABSCOEF(I) ADD00103 C ADD00104 C DOUBLE TRANSITIONS (2,2,3,3) ADD00105 C ADD00106 SCAL=3.3810D0 ADD00107 IF (TEMP-145.D0) 110,110,120 ADD00108 110 S=Y(X,3.3810*5.9340D-62,-1.5429D0,.15677D0) ADD00109 GO TO 130 ADD00110 120 S=Y(X,3.3810*5.5440D-63,-0.58016D0,.05844D0) ADD00111 130 CONTINUE ADD00112 CALL ADDSPEC (S,E,T1,T2,T3,T4,TEMP,NF,FREQ,ABSCOEF,0,0,2,2,3,3) ADD00113 DO 140 I=1,NF ADD00114 140 ALFATOT(I)=ALFATOT(I)+ABSCOEF(I) ADD00115 CALL BOUND54 (TEMP,RSI,NSOL) ADD00116 DO 150 I=1,NSOL ADD00117 RSILO(I)=DLOG(RSI(I)*1.D80) ADD00118 150 OMEG(I)=DFLOAT(I-1)*DX ADD00119 CALL SPLINE (NSOL,1,0,EPS,OMEG,RSILO,TT,SS,SI,NR,RSIGG) ADD00120 C ADD00121 C PARAMETERS FOR 4045 (HYDROGEN HEXADECAPOLE) COMPONENTS ADD00122 C ADD00123 S=Y(X,17.7052*9.7163D-65,-1.8937D0,2.01078D-1) ADD00124 E=-Y(X,9.8311D-7,4.9861D0,-4.61612D-1) ADD00125 T1=Y(X,2.5235D-15,1.6788D0,-1.86225D-1) ADD00126 T2=Y(X,3.27189D-13,-0.2089D0,-3.3745D-2) ADD00127 T3=Y(X,8.5546D-9,-3.8389D0,3.3309D-1) ADD00128 T4=Y(X,1.33844D-13,0.3935D0,-7.49592D-2) ADD00129 SCAL=17.7052d0 ADD00130 CALL ADDSPEC (S,E,T1,T2,T3,T4,TEMP,NF,FREQ,ABSCOEF,0,0,4,0,4,5) ADD00131 DO 160 I=1,NF ADD00132 160 ALFATOT(I)=ALFATOT(I)+ABSCOEF(I) ADD00133 C ADD00134 C PARAMETERS FOR 0445 - NITROGEN'S HEXADECAPOLE ADD00135 C ALL AS FOR 4045 EXCEPT PARAMETER S ADD00136 C ADD00137 S=Y(X,3181.456*9.7163D-65,-1.8937D0,2.01078D-1) ADD00138 SCAL=3181.456D0 ADD00139 CALL ADDSPEC (S,E,T1,T2,T3,T4,TEMP,NF,FREQ,ABSCOEF,0,0,0,4,4,5) ADD00140 DO 170 I=1,NF ADD00141 170 ALFATOT(I)=ALFATOT(I)+ABSCOEF(I) ADD00142 C ADD00143 WRITE(6,190) FREQ(1),FREQ(NF),FREQ(2)-FREQ(1),TEMP ADD00146 WRITE(6,200) (ALFATOT(I),I=1,NF) ADD00147 C ADD00148 180 FORMAT ( 34H1ABSORPTION SPECTRA OF H2 - N2 AT,F8.1, 2H K,/1X,43(ADD00149 11H=),/, 11H MIN.FREQ.=,F8.1, 5H CM-1,10X, 10HMAX.FREQ.=,F8.1, 5HADD00150 2 CM-1,10X, 15HFREQ.INCREMENT=,F8.2, 5H CM-1,5X, 2HIN,I5, 6H STEADD00151 3PS,//) ADD00152 190 FORMAT (/, 'ABSORPTION COEFFICIENT ALPHA(FNU), FROM', 1 F7.1, ' CM-1 TO', F7.1, ' CM-1', /, ' AT', F7.2, 1 ' CM-1 INCREMENTS, TEMPERATURE=', F7.2, 1 ' K', /, ' UNITS [CM-1 AMAGAT-2]',/) ADD00155 200 FORMAT (101(6E12.4,/)) END ADD00158 SUBROUTINE ADDSPEC (G0,EP,TAU1,TAU2,TAU5,TAU6,TEMP,NF,FREQ,ABSCOEFADD00159 1,MP,LIKE,LAMBDA1,LAMBDA2,LAMBDA,LVALUE) ADD00160 C ADD00161 C THIS PROGRAM GENERATES LISTING OF CIA TR ALFA(OMEGA) ADD00162 C IF EITHER LAMBDA1 OR LAMBDA2 EQUAL TO ZERO - SINGLE TRANSITIONS; ADD00163 C LAMBDA1 CORRESPONDS TO H2, LAMBDA2 TO N2 ADD00164 C DOUBLE TRANSITIONS ARE ASSUMED OTHERWISE. ADD00165 C LIKE=1 FOR LIKE SYSTEMS (AS H2-H2), SET LIKE=0 ELSE. ADD00167 C ADD00168 IMPLICIT DOUBLE PRECISION (A-h,o-z) LOGICAL COND ADD00169 COMMON /BB/ OMEG(201),RSI(201),RSIGG(201),NSOL,BETA,SCAL ADD00170 COMMON /RSILO/ RSILO(201) ADD00171 COMMON /APP3/ SLIT,DX,NSRI,WNRMAX3,NS,NSRIUP,JSLIT ADD00172 COMMON /H2PART/ Q,WH2(2),B0,D0,JRANGE1 ADD00173 COMMON /N2PART/ Q1,WN2(2),B01,D01,JRANGE2 ADD00174 DIMENSION ABSCOEF(NF), FREQ(NF) ADD00175 DATA CLOSCHM,BOLTZWN/2.68675484D19,.6950304/ ADD00176 DATA HBAR,PI,CLIGHT/1.054588757D-27,3.1415926535898,2.9979250D10/ ADD00177 EH2(I)=(B0-DFLOAT(I)*D0)*DFLOAT(I) ADD00178 PH2(J,T)=DFLOAT(2*J+1)*WH2(1+MOD(J,2))*DEXP(-1.4387859/T* 1 EH2(J*(J+1))) EN2(I)=(B01-DFLOAT(I)*D01)*DFLOAT(I) ADD00181 PN2(J,T)=DFLOAT(2*J+1)*WN2(1+MOD(J,2))*DEXP(-1.4387859/T* 1 EN2(J*(J+1))) TWOPIC=2.*PI*CLIGHT ADD00184 C IF (LIKE.NE.1) LIKE=0 ADD00186 CALIB=TWOPIC*((4.*PI**2)/(3.*HBAR*CLIGHT))*CLOSCHM**2 ADD00187 CALIB=CALIB/DFLOAT(1+LIKE) ADD00188 BETA=1./(BOLTZWN*TEMP) ADD00189 LIST=NF ADD00190 DO 10 I=1,LIST ADD00191 10 ABSCOEF(I)=0.0 ADD00192 C ADD00193 C ROTATIONAL SPECTRUM; DETAILED LISTING ******************* ADD00194 C ADD00195 WRITE (6,160) LAMBDA1,LAMBDA2,LAMBDA,LVALUE,G0,EP,TAU1,TAU2, 1 TAU5,TAU6 ADD00198 1,G0*BGAMA(0.,TAU1,TAU2,EP,TAU5,TAU6,TEMP) ADD00199 C ADD00200 IF ((LAMBDA1.EQ.0).OR.(LAMBDA2.EQ.0)) GO TO 70 ADD00201 JPLUSL=JRANGE1+MAX0(LAMBDA1,LAMBDA2) ADD00202 JPLU2=JRANGE2+MAX0(LAMBDA1,LAMBDA2) ADD00203 DO 60 I1=1,JRANGE1 ADD00204 J1=I1-1 ADD00205 DO 60 IP1=1,JPLUSL ADD00206 JP1=IP1-1 ADD00207 CG1S=CLEBSQR(J1,LAMBDA1,JP1) ADD00208 IF (CG1S) 60,60,20 ADD00209 20 P1=PH2(J1,TEMP)/Q ADD00210 OMEGA1=EH2(JP1*IP1)-EH2(J1*I1) ADD00211 DO 50 I2=1,JRANGE2 ADD00212 J2=I2-1 ADD00213 DO 50 IP2=1,JPLU2 ADD00214 JP2=IP2-1 ADD00215 CG2S=CLEBSQR(J2,LAMBDA2,JP2) ADD00216 IF (CG2S) 50,50,30 ADD00217 30 P2=PN2(J2,TEMP)/Q1 ADD00218 OMEGA2=EN2(JP2*IP2)-EN2(J2*I2) ADD00219 FAC=CALIB*P1*P2*CG1S*CG2S ADD00220 DO 40 I=1,LIST ADD00221 FRQ=FREQ(I)-OMEGA1-OMEGA2 ADD00222 WKI=FREQ(I)*(1.-DEXP(-BETA*FREQ(I))) ADD00223 WKF=WKI*FAC ADD00224 XBG=G0*BGAMA(FRQ,TAU1,TAU2,EP,TAU5,TAU6,TEMP) ADD00225 COND=ABS(FRQ).LE.WNRMAX3 ADD00226 IF (COND) XBG=XBG+SCAL*SPECFCT(FRQ,OMEG,RSILO,RSIGG,NSOL,ADD00227 1 BETA) ADD00228 ABSCOEF(I)=ABSCOEF(I)+XBG*WKF ADD00229 40 CONTINUE ADD00230 50 CONTINUE ADD00231 60 CONTINUE ADD00232 GO TO 150 ADD00233 70 IF (LAMBDA1.EQ.0) GO TO 110 ADD00234 C ADD00235 C SINGLE TRANSITIONS ON HYDROGEN'S ROTATIONAL FREQUENCIES. ADD00236 C ADD00237 JPLUSL=JRANGE1+LAMBDA ADD00238 DO 100 I=1,JRANGE1 ADD00239 J=I-1 ADD00240 DO 100 IP=1,JPLUSL ADD00241 JP=IP-1 ADD00242 CGS=CLEBSQR(J,LAMBDA,JP) ADD00243 c if(lambda1.eq.2) print 8877, j1,jp1,j2,jp2 8877 format(' J1,J1p,j2,j2p=', 4i5) IF (CGS) 100,100,80 ADD00244 80 P=PH2(J,TEMP)/Q ADD00245 OMEGA1=EH2(JP*IP)-EH2(J*I) ADD00246 FAC=CALIB*P*CGS ADD00247 DO 90 IQ=1,LIST ADD00248 FRQ=FREQ(IQ)-OMEGA1 ADD00249 WKI=FREQ(IQ)*(1.-DEXP(-BETA*FREQ(IQ))) ADD00250 WKF=WKI*FAC ADD00251 XBG=G0*BGAMA(FRQ,TAU1,TAU2,EP,TAU5,TAU6,TEMP) ADD00252 c if(j1.ne.0) go to 978 c if(lambda1.eq.2) print 555, freq(i), frq, xbg c555 format(' freq(i), frq', 2f10.2, c 1 ' g0*bgama(frq...) 32 H2', e12.5) c978 continue COND=ABS(FRQ).LE.WNRMAX3 ADD00253 IF (COND) XBG=XBG+SCAL*SPECFCT(FRQ,OMEG,RSILO,RSIGG,NSOL,BETADD00254 1 A) ADD00255 c if(cond ) go to 9999 c go to 7777 c9999 if(lambda1.eq.2 .and. j1.eq.0) print 3333, c 1 SCAL*SPECFCT(FRQ,OMEG,RSILO,RSIGG,NSOl,beta) c3333 format(' scal*specfct at FRQ ( 32 H2 b-b)', e12.5) c7777 continue ABSCOEF(IQ)=ABSCOEF(IQ)+XBG*WKF ADD00256 90 CONTINUE ADD00257 100 CONTINUE ADD00258 GO TO 150 ADD00259 C ADD00260 C SINGLE TRANSITIONS ON NITROGEN'S ROTATIONAL FREQUENCIES. ADD00261 C ADD00262 110 JPLU2=JRANGE2+LAMBDA ADD00263 DO 140 I=1,JRANGE2 ADD00264 J=I-1 ADD00265 DO 140 IP=1,JPLU2 ADD00266 JP=IP-1 ADD00267 CGS=CLEBSQR(J,LAMBDA,JP) ADD00268 IF (CGS) 140,140,120 ADD00269 120 P=PN2(J,TEMP)/Q1 ADD00270 OMEGA1=EN2(JP*IP)-EN2(J*I) ADD00271 FAC=CALIB*P*CGS ADD00272 DO 130 IQ=1,LIST ADD00273 FRQ=FREQ(IQ)-OMEGA1 ADD00274 WKI=FREQ(IQ)*(1.-DEXP(-BETA*FREQ(IQ))) ADD00275 WKF=WKI*FAC ADD00276 XBG=G0*BGAMA(FRQ,TAU1,TAU2,EP,TAU5,TAU6,TEMP) ADD00277 COND=ABS(FRQ).LE.WNRMAX3 ADD00278 IF (COND) XBG=XBG+SCAL*SPECFCT(FRQ,OMEG,RSILO,RSIGG,NSOL,BETADD00279 1 A) ADD00280 ABSCOEF(IQ)=ABSCOEF(IQ)+XBG*WKF ADD00281 130 CONTINUE ADD00282 140 CONTINUE ADD00283 150 CONTINUE ADD00284 WRITE(6,170) (ABSCOEF(I),I=1,LIST) ADD00285 RETURN ADD00287 C ADD00288 160 FORMAT (/, 32H LAMBDA1,LAMBDA2, LAMBDA,LVALUE=,2I3,2X,2I3, 12H COMADD00289 1PONENT .,/15X, 22HLINE SHAPE PARAMETERS:,6E12.3,5X, 5HG(0)=,E12.3ADD00290 2/) ADD00291 170 FORMAT ((1X,10E12.4,/)) ADD00292 C ADD00293 END ADD00294 FUNCTION CLEBSQR (L,LAMBDA,LP) CLE00295 C CLE00296 C SQUARE OF CLEBSCH-GORDAN COEFFICIENT (L,LAMBDA,0,0;LP,0) CLE00297 C FOR INTEGER ARGUMENTS ONLY CLE00298 C NOTE THAT LAMBDA SHOULD BE SMALL, MAYBE @10 OR SO. CLE00299 C CLE00300 IMPLICIT DOUBLE PRECISION (A-h,o-z) FC=DFLOAT(2*LP+1) CLE00301 GO TO 10 CLE00302 C CLE00303 ENTRY THREEJ2 CLE00304 C CLE00305 C THIS ENTRY RETURNS THE SQUARED 3-J SYMBOL L LAMBDA LP CLE00306 C 0 0 0 CLE00307 C INSTEAD OF THE CLEBSCH-GORDAN COEFFICIENT CLE00308 C (LIMITATION TO INTEGER ARGUMENTS ONLY) CLE00309 C CLE00310 C NOTE THAT THE THREE-J SYMBOLS ARE COMPLETELY SYMMETRIC IN THE CLE00311 C ARGUMENTS. IT WOULD BE ADVANTAGEOUS TO REORDER THE INPUT ARGUMENT CLE00312 C LIST SO THAT LAMBDA BECOMES THE SMALLEST OF THE 3 ARGUMENTS. CLE00313 C CLE00314 FC=1. CLE00315 10 CLEBSQR=0. CLE00316 IF (((L+LAMBDA).LT.LP).OR.((LAMBDA+LP).LT.L).OR.((L+LP).LT.LAMBDA)CLE00317 1) RETURN CLE00318 IF (MOD(L+LP+LAMBDA,2).NE.0) RETURN CLE00319 IF ((L.LT.0).OR.(LP.LT.0).OR.(LAMBDA.LT.0)) RETURN CLE00320 F=1./DFLOAT(L+LP+1-LAMBDA) CLE00321 IF (LAMBDA.EQ.0) GO TO 30 CLE00322 I1=(L+LP+LAMBDA)/2 CLE00323 I0=(L+LP-LAMBDA)/2+1 CLE00324 DO 20 I=I0,I1 CLE00325 20 F=F*DFLOAT(I)/DFLOAT(2*(2*I+1)) CLE00326 30 P=FC*F*FCTL(LAMBDA+L-LP)*FCTL(LAMBDA+LP-L) CLE00327 CLEBSQR=P/(FCTL((LAMBDA+L-LP)/2)*FCTL((LAMBDA+LP-L)/2))**2 CLE00328 RETURN CLE00329 C CLE00330 END CLE00331 FUNCTION FCTL (N) FCT00332 IMPLICIT DOUBLE PRECISION (A-h,o-z) P(Z)=((((-2.294720936D-4)/Z-(2.681327160D-3))/Z+(3.472222222D-3))/FCT00333 1Z+(8.333333333D-2))/Z+1. FCT00334 FCTL=1. FCT00335 IF (N.LE.1) RETURN FCT00336 IF (N.GT.15) GO TO 20 FCT00337 J=1 FCT00338 DO 10 I=2,N FCT00339 10 J=J*I FCT00340 FCTL=DFLOAT(J) FCT00341 RETURN FCT00342 20 Z=DFLOAT(N+1) FCT00343 FCTL=(DEXP(-Z)*(Z**(Z-0.5))*P(Z)*2.5066282746310) FCT00344 RETURN FCT00345 C FCT00346 END FCT00347 FUNCTION BGAMA (FNU,T1,T2,EPS,T3,T4,TEMP) BGA00348 C BGA00349 C EQUATION 13, SO-CALLED EBC MODEL, OF BORYSOW,TRAFTON,FROMMHOLD, BGA00350 C AND BIRNBAUM, ASTROPHYS.J., TO BE PUBLISHED (1985) BGA00351 C BGA00352 IMPLICIT DOUBLE PRECISION (A-h,o-z) REAL*8 K0 BGA00353 DATA PI,CLIGHT/3.1415926535898,2.99792458D10/ BGA00354 DATA HBAR,BOLTZ/1.0545887D-27,1.380662D-16/ BGA00355 P1(X)=((((((.0045813*X+.0360768)*X+.2659732)*X+1.2067492)*X+3.0899BGA00356 1424)*X+3.5156229)*X+1.) BGA00357 P2(X)=((((((.00000740*X+.00010750)*X+.00262698)*X+.03488590)*X+.23BGA00358 1069756)*X+.42278420)*X-.57721566) BGA00359 P3(X)=((((((.00032411*X+.00301532)*X+.02658733)*X+.15084934)*X+.51BGA00360 1498869)*X+.87890594)*X+.5) BGA00361 P4(X)=((((((-.00004686*X-.00110404)*X-.01919402)*X-.18156897)*X-.6BGA00362 17278579)*X+.15443144)*X+1.) BGA00363 P5(X)=((((((.00053208*X-.00251540)*X+.00587872)*X-.01062446)*X+.02BGA00364 1189568)*X-.07832358)*X+1.25331414) BGA00365 P6(X)=((((((-.00068245*X+.00325614)*X-.00780353)*X+.01504268)*X-.0BGA00366 13655620)*X+.23498619)*X+1.25331414) BGA00367 C BGA00368 OMEGA=2.*PI*CLIGHT*FNU BGA00369 T0=HBAR/(2.*BOLTZ*TEMP) BGA00370 Z=DSQRT((1.+(OMEGA*T1)**2)*(T2*T2+T0*T0))/T1 BGA00371 IF (Z-2.) 10,10,20 BGA00372 10 XK1=Z*Z*DLOG(Z/2.)*P3((Z/3.75)**2)+P4((Z/2.)**2) BGA00373 GO TO 30 BGA00374 20 XK1=dSQRT(Z)*DEXP(-Z)*P6(2./Z) BGA00375 30 ZP=DSQRT((1.+(OMEGA*T4)**2)*(T3*T3+T0*T0))/T4 BGA00376 IF (ZP-2.) 40,40,50 BGA00377 40 K0=-DLOG(ZP/2.)*P1((ZP/3.75)**2)+P2((ZP/2.)**2) BGA00378 GO TO 60 BGA00379 50 K0=DEXP(-ZP)*P5(2./ZP)/DSQRT(ZP) BGA00380 60 BGAMA=((T1/PI)*DEXP(T2/T1+T0*OMEGA)*XK1/(1.+(T1*OMEGA)**2)+ BGA00381 1 EPS*(T3/PI)*DEXP(T3/T4+T0*OMEGA)*K0)/(1.+EPS) BGA00382 RETURN BGA00383 C BGA00384 END BGA00385 SUBROUTINE PARTSUM (TEMP) PAR00386 C PAR00387 C H2 ROTATIONAL PARTITION SUM Q = Q(T). PAR00388 C PAR00389 IMPLICIT DOUBLE PRECISION (A-h,o-z) COMMON /H2PART/ Q,WH2(2),B0,D0,JRANGE1 PAR00390 COMMON /N2PART/ Q1,WN2(2),B01,D01,JRANGE2 PAR00391 DATA B0,D0,WH2(1),WH2(2)/59.3392,0.04599,1.,3./ PAR00392 DATA B01,D01,WN2(1),WN2(2)/1.98957,.58D-5,2.,1./ PAR00393 EH2(I)=(B0-DFLOAT(I)*D0)*FLOAT(I) PAR00394 PH2(J,T)=DFLOAT(2*J+1)*WH2(1+MOD(J,2))*DEXP(-1.4387859* 1 EH2(J*(J+1))/ T) PAR00396 EN2(I)=(B01-DFLOAT(I)*D01)*FLOAT(I) PAR00397 PN2(J,T)=DFLOAT(2*J+1)*WN2(1+MOD(J,2))*DEXP(-1.4387859* 1 EN2(J*(J+1))/T) PAR00398 PAR00399 C PAR00400 C Q,B01,D01,WN2 - PARTITION FCT., ROT.CONSTANTS, WEIGHTS FOR N2 PAR00401 C Q,B0,D0,WH2 - PARTITION FCT., ROT.CONSTANTS, WEIGHTS FOR H2 PAR00402 C PAR00403 Q=0. PAR00404 J=0 PAR00405 10 DQ=PH2(J,TEMP) PAR00406 Q=Q+DQ PAR00407 J=J+1 PAR00408 IF (DQ.GT.Q/900.) GO TO 10 PAR00409 JRANGE1=J PAR00410 Q1=0. PAR00412 J=0 PAR00413 20 DQ1=PN2(J,TEMP) PAR00414 Q1=Q1+DQ1 PAR00415 J=J+1 PAR00416 IF (DQ1.GT.Q1/900.) GO TO 20 PAR00417 JRANGE2=J PAR00418 RETURN PAR00420 C PAR00421 30 FORMAT (/, 40H ROTATIONAL PARTITION FUNCTION OF H2: Q=,F8.2,10X, PAR00422 17HJ MAX =,I3/) PAR00423 40 FORMAT (/, 40H ROTATIONAL PARTITION FUNCTION OF N2: Q=,F8.2,10X, PAR00424 17HJ MAX =,I3/) PAR00425 C PAR00426 END PAR00427 SUBROUTINE PROFILE (X,Y) PRO00428 IMPLICIT DOUBLE PRECISION (A-h,o-z) COMMON /BL3/ RSI(401) PRO00429 C PRO00430 C ATRIANGULAR SLIT FUNCTION IS USED. PRO00431 C PRO00432 COMMON /APP3/ SLIT,DX,NSRI,WNRMAX3,NS,NSRIUP,JSLIT PRO00433 IF (Y) 50,60,10 PRO00434 10 CONTINUE PRO00435 X0=(NSRI+1.)+X/DX PRO00436 NC=X0 PRO00437 N1=NC+1 PRO00438 SLOPE=Y/SLIT PRO00439 NU=X0-NS PRO00440 IF (NU.LT.1) NU=1 PRO00441 IF (NU.GT.NSRIUP) RETURN PRO00442 NO=X0+NS PRO00443 IF (NO.GT.NSRIUP) NO=NSRIUP PRO00444 IF (NO.LT.1) RETURN PRO00445 IF (NC.GT.NSRIUP) NC=NSRIUP PRO00446 IF (NC.LE.1) GO TO 30 PRO00447 DO 20 I=NU,NC PRO00448 XI=(I-1.)*DX-WNRMAX3 PRO00449 DR=SLOPE*(XI-(X-SLIT)) PRO00450 IF (DR.LE.0.) GO TO 20 PRO00451 RSI(I)=RSI(I)+DR PRO00452 20 CONTINUE PRO00453 30 IF (NC.GE.NSRIUP) RETURN PRO00454 IF (N1.LT.1) N1=1 PRO00455 DO 40 I=N1,NO PRO00456 XI=(I-1.)*DX-WNRMAX3 PRO00457 DR=Y-SLOPE*(XI-X) PRO00458 IF (DR.LE.0.) GO TO 40 PRO00459 RSI(I)=RSI(I)+DR PRO00460 40 CONTINUE PRO00461 RETURN PRO00462 50 WRITE(6,70) SLIT PRO00463 60 CONTINUE PRO00464 RETURN PRO00465 C PRO00466 70 FORMAT (/, 30H A TRIANGULAR SLIT FUNCTION OF,F6.3, 23H CM-1 HALFWIPRO00467 1DTH IS USED,/) PRO00468 C PRO00469 END PRO00470 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) DATA SCALEF/1.D-80/ DELY(I)=(PHI(I+1)-PHI(I))/(OMEGA(I+1)-OMEGA(I)) c 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 INTERPOLATION IS BEING MADE BELOW: 30 I = 2 90 IF (f(1)-OMEGA(I))200,200,100 100 I = I+1 GOTO 90 200 I = I-1 HT1=f(1)-OMEGA(I) HT2=f(1)-OMEGA(I+1) PROD=HT1*HT2 S3=(PHI2(I+1)-PHI2(I))/(OMEGA(I+1)-OMEGA(I)) SS2=PHI2(I)+HT1*S3 DELSQS=(PHI2(I)+PHI2(I+1)+SS2)/6. gp(1)=PHI(I)+HT1*DELY(I)+PROD*DELSQS SPECFCT=DEXP( tfac + gp(1) ) * SCALEF RETURN SPE00491 END SPE00493 SUBROUTINE SPLINE (L,M,K,EPS,X,Y,T,SS,SI,NR,S2) SPL00494 C SPL00495 C SPLINE INTERPOLATION AND QUADRATURE, THIRD ORDER AFTER GREVILLE. SPL00496 C INPUT ARGUMENTS L...Y, OUTPUT SS...NR. SPL00497 C L DATA POINTS X(1), Y(1) ... X(L),Y(L) SPL00498 C EPS=ERROR CRITERION, TYPICALLY EPS=1.E-5 FOR 5 DECI. PLACES ACCURASPL00499 C M ARGUMENTS T(1)..T(M) FOR WHICH FUNCTION VALUES SS(1)..SS(M), FORSPL00500 C K=0; OR FIRST OR SECOND DERIVATIVE FOR K=1 OR -1, RESPECTIVELY. SPL00501 C NOTE THAT M HAS TO BE AT LEAST EQUAL TO 1. SPL00502 C SI=INTEGRAL (OVER WHOLE INTERVAL) FOR K=2 ONLY. SPL00503 C FOR 'NATURAL' SPLINE FUNCTIONS, S2(1)=S2(L)=0. MUST BE INPUT*NOTE*SPL00504 C N0 INDICATES THE NUMBER OF OUT-OF-RANGE CALLS. X(1)@T(I)@X(L) SPL00505 C EXTRAPOLATE WITH CAUTION. (ASSUMPTION D2Y/DX2 = 0.) SPL00506 C S2(I) IS THE 2ND DERIVATIVE AT X=X(I) AND IS COMPUTED INTERNALLY. SPL00507 C SPL00508 IMPLICIT DOUBLE PRECISION (A-h,o-z) DIMENSION X(L), Y(L), T(M), SS(M), S2(L) SPL00509 DELY(I)=(Y(I+1)-Y(I))/(X(I+1)-X(I)) SPL00510 B(I)=(X(I)-X(I-1))*0.5/(X(I+1)-X(I-1)) SPL00511 C(I)=3.*(DELY(I)-DELY(I-1))/(X(I+1)-X(I-1)) SPL00512 N=L SPL00513 N1=N-1 SPL00514 NR=0 SPL00515 DO 10 I=2,N1 SPL00516 10 S2(I)=C(I)/1.5 SPL00517 OMEGA=1.0717968 SPL00518 IC=0 SPL00519 C SPL00520 C 'NATURAL' SPLINE FUNCTIONS OF THIRD ORDER. SPL00521 C SPL00522 S2(N)=0. SPL00523 S2(1)=S2(N) SPL00524 20 ETA=0. SPL00525 IC=IC+1 SPL00526 SM=ABS(S2(1)) SPL00527 DO 30 I=2,N SPL00528 IF (ABS(S2(I)).GT.SM) SM=ABS(S2(I)) SPL00529 30 CONTINUE SPL00530 EPSI=EPS*SM SPL00531 DO 50 I=2,N1 SPL00532 W=(C(I)-B(I)*S2(I-1)-(0.5-B(I))*S2(I+1)-S2(I))*OMEGA SPL00533 IF (ABS(W)-ETA) 50,50,40 SPL00534 40 ETA=ABS(W) SPL00535 50 S2(I)=S2(I)+W SPL00536 IF (ETA-EPSI) 60,20,20 SPL00537 ENTRY IXPOLAT SPL00538 C SPL00539 C THIS ENTRY USEFUL WHEN ITERATION PREVIOUSLY COMPLETED SPL00540 C SPL00541 N=L SPL00542 N1=N-1 SPL00543 NR=0 SPL00544 IC=-1 SPL00545 60 IF (K.EQ.2) GO TO 260 SPL00546 GO TO 70 SPL00547 70 DO 250 J=1,M SPL00548 I=1 SPL00549 IF (T(J)-X(1)) 110,210,80 SPL00550 80 IF (T(J)-X(N)) 100,190,150 SPL00551 90 IF (T(J)-X(I)) 200,210,100 SPL00552 100 I=I+1 SPL00553 GO TO 90 SPL00554 110 NR=NR+1 SPL00555 HT1=T(J)-X(1) SPL00556 HT2=T(J)-X(2) SPL00557 YP1=DELY(1)+(X(1)-X(2))*(2.*S2(1)+S2(2))/6. SPL00558 IF (K) 140,130,120 SPL00559 120 SS(J)=YP1+HT1*S2(1) SPL00560 GO TO 250 SPL00561 130 SS(J)=Y(1)+YP1*HT1+S2(1)*HT1*HT1/2. SPL00562 GO TO 250 SPL00563 140 SS(J)=S2(I) SPL00564 GO TO 250 SPL00565 150 HT2=T(J)-X(N) SPL00566 HT1=T(J)-X(N1) SPL00567 NR=NR+1 SPL00568 YPN=DELY(N1)+(X(N)-X(N1))*(S2(N1)+2.*S2(N))/6. SPL00569 IF (K) 180,170,160 SPL00570 160 SS(J)=YPN+HT2*S2(N) SPL00571 GO TO 250 SPL00572 170 SS(J)=Y(N)+YPN*HT2+S2(N)*HT2*HT2/2. SPL00573 GO TO 250 SPL00574 180 SS(J)=S2(N) SPL00575 GO TO 250 SPL00576 190 I=N SPL00577 200 I=I-1 SPL00578 210 HT1=T(J)-X(I) SPL00579 HT2=T(J)-X(I+1) SPL00580 PROD=HT1*HT2 SPL00581 S3=(S2(I+1)-S2(I))/(X(I+1)-X(I)) SPL00582 SS2=S2(I)+HT1*S3 SPL00583 DELSQS=(S2(I)+S2(I+1)+SS2)/6. SPL00584 IF (K) 240,220,230 SPL00585 220 SS(J)=Y(I)+HT1*DELY(I)+PROD*DELSQS SPL00586 GO TO 250 SPL00587 230 SS(J)=DELY(I)+(HT1+HT2)*DELSQS+PROD*S3/6. SPL00588 GO TO 250 SPL00589 240 SS(J)=SS2 SPL00590 250 CONTINUE SPL00591 260 SI=0. SPL00592 DO 270 I=1,N1 SPL00593 H=X(I+1)-X(I) SPL00594 270 SI=SI+0.5*H*(Y(I)+Y(I+1))-H**3*(S2(I)+S2(I+1))/24. SPL00595 IF (K.EQ.2) NR=IC SPL00596 RETURN SPL00597 C SPL00598 END SPL00599 SUBROUTINE BOUND32 (TEMP,RSI,NSOL) BOU00600 IMPLICIT DOUBLE PRECISION (A-h,o-z) DIMENSION EB(2,8), A(8,8), RSI(201) BOU00601 C BOU00602 C EB(I,K) - BOUND ENERGIES FOR VIB.NR.V=I-1, ROT.NR J=K-1 BOU00603 C A(I,J)=( C(I,L,J)**2) * (MTX.EL. (I,BETA,J)**2 ) BOU00604 C I,J-DENOTE ROTATIONAL STATES FOR VIB. LEVEL V=0 BOU00605 C SEE TABLE 10 OF AP. J. VOL.303, P.508 (1986) BOU00606 C BOU00607 COMMON /APP3/ SLIT,DX,NSRI,WNRMAX3,NS,NSRIUP,JSLIT BOU00608 COMMON /BL3/ RSIBB(401) BOU00609 DATA (EB(1,J),J=1,8)/-17.0467,-15.9605,-13.8056,-10.6201,-6.4709,-BOU00610 11.4752,4.0588,10.1923/ BOU00611 DATA (EB(2,J),J=1,8)/-0.2735,0.,0.,0.,0.,0.,0.,0./ BOU00612 DATA TWOPIC/1.88365183D11/ BOU00613 C BOU00614 C NSRI - HOW MANY POINTS FOR B-B SPECTRAL FUNCTION TO BE GIVEN BOU00615 C WNRMAX3 - THE FREQUENCY RANGE OF B-B CONTRIBUTION BOU00616 C SLIT- THE HALFWIDTH OF THE SPECTRAL PROFILE CONVOLUTED WITH BOU00617 C B-B SPECTRUM, IN [CM-1]. BOU00618 C BOU00619 NSRI=129 BOU00620 WNRMAX3=20.64 BOU00621 c This is an estimated range of frequencies (from 0 to WNRMAX) c from the line center, in cm-1, where bound-bound intensities matter c at slitwidth approx. 4.3 cm-1 NSRIUP=2*NSRI+1 BOU00622 DX=WNRMAX3/DFLOAT(NSRI) BOU00623 DO 10 I=1,401 BOU00624 10 RSIBB(I)=0.0 BOU00625 NS=INT(SLIT/DX) BOU00626 DO 20 I=1,8 BOU00627 DO 20 J=1,8 BOU00628 C BOU00629 C A(I,J) FOR L=3 CONTRIBUTION (REDUCED BETA FUNCTION) BOU00630 C BOU00631 20 A(I,J)=0.0 BOU00632 A(1,4)=.1722 BOU00633 A(2,3)=.2272 BOU00634 A(2,5)=.2767 BOU00635 A(3,4)=.2234 BOU00636 A(3,6)=.3472 BOU00637 A(4,5)=.2519 BOU00638 A(4,7)=.3784 BOU00639 A(5,6)=.2686 BOU00640 A(5,8)=.4007 BOU00641 A(6,7)=.2639 BOU00642 A(7,8)=.2550 BOU00643 DO 30 I=1,8 BOU00644 DO 30 J=1,8 BOU00645 30 A(I,J)=A(I,J)*1.D-41 BOU00646 ALFA=1./(0.69519*TEMP) BOU00647 RM=1.8657*1.672649D-24 BOU00648 PI=3.141592654 BOU00649 PF=((RM*(1.380662D-16)*TEMP*2.*PI/(6.626176D-27)**2)**1.5) BOU00650 DO 50 L=1,7 BOU00651 STOKE1=EB(1,L+1)-EB(1,L) BOU00652 C BOU00653 C FREQUENCY SHIFT FOR DELTA J=+1,DELTA V=0 BOU00654 C BOU00655 IF (L.GT.5) GO TO 40 BOU00656 STOKE3=EB(1,L+3)-EB(1,L) BOU00657 C BOU00658 C STOKE3- FREQUENCY SHIFT FOR DELTA J=+3, DELTA V=0 BOU00659 C BOU00660 STOKI=A(L,L+3)*DEXP(-ALFA*EB(1,L))/PF BOU00661 CALL PROFILE (STOKE3,STOKI) BOU00662 STOKIP=A(L,L+3)*DEXP(-ALFA*EB(1,L+3))/PF BOU00663 CALL PROFILE (-STOKE3,STOKIP) BOU00664 40 STOKI=A(L,L+1)*DEXP(-ALFA*EB(1,L))/PF BOU00665 CALL PROFILE (STOKE1,STOKI) BOU00666 STOKIP=A(L,L+1)*DEXP(-ALFA*EB(1,L+1))/PF BOU00667 CALL PROFILE (-STOKE1,STOKIP) BOU00668 50 CONTINUE BOU00669 AV1=0.3255D-43 BOU00670 STOKE=EB(2,1)-EB(1,4) BOU00671 STOKI=AV1*DEXP(-ALFA*EB(1,4))/PF BOU00672 CALL PROFILE (STOKE,STOKI) BOU00673 STOKIP=AV1*DEXP(-ALFA*EB(2,1))/PF BOU00674 CALL PROFILE (-STOKE,STOKIP) BOU00675 DO 60 N=1,NSRIUP BOU00676 60 RSIBB(N)=RSIBB(N)/TWOPIC/SLIT BOU00677 NSOL=NSRI+1 BOU00678 DO 70 I=1,NSOL BOU00679 K=(NSRI+1)+(I-1) BOU00680 70 RSI(I)=RSIBB(K) BOU00681 C BOU00682 C RSI - BOUND-BOUND CONTRIBUTION FOR POSITIVE FREQUENCY SHIFTS BOU00683 C BOU00684 WRITE(6,80) SLIT BOU00685 80 FORMAT ( 10H SLITWIDTH,F10.4, 13H CM-1 ASSUMED,/) BOU00688 RETURN BOU00686 END BOU00690 SUBROUTINE BOUND54 (TEMP,RSI,NSOL) BOU00691 IMPLICIT DOUBLE PRECISION (A-h,o-z) DIMENSION EB(2,8), A(8,8), RSI(201) BOU00692 C BOU00693 C EB(I,K) - BOUND ENERGIES FOR VIB.NR.V=I-1, ROT.NR J=K-1 BOU00694 C A(I,J)=( C(I,L,J)**2) * (MTX.EL. (I,BETA,J)**2 ) BOU00695 C I,J-DENOTE ROTATIONAL STATES FOR VIB. LEVEL V=0 BOU00696 C SEE TABLE 10 OF AP. J. VOL.303, P.508 (1986) BOU00697 C BOU00698 COMMON /APP3/ SLIT,DX,NSRI,WNRMAX,NS,NSRIUP,JSLIT BOU00699 COMMON /BL3/ RSIBB(401) BOU00700 DATA (EB(1,J),J=1,8)/-17.0467,-15.9605,-13.8056,-10.6201,-6.4709,-BOU00701 11.4752,4.0588,10.1923/ BOU00702 DATA (EB(2,J),J=1,8)/-0.2735,0.,0.,0.,0.,0.,0.,0./ BOU00703 DATA TWOPIC/1.88365183D11/ BOU00704 C BOU00705 C NSRI - HOW MANY POINTS FOR B-B SPECTRAL FUNCTION TO BE GIVEN BOU00706 C WNRMAX - THE FREQUENCY RANGE OF B-B CONTRIBUTION BOU00707 C SLIT- THE HALFWIDTH OF THE SPECTRAL PROFILE CONVOLUTED WITH BOU00708 C B-B SPECTRUM, IN [CM-1]. BOU00709 C BOU00710 NSRI=175 BOU00711 WNRMAX=28.16 BOU00712 NSRIUP=2*NSRI+1 BOU00713 DX=WNRMAX/DFLOAT(NSRI) BOU00714 c This is an estimated range of frequencies (from 0 to WNRMAX) c from the line center, in cm-1, where bound-bound intensities matter c at slitwidth approx. 4.3 cm-1 DO 10 I=1,401 BOU00715 10 RSIBB(I)=0.0 BOU00716 NS=INT(SLIT/DX) BOU00717 DO 20 I=1,8 BOU00718 DO 20 J=1,8 BOU00719 C BOU00720 C A(I,J) FOR L=5 CONTRIBUTION (REDUCED BETA FUNCTION) BOU00721 C BOU00722 20 A(I,J)=0.0 BOU00723 A(1,6)=0.08605 BOU00724 A(2,5)=0.1294 BOU00725 A(2,7)=0.1178 BOU00726 A(3,4)=0.1506 BOU00727 A(3,6)=0.1065 BOU00728 A(3,8)=0.1360 BOU00729 A(4,5)=0.1116 BOU00730 A(4,7)=0.1008 BOU00731 A(5,6)=0.1037 BOU00732 A(5,8)=0.09753 BOU00733 A(6,7)=0.09175 BOU00734 A(7,8)=0.08088 BOU00735 DO 30 I=1,8 BOU00736 DO 30 J=1,8 BOU00737 30 A(I,J)=A(I,J)*1.D-44 BOU00738 ALFA=1./(0.69519*TEMP) BOU00739 RM=1.8657*1.672649D-24 BOU00740 PI=3.141592654 BOU00741 PF=((RM*(1.380662D-16)*TEMP*2.*PI/(6.626176D-27)**2)**1.5) BOU00742 DO 60 L=1,7 BOU00743 STOKE1=EB(1,L+1)-EB(1,L) BOU00744 C BOU00745 C FREQUENCY SHIFT FOR DELTA J=+1,DELTA V=0 BOU00746 C BOU00747 IF (L.GT.5) GO TO 50 BOU00748 STOKE3=EB(1,L+3)-EB(1,L) BOU00749 C BOU00750 C STOKE3- FREQUENCY SHIFT FOR DELTA J=+3, DELTA V=0 BOU00751 C BOU00752 IF (L.GT.3) GO TO 40 BOU00753 STOKE5=EB(1,L+5)-EB(1,L) BOU00754 STOKI=A(L,L+5)*DEXP(-ALFA*EB(1,L))/PF BOU00755 CALL PROFILE (STOKE5,STOKI) BOU00756 STOKIP=A(L,L+5)*DEXP(-ALFA*EB(1,L+5))/PF BOU00757 CALL PROFILE (-STOKE5,STOKIP) BOU00758 40 STOKI=A(L,L+3)*DEXP(-ALFA*EB(1,L))/PF BOU00759 CALL PROFILE (STOKE3,STOKI) BOU00760 STOKIP=A(L,L+3)*DEXP(-ALFA*EB(1,L+3))/PF BOU00761 CALL PROFILE (-STOKE3,STOKIP) BOU00762 50 STOKI=A(L,L+1)*DEXP(-ALFA*EB(1,L))/PF BOU00763 CALL PROFILE (STOKE1,STOKI) BOU00764 STOKIP=A(L,L+1)*DEXP(-ALFA*EB(1,L+1))/PF BOU00765 CALL PROFILE (-STOKE1,STOKIP) BOU00766 60 CONTINUE BOU00767 DO 70 N=1,NSRIUP BOU00768 70 RSIBB(N)=RSIBB(N)/TWOPIC/SLIT BOU00769 NSOL=NSRI+1 BOU00770 DO 80 I=1,NSOL BOU00771 K=(NSRI+1)+(I-1) BOU00772 80 RSI(I)=RSIBB(K) BOU00773 c print 333, (rsi(i), i=1, nsol) 333 format(' RSI(i)',/,101(6e12.3,/)) C BOU00774 C RSI - BOUND-BOUND CONTRIBUTION FOR POSITIVE FREQUENCY SHIFTS BOU00775 C BOU00776 RETURN BOU00777 END BOU00779