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 * * OMEGA**2 * I8(R/SIGMA) / K C (UNITS: KELVIN*ANGSTROM**6) C C 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 )**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