C PROGRAM TO GENERATE EXCITATION RATES FOR CO IN COLLISIONS W/ H2O C USING 'ECS-EP' RATE LAW; MILLOT, J. CHEM. PHYS. 93, 8001 (1990) C SEE S. GREEN, AP. J. 412, 436 (1993) C C FOLLOWING VALUES CONTROL CALCULATION -- CHANGE AS NECESSARY. C TEMP IS THE KINETIC TEMPERATURE C JMAX IS THE HIGHEST ROTATIONAL LEVEL, I.E., C RATE(JI->JF) IS CALCULATED FOR JI=0,JMAX AND JF=0,JMAX C C SYSTEM DEPENDENT PARAMETERS ARE SET IN SUBROUTINE MILLOT: C A0,T0,XN; BETA,GAMMA; XLC; LABEL; BE ALSO SET THERE C PARAMETERS FOR N2-H2O ARE USED; HOWEVER, THE FUNDAMENTAL RATES C ARE MODIFIED TO APPROXIMATE VALUES FOR CO-H2O COLLISIONS: C (1) AN OVERALL SCALING FACTOR IS USED TO INCREASE C TOTAL RATES TO AGREE WITH ROOM TEMP CARS LINEWIDTH VALUES C (TO APPROX. 10%). (2) RATES ARE APPORTIONED TO ODD DELTA-J C (3/7) AS WELL AS EVEN DELTA-J (4/7). NOTE THAT THE RELATIVE C WEIGHTS IN THIS SPLIT ARE ONLY AN EDUCATED GUESS. C C NEEDS ANGULAR MOMENTUM COUPLING FUNCTION SUBROUTINES C THREEJ, BINOM, PARITY C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DOUBLE PRECISION KBOLZ LOGICAL LPRT DIMENSION W(90),OMEGA(90) COMMON/CONSTS/BE,BET,TEMP,URED C N.B. LMAX HAS TO ACCOMMODATE HIGHEST L-VALUE NEEDED IN IOS/ECS C SUMMATION. LMAX=80 SHOULD ACCOMODATE TO TEMPS OF 1000 K. C W AND OMEGA MUST BE DIMENSIONED LMAX OR LARGER DATA LMAX/80/ C JMAX IS HIGHEST J-VALUE DESIRED; VALUES HIGHER THAN ABOUT 30-40 C MAY STRETCH VALIDITY OF THE MODEL DATA JMAX/20/ C C LOGICAL LPRT CONTROLS PRINT OUTPUT C LPRT=.FALSE. SUPPRESSES ALL OUTPUT EXCEPT FINAL RATE VALUES. DATA LPRT/.TRUE./ C C SOME PHYSICAL CONSTANTS... DATA PI/3.14159 26535 89793 D0/,AVOGAD/6.024D23/,KBOLZ/1.3805D-16/ C STATEMENT FUNCTION FOR 2*L+1 ... F(L)=2.D0*L+1.D0 C C SUPPRESS UNDERFLOWS VIA FORTRAN H EXTENDED ERROR FACILITY C CALL ERRSET(208,256,-1,1,1,1) C C SET TEMPERATURE -- PASSED IN COMMON TEMP=200. C MX=LMAX CALL MILLOT(MX,W,OMEGA,LPRT) C C CALCULATE AVERAGE VELOCITY = SQRT(8*KT/PI*URED) IN CM/SEC VBAR=SQRT(8.D0*KBOLZ*TEMP/(PI*(URED/AVOGAD))) C C CONV * SIG(ANG**2) -> FWHH IN 1.E-3 CM**-1/ATM (I.E., MK/ATM) CONV=8.75*SQRT(URED*TEMP) CONV=1.D3/CONV C C LOOP OVER JI (INITIAL) AND JF (FINAL) ROTOR LEVELS DO 1000 JI=0,JMAX SUM=0. DO 1100 JF=0,JMAX IF (JI.EQ.JF) GO TO 1100 C CALCULATE RATE JI->JF JX=MAX(JI,JF) EX=BET*(JI*(JI+1)-JX*(JX+1)) EX=EXP(EX/TEMP) AJ=OMEGA(JX ) LTOP=MIN(LMAX,JI+JF) LMIN=IABS(JI-JF) RATE=0. DO 2000 L=LMIN,LTOP TJ=THREEJ(JI,JF,L) WL=W(L ) AL=OMEGA(L ) 2000 RATE=RATE+F(L)*TJ*TJ*WL/(AL*AL) RATE=RATE*F(JF)*(EX*AJ*AJ) C N.B. RATE HERE IS HWHH IN CM**-1/ATM (MK/ATM); CONVERT TO ANG**2 SIG=RATE*(2./CONV) C AND THEN MULT BY AVG VELOCITY TO CONVERT CROSS SECTION TO RATE ROUT=VBAR*1.D-16*SIG C IF (LPRT.AND.RATE.GE.1.D-3) WRITE(6,600) JI,JF,SIG,ROUT C 600 FORMAT(' RATE(',I3,'->',I3,') =',F12.2,' ANG**2 =',1P,D15.3, C 1 ' CM**3/SEC ') WRITE(6,600) JI,JF,TEMP,ROUT 600 FORMAT(' RATE(',I2,'->',I2,',T=',F6.1,') =',1P,D14.3, 1 ' CM**3/SEC') SUM=SUM+RATE 1100 CONTINUE C CONVERSION TO MK/AMAGAT TO CF SOME DATA SX=SUM*(TEMP/273.) 1000 IF (LPRT) WRITE(6,601) SUM,SUM*(2./CONV) 601 FORMAT(5X,'GAMMA(HWHH) =',F12.2,' MK/ATM =',F12.2,' ANG**2 ') STOP END SUBROUTINE MILLOT(MX,W,OMEGA,LPRT) IMPLICIT DOUBLE PRECISION (A-H,O-Z) LOGICAL LPRT DIMENSION W(MX),OMEGA(MX) CHARACTER*8 LABEL(2) COMMON/CONSTS/BE,BET,TEMP,URED C C ECS-EP CALCULATION (SEE MILLOT JCP 93, 8001 (1990) C IF (LPRT) WRITE(6,601) 601 FORMAT(' ECS-EP (MILLOT) CALCULATION OF W-MATRIX'/) C C -------------------------------------------------------- C SET SYSTEM DEPENDENT PARAMETERS C BELOW FOR N2-H2O; OBTAINED BY FITTING RAMAN LINEWIDTHS. C WITH CORRECT ROTATION CONSTANT FOR CO, AND C XN FOR CO-H2O FROM HARTMANN, ET AL., APPL. OPT 27, 3063 (1988) C N.B. TEMPERATURE (TEMP=XXXX) IS SET IN MAIN PROGRAM. LABEL(1)='N2-H2O: ' LABEL(2)=' MILLOT ' C BE=2.01 BE=1.92265 A0=23.71 T0=300. C XN= 1.155 XN=0.75 GAMMA=0.7158 BETA=0.05416 AMU=10.96 XLC=11.40 C THIS ENDS SYSTEM DEPENDENT CONSTANT INITIALIZATION C -------------------------------------------------------- C C N.B. TEMP, REDUCED MASS, AND ROTATION CONSTANT PASSED VIA COMMON C ROTATION CONSTANT GIVEN IN CM**-1 AND CONVERTED TO KELVINS BET=BE*(300./208.5) AT=A0*(T0/TEMP)**XN URED=AMU C CALCULATE ECS ADIABATICITY FACTOR ... OFACT=.06475*BE*XLC*SQRT(AMU/TEMP) C IF (LPRT) WRITE(6,602) LABEL,A0,T0,XN,BETA,GAMMA,XLC 602 FORMAT(1X,'PARAMETERS FOR ',2A8/ 1 /' A0, T0, XN =',F9.3,F8.1,F9.4/ 2 ' BETA, GAMMA =',2F12.6/ 3 ' SCALING LENGTH, LC =',F6.2,' ANGSTROMS') IF (LPRT) WRITE(6,603) BE,BET,TEMP,AT,AMU 603 FORMAT(' ROTATION CONSTANT =',F8.4,' CM**-1',' = ',F9.4,' KELVIN'/ 1 ' REQUESTED TEMPERATURE =',F8.1,' A(T) = ',F8.2,' MK/ATM'/ 2 ' REDUCED MASS =',F9.4,' A.M.U'/) IF (LPRT) WRITE(6,604) 604 FORMAT(' MODIFIED FOR CO-H2O BY S. GREEN NOV 92'//) C DO 1000 I=1,MX L= I C CALCULATE BASE RATE W(0->L) BY 'EXPONENTIAL-POWER LAW' FORMULA EL=BET*L*(L+1)/TEMP EX=EXP(-BETA*EL) XL=1.D0/(L*(L+1))**GAMMA W(I)=AT*EX*XL C SCALE UP BY 1.15 TO APPROX MATCH DIFFERENCE IN LINEWIDTHS C BETWEEN N2 AND CO, AND ALLOCATE BETWEEN EVEN AND ODD. IF (L-2*(L/2).EQ.0) THEN W(I)=W(I)*(1.15 )*(4./7.) ELSE W(I)=W(I)*(1.15 )*(3./7.) ENDIF C CALCULATE ECS OMEGA(L) USING E(L)-E(L-2) TAU=OFACT*(4*L-2) OMEGA(I)=1./(1.+TAU*TAU/6.) IF (LPRT) WRITE(6,600) L,W(I),OMEGA(I) 600 FORMAT(' L, W(0->L), A(L)',I4,F8.2,F8.4) 1000 CONTINUE RETURN END FUNCTION THREEJ (J1,J2,J3) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C C COMPUTATION OF SPECIAL WIGNER 3J COEFFICIENT WITH C VANISHING PROJECTIONS. SEE EDMUNDS, P. 50. C C STATEMENT FUNCTION FOR DELTA ASSOCIATED W/ RACAH AND SIXJ SYMBOLS C SEE E.G. EDMUNDS P. 99. DELTA(I,J,K)= DSQRT(1.D0/ ( BINOM(I+J+K+1,I+J-K) * 1 BINOM(K+K+1,I-J+K) * DFLOAT(K+J-I+1) ) ) C I1=J1+J2+J3 IF (PARITY(I1).LT.0.D0) GO TO 8 1 I2=J1-J2+J3 IF (I2) 8,2,2 2 I3=J1+J2-J3 IF (I3) 8,3,3 3 I4=-J1+J2+J3 IF (I4) 8,4,4 4 I5=I1/2 I6=I2/2 SIGN=PARITY(I5) 7 THREEJ=SIGN*DELTA(J1,J2,J3)*BINOM(I5,J1)*BINOM(J1,I6) RETURN 8 THREEJ=0D0 RETURN END FUNCTION BINOM(N,M) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C C *** COMPUTATION OF BINOMIAL COEFFICIENTS *** C FN = N+1 NM = N-M MNM = MIN0(NM,M) F = 0.D0 B = 1.D0 IF(MNM) 3,3,1 1 DO 2 I = 1,MNM F = F+1.D0 C = (FN-F)*B 2 B = C/F 3 BINOM = B RETURN END FUNCTION PARITY(I) IMPLICIT DOUBLE PRECISION (A-H,O-Z) PARITY=1.D0 IF (I-2*(I/2).EQ.0) RETURN PARITY=-PARITY RETURN END