PROGRAM H2H2model03_LT
c
C	This model allows for choice:
c	with, or without J-dependence,
c	and choice of normal and equilibrium H2

C       This program reads file 'par.beta'
C
C =========================================================================
C       Copyright (C) 1998 by Yi Fu and A. Borysow
C =========================================================================
C Copyright Notice:
C You may use this program for your scientific applications,
C but please do not distribute it yourself.
C Direct your inquires to: aborysow@nbi.dk
C Final request: If you publish  your work which benefited from this
C program, please acknowledge using this program 
C by referring to:
C       Y. Fu, C. Zheng and A. Borysow, JQSRT, 1999, (submitted), where original
c	quantum mechanical computations equivalent to this model have been described. 
c	The detailed description of this model is  available only on WWW (unpublished model)

c	For more information, refer also to the paper by
c	C. Brodbeck and J.-P. Bouanich and  Nguyen-Van-Thanh  and
c       Y. Fu and A. Borysow, J. Chem. Phys., vol. 110, p. 4750, 1999.

c	This model is IDENTICAL with the theoretical results presented in
c	Fu et al. 1999 and Brodbeck et al. 1999

C ============================================================================
C
        IMPLICIT double precision (A-H,O-Z)
        CHARACTER*10  MARKER
        dimension ibgama(6),LA1(6),LA2(6),LA(6),LL(6)
	dimension NVIB(10, 4), w(2),  Npar(6)
        Dimension BM0(12,4), BTAU1(12,4), BTAU2(12,4)
 	dimension alfatot(0:12, 600)
	common/ integ / rlo, rup, frac, Ninteg
	common/temp/temp
        COMMON /PRTSUM/ q1, w, NTYPE, JRANGE1, MARKER
        COMMON /SPECIT/ LIST
	common/ specit1/ FREQ(600), ABSCOEF(600)
	common/ read/ go
	common/om/omega_0
	COMMON/VIBIF/ NVIB1, NVIB1p, NVIB2, NVIB2p
	COMMON /bpar/ Bpar(9, 4), NB
	COMMON /Jdep / JDEP

	DATA  LIKE/ 1 /
 	DATA W(1),  W(2)/1.0, 3.0/   ! rotational weights;
	data nf/ 600 /
	data NVIB/0, 0, 0, 1, 1, 1, 2, 2, 2, 3,       ! == V1
     &            0, 1, 1, 1, 2, 2, 2, 3, 3, 3,       ! == V1'
     &            0, 0, 1, 0, 0, 1, 0, 0, 1, 0,       ! == V2
     &            3, 2, 3, 3, 2, 3, 3, 2, 3, 3/       ! == V2'

C       IBGAMA = 0 uses K0, IBGAMA = 1 uses BC, BGAMA = 2 uses K1
        data ibgama/0,1,1,2,0,1/    ! IBGAMA = 0 uses K0, IBGAMA = 1 uses BC
C       data ibgama/0,0,1,2,0,0,1,1,2,0/
        data LA1   /0,2,0,2,0,2/
        data LA2   /0,0,2,0,2,2/
        data LA    /0,2,2,2,2,3/
        data LL    /1,3,3,1,1,3/

         data Npar  /7,4,4,7,7,4/
      DATA BM0/-64.6632, -64.7965, -64.5628, -70.1186, -64.5327,0.d0,
     &        0.0d0,-64.5732,-64.5689,-65.8634,-65.1641,-65.7947, !B0
     &      -3.79051, -3.18043, -2.90725, 1.67553, -3.42757, 0.d0,
     &      0.0d0,-2.93532,-2.97008,-3.26495,-4.97074,-2.84093,     !B1
     &      1.69722, 1.22891, 1.1752,  -1.51323,  1.39377, 0.d0,
     &      0.0d0, 1.13392, 1.15348, 1.30463, 2.15157,1.09193,     !B2
     &    -0.237117, -0.120397, -0.156739,  0.37531, -0.156109, 0.d0,
     &     0.0d0,-0.123712,-0.127921,-0.145022,-0.244661,-0.122333/ !B3

      DATA BTAU1/-11.8974, -12.1305, -10.57, -14.4108, -12.3188,0.d0,
     &          0.0d0,-11.9119,-11.813,-12.0872,-13.5173,-11.7945,  !B0
     &          -1.19800, -0.469693, -2.77937, 3.35085,-0.386058,0.d0,
     &        0.0d0,-0.595612,-0.741825,-0.251777,1.39961,-0.712328, !B1
     &        0.604604, -0.0484486, 1.12265, -2.03575, -0.0590297,0.d0,
     &        0.0d0,0.0174677,0.115793,-0.140516,-1.01606,0.0745054, !B2
     &      -0.140778, 0.00534555, -0.16634, 0.327846, 0.0152675,0.d0,
     &     0.0d0,-0.00576595,-0.0310663,0.0232436,0.183536,-0.0117238/ !B3

      DATA BTAU2/-19.6454,-12.6634,-12.6675,0.0d0,-7.49955,0.d0,
     &           0.0d0,-12.5616,-12.6915,0.0d0,-9.8152,-12.5353,  !B0
     &           12.8185, -0.1798,0.297043,  0.0d0, -4.75913,0.d0,
     &           0.d0,-0.179975,0.0680707,0.d0,-4.75287,-0.159049,  !B1
     &           -7.74865,-0.134118,-0.267489,0.d0,1.50426,0.d0,
     &           0.d0,-0.129425,-0.291291,0.d0,2.21046,-0.129466,    !B2
     &           1.47426,0.0170518,0.00264082, 0.d0,-0.201579,0.d0,
     &           0.d0,0.0152609,0.0527647,0.d0,-0.376787,0.0104684/  !B3


       open(unit=10,file='data',form='formatted', status='unknown')      !ouput
       open(unit=8, file='par.beta', form='formatted',status='old')   
c						        !input par for Beta_(R) | Bpar

	write (*,*) 
     1 ' This is a model for the RV CIA: 2nd overtone band'
	write (*,*) 
     2 ' You can choose to run calculations WITH (JDEP=1) '
	write (*,*) 
     3 ' or WITHOUT (JDEP=0)  <<J-dependence>> (see paper)'
	write (*,*) 
     4 ' Refer to Borysow and Fu, Icarus, 1998, submitted'
	write (*,*)  ' JDEP=?'
	read (*,*) JDEP

	write(*,*)  'Temperature? (min. 20., max. 500.) '
	read (*,*)   temp
	if(temp.gt.500.D0) 
     1	write(*,*) 
     2 'Temperature beyond allowed upper limit: proceed on your own!!!' 

	write(*,*) ' Normal (NTYPE=1) or Equilibrium (Ntype=0) H2?'
	write(*,*) ' NTYPE=?  '
	read (*,*) NTYPE
	if (NTYPE .eq. 1) write (*,*) ' Normal H2 assumed'
	if (NTYPE .eq. 0) write (*,*) ' Equilibrium H2 assumed'
	if ((NTYPE .eq.1) .or. (NTYPE.eq.0))  then  ! OK then
	else
	stop 1001
	endif

	write(*,*) 'Input minimum, maximum, frequency step (in 1/cm):'
	read(*,*) FREQLO,FREQMAX,DELTASTEP
	list = INT (FREQMAX - FREQLO)/DELTASTEP
	write(*,*) 'You will find the results in file: data' 

      DO 10 I=1,list
10    FREQ(I)=FREQLO+DFLOAT(I-1)*DELTASTEP

	if(JDEP.eq.0) write (10,*) ' No J-dependence assumed'
	if(JDEP.eq.1) write (10,*) ' J-dependence assumed'
	write(10,*) ' Refer to Borysow & Fu, Icarus, 1998'  
	if (NTYPE .eq. 1) write (10,*) '# Normal H2 assumed'
	if (NTYPE .eq. 0) write (10,*) '# Equilibrium H2 assumed'
	write(10 , 1212) temp, freqlo, freqmax, deltastep
1212	format('##   2nd Overtone Band RV CIA SPECTRA of H2-H2 ',/,
     &         '##         Temperature = ', f9.2, ' K', /,
     &         '##   Frequency Min., Max., Step:',3F9.1, ' cm-1',/
     &         '##   freq(cm-1)   alfa*1E+8 (cm-1amagat-2)',/)

C +++++++++++++	Compute the partition function Q1: +++++++++++++++++++

      CALL PARTFCT (TEMP, Q1, NTYPE, JRANGE1, MARKER, W )

        ALFA =1./(0.69519*TEMP)
        hbar = 1.054588757D-27
	boltzk = 1.38054d-16
	tau0 = hbar/(2.*boltzk*temp)
C *****************************************************************

        TLG = DLOG10(TEMP)
	do 11 N_VIB = 1, 2           ! LOOP for Vibrational Transitions
	NVIB1  =  NVIB(N_VIB, 1)
	NVIB1p =  NVIB(N_VIB, 2)
	NVIB2  =  NVIB(N_VIB, 3)
	NVIB2p =  NVIB(N_VIB, 4)

	do 11 N_term = 1, 6           ! LOOP for La1La2LaL terms
	NB = Npar(N_term)

	   read(8,*)
         do N_B = 1, 4                 ! LOOP for reading B_par of J_Dependance
        read(8,*) (Bpar(i, N_B), i = 1, 9)
        end do

       if ( Bpar(1, 1) .eq. 0.0 ) goto 11
       IF(N_VIB.EQ.2 .AND. N_term.EQ.1) GOTO 11
          GM0 = 0.0d0
          GTau_1 = 0.0d0
          GTau_2 = 0.0d0
        Np_term = N_term+ (N_VIB-1)*6
          DO J=1,4
          GM0 = GM0 + BM0(Np_term,J)* TLG **(J-1)
          GTau_1 = GTau_1 + BTAU1(Np_term,J)*TLG **(J-1)
          GTau_2 = GTau_2 + BTAU2(Np_term,J)*TLG **(J-1)
          ENDDO

          XM0 = 10**GM0
          Tau1 = 10.**GTau_1
          Tau2 = 10.**GTau_2
        

c	* * * * ** * * * * * * * * * * ** * * * * * * * * * * ** * * * * *

        CALL ADDSPEC(TEMP,LIKE, LA1(N_term), LA2(N_term),
     &                 LA(N_term), LL(N_term), freqcut,
     &                 ibgama(N_term),XM0,tau1,tau2)

	do 198 l=1,list
	alfatot(0, l) = alfatot(0, l) + ABSCOEF(l) * 2. ! total
       alfatot((N_VIB-1)*6. + N_term, l) = ABSCOEF(l) * 2. ! separate terms
198     continue

11     continue
       write (10,9997) (freq(l), alfatot(0,l)*1.D8, l=1, list)
9997	 format(3x,f10.2, 2x, e12.4)
       END
       
      SUBROUTINE ADDSPEC(TEMP, LIKE,LAMBDA1,LAMBDA2,LAMBDA,LVALUE,
     &                      CUT, ibgama,  G0, tau1,tau2)
C

C     This subroutine generates a detailed listing of CIA ALFA(OMEGA)
c	IBGAMA=0 means uses K0, IBGAMA=1 (uses BC), IBGAMA=2 (uses K1)
C     LIKE=1 for like systems (as H2-H2) 
C
      IMPLICIT double precision (A-H,O-Z)
      CHARACTER*10 MESSAGE, MARKER
      COMMON /PRTSUM/ q1, w1(2), NTYPE, JRANGE1, MARKER
      COMMON /SPECIT/ LIST
	common/specit1/ FREQ(600),ABSCOEF(600)
	common/ integ / rlo, rup, frac, ninteg
      COMMON/VIBIF/ NVIB1, NVIB1p, NVIB2, NVIB2p
      COMMON/ROTIF/ J1, JP1, J2, JP2
	common /Jdep / JDEP
	data nf /600 /

      DATA  SCALEF/1.D80/,  BOLTZWN/ 0.6950304D0/
      DATA HBAR,PI,CLIGHT/1.054588757D-27, 3.141592653589797D0,
     1                    2.9979250D10/
      DATA Boltzk / 1.380658d-16 /

      DO 10 I=1, nf
   10 ABSCOEF(I)=0.D0
      TWOPIC=2.*PI*CLIGHT
      MESSAGE='NON-COMUL.'
      CALIB=TWOPIC*((4.*PI**2)/(3.*HBAR*CLIGHT))*(2.68676D19)**2
      CALIB=CALIB/DFLOAT(1+LIKE)
      FREQCUT = CUT		! do not alter the outside parameter CUT
      BETA1 = 1.D0/(BOLTZWN*TEMP)

      JPLUSL=JRANGE1+MAX0(LAMBDA1,LAMBDA2)

      DO 160 I1=1,JRANGE1
         J1=I1-1
      DO 160 IP1=1,JPLUSL
         JP1=IP1-1
         CG1S=CLEBSQR(J1,LAMBDA1,JP1)
         IF (CG1S) 160,160,110
  110    P1=H2POPL(NVIB1,J1,TEMP,W1)/Q1
         OMEGA1=H2ELEV(NVIB1p,JP1)-H2ELEV(NVIB1,J1)

      DO 150 I2=1,JRANGE1
         J2=I2-1
      DO 150 IP2=1,JPLUSL
         JP2=IP2-1
         CG2S=CLEBSQR(J2,LAMBDA2,JP2)
         IF (CG2S) 150,150,120
  120    P2=H2POPL(NVIB2,J2,TEMP,W1)/Q1
         OMEGA2=H2ELEV(NVIB2p,JP2)-H2ELEV(NVIB2,J2)

        FAC=CALIB*P1*P2*CG1S*CG2S   
	  if(JDEP.eq.1) FAC=FAC*xscale() 
	  if(JDEP.eq.0) FAC=FAC
	
      DO 140 I=1,LIST
      FRQ=FREQ(I)-OMEGA1-OMEGA2


  130   WKF=FREQ(I) * (1.-DEXP(- BETA1 * FREQ(I))) * FAC

	if (ibgama.eq.0) xx = g0 * bgama0(frq,tau1,tau2,temp)
	if (ibgama.eq.1) xx = g0 * bgama1(frq,tau1,tau2,temp)
      if (ibgama.eq.2) xx = g0 * bgama_H(frq,tau1,temp)
        ABSCOEF(I) = ABSCOEF(I) + xx * WKF
  140   CONTINUE

  150   CONTINUE
  160   CONTINUE

      RETURN
      END

      SUBROUTINE PARTFCT (TEMP, Q, NTYPE, JRANGE, MARKER, W)
C
C     NTYPE = 0 for equilibrium hydrogen 
C     NTYPE = 1 redefines the weights W for "normal" hydrogen with a
C               ratio of ortho:para H2 of 3:1. 
C
      IMPLICIT double precision (A-H,O-Z)
      CHARACTER*10  MARKER
      DIMENSION W(2)
      DATA  BOLTZWN /.6950304D0/

      Q = 0.0
      J = -1

      IF (NTYPE.EQ.1) GO TO 100
40    J=J+1
      DDD = H2POPL(0,J,TEMP,W) + h2popl(1,j,temp,w) + 
     1      h2popl(2,j,temp,w) + h2popl(3,j,temp,w)

      Q = Q + DDD
      IF (DDD.GE.0.002*Q) GO TO 40
      JRANGE=J+2
      MARKER='EQUILIBRM '
      GO TO 200

C     "NORMAL" HYDROGEN REQUIRES REDEFINITION OF W(2):

100   continue
      Qeven = 0.0
70    J=J+1
      DDD = H2POPL(0,J,TEMP,W) + h2popl(1,j,temp,w) + 
     1      h2popl(2,j,temp,w) + h2popl(3,j,temp,w)

      Q = Q + DDD
      IF (MOD(J,2).eq.0) Qeven = Qeven + DDD
      IF (DDD.GE.0.002*Q) GO TO 70
      JRANGE=J+2
      MARKER='NORMAL    '

C
C     DEFINITION OF "NORMAL" HYDROGEN: 3*S(EVEN) = S(ODD)
C
      W(2) = W(2) * (3.0 * Qeven) / (Q - Qeven)  !normaL
      Q = 4.0 * Qeven
 
C     W(2) = W(2) * (0.7544 * Qeven) / (Q - Qeven) !c_para=.57
C     Q = 1.7544 * Qeven

200   continue
      RETURN
      END


      FUNCTION H2POPL (NVIB, J, T, W)
C     UNNORMALIZED (!) Boltzmann factor of the roto- vibrational levels
C          of the H2 molecule.  (To be normalized by Q = partition sum)
C
      IMPLICIT double precision (A-H,O-Z)
      DIMENSION W(2)
      DATA  BOLTZWN /.6950304D0/
C
20	DDD = H2ELEV(nvib,J)
        H2POPL = DFLOAT(2*J+1)*W(1+MOD(J,2))*DEXP(-DDD/(BOLTZWN*T))
      RETURN
      END

      FUNCTION H2ELEV (NVIB, JROT)
C     Roto- vibrational energy levels of the H2 molecule
C     Needed for the computation of the partition sum Q of H2 ******
C     This function is for low temperatures.
C
C     All D, B, H and wavenumber are from the reference By 
C     Uwe Fink T.A. Wiggins and D.H. Rank, JMS 18, 384-395 (1965)
C     Checked by Yi Fu in Jan 29, 1998
 
      IMPLICIT double precision (A-H,O-Z)
 
      EH2(I,B, D, H)=B*DFLOAT(I) -D*DFLOAT(I)**2 + H *DFLOAT(I)**3
 
        if(nvib.eq.0) h2elev=eh2(jrot*(jrot+1),59.3362D0, 
     &           4.583D-2, 4.87D-5)
 
        if(nvib.eq.1) h2elev=eh2(jrot*(jrot+1),56.3722D0, 
     &         4.3880D-2, 3.97D-5) + 4161.1782d0
C       4161.1782cm^-1 is the wavenumber of the vibrational
C       transition of v= 0-> 1

        if(nvib.eq.2) h2elev=eh2(jrot*(jrot+1),53.4823D0,
     &        4.272D-2,3.97D-5) + 8087.d0
C       8087.000d0 cm^-1 is the wavenumber of the vibrational
C       transition of v= 0 -> 2

        if(nvib.eq.3) h2elev=eh2(jrot*(jrot+1),50.6245D0,
     &       4.08d-2, 3.68d-5) + 11782.355d0
C     11782.355d0cm^-1 is the wavenumber of the vibrational
C     transition of v= 0-> 3

        if(nvib.gt.3) stop 888
C
      RETURN
      END


      FUNCTION XK0(X)
C     MODIFIED BESSEL FUNCTION K0(X)
C     ABRAMOWITZ AND STEGUN P.379
	implicit double precision (a-h,o-z)
      IF(X-2.) 10,10,20
   10 T=(X/3.75d0)**2
      FI0=(((((.0045813*T+.0360768)*T+.2659732)*T
     1 +1.2067492)*T+3.0899424)*T+3.5156229)*T+1.
      T=(X/2.)**2
      P=(((((.00000740*T+.00010750)*T+.00262698)*T
     1 +.03488590)*T+.23069756)*T+.42278420)*T+
     2 (-.57721566)
      X=dABS(X)
      XK0=-dLOG(X/2.)*FI0+P
      RETURN
   20 T=(2./X)
      P=(((((.00053208*T-.00251540)*T+.00587872)*T
     1 -.01062446)*T+.02189568)*T-.07832358)*T+
     2 1.25331414
      X=dMIN1(X,330.d0)
      XK0=dEXP(-X)*P/dSQRT(X)
      RETURN
      END

      FUNCTION XK1(X)
C     MODIFIED BESSEL FUNCTION K1(X) TIMES X
C     PRECISION IS BETTER THAN 2.2E-7 EVERYWHERE.
C     ABRAMOWITZ AND S,TEGUN, P.379; TABLES P.417.
	implicit double precision (a-h,o-z)
      IF(X-2.) 10,10,20
   10 T=(X/3.75)**2
      FI1=X*((((((.00032411*T+.00301532)*T+.02658733)*T+.15084934)
     1 *T+.51498869)*T+.87890594)*T+.5)
      T=(X/2.)**2
      P=(((((-.00004686*T-.00110404)*T-.01919402)*T-.18156897)*T-
     1 .67278579)*T+.15443144)*T+1.
      XK1=X*dLOG(X/2)*FI1+P
      RETURN
   20  T=2./X
      P=(((((-.00068245*T+.00325614)*T-.00780353)*T+.01504268)*T-
     1 .03655620)*T+.23498619)*T+1.25331414
      X=dMIN1(X,330.d0)
      XK1=dSQRT(X)*dEXP(-X)*P
      RETURN
      END


      FUNCTION BGAMA0(FNU,TAU1,TAU2,TEMP)
C     K0 LINE SHAPE MODEL
C     NORMALIZATION SUCH THAT 0TH MOMENT EQUALS UNITY
C     FNU IS THE FREQUENCY IN CM-1; TEMP IN KELVIN.
	implicit double precision (a-h,o-z)
      DATA TWOPIC,BKW,HBOK,PI/1.883651568d11,.6950304256d0,
     1 7.638280918d-12, 3.141592654d0/
	TAU0 = HBOK /(2.*TEMP)
      TAU4=DSQRT(TAU1*TAU1+(TAU0)**2)
      OMEGA=TWOPIC*FNU
      XNOM=1.d0/(TAU2*TAU2)+OMEGA*OMEGA
      X=TAU4*dSQRT(XNOM)
      TAU12=TAU1/TAU2
      TAU12=DMIN1(TAU12,430.d0)
      BGAMA0=(TAU1/PI)*DEXP(TAU12+FNU/(2.d0*BKW*TEMP))*
     1 XK0(X)
	end


      FUNCTION CLEBSQR (L,LAMBDA,LP)
C     SQUARE OF CLEBSCH-GORDAN COEFFICIENT (L,LAMBDA,0,0;LP,0)
C     FOR INTEGER ARGUMENTS ONLY
C     NOTE THAT LAMBDA SHOULD BE SMALL, MAYBE @10 OR SO.
C
      IMPLICIT REAL*8 (A-H,O-Z)
      FC=DFLOAT(2*LP+1)
      GO TO 10
C
      ENTRY THREEJ2
C
C     THIS ENTRY RETURNS THE SQUARED 3-J SYMBOL   L LAMBDA LP
C                                                 0    0    0
C     INSTEAD OF THE CLEBSCH-GORDAN COEFFICIENT
C     (LIMITATION TO INTEGER ARGUMENTS ONLY)
C
C     NOTE THAT THE THREE-J SYMBOLS ARE COMPLETELY SYMMETRIC IN THE
C     ARGUMENTS. IT WOULD BE ADVANTAGEOUS TO REORDER THE INPUT ARGUMENT
C     LIST SO THAT LAMBDA BECOMES THE SMALLEST OF THE 3 ARGUMENTS.
C
      FC=1.
   10 CLEBSQR=0.
      IF (((L+LAMBDA).LT.LP).OR.((LAMBDA+LP).LT.L).OR.((L+LP).LT.LAMBDA)
     1) RETURN
      IF (MOD(L+LP+LAMBDA,2).NE.0) RETURN
      IF ((L.LT.0).OR.(LP.LT.0).OR.(LAMBDA.LT.0)) RETURN
      F=1./DFLOAT(L+LP+1-LAMBDA)
      IF (LAMBDA.EQ.0) GO TO 30 
      I1=(L+LP+LAMBDA)/2        
      I0=(L+LP-LAMBDA)/2+1
      DO 20 I=I0,I1       
   20 F=F*DFLOAT(I)/DFLOAT(2*(2*I+1))
   30 P=FC*F*FCTL(LAMBDA+L-LP)*FCTL(LAMBDA+LP-L)
      CLEBSQR=P/(FCTL((LAMBDA+L-LP)/2)*FCTL((LAMBDA+LP-L)/2))**2
      RETURN
      END

      FUNCTION FCTL (N)
C     Simple FACTORIALS program
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
      P(Z)=((((-2.294720936D-4)/Z-(2.681327160D-3))/Z+(3.472222222D-3))/
     1Z+(8.3333333333D-2))/Z+1.
      FCTL=1            
      IF (N.LE.1) RETURN
      IF (N.GT.15) GO TO 20
      J=1           
      DO 10 I=2,N   
   10 J=J*I         
      FCTL=DFLOAT(J)
      RETURN        
   20 Z=DFLOAT(N+1) 
      FCTL=(DEXP(-Z)*(Z**(Z-0.5))*P(Z)*2.5066282746310)
      RETURN
      END

      FUNCTION BGAMA1(FNU,TAU1,TAU2,TEMP)
C     BIRNBAUM's CIA LINE SHAPE MODEL (K1)
C     NORMALIZATION SUCH THAT 0TH MOMENT EQUALS UNITY
C     FNU IS THE FREQUENCY IN CM-1; TEMP IN KELVIN.
      implicit double precision (a-h,o-z)
      DATA TWOPIC,BKW,HBOK,PI/1.883651568d11,.6950304256d0,
     1 7.638280918d-12, 3.141592654d0/

      TAU3=dSQRT(TAU2*TAU2+(HBOK/(TEMP*2.))**2)
      OMEGA=TWOPIC*FNU
      DENOM=1.d0+(OMEGA*TAU1)**2
      X=(TAU3/TAU1)*DSQRT(DENOM)
      AAA=TAU2/TAU1
      AAA=dMIN1(AAA,430.d0)
      BGAMA1=(TAU1/PI)*dEXP(AAA+FNU/(2.d0*BKW*TEMP))*
     1 XK1(X)/DENOM
      RETURN
      END

      FUNCTION BGAMA_H(FNU,tau ,TEMP)
c	"K1" one parameter model lineshape (14-Jul-1989 )
c	with Egelstaff time:
c	Authorship: J. Borysow & A. Borysow (1989)
C     NORMALIZATION SUCH THAT 0TH MOMENT EQUALS UNITY
C     FNU IS THE FREQUENCY IN CM-1; TEMP IN KELVIN.
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DATA TWOPIC,BKW,HBOK,PI/1.883651568E11,.6950304256,
     1 7.638280918D-12, 3.141592654/

      TAU3=DSQRT( 3./2.*TAU**2 + (HBOK/(TEMP*2.))**2 )
      OMEGA=TWOPIC*FNU
      ARG = DABS(OMEGA)*TAU3
	IF(ARG.EQ.0.D0) GO TO 5
      BGAMA_H =1./PI * DEXP(FNU/(2.*BKW*TEMP))* XK1(ARG)
     1      * (TAU/TAU3)**2 * TAU * (3./2.)**(3./2.)
	RETURN
5	BGAMA_H =1./PI * DEXP(FNU/(2.*BKW*TEMP))
     1      * (TAU/TAU3)**2 * TAU * (3./2.)**(3./2.)
C	LIMIT K1(X) (X-->0) = 1/X; SO 
C	LIMIT X*K1(X) (X-->0) = X/X = 1. EXACT
c	G1 = (beta *hbar)/(tau**2)
c	G2 = (5./3.)*beta*beta*hbar*hbar/(tau**4) + 2./(tau**2)
c	where 	beta = 1./(boltzk * temp) in cgs units
      RETURN
      END

      SUBROUTINE DERIVAS(FCT,X0,DX1, Y,Y1,Y2,Y3,Y4,IFLAG)
C     COMPUTES OTH..4TH DERIVATIVES OF Y=FCT(X) AT X=X0 USING FIVE
C     ABSCISSAE AND CENTRAL DIFFERENCES. FCT MUST BE DECLARED EXTERNAL
C     IN CALLING ROUTINE. NOTE THAT FCT(X) MUST BE DEFINED FOR
C     X0-2DX<=X<=X0+2DX
C     A TEST OF THE DEFINITION OF THE FIRST DERIVATIVE: IFLAG=1 OUTPUT
C     SEE ABRAMOWITZ AND STEGUN, P. 914
	IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION  F(7)
9     X=X0 
      H=DX1 
      HSQ=H*H 
      IFLAG=0
      DO 10 I=1,7
10    F(I)=FCT(X-DFLOAT(4-I)*H)
      Y=F(4)
      A=F(1)+F(7)   
      A1=F(1)-F(7)
      B=F(2)+F(6)   
      B1=F(2)-F(6)
      C=F(3)+F(5)   
      C1=F(3)-F(5)

C     THESE ARE HIGHEST ORDER DERIVATIVES
      Y1 = (-A1+9.*B1-45.*C1)/(60.*H)
      Y2 = (A-13.5*B+135.*C-245.*Y)/(90.*HSQ)
      Y3 = (A1-8.*B1+13.*C1)/(8.*HSQ*H)
      Y4 = (-A+12.*B-39.*C+56.*Y)/(6.*HSQ*HSQ)
      RETURN
      END

        function xscale()
        IMPLICIT double precision (A-H,O-Z)
        external va
        common /temp/temp
        DATA ANGAU, FOURPI, HBAR/ 0.52917706, 12.566370614359,
     &                            1.05458876D-34/
        T= TEMP * 1.380662D-23
        FMP=2.0*1.67265D-27
        F2=HBAR*HBAR/(12.*FMP*T*T*ANGAU**2*1.D-20)
        h=0.05
        g0=0.0
        g1=0.0
        do 100 i=30,160
        r=h*dfloat(i)       ! in A.U.
        dr=r/100.0
	 CALL BETAJ(R, B0, BJ)
        CALL DERIVAS(VA,R,DR,V,DV,DDV,Z3,Z4,JFLAG)
        G = DEXP(-V/T)*R*R
        Gp=G*(1.+F2*((DV/T-4./R)*DV-2.*DDV))
        g0=g0+gp*B0**2
        g1=g1+gp*BJ**2
 100    continue
        g0=FOURPI*G0*h*(ANGAU*1.D-8)**3
        g1=FOURPI*G1*h*(ANGAU*1.D-8)**3
        xscale=g1/g0
        return
        end

        SUBROUTINE BETAJ(R, BETA0, BETA)

      IMPLICIT double precision (A-H,O-Z)
      CHARACTER*10 NOTE(2),LGAS,LABEL
      dimension BJ(9)
      COMMON /bpar/ Bpar(9, 4), n
      COMMON /COMUNIC/ NOTE,LGAS,LABEL,FMP,EPSP(2),RMP(2),LVALUE
      COMMON/ROTIF/ J1, JP1, J2, JP2

      DATA DEBYE / 2.541769713d-18 /

      x = R-6.0
	
	do i = 1, 9
	BJ(i) = Bpar(i, 1) / R**N + Bpar(i, 2) *
     &          DEXP( Bpar(i, 3) * x + Bpar(i, 4) * x**2 )
	end do

      Beta = BJ(1) + BJ(2) * J1*(J1+1) + BJ(3) *(J1*(J1+1))**2 + 
     &       BJ(4) * J2*(J2+1) + BJ(5) * (J2*(J2+1))**2 +
     &       BJ(6) * JP1*(JP1+1) + BJ(7) * (JP1*(JP1+1))**2 + 
     &       BJ(8) * JP2*(JP2+1) + BJ(9) * (JP2*(JP2+1))**2

      BETA = Beta * debye*1.0d-6
      BETA0= BJ(1)* debye*1.0d-6

      RETURN
      END


      FUNCTION VA(R)
C
C	double potential H2(v=0) - H2(v=0) initial ---> 
C	H2(v=1) - H2(v=1) as  final
C     H2(v=0) - H2(v=0) initial as the initial and final for 2nd
C     overtone
C     ==============================================================
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)

       FACTOR = (4.803250D-10)**2/(0.52916607D-8) *(1.D-7)

C	POTENTIAL H2(v=0) - H2(v=0)
c        NOTE='00_-_00_H2'
        RM  = 6.5*0.52917706D0
	RMIN1= 6.5D0
	E=    0.11163D-3
        EPSP= 0.11163D-3 * FACTOR
	A=    0.37996D7
	GAMA= 2.39588D0
	ALPHA=14.84169D0
	D =   1.02354D0
	C6=  -12.134D0
	C8=  -221.26D0
	C10= -0.50358D4

	C6S = C6 /( E * RMIN1**6 )
	C8S = C8 /( E * RMIN1**8 )
     	C10S = C10 /( E * RMIN1**10 )

	X = R/RM

	V_REP = A * X**GAMA  * DEXP(- ALPHA*X)
	F = 1.D0
	IF(X .LE. D ) F = DEXP(- ( D/X-1.)**2)
	V_DISP = F * (C6S/X**6 + C8S/X**8 + C10S/X**10)
	FF = V_REP + V_DISP
        VA= FF * EPSP
        RETURN
        END