        BLOCK DATA INPUT                                                  DAT00001

c	=================================================================
c	Copyright, Aleksandra Borysow, 1988
c	=================================================================
c
c	Clean F77 version of an old CYBER program
c	Runs on VAX ( tested 17-NOV-1988) OK
C                                                                       DAT00002
c	==================================================
c	If you use this program, please cite original work:
c       A. Borysow and L. Frommhold,
c       "Collision induced rototranslational absorption spectra of
c       CH4-CH4 pairs at temperatures from 50 to 300K",
c	Astrophysical Journal, vol.318, p.940-943, (1987). 
c	==================================================

C     THESE DATA SERVE AS THE INPUT TO THE MAIN PROGRAM ADDEM.          DAT00003
C     TEMP = TEMPERATURE IN KELVIN, SHOULD BE BETWEEN 40. AND 300.      DAT00004
C     FNUMIN = LOWEST FREQUENCY IN CM-1, FOR LISTING OF ALPHA(FNU)      DAT00005
C     FNUMAX = HIGHEST FREQUENCY IN CM-1, FOR LISTING OF ALPHA(FNU)     DAT00006
C     DNU = FREQUENCY INCREMENT IN CM-1. DNU SHOULD BE CHOSEN SO        DAT00007
C           THAT NOT MORE THAN 600 STEPS ARE NEEDED TO GO FROM          DAT00008
C           FNUMIN TO FNUMAX (ELSE ARRAY DIMENSIONS OF FREQ,ABSCOEF     DAT00009
C           MUST BE ADJUSTED IN ADDEM).                                 DAT00010
C	LIKE = 1 (fixed)
C                                                                       DAT00011
      implicit double precision (a-h,o-z)
      CHARACTER*10 LGAS                                                 DAT00012
      COMMON /BLOCKIN/ TEMP,FNUMIN,FNUMAX,DNU                           DAT00013
      COMMON /LIKE/ LIKE,LGAS                                           DAT00014

c     USER: You may change data only here:

      DATA TEMP/163./                                                   DAT00015
      DATA FNUMIN/0./                                                   DAT00016
      DATA FNUMAX/600./                                                 DAT00017
      DATA DNU/10./                                                     DAT00018
      DATA LIKE/1/                                                      DAT00019
      DATA LGAS/' CH4 - CH4' /                                          DAT00020
      END                                                               DAT00022

        PROGRAM ADDEM 
C                                                                       ADD00024
C       PROGRAM GENERATES THE CH4-CH4 CIA SPECTRA at temperatures 50-300K

C       WRITTEN BY ALEKSANDRA BORYSOW, University of Texas  at AUSTIN 
c	(MARCH 25, 1987); now at Joint Institute of Laboratory Astrophysics, 
c	University of Colorado, Boulder.
C       Program compatible with work described in the paper by
c	A. Borysow & L. Frommhold;
c	ASTROPHYSICAL JOURNAL, VOL.318, pp. 940-943, 1987.
C                                                                       ADD00029
C       File CORRCH4. contains the STATISTICAL (NUCLEAR) CORRECTIONS    ADD00030
C                                                                       ADD00031
	implicit double precision (a-h,o-z)
	character*10 LGAS
      COMMON /BLOCKIN/ TEMP,FNUMIN,FNUMAX,DNU                           ADD00032
      COMMON /ROT/ AMUL0,DELTA                                          ADD00033
      COMMON /STAT/ QW3(51,3),QW4(51,4)                                 ADD00034
      COMMON /LIKE/ LIKE,LGAS                                           ADD00035
      DIMENSION FREQ(601), ABSCOEF(601), ALFATOT(601)                   ADD00036
      DATA OMEGA0,DELOM/2.599,0.0071/                                   ADD00037
      DATA PHI,DELPHI/-8.19,0.033/                                      ADD00038
      Y(X,A,B,C)=A*DEXP((C*X+B)*X)                                      ADD00039
C                                                                       ADD00040
	OPEN(UNIT=33, FILE='corrch4', STATUS='OLD')
	OPEN(UNIT=6, FILE='output',   STATUS='unknown')

      NF=INT((FNUMAX-FNUMIN)/DNU+0.5)+1                                 ADD00041
      IF (NF.GT.601) NF=601                                             ADD00042
      FNUMAX=FNUMIN+FLOAT(NF-1)*DNU                                     ADD00043
      WRITE (6, 80) LGAS,TEMP,FNUMIN,FNUMAX,DNU,NF-1                        
      CALL PARTSUM (TEMP)                                               ADD00045

C     THE CH4 - CH4 SPECTRA   FOR 50-300K                               ADD00047
C     =================                                                 ADD00048
C                                                                       ADD00049
      X=DLOG(TEMP)                                                      ADD00050
      DO 10 I=1,NF                                                      ADD00051
         FREQ(I)=FNUMIN+FLOAT(I-1)*DNU                                  ADD00052
         ALFATOT(I)=0.0                                                 ADD00053
   10 ABSCOEF(I)=0.                                                     ADD00054
C                                                                       ADD00055
C      READ TAPE 33 THE QUANTUM CORRECTIONS (12 Y(J,JP,LAM)-1)          ADD00056
C                                                                       ADD00057
      REWIND 33                                                         ADD00058
      DO 30 N=1,3                                                       ADD00059
         DO 20 I=1,5                                                    ADD00060
            K=10*(I-1)+1                                                ADD00061
            K1=10*I                                                     ADD00062
            READ (33,100) (QW3(J,N),J=K,K1)                             ADD00063
   20    CONTINUE                                                       ADD00064
         READ (33,90) QW3(51,N)                                         ADD00065
   30 CONTINUE                                                          ADD00066
      DO 50 N=1,4                                                       ADD00067
         DO 40 I=1,5                                                    ADD00068
            K=(I-1)*10+1                                                ADD00069
            K1=10*I                                                     ADD00070
            READ (33,100) (QW4(J,N),J=K,K1)                             ADD00071
   40    CONTINUE                                                       ADD00072
         READ (33,100) QW4(51,N)                                        ADD00073
   50 CONTINUE                                                          ADD00074
	close(33)
C                                                                       ADD00075
C     THESE STATISTICAL FACTORS ARE VALID ONLY FOR DELTA J>0            ADD00076
C     FOR DELTA J<0 THESE WILL BE EQUAL TO 1.                           ADD00077
C     METHANE INDUCED COMPONENTS                                        ADD00078
C     ===============================                                   ADD00079
C     OCTOPOLE- INDUCED TERM (43)                                       ADD00080
C                                                                       ADD00081
      AMUL0=OMEGA0                                                      ADD00082
      DELTA=DELOM                                                       ADD00083
      S=Y(X,0.31933D-51/7.,-8.24624D0,0.73213D0)                            
      E=Y(X,-0.61214D0,0.18496D0,-0.04047D0)
      T1=Y(X,0.36548D-11,-0.63873D0,0.00818D0) 
      T2=Y(X,0.16277D-11,-0.52965D0,0.00404D0)                           
      T3=Y(X,0.46747D-11,-0.73107D0,0.01892D0)                            
      T4=Y(X,0.20588D-11,-0.04288D0,-0.04548D0)                             
      CALL ADSPEC1 (S,E,T1,T2,T3,T4,TEMP,NF,FREQ,ABSCOEF,0,LIKE,0,3,3,4)
      DO 60 I=1,NF                                                      ADD00091
   60 ALFATOT(I)=ABSCOEF(I)+ALFATOT(I)                                  ADD00092
C                                                                       ADD00093
C     HEXADECAPOLE-INDUCED TERM                                         ADD00094
C     =============================                                     ADD00095
C                                                                       ADD00096
      AMUL0=PHI                                                         ADD00097
      DELTA=DELPHI                                                      ADD00098
      S=Y(X,0.12504D-53/9.,-7.59574D0,0.66118D0)                            
      E=Y(X,-0.5D0,0.0D0,0.0D0)                                               
      T1=Y(X,0.40711D-11,-1.20681D0,0.099D0)                                
      T2=Y(X,0.1743D-12,0.37304D0,-0.08709D0)                               
      T3=Y(X,0.42016D-12,0.09897D0,-0.05468D0)                              
      T4=Y(X,0.15957D-7,-5.10481D0,0.55655D0)                               
      CALL ADSPEC1 (S,E,T1,T2,T3,T4,TEMP,NF,FREQ,ABSCOEF,0,LIKE,0,4,4,5)ADD00105
      DO 70 I=1,NF                                                      ADD00106
   70 ALFATOT(I)=ABSCOEF(I)+ALFATOT(I)                                  ADD00107
      write(6,110) FREQ(1),FREQ(NF),FREQ(2)-FREQ(1),TEMP                  
      write(6,120) (ALFATOT(I),I=1,NF)                                    
C                                                                       ADD00110
   80 FORMAT ( 23H1ABSORPTION SPECTRA OF ,A10,  3H AT,F8.1,  2H K,/1X,43ADD00111
     1(1H=),/, 11H MIN.FREQ.=,F8.1,  5H CM-1,10X, 10HMAX.FREQ.=,F8.1,  5ADD00112
     2H CM-1,10X, 15HFREQ.INCREMENT=,F8.2,  5H CM-1,5X,  2HIN,I5,  6H STADD00113
     3EPS,//)                                                           ADD00114
   90 FORMAT (F8.5)                                                     ADD00115
  100 FORMAT (10F8.5)                                                   ADD00116
  110 FORMAT (/,' ABSORPTION COEFFICIENT ALPHA(FNU), FROM', F10.2,      ADD00117
     1' CM-1 TO', F10.2,' CM-1',/,' AT', F6.2,' CM-1 INCREMENTS',
     1 ' TEMPERATURE =',F7.2, ' K,  IN UNITS OF [CM-1 AMAGAT-2]',/)       
  120 FORMAT (  101 (6E13.5,/))                                         ADD00120
      END                                                               ADD00122

      SUBROUTINE ADSPEC1 (G0,EP,TAU1,TAU2,TAU5,TAU6,TEMP,NF,FREQ,ABSCOEFADS00123
     1,MP,LIKE,LAMBDA1,LAMBDA2,LAMBDA,LVALUE)                           ADS00124
C                                                                       ADS00125
C     FOR INDUCTION BY  A  T E T R A H E D R A L M O L E C U L E        ADS00126
C     THIS PROGRAM GENERATES LISTING OF CIA TR ALFA(OMEGA)              ADS00127
C                                                                       ADS00128
      implicit double precision (a-h,o-z)
      COMMON /CHPART/ Q1,WCH(2),B01,D01,JRANGE2                         ADS00130
      COMMON /STAT/ QQW3(51,3),QQW4(51,4)                               ADS00131
      COMMON /ROT/ AMUL0,DELTA                                          ADS00132
      DIMENSION ABSCOEF(NF), FREQ(NF)                                   ADS00133
      DATA CLOSCHM,BOLTZWN/2.68675484D19,.6950304/                      ADS00134
      DATA HBAR,PI,CLIGHT/1.054588757D-27,3.1415926535898,2.9979250D10/ ADS00135
      ECH4(I)=B01*FLOAT(I)                                              ADS00136
      PCH4(J,T)=DEXP(-1.4387859/T*ECH4(J*(J+1)))                        ADS00137
      TWOPIC=2.*PI*CLIGHT                                               ADS00138
      CALIB=TWOPIC*((4.*PI**2)/(3.*HBAR*CLIGHT))*CLOSCHM**2             ADS00139
      CALIB=CALIB/FLOAT(1+LIKE)                                         ADS00140
      BETA=1./(BOLTZWN*TEMP)                                            ADS00141
      LIST=NF                                                           ADS00142
      DO 10 I=1,LIST                                                    ADS00143
   10 ABSCOEF(I)=0.0                                                    ADS00144
C                                                                       ADS00145
      WRITE (6,70) LAMBDA1,LAMBDA2,LAMBDA,LVALUE,G0,EP,TAU1,TAU2,
     1  TAU5,TAU6, G0*BGAMA(0.,TAU1,TAU2,EP,TAU5,TAU6,TEMP)              
C                                                                       ADS00148
      DO 60 I=1,JRANGE2                                                 ADS00149
         J=I-1                                                          ADS00150
         P=FLOAT(2*J+1)*PCH4(J,TEMP)/Q1                                 ADS00151
         AMUL2=(AMUL0+DELTA*FLOAT(J*(J+1)))**2                          ADS00152
         DO 40 IP=1,LAMBDA                                              ADS00153
C                                                                       ADS00154
C         POSITIVE DELTA J                                              ADS00155
C                                                                       ADS00156
            JP=J+IP                                                     ADS00157
            OMEGA1=ECH4(JP*(JP+1))-ECH4(J*I)                            ADS00158
            IF (LAMBDA.EQ.3) CC=(1.+QQW3(I,IP)/4.)                      ADS00159
            IF (LAMBDA.EQ.4) CC=(1.+QQW4(I,IP)/4.)                      ADS00160
            FAC=CALIB*P*CC                                              ADS00161
            DO 20 IQ=1,LIST                                             ADS00162
               FRQ=FREQ(IQ)-OMEGA1                                      ADS00163
               WKF=FREQ(IQ)*(1.-DEXP(-BETA*FREQ(IQ)))*FAC*
     1         (2.*DFLOAT(JP)+1.) 
               XBG=G0*BGAMA(FRQ,TAU1,TAU2,EP,TAU5,TAU6,TEMP)*AMUL2      ADS00166
               ABSCOEF(IQ)=ABSCOEF(IQ)+XBG*WKF                          ADS00167
   20       CONTINUE                                                    ADS00168
C                                                                       ADS00169
C     NEGATIVE DELTA J                                                  ADS00170
C                                                                       ADS00171
            JP=J-IP                                                     ADS00172
            IF (JP.LT.0) GO TO 40                                       ADS00173
            OMEGA1=ECH4(JP*(JP+1))-ECH4(J*I)                            ADS00174
            FAC=CALIB*P                                                 ADS00175
            DO 30 IQ=1,LIST                                             ADS00176
               FRQ=FREQ(IQ)-OMEGA1                                      ADS00177
               WKF=FREQ(IQ)*(1.-DEXP(-BETA*FREQ(IQ)))*FAC*
     1         (2.*DFLOAT(JP)+1.)                                       ADS00179
               XBG=G0*BGAMA(FRQ,TAU1,TAU2,EP,TAU5,TAU6,TEMP)*AMUL2      ADS00180
               ABSCOEF(IQ)=ABSCOEF(IQ)+XBG*WKF                          ADS00181
   30       CONTINUE                                                    ADS00182
   40    CONTINUE                                                       ADS00183
C                                                                       ADS00184
C     DELTA J=0                                                         ADS00185
C                                                                       ADS00186
         JP=J                                                           ADS00187
         FAC=CALIB*P                                                    ADS00188
         DO 50 IQ=1,LIST                                                ADS00189
            FRQ=FREQ(IQ)                                                ADS00190
            WKF=FREQ(IQ)*(1.-DEXP(-BETA*FREQ(IQ)))*FAC*
     1      (2.*FLOAT(JP)+1.)                                           ADS00191
            XBG=G0*BGAMA(FRQ,TAU1,TAU2,EP,TAU5,TAU6,TEMP)*AMUL2         ADS00192
            ABSCOEF(IQ)=ABSCOEF(IQ)+XBG*WKF                             ADS00193
   50    CONTINUE                                                       ADS00194
   60 CONTINUE                                                          ADS00195
      write(6,80) (ABSCOEF(I),I=1,LIST)                                 ADS00196
      RETURN                                                            ADS00197
C                                                                       ADS00198
   70 FORMAT (/, 32H LAMBDA1,LAMBDA2, LAMBDA,LVALUE=,2I3,2X,2I3, 12H COMADS00199
     1PONENT .,/15X, 22HLINE SHAPE PARAMETERS:,6E12.3,5X,  5HG(0)=,E12.3ADS00200
     2/)                                                                ADS00201
   80 FORMAT ((1X,10E12.4,/))                                           ADS00202
      END                                                               ADS00204

      FUNCTION BGAMA (FNU,T1,T2,EPS,T3,T4,TEMP)                         BGA00205
C                                                                       BGA00206
	implicit double precision (a-h,o-z)
      REAL*8 K0                                                         BGA00207
      DATA PI,CLIGHT/3.1415926535898,2.99792458D10/                     BGA00208
      DATA HBAR,BOLTZ/1.0545887D-27,1.380662D-16/                       BGA00209
      P1(X)=((((((.0045813*X+.0360768)*X+.2659732)*X+1.2067492)*X+3.0899BGA00210
     1424)*X+3.5156229)*X+1.)                                           BGA00211
      P2(X)=((((((.00000740*X+.00010750)*X+.00262698)*X+.03488590)*X+.23BGA00212
     1069756)*X+.42278420)*X-.57721566)                                 BGA00213
      P3(X)=((((((.00032411*X+.00301532)*X+.02658733)*X+.15084934)*X+.51BGA00214
     1498869)*X+.87890594)*X+.5)                                        BGA00215
      P4(X)=((((((-.00004686*X-.00110404)*X-.01919402)*X-.18156897)*X-.6BGA00216
     17278579)*X+.15443144)*X+1.)                                       BGA00217
      P5(X)=((((((.00053208*X-.00251540)*X+.00587872)*X-.01062446)*X+.02BGA00218
     1189568)*X-.07832358)*X+1.25331414)                                BGA00219
      P6(X)=((((((-.00068245*X+.00325614)*X-.00780353)*X+.01504268)*X-.0BGA00220
     13655620)*X+.23498619)*X+1.25331414)                               BGA00221
C                                                                       BGA00222
      OMEGA=2.*PI*CLIGHT*FNU                                            BGA00223
      T0=HBAR/(2.*BOLTZ*TEMP)                                           BGA00224
      Z=DSQRT((1.+(OMEGA*T1)**2)*(T2*T2+T0*T0))/T1                      BGA00225
      IF (Z-2.) 10,10,20                                                BGA00226
   10 XK1=Z*Z*DLOG(Z/2.)*P3((Z/3.75)**2)+P4((Z/2.)**2)                  BGA00227
      GO TO 30                                                          BGA00228
   20 XK1=DSQRT(Z)*DEXP(-Z)*P6(2./Z)                                    BGA00229
   30 ZP=DSQRT((1.+(OMEGA*T4)**2)*(T3*T3+T0*T0))/T4                     BGA00230
      IF (ZP-2.) 40,40,50                                               BGA00231
   40 K0=-DLOG(ZP/2.)*P1((ZP/3.75)**2)+P2((ZP/2.)**2)                   BGA00232
      GO TO 60                                                          BGA00233
   50 K0=DEXP(-ZP)*P5(2./ZP)/DSQRT(ZP)                                  BGA00234
   60 BGAMA=((T1/PI)*DEXP(T2/T1+T0*OMEGA)*XK1/
     1  (1.+(T1*OMEGA)**2)+EPS*(T3/
     1PI)*DEXP(T3/T4+T0*OMEGA)*K0)/(1.+EPS)                             BGA00236
      RETURN                                                            BGA00237
      END                                                               BGA00239

      SUBROUTINE PARTSUM (TEMP)                                         PAR00240
	implicit double precision (a-h,o-z)
      COMMON /CHPART/ Q1,WCH(2),B01,D01,JRANGE2                         PAR00241
      DATA B01,D01,WCH(1),WCH(2)/5.24,0.,1.,1./                         PAR00242
      ECH4(I)=B01*FLOAT(I)                                              PAR00243
      PCH4(J,T)=DEXP(-1.4387859/T*ECH4(J*(J+1)))                        PAR00244
C                                                                       PAR00245
C     Q1,B01,D01,WCH - PARTITION FCT., ROT.CONSTANTS, WEIGHTS FOR CH4   PAR00246
C     *** PARTITION FUNCTION FOR CH4 *******************                PAR00247
C                                                                       PAR00248
      Q1=0.                                                             PAR00249
      J=0                                                               PAR00250
   10 DQ=PCH4(J,TEMP)*FLOAT((2*J+1)**2)                                 PAR00251
      Q1=Q1+DQ                                                          PAR00252
      J=J+1                                                             PAR00253
      IF (DQ.GT.Q1/1000.) GO TO 10                                      PAR00254
      JRANGE2=J+3                                                       PAR00255
      RETURN                                                            PAR00256
      END                                                               PAR00258
