        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
