      SUBROUTINE PIAH2(NU,NT,XROT,XTRA,T)
C
C     XKH2 IS CALLED BY TH2HE TO COMPUT XK1, A FACTOR IN THE
C     H2-H2 COLLISION-INDUCED ABSORPTION.
C
C     THIS VERSION INCLUDES RESULTS OF RECENT LOW-TEMPERATURE FITS
C     BY G. BIRNBAUM ET AL (NBS), AS IMPLEMENTED BY VIRGIL KUNDE
C     AND GORDON BJORAKER, AUGUST 1982.
C
C     THIS VERSION ALLOWS FOR A NON-EQUILIBRIUM PARA-HYDROGEN FRACTION.
C
C     LAST CHANGED BY JOHN HORNSTEIN   CSC   FEB 28,1983
C     THIS CHANGE DELETED THE DOUBLE-TRANSITIONS, WHICH HAD BEEN
C     INCORRECTLY FORMULATED IN BACHET ET AL. THE DELETION IS BY
C     COMMENTING OUT WITH A "CD".
C
C  THREE PARAMETERS ARE USED TO FIT THE H2-H2 OPACITY FROM
C  0 TO 2000 CM-1:  STREN, TAU1, & TAU2.  THE TEMPERATURE
C  DEPENDENCE OF EACH PARAMETER IS HANDLED DIFFERENTLY.
C  STREN WAS FIT TO EXPERIMENTAL VALUES IN DORE ET AL AND
C  BACHET ET AL (1982) USING A LINEAR TEMPERATURE RELATION.
C  TAU1 & TAU2 FOLLOW A POWER LAW OF THE FORM
C  TAU = TAU0 * (T/T0)**BT AS SUGGESTED BY DORE ET AL.
C  CONSTANTS FOR TAU1 ARE FROM DORE ET AL
C  TAU2 = 1.23 * (TAU2 EVALUATED USING DORE POWER LAW). THIS MATCHES
C  BACHET ET AL VALUES AT 195K & 297K (TAIL WT SINGLE +
C  DOUBLE TRANSITION COEFFICIENTS).
C
C  PARAMETER DESCRIPTION
C
C  STRENGTH PARAMETER:  STREN     (SEE BACHET ET AL EQUATION 4B)
C  STREN IS REALLY S/KB, TO AVOID LARGE POWERS OF 10.
C
C  STREN = 6 * <A**2> * OMEGA**2 * I8(R/SIGMA) / K
C  (UNITS:  KELVIN*ANGSTROM**6)
C
C  <A**2> IS MEAN SQUARE POLARIZABILITY   UNIT:  ANGSTROM**6
C  OMEGA  IS QUADRUPOLE MOMENT            UNIT:  ESU*CM**2 
C  R IS DISTANCE, SIGMA IS MOLECULAR SEPARATION, X = R/SIGMA
C  I8(X) = 4*PI*INT(X**-8*EXP(-4.*E/(K*T)*(X**-12-X**-6))*X**2*DX)      00040000
C  E & SIGMA ARE PARAMETERS OF LENNARD-JONES INTERMOLECULAR POTENTIAL
C  INT(...DX) IS INTEGRAL 0 TO INFINITY
C
C  TIME PARAMETERS:  TAU1 & TAU2
C
C  TAU1    CONTROLS HALF WIDTH NEAR LINE CENTER
C  TAU2    CONTROLS EXPONENTIAL DECAY OF WINGS
C
C
C  DOUBLE TRANSITIONS ARE INCLUDED.
C  THE ABSORPTION DUE TO DOUBLE TRANSITIONS IS PROPORTIONAL TO C2
C  C2 = 4./45. * KAPPA**2
C  KAPPA = (APL-APP)/(1./3. * (APL+2.*APP))
C  APL & APP ARE POLARIZABILITIES PARALLEL
C  AND PERPENDICULAR TO INTERNUCLEAR AXIS
C  KAPPA = .375 FOR GROUND VIB STATE & J=0,1 OR 2 SEE KOLOS & WOLNIEWICZ
C  JOURNAL OF CHEMICAL PHYSICS 46, 1426 (1967) TABLE 3
      IMPLICIT REAL (A-H,O-Z)
      REAL NU,NSQR,K
      LOGICAL ILIST
      DIMENSION XROT(251),XTRA(251),XJ(10),G(10),A2(250),RHO(250,10),
     A E(10),WN(10),WW(10),CG(10),BH(250),STREN(250),TAU1R(250),
     B TAU2R(250),TAU1T(250),TAU2T(250),COREG0(2),COREG(10,2)
      DIMENSION T(251)
      DATA JMAX1/8/,JMAX2/10/,
     .   PI/3.141593/,HBAR/1.05450E-27/,C/2.997925E+10/,K/1.38054E-16/,
     .   HCK/1.43879/,XNLOS/2.687E19/,S0/178./,DSDT/.4091/,
     .   TAU1R0/7.25/,TAU2R0/3.45/,TAU1T0/7.25/,TAU2T0/3.45/,
     .   BT1/-.605/,BT2/-.607/,T0/273.15/,C2/0.0125/
C
      DIMENSION  FPARA(99)
      save pi,c,a2,ww,tau1t,tau2t,coreg0,rho,tau1r,tau2r,c1,k,bh,stren
      save cg,coreg
C
C**********************************************************************
C
C
C  G(J) ARE STATISTICAL WTS:
C  G(J)=1 EVEN ROTATIONAL STATES,
C  G(J)=3 ODD  ROTATIONAL STATES.
      ILIST=.FALSE.
C SET THE ORTHO/PARA TO EQUILIBRIUM -CPM
      NPARA=-1.
      DO 10  J = 1,JMAX2
        COREG(J,1)=1.
        COREG(J,2)=0.
        XJ(J) = FLOAT(J-1)
        G(J) = 1.0
        IF (MOD(J,2) .EQ. 0)  G(J) = 3.0
 10     E(J) = H2ENER(0.0,XJ(J))
C  ******************   WARNING   ***********************
C  CHANGE MADE BY REGIS COURTIN (APRIL 1986).
C  THE FOLLOWING COEFFICIENTS ARE INTRODUCED FOR THE CASE
C  OF COLLISIONS BETWEEN H2 AND N2 MOLECULES.
C  THE CORRECTED ABSORPTION COEFFICIENTS FIT THE DATA OF
C  DORE ET AL (1986), PREPRINT.
C  ******************************************************
      COREG0(1)=22.35
      COREG0(2)=-0.56
      COREG(1,1)=0.36
      COREG(1,2)=0.18
      COREG(2,1)=10.91
      COREG(2,2)=-0.433
C  ******************************************************
C CG(J) = (2*J+1) * (CLEBSCH-GORDAN COEFF <J 2 J'> )**2
      DO 20  J = 1,JMAX1
        CG(J) = (2.0*XJ(J)+1.0) * (3.0*(XJ(J)+1.0)*(XJ(J)+2.0)) /
     .          (2.0*(2.0*XJ(J)+1.0)*(2.0*XJ(J)+3.0))
        WN(J) = E(J+2) - E(J)
 20     WW(J) = 2.0 * PI * C * WN(J)
      NSQR = (XNLOS*1.0E-24)**2
      C1 = 2. * PI**2 * NSQR/(3. * HBAR * C)
C     LIST HEADING FOR ORTHO AND PARA PROFILES:
      IF ( (NPARA .LE. 0) .OR. (.NOT. ILIST))  GO TO 7000
      WRITE(6,5000)
5000  FORMAT(//,' EQUILIBRIUM AND ACTUAL FRACTIONS OF PARA-',
     .      ' AND ORTHO-HYDROGEN:',/,' LAYER',7X,'FPARA(EQU)',
     .      7X,'FPARA(HERE)',7X,'FORTHO(EQU)',7X,'FORTHO(HERE)')
7000  DO 90 IT = 1,NT
        BH(IT) = HBAR / (K*T(IT))
C     EVALUATE    STREN, TAU1R, TAU2R AT DESIRED T
        STREN(IT) = (S0 + DSDT*T(IT))
        TAU1R(IT) = TAU1R0 * (T(IT)/T0)**BT1 * 1.0E-14
        TAU2R(IT) = TAU2R0 * (T(IT)/T0)**BT2 * 1.0E-14
        TAU1T(IT) = TAU1T0 * (T(IT)/T0)**BT1 * 1.0E-14
        TAU2T(IT) = TAU2T0 * (T(IT)/T0)**BT2 * 1.0E-14
C     COMPUTE THE FULL PARTION FUNCTION Z, USED IN EQUILIBRIUM,
C     WHERE ORTHO AND PARA ARE CONVENIENTLY TREATED AS DIFFERENT
C     STATES OF THE SAME SPECIES.
        Z = 0.0
        DO 50 J = 1,JMAX2
          RHO(IT,J) = EXP(-1.*HCK*E(J)/T(IT))
 50       Z = Z + (2.0*XJ(J)+1.0)*G(J)*RHO(IT,J)
        IF (NPARA .LE. 0)  GO TO 54
C     COMPUTE THE PARTITION FUNCTIONS ZPARA AND ZORTHO USED FOR
C     NON-EQUILIBRIUM RATIOS, WHERE IT IS CONVENIENT TO TREAT
C     ORTHO AND PARA AS DISTINCT SPECIES. THE NUCLEAR SPIN
C     FACTORS G(J) CANCEL OUT IN THIS CASE.
        ZPARA = 0.
        ZORTHO = 0.
        DO 1000  J=1,JMAX2,2
          ZPARA = ZPARA + (2.*XJ(J) + 1.)*RHO(IT,J)
          JJ = J + 1
1000      ZORTHO = ZORTHO + (2.*XJ(JJ) + 1.)*RHO(IT,JJ)
C     COMPUTE AND LIST THE EQUILIBRIUM AND ACTUAL PROFILES:
        FEPARA = ZPARA/Z
        FEORTH = 3.*ZORTHO/Z
        FORTHO = 1. - FPARA(IT)
        IF (ILIST)  WRITE(6,2000)  IT, FEPARA, FPARA(IT), FEORTH, FORTHO
2000    FORMAT(' ',I4,2F15.5,10X,2F15.5)
C     FORM A NEW RHO WHICH EQUALS RHO(PARA) WHEN XJ IS EVEN
C     (INDEX J IS ODD) AND EQUALS RHO(ORTHO) WHEN XJ IS ODD.
      DO 3000  J=1,JMAX2,2
        RHO(IT,J) = FPARA(IT)*RHO(IT,J)/ZPARA
        JJ = J + 1
3000    RHO(IT,JJ) = FORTHO*RHO(IT,JJ)/ZORTHO
      GO TO 57
C  EQUILIBRIUM HYDROGEN STATISTICAL WEIGHTS
54      DO 55  J = 1,JMAX2
55        RHO(IT,J) = G(J)*RHO(IT,J) / Z
57      A2(IT) = 0.0
        DO 60  J = 1,JMAX1
60        A2(IT) = A2(IT) + XJ(J)*(XJ(J)+1.0)*(2.0*XJ(J)+1.0)*RHO(IT,J)/
     .                  ((2.0*XJ(J)-1.0)*(2.0*XJ(J)+3.0))
 90   CONTINUE
      RETURN
      ENTRY OPAH2(NU,NT,XROT,XTRA,T)
      W = 2.0 * PI * C * NU
      DO 200  IT = 1,NT
        DBL = 0.
        FR = 0.
C    EVALUATE TRANSLATIONAL PART OF SHAPE FACTOR F
        FT = A2(IT) * GAMFCN(W,T(IT),TAU1T(IT),TAU2T(IT))
     .* COREG0(1)*T(IT)**COREG0(2)
        DO 125 J = 1,JMAX1
C     EVALUATE THE ROTATIONAL PART OF F
          FR = FR + CG(J) * (RHO(IT,J)*GAMFCN(W-WW(J),T(IT),TAU1R(IT),
     .         TAU2R(IT))+RHO(IT,J+2)*GAMFCN(W+WW(J),T(IT),TAU1R(IT),
     .         TAU2R(IT)))
     .                     * COREG(J,1)*T(IT)**COREG(J,2)
C         DBL = DBL + CG(J) * (RHO(IT,J) + RHO(IT,J+2))
 125  CONTINUE
C  ADD ON PART OF DOUBLE TRANSITION OPACITY (BACHET ET AL EQN 11)
CD    FR = FR * (1.0 + C2 * A2(IT))
C  MORE DOUBLE TRANS  (BACHET ET AL EQN 12)
CD    FT = FT * (1.0 + C2 * A2(IT))
C  AND STILL MORE DOUBLE TRANS (BACHET ET AL EQN 13)
CD    FT = FT + C2 * DBL*DBL*GAMFCN(W,T(IT),TAU1T(IT),TAU2T(IT))
      F = FT + FR
C  EVALUATE & ADD ON DOUBLE TRANSITIONS (BACHET ET AL EQN 14)
CD    DO 130 J1=1,4
CD    DO 130 J2=1,4
CD    IF (J1.NE.(J1+2).OR.J2.NE.(J1+2))
CD   A  F = F + C2*RHO(IT,J1)*CG(J1)*RHO(IT,J2)*CG(J2)
CD   B  * GAMFCN(W-WW(J1)-WW(J2),T(IT),TAU1R(IT),TAU2R(IT))
CD    IF ((J1+2).NE.(J2+2).OR.J2.NE.J1)
CD   A  F = F + C2*RHO(IT,J1+2)*CG(J1)*RHO(IT,J2)*CG(J2)
CD   B  * GAMFCN(W+WW(J1)-WW(J2),T(IT),TAU1R(IT),TAU2R(IT))
CD    IF ((J1+2).NE.J2.OR.(J2+2).NE.J1)
CD   A  F = F + C2*RHO(IT,J1+2)*CG(J1)*RHO(IT,J2+2)*CG(J2)
CD   B  * GAMFCN(W+WW(J1)+WW(J2),T(IT),TAU1R(IT),TAU2R(IT))
CD    IF (J1.NE.J2.OR.(J2+2).NE.(J1+2))
CD   A  F = F + C2*RHO(IT,J1)*CG(J1)*RHO(IT,J2+2)*CG(J2)
CD   B  * GAMFCN(W-WW(J1)+WW(J2),T(IT),TAU1R(IT),TAU2R(IT)
CD130 CONTINUE
      XROT(IT) = C1 * K * W * (1.0-EXP(-BH(IT)*W)) * STREN(IT)*FR
      XTRA(IT) = XROT(IT)*FT/FR
200   CONTINUE
      RETURN
      END
