| 1 | SUBROUTINE PIAH2(NU,NT,XROT,XTRA,T) |
|---|
| 2 | C |
|---|
| 3 | C XKH2 IS CALLED BY TH2HE TO COMPUT XK1, A FACTOR IN THE |
|---|
| 4 | C H2-H2 COLLISION-INDUCED ABSORPTION. |
|---|
| 5 | C |
|---|
| 6 | C THIS VERSION INCLUDES RESULTS OF RECENT LOW-TEMPERATURE FITS |
|---|
| 7 | C BY G. BIRNBAUM ET AL (NBS), AS IMPLEMENTED BY VIRGIL KUNDE |
|---|
| 8 | C AND GORDON BJORAKER, AUGUST 1982. |
|---|
| 9 | C |
|---|
| 10 | C THIS VERSION ALLOWS FOR A NON-EQUILIBRIUM PARA-HYDROGEN FRACTION. |
|---|
| 11 | C |
|---|
| 12 | C LAST CHANGED BY JOHN HORNSTEIN CSC FEB 28,1983 |
|---|
| 13 | C THIS CHANGE DELETED THE DOUBLE-TRANSITIONS, WHICH HAD BEEN |
|---|
| 14 | C INCORRECTLY FORMULATED IN BACHET ET AL. THE DELETION IS BY |
|---|
| 15 | C COMMENTING OUT WITH A "CD". |
|---|
| 16 | C |
|---|
| 17 | C THREE PARAMETERS ARE USED TO FIT THE H2-H2 OPACITY FROM |
|---|
| 18 | C 0 TO 2000 CM-1: STREN, TAU1, & TAU2. THE TEMPERATURE |
|---|
| 19 | C DEPENDENCE OF EACH PARAMETER IS HANDLED DIFFERENTLY. |
|---|
| 20 | C STREN WAS FIT TO EXPERIMENTAL VALUES IN DORE ET AL AND |
|---|
| 21 | C BACHET ET AL (1982) USING A LINEAR TEMPERATURE RELATION. |
|---|
| 22 | C TAU1 & TAU2 FOLLOW A POWER LAW OF THE FORM |
|---|
| 23 | C TAU = TAU0 * (T/T0)**BT AS SUGGESTED BY DORE ET AL. |
|---|
| 24 | C CONSTANTS FOR TAU1 ARE FROM DORE ET AL |
|---|
| 25 | C TAU2 = 1.23 * (TAU2 EVALUATED USING DORE POWER LAW). THIS MATCHES |
|---|
| 26 | C BACHET ET AL VALUES AT 195K & 297K (TAIL WT SINGLE + |
|---|
| 27 | C DOUBLE TRANSITION COEFFICIENTS). |
|---|
| 28 | C |
|---|
| 29 | C PARAMETER DESCRIPTION |
|---|
| 30 | C |
|---|
| 31 | C STRENGTH PARAMETER: STREN (SEE BACHET ET AL EQUATION 4B) |
|---|
| 32 | C STREN IS REALLY S/KB, TO AVOID LARGE POWERS OF 10. |
|---|
| 33 | C |
|---|
| 34 | C STREN = 6 * <A**2> * OMEGA**2 * I8(R/SIGMA) / K |
|---|
| 35 | C (UNITS: KELVIN*ANGSTROM**6) |
|---|
| 36 | C |
|---|
| 37 | C <A**2> IS MEAN SQUARE POLARIZABILITY UNIT: ANGSTROM**6 |
|---|
| 38 | C OMEGA IS QUADRUPOLE MOMENT UNIT: ESU*CM**2 |
|---|
| 39 | C R IS DISTANCE, SIGMA IS MOLECULAR SEPARATION, X = R/SIGMA |
|---|
| 40 | C I8(X) = 4*PI*INT(X**-8*EXP(-4.*E/(K*T)*(X**-12-X**-6))*X**2*DX) 00040000 |
|---|
| 41 | C E & SIGMA ARE PARAMETERS OF LENNARD-JONES INTERMOLECULAR POTENTIAL |
|---|
| 42 | C INT(...DX) IS INTEGRAL 0 TO INFINITY |
|---|
| 43 | C |
|---|
| 44 | C TIME PARAMETERS: TAU1 & TAU2 |
|---|
| 45 | C |
|---|
| 46 | C TAU1 CONTROLS HALF WIDTH NEAR LINE CENTER |
|---|
| 47 | C TAU2 CONTROLS EXPONENTIAL DECAY OF WINGS |
|---|
| 48 | C |
|---|
| 49 | C |
|---|
| 50 | C DOUBLE TRANSITIONS ARE INCLUDED. |
|---|
| 51 | C THE ABSORPTION DUE TO DOUBLE TRANSITIONS IS PROPORTIONAL TO C2 |
|---|
| 52 | C C2 = 4./45. * KAPPA**2 |
|---|
| 53 | C KAPPA = (APL-APP)/(1./3. * (APL+2.*APP)) |
|---|
| 54 | C APL & APP ARE POLARIZABILITIES PARALLEL |
|---|
| 55 | C AND PERPENDICULAR TO INTERNUCLEAR AXIS |
|---|
| 56 | C KAPPA = .375 FOR GROUND VIB STATE & J=0,1 OR 2 SEE KOLOS & WOLNIEWICZ |
|---|
| 57 | C JOURNAL OF CHEMICAL PHYSICS 46, 1426 (1967) TABLE 3 |
|---|
| 58 | IMPLICIT REAL (A-H,O-Z) |
|---|
| 59 | REAL NU,NSQR,K |
|---|
| 60 | LOGICAL ILIST |
|---|
| 61 | DIMENSION XROT(251),XTRA(251),XJ(10),G(10),A2(250),RHO(250,10), |
|---|
| 62 | A E(10),WN(10),WW(10),CG(10),BH(250),STREN(250),TAU1R(250), |
|---|
| 63 | B TAU2R(250),TAU1T(250),TAU2T(250),COREG0(2),COREG(10,2) |
|---|
| 64 | DIMENSION T(251) |
|---|
| 65 | DATA JMAX1/8/,JMAX2/10/, |
|---|
| 66 | . PI/3.141593/,HBAR/1.05450E-27/,C/2.997925E+10/,K/1.38054E-16/, |
|---|
| 67 | . HCK/1.43879/,XNLOS/2.687E19/,S0/178./,DSDT/.4091/, |
|---|
| 68 | . TAU1R0/7.25/,TAU2R0/3.45/,TAU1T0/7.25/,TAU2T0/3.45/, |
|---|
| 69 | . BT1/-.605/,BT2/-.607/,T0/273.15/,C2/0.0125/ |
|---|
| 70 | C |
|---|
| 71 | DIMENSION FPARA(99) |
|---|
| 72 | save pi,c,a2,ww,tau1t,tau2t,coreg0,rho,tau1r,tau2r,c1,k,bh,stren |
|---|
| 73 | save cg,coreg |
|---|
| 74 | C |
|---|
| 75 | C********************************************************************** |
|---|
| 76 | C |
|---|
| 77 | C |
|---|
| 78 | C G(J) ARE STATISTICAL WTS: |
|---|
| 79 | C G(J)=1 EVEN ROTATIONAL STATES, |
|---|
| 80 | C G(J)=3 ODD ROTATIONAL STATES. |
|---|
| 81 | ILIST=.FALSE. |
|---|
| 82 | C SET THE ORTHO/PARA TO EQUILIBRIUM -CPM |
|---|
| 83 | NPARA=-1. |
|---|
| 84 | DO 10 J = 1,JMAX2 |
|---|
| 85 | COREG(J,1)=1. |
|---|
| 86 | COREG(J,2)=0. |
|---|
| 87 | XJ(J) = FLOAT(J-1) |
|---|
| 88 | G(J) = 1.0 |
|---|
| 89 | IF (MOD(J,2) .EQ. 0) G(J) = 3.0 |
|---|
| 90 | 10 E(J) = H2ENER(0.0,XJ(J)) |
|---|
| 91 | C ****************** WARNING *********************** |
|---|
| 92 | C CHANGE MADE BY REGIS COURTIN (APRIL 1986). |
|---|
| 93 | C THE FOLLOWING COEFFICIENTS ARE INTRODUCED FOR THE CASE |
|---|
| 94 | C OF COLLISIONS BETWEEN H2 AND N2 MOLECULES. |
|---|
| 95 | C THE CORRECTED ABSORPTION COEFFICIENTS FIT THE DATA OF |
|---|
| 96 | C DORE ET AL (1986), PREPRINT. |
|---|
| 97 | C ****************************************************** |
|---|
| 98 | COREG0(1)=22.35 |
|---|
| 99 | COREG0(2)=-0.56 |
|---|
| 100 | COREG(1,1)=0.36 |
|---|
| 101 | COREG(1,2)=0.18 |
|---|
| 102 | COREG(2,1)=10.91 |
|---|
| 103 | COREG(2,2)=-0.433 |
|---|
| 104 | C ****************************************************** |
|---|
| 105 | C CG(J) = (2*J+1) * (CLEBSCH-GORDAN COEFF <J 2 J'> )**2 |
|---|
| 106 | DO 20 J = 1,JMAX1 |
|---|
| 107 | CG(J) = (2.0*XJ(J)+1.0) * (3.0*(XJ(J)+1.0)*(XJ(J)+2.0)) / |
|---|
| 108 | . (2.0*(2.0*XJ(J)+1.0)*(2.0*XJ(J)+3.0)) |
|---|
| 109 | WN(J) = E(J+2) - E(J) |
|---|
| 110 | 20 WW(J) = 2.0 * PI * C * WN(J) |
|---|
| 111 | NSQR = (XNLOS*1.0E-24)**2 |
|---|
| 112 | C1 = 2. * PI**2 * NSQR/(3. * HBAR * C) |
|---|
| 113 | C LIST HEADING FOR ORTHO AND PARA PROFILES: |
|---|
| 114 | IF ( (NPARA .LE. 0) .OR. (.NOT. ILIST)) GO TO 7000 |
|---|
| 115 | WRITE(6,5000) |
|---|
| 116 | 5000 FORMAT(//,' EQUILIBRIUM AND ACTUAL FRACTIONS OF PARA-', |
|---|
| 117 | . ' AND ORTHO-HYDROGEN:',/,' LAYER',7X,'FPARA(EQU)', |
|---|
| 118 | . 7X,'FPARA(HERE)',7X,'FORTHO(EQU)',7X,'FORTHO(HERE)') |
|---|
| 119 | 7000 DO 90 IT = 1,NT |
|---|
| 120 | BH(IT) = HBAR / (K*T(IT)) |
|---|
| 121 | C EVALUATE STREN, TAU1R, TAU2R AT DESIRED T |
|---|
| 122 | STREN(IT) = (S0 + DSDT*T(IT)) |
|---|
| 123 | TAU1R(IT) = TAU1R0 * (T(IT)/T0)**BT1 * 1.0E-14 |
|---|
| 124 | TAU2R(IT) = TAU2R0 * (T(IT)/T0)**BT2 * 1.0E-14 |
|---|
| 125 | TAU1T(IT) = TAU1T0 * (T(IT)/T0)**BT1 * 1.0E-14 |
|---|
| 126 | TAU2T(IT) = TAU2T0 * (T(IT)/T0)**BT2 * 1.0E-14 |
|---|
| 127 | C COMPUTE THE FULL PARTION FUNCTION Z, USED IN EQUILIBRIUM, |
|---|
| 128 | C WHERE ORTHO AND PARA ARE CONVENIENTLY TREATED AS DIFFERENT |
|---|
| 129 | C STATES OF THE SAME SPECIES. |
|---|
| 130 | Z = 0.0 |
|---|
| 131 | DO 50 J = 1,JMAX2 |
|---|
| 132 | RHO(IT,J) = EXP(-1.*HCK*E(J)/T(IT)) |
|---|
| 133 | 50 Z = Z + (2.0*XJ(J)+1.0)*G(J)*RHO(IT,J) |
|---|
| 134 | IF (NPARA .LE. 0) GO TO 54 |
|---|
| 135 | C COMPUTE THE PARTITION FUNCTIONS ZPARA AND ZORTHO USED FOR |
|---|
| 136 | C NON-EQUILIBRIUM RATIOS, WHERE IT IS CONVENIENT TO TREAT |
|---|
| 137 | C ORTHO AND PARA AS DISTINCT SPECIES. THE NUCLEAR SPIN |
|---|
| 138 | C FACTORS G(J) CANCEL OUT IN THIS CASE. |
|---|
| 139 | ZPARA = 0. |
|---|
| 140 | ZORTHO = 0. |
|---|
| 141 | DO 1000 J=1,JMAX2,2 |
|---|
| 142 | ZPARA = ZPARA + (2.*XJ(J) + 1.)*RHO(IT,J) |
|---|
| 143 | JJ = J + 1 |
|---|
| 144 | 1000 ZORTHO = ZORTHO + (2.*XJ(JJ) + 1.)*RHO(IT,JJ) |
|---|
| 145 | C COMPUTE AND LIST THE EQUILIBRIUM AND ACTUAL PROFILES: |
|---|
| 146 | FEPARA = ZPARA/Z |
|---|
| 147 | FEORTH = 3.*ZORTHO/Z |
|---|
| 148 | FORTHO = 1. - FPARA(IT) |
|---|
| 149 | IF (ILIST) WRITE(6,2000) IT, FEPARA, FPARA(IT), FEORTH, FORTHO |
|---|
| 150 | 2000 FORMAT(' ',I4,2F15.5,10X,2F15.5) |
|---|
| 151 | C FORM A NEW RHO WHICH EQUALS RHO(PARA) WHEN XJ IS EVEN |
|---|
| 152 | C (INDEX J IS ODD) AND EQUALS RHO(ORTHO) WHEN XJ IS ODD. |
|---|
| 153 | DO 3000 J=1,JMAX2,2 |
|---|
| 154 | RHO(IT,J) = FPARA(IT)*RHO(IT,J)/ZPARA |
|---|
| 155 | JJ = J + 1 |
|---|
| 156 | 3000 RHO(IT,JJ) = FORTHO*RHO(IT,JJ)/ZORTHO |
|---|
| 157 | GO TO 57 |
|---|
| 158 | C EQUILIBRIUM HYDROGEN STATISTICAL WEIGHTS |
|---|
| 159 | 54 DO 55 J = 1,JMAX2 |
|---|
| 160 | 55 RHO(IT,J) = G(J)*RHO(IT,J) / Z |
|---|
| 161 | 57 A2(IT) = 0.0 |
|---|
| 162 | DO 60 J = 1,JMAX1 |
|---|
| 163 | 60 A2(IT) = A2(IT) + XJ(J)*(XJ(J)+1.0)*(2.0*XJ(J)+1.0)*RHO(IT,J)/ |
|---|
| 164 | . ((2.0*XJ(J)-1.0)*(2.0*XJ(J)+3.0)) |
|---|
| 165 | 90 CONTINUE |
|---|
| 166 | RETURN |
|---|
| 167 | ENTRY OPAH2(NU,NT,XROT,XTRA,T) |
|---|
| 168 | W = 2.0 * PI * C * NU |
|---|
| 169 | DO 200 IT = 1,NT |
|---|
| 170 | DBL = 0. |
|---|
| 171 | FR = 0. |
|---|
| 172 | C EVALUATE TRANSLATIONAL PART OF SHAPE FACTOR F |
|---|
| 173 | FT = A2(IT) * GAMFCN(W,T(IT),TAU1T(IT),TAU2T(IT)) |
|---|
| 174 | .* COREG0(1)*T(IT)**COREG0(2) |
|---|
| 175 | DO 125 J = 1,JMAX1 |
|---|
| 176 | C EVALUATE THE ROTATIONAL PART OF F |
|---|
| 177 | FR = FR + CG(J) * (RHO(IT,J)*GAMFCN(W-WW(J),T(IT),TAU1R(IT), |
|---|
| 178 | . TAU2R(IT))+RHO(IT,J+2)*GAMFCN(W+WW(J),T(IT),TAU1R(IT), |
|---|
| 179 | . TAU2R(IT))) |
|---|
| 180 | . * COREG(J,1)*T(IT)**COREG(J,2) |
|---|
| 181 | C DBL = DBL + CG(J) * (RHO(IT,J) + RHO(IT,J+2)) |
|---|
| 182 | 125 CONTINUE |
|---|
| 183 | C ADD ON PART OF DOUBLE TRANSITION OPACITY (BACHET ET AL EQN 11) |
|---|
| 184 | CD FR = FR * (1.0 + C2 * A2(IT)) |
|---|
| 185 | C MORE DOUBLE TRANS (BACHET ET AL EQN 12) |
|---|
| 186 | CD FT = FT * (1.0 + C2 * A2(IT)) |
|---|
| 187 | C AND STILL MORE DOUBLE TRANS (BACHET ET AL EQN 13) |
|---|
| 188 | CD FT = FT + C2 * DBL*DBL*GAMFCN(W,T(IT),TAU1T(IT),TAU2T(IT)) |
|---|
| 189 | F = FT + FR |
|---|
| 190 | C EVALUATE & ADD ON DOUBLE TRANSITIONS (BACHET ET AL EQN 14) |
|---|
| 191 | CD DO 130 J1=1,4 |
|---|
| 192 | CD DO 130 J2=1,4 |
|---|
| 193 | CD IF (J1.NE.(J1+2).OR.J2.NE.(J1+2)) |
|---|
| 194 | CD A F = F + C2*RHO(IT,J1)*CG(J1)*RHO(IT,J2)*CG(J2) |
|---|
| 195 | CD B * GAMFCN(W-WW(J1)-WW(J2),T(IT),TAU1R(IT),TAU2R(IT)) |
|---|
| 196 | CD IF ((J1+2).NE.(J2+2).OR.J2.NE.J1) |
|---|
| 197 | CD A F = F + C2*RHO(IT,J1+2)*CG(J1)*RHO(IT,J2)*CG(J2) |
|---|
| 198 | CD B * GAMFCN(W+WW(J1)-WW(J2),T(IT),TAU1R(IT),TAU2R(IT)) |
|---|
| 199 | CD IF ((J1+2).NE.J2.OR.(J2+2).NE.J1) |
|---|
| 200 | CD A F = F + C2*RHO(IT,J1+2)*CG(J1)*RHO(IT,J2+2)*CG(J2) |
|---|
| 201 | CD B * GAMFCN(W+WW(J1)+WW(J2),T(IT),TAU1R(IT),TAU2R(IT)) |
|---|
| 202 | CD IF (J1.NE.J2.OR.(J2+2).NE.(J1+2)) |
|---|
| 203 | CD A F = F + C2*RHO(IT,J1)*CG(J1)*RHO(IT,J2+2)*CG(J2) |
|---|
| 204 | CD B * GAMFCN(W-WW(J1)+WW(J2),T(IT),TAU1R(IT),TAU2R(IT) |
|---|
| 205 | CD130 CONTINUE |
|---|
| 206 | XROT(IT) = C1 * K * W * (1.0-EXP(-BH(IT)*W)) * STREN(IT)*FR |
|---|
| 207 | XTRA(IT) = XROT(IT)*FT/FR |
|---|
| 208 | 200 CONTINUE |
|---|
| 209 | RETURN |
|---|
| 210 | END |
|---|
