C     THIS IS A LITTLE DRIVER FOR SUBROUTINE HNCO TO DEMONSTRATE HOW TO 
C       USE IT TO CALCULATE RATES AMONG VARIOUS LEVELS FOR A GIVEN TEMP.
C     NLEVEL IS THE NUMBER OF LEVELS; JLEVEL(2,NLEVEL) DESCRIBES THE
C       LEVELS AS JLEVEL(1,I)=J, JLEVEL(2,I)=K
C       FOR K.NE.O -/+K ARE LOWER/UPPER MEMBER OF DOUBLET
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NLEVEL=27)
      DIMENSION JLEVEL(2,NLEVEL),SIG(NLEVEL,NLEVEL)
      DATA  JLEVEL/0,0, 1,0, 2,0, 3,0, 4,0, 5,0, 6,0, 7,0, 8,0, 9,0,
     1  10,0, 11,0, 12,0, 1,-1, 1,1, 2,-1, 2,1, 3,-1, 3,1, 4,-1, 4,1,
     2         5,-1, 5,1, 6,-1, 6,1, 7,-1, 7,1/
      TEMP=80.
      WRITE(6,*)' ***** BELOW FOR TEMPERATURE',TEMP
      WRITE(6,*)
      CALL HNCO(NLEVEL,JLEVEL,SIG,TEMP)
      DO 5100 I=1,NLEVEL
      WRITE(6,654) I,JLEVEL(1,I),JLEVEL(2,I)
  654 FORMAT('0 INITIAL LEVEL =',I4,'      J, K =',2I4)
 5100 WRITE(6,644) (IF,SIG(I,IF),IF=1,NLEVEL)
  644 FORMAT(6(4X,I3,1P,D12.4,1X))
      STOP
      END
      SUBROUTINE HNCO(NLEV,JLEV,SIG,TEMP)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION INTEG,KBOLTZ
      DIMENSION JLEV(2,NLEV),SIG(NLEV,NLEV)
C
C     PROGRAM TO COMPUTE MATRIX OF RATE CONSTANTS FOR HNCO-HE
C        USING IOS-SCALING.
C
C     PROGRAM CORRESPONDS TO NASA TECH. MEMO. TM 87791
C
C     PROGRAM DEVELOPED  FEBRUARY 1986 BY:
C        SHELDON GREEN,
C        NASA GODDARD INSTITUTE FOR SPACE STUDIES,
C        2880 BROADWAY,
C        NEW YORK, NY 10025
C        AGXSG@GISS.NASA.GOV
C     PLEASE CONTACT AUTHOR FOR QUESTIONS OR PROBLEMS.
C
C     PROGRAM INPUTS L,M1,M2,(QL(IE),IE=1,NNRG) FROM DATA SET HNCO.DAT;
C     THESE ARE QL IN ANG**2 FROM IOS CALCS AT VARIOUS ENERGIES,
C         ENERGY(IE), IE=1,NNRG.
C     DATA CARDS ARE READ ONLY ON THE FIRST CALL TO SUBROUTINE.
C
C     STATE-TO-STATE RATES IN CM**3/SEC RETURNED IN SIG(I,F) FOR
C       LEVELS DESCRIBED IN INPUT JLEV(1,I)=J
C                                 JLEV(2,I)=K,  I=1,NLEV
C       NOTE:        (K.LT.0 FOR (K)-(-K), K.GT.0 FOR (K)+(-K))
C
C     DIMENSIONS FOR UP TO MXQL L,M1,M2-VALUES AND MXNRG ENERGIES
C       THESE SHOULD NOT REQUIRE CHANGING:
      PARAMETER (MXQL=165, MXNRG=9)
      DIMENSION RATES(MXQL),QL(MXNRG,MXQL),ENERGY(MXNRG)
     1     ,RIU(MXQL),RID(MXQL),RIA(MXQL),EVAL(MXQL)
      DIMENSION LMM(3,MXQL)
C
C     FOLLOWING DIMENSIONS FOR STATE-TO-STATE RATES MUST BE CHANGED
C       IF NUMBER OF LEVELS (NLEVEL) .GT. 100 IS DESIRED.
      PARAMETER (MXLVL=100)
      DIMENSION ELEV(MXLVL)
      CHARACTER*1 S(MXLVL)
C
C     ------------------------------------------------------------------
C     THE LOGICAL VARIABLES BELOW CONTROL INTERPRETATION OF ENERGY
C       IN THE IOS METHOD.  CURRENT SETTINGS APPEAR BEST CF CS CALCS.
C
C     LOGICAL LAVG = .TRUE. GETS (GEOMETRIC) AVERAGE OF UP/DOWN IOS
      LOGICAL LAVG
C     LOGICAL LDOWN = .TRUE. GETS UPWARD RATES VIA DOWNWARD PLUS D.B.
C             LDOWN=.FALSE. DOES JUST THE OPPOSITE.
      LOGICAL LDOWN
      DATA LAVG/.FALSE./, LDOWN/.FALSE./
C     THE FOLLOWING CONTROLS PRINT OUTPUT FROM THIS ROUTINE
C         LPRT/.TRUE./  PRINTS ALL RELEVANT INFO. DURING RUN
C         LPRT/.FALSE./ PRINTS ONLY WARNING/ERROR MESSAGES.
      LOGICAL LPRT
      DATA LPRT/.FALSE./
C     ------------------------------------------------------------------
C     DATA FOR HNCO-HE CONSTANTS, AND TO LIMIT INPUT TO 1ST CALL.
      DATA A,B,C/.3693,.3639,30.6381/, URED/3.66/
      DATA INIT/-1/
C     DATA FOR INPUT ENERGIES OF QLM
      DATA ENERGY/15.,30.,50.,75.,100.,150.,200.,400.,600./
C     DATA FOR RECOMMENDED MINIMUM/MAXIMUM TEMPERATURES
      DATA TMPMIN/30./,TMPMAX/250./
C     STATEMENT FUNCTION  . . .
      F(I)=FLOAT(2*I+1)
C
C     SET DEFAULT VALUES BEFORE READ
      NNRG=MXNRG
      NQL=MXQL
C
      IF (LPRT) WRITE(6,600) URED,TEMP
  600 FORMAT(' PROGRAM TO COMPUTE ADIABATICALLY-CORRECTED-SUDDEN',
     1       ' MATRIX OF RATE COEFFICIENTS FROM IOS Q(L,M1,M2)'/
     2       /' HNCO - HE ELECTRON GAS POTENTIAL'/
     3       /' URED =',F12.5,' A.M.U.'/
     6       /' TEMPERATURE REQUESTED =',F8.2)
      IF (LPRT.AND.LAVG) WRITE(6,639)
  639 FORMAT(/' * * * GEOMETRIC AVERAGE OF INITIAL UP/DOWN RATES',
     &       ' WILL BE USED IN IOS SCALING')
      IF (LPRT.AND.LDOWN) WRITE(6,638)
  638 FORMAT(/' * * * SCALING USED FOR DOWNWARD RATES; UP RATES VIA',
     &       ' DETAILED BALANCE.')
      IF (LPRT.AND.(.NOT.LDOWN)) WRITE(6,648)
  648 FORMAT(/' * * * NOTE.  LDOWN=FALSE.  SCALING FOR UP RATES; DOWN',
     &       ' RATES VIA DETAILED BALANCE.')
C
C     INPUT L,M1,M2,(QL(IE),IE=1,NNRG), IN ANGSTROM**2,  FROM CARDS
      IF (INIT.GE.0) GO TO 1999
      DO 1033 I=1,NQL
 1033 READ(5,599,END=1998) (LMM(J,I),J=1,3),(QL(IE,I),IE=1,NNRG)
  599 FORMAT(3I2,9E8.2)
      INIT=1
      GO TO 1999
 1998 WRITE(6,637) NQL,I
  637 FORMAT('0 END OF DATA BEFORE REQUESTED NQL',2I4)
      STOP
C
 1999 IF (LPRT) WRITE(6,658) NQL
  658 FORMAT(/' INPUT COMPLETED FOR',I4,'  QL(L,M1,M2)')
      IF (NLEV.LE.MXLVL) GO TO 1997
      WRITE(6,642) NLEV,MXLVL
  642 FORMAT('0 *** ERROR. NLEV REDUCED TO MXLVL',2I8)
      NLEV=MXLVL
C     GET 'ENERGY' EVAL ASSOCIATED WITH L,M1,M2
 1997 DO 9100 I=1,NQL
      EVAL(I)=0.25*(ENRG(LMM(1,I),LMM(2,I),A,B,C)+
     1    ENRG(LMM(1,I),-LMM(2,I),A,B,C)+ENRG(LMM(1,I),LMM(3,I),A,B,C)
     2           +ENRG(LMM(1,I),-LMM(3,I),A,B,C) )
 9100 CONTINUE
      IF (LPRT) WRITE(6,623) A,B,C,
     &            (I,(LMM(J,I),J=1,3),EVAL(I),I=1,NQL)
  623 FORMAT(/' ***'/' *** NEAR SYMMETRIC TOP ENERGIES FROM A, B, C=',
     1 3F10.4   /   ' ***'/
     2       /'  LEVEL   L  M1  M2     ENERGY'/
     3       (I7,3I4,F10.4))
C
C     GET RATES AS BOLTZMANN AVERAGE OF QL INPUT
C     CONVERSION FACTORS TO GIVE RATES IN CM**3/SEC
      KBOLTZ=1.3805E-16
      AVOGAD=6.024E23
      PI=3.1415926
      ERGPCM=KBOLTZ*300./208.4
      ERGPCM=ERGPCM*ERGPCM
      FACT=(PI*URED/AVOGAD)*8.*KBOLTZ**3
      TFACT=8.E-16*ERGPCM/SQRT(FACT*TEMP**3)
      SFACT=(208.4/300.)*TEMP
      FKT=300./(208.4*TEMP)
      IF (TEMP.GE.TMPMIN.OR.TEMP.LE.TMPMAX)  GO TO 2001
      WRITE(6,603) TMPMIN,TMPMAX
  603 FORMAT(/' * * * ERROR.  TEMP NOT IN RECOMMENDED RANGE',2F10.2/
     1        ' * * *        PROGRAM CAN BE FORCED BY OVERRIDING',
     2        ' TMPMIN OR TMPMAX')
      STOP
 2001 IF (LPRT) WRITE(6,604) NQL,NNRG,TEMP
  604 FORMAT(//' PROCESSING OF',I4,' QL AND',I4,' ENERGIES FOR'
     &       ,' BOLTZMANN AVERAGE AT TEMPERATURE =',F8.2)
C     ORDER QL ON INCREASING ENERGY
      CALL EORDR(MXNRG,NNRG,ENERGY,NQL,QL)
C
      DO 2000 L=1,NQL
      IF (LPRT)
     & WRITE(6,605) (LMM(J,L),J=1,3),(ENERGY(IE),QL(IE,L),IE=1,NNRG)
  605 FORMAT(/'  FOR L, M1, M2 =',3I3,'   ENERGY   QL(ANGSTROMS)',
     &' ARE'/(25X,F10.2,F10.4)  )
      THRESH=0.
      RAVG=INTEG(NNRG,ENERGY,QL(1,L),THRESH,TEMP,ERR)
      RATES(L)=TFACT*RAVG
      RAVG=RAVG/(SFACT*SFACT)
      IF (LPRT)  WRITE(6,606)  RAVG,RATES(L)
  606 FORMAT(/'  AVERAGE IS',F12.4,'  ANG**2  =',1PE12.4,'  CM**3/SEC')
      IF ((ERR.GT.0.1).AND. LPRT)  WRITE(6,696)  ERR
  696 FORMAT('  * * * WARNING.  FRACTION IN HIGH ENERGY',
     &       ' EXTRAPOLATION =',2F8.3)
C ***
C     GET RIU AND RID (INITIAL ENERGY UPWARD, INITIAL ENERGY DOWNWARD)
C       AS DEFINED BY CHAPMAN AND GREEN CHEM PHYS LETT 112, 436 (1984).
      RIU(L)=RATES(L)
      RID(L)=RIU(L)*EXP(-EVAL(L)*FKT)
      RIA(L)=SQRT(RIU(L)*RID(L))
      IF (RIU(L).LT.0.) RIA(L)=-RIA(L)
      IF (LAVG) RATES(L)=RIA(L)
 2000 CONTINUE
C
C     OUTPUT RIU, RID
      IF (LPRT)
     &WRITE(6,677) TEMP,((LMM(J,I),J=1,3),RIU(I),RIA(I),RID(I),I=1,NQL)
  677 FORMAT(//' CHAPMAN-GREEN INITIAL UP, AVERAGE, AND DOWN  RATES IN',
     1       ' CM**3/SEC'//' TEMPERATURE =',F10.2/
     2       /'   L  M1  M2     R(IU)         R(IA)         R(ID)'/
     2       (3I4,1P,3D14.4) )
C
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C                                                                   *
C     AT THIS POINT  CALCULATE MATRIX OF STATE-TO-STATE RATES       *
C             WITH OR WITHOUT ADIABATIC CORRECTION.                 *
C     N.B. RATES() IS IN SAME ORDER AS QL, SO USE IXQL ON LMM TO GET*
C                                                                   *
      IF (LPRT) WRITE(6,652)  TEMP
  652 FORMAT(//' RATES WILL BE COMPUTED AMONG FOLLOWING LEVELS FOR TEMP'
     &      ,' =',F8.2/' LEVEL J   K   ENERGY')
      IMSG=0
      IMSGX=0
C      CALL IXQLFI(LMM,NQL)
      DO 5002 I=1,NLEV
      ELEV(I)=ENRG(JLEV(1,I),JLEV(2,I),A,B,C)
 5002 IF (LPRT) WRITE(6,653) I,JLEV(1,I),JLEV(2,I),ELEV(I)
  653 FORMAT(3I4,F12.4)
      IF (LPRT.AND.LAVG) WRITE(6,639)
      IF (LPRT.AND.LDOWN) WRITE(6,638)
      IF (LPRT.AND.(.NOT.LDOWN)) WRITE(6,648)
C
      DO 5100 I=1,NLEV
      JI=JLEV(1,I)
      KI=IABS(JLEV(2,I))
      EPSI=1.
      IF (JLEV(2,I).LT.0) EPSI=-1.
      IF (KI.EQ.0)  EPSI=0.
      IF (LPRT) WRITE(6,654) I,JI,KI,EPSI
  654 FORMAT(/' INITIAL LEVEL =',I4,'      J, K, PARITY =',2I4,F6.1)
      DO 5101 IF=1,NLEV
      SIG(I,IF)=0.
      IF (IF.EQ.I) GO TO 5101
      JF=JLEV(1,IF)
      KF=IABS(JLEV(2,IF))
      EPSF=1.
      IF (JLEV(2,IF).LT.0) EPSF=-1.
      IF (KF.EQ.0)  EPSF=0.
      DELTAE=ELEV(IF)-ELEV(I)
      IF (.NOT.LDOWN) GO TO 5191
C     LDOWN = .TRUE. CASE; USE IOS SCALING ONLY ON DOWNWARD TRANSITIONS
      IF (DELTAE.LE.0.)  GO TO 5291
      CALL SIG5(JF,KF,EPSF,JI,KI,EPSI,SIG(I,IF),S(IF),IMSG,RATES,
     1          LMM,NQL)
      SIG(I,IF)=SIG(I,IF)*F(JF)*EXP(-DELTAE*FKT)/F(JI)
      GO TO 5102
 5291 CALL SIG5(JI,KI,EPSI,JF,KF,EPSF,SIG(I,IF),S(IF),IMSG,RATES,
     1          LMM,NQL)
      GO TO 5102
C     BELOW FOR LDOWN=.FALSE.; CALC UP RATES DIRECTLY, DOWN VIA D.B.
 5191 IF (DELTAE.GE.0.)  GO TO 5391
      CALL SIG5(JF,KF,EPSF,JI,KI,EPSI,SIG(I,IF),S(IF),IMSG,RATES,
     1          LMM,NQL)
      SIG(I,IF)=SIG(I,IF)*F(JF)*EXP(-DELTAE*FKT)/F(JI)
      GO TO 5102
 5391 CALL SIG5(JI,KI,EPSI,JF,KF,EPSF,SIG(I,IF),S(IF),IMSG,RATES,
     1          LMM,NQL)
 5102 IF (IMSG.LE.0) GO TO 5101
      IMSGX=IMSGX+1
      WRITE(6,692) I,IF,IMSG
  692 FORMAT(' WARNING. FOR LEVEL',I3,' TO',I3,',  MISSING',I4,
     1      '  L,M1,M2 CONTRIBUTIONS TO RATE.')
      IMSG=0
 5101 CONTINUE
 5100 IF (LPRT) WRITE(6,644) (IF,SIG(I,IF),S(IF),IF=1,NLEV)
  644 FORMAT(6(4X,I3,1PE12.4,A1))
      IF (LPRT.AND.IMSGX.GT.0)  WRITE(6,699)
  699 FORMAT(    '0 ***** NOTE.  FOR CROSS SECTIONS MARKED WITH A STAR,
     1ALL CONTRIBUTING Q(L) ARE NOT AVAILABLE.')
C * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
 9999 RETURN
      END
      SUBROUTINE SIG5(JI,KI,EPSI,JF,KF,EPSF,SIG,S,IMSG,RATES,LMM,NQL)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION RATES(2),LMM(3,NQL)
      CHARACTER*1 S
      CHARACTER*1 BLANK/' '/,STAR/'*'/
      DATA ZTOL/1.E-8/
C
C     STATEMENT FUNCTIONS. . .
      F(J)=FLOAT(2*J+1)
      XNORM(EPS)=1./(1.+ABS(EPS))
C
      XJI=JI
      XKI=KI
      XJF=JF
      XKF=KF
      XNI=XNORM(EPSI)
      XNF=XNORM(EPSF)
      LLO=IABS(JI-JF)
      LHI=JI+JF
      PJK=PARITY(JI+JF+KI+KF)
      MPLS=KI+KF
      MMIN=IABS(KI-KF)
      P2=1.
      IF (KI-KF.LT.0)  P2=PARITY(MMIN)
      SIG=0.
      S=BLANK
      TMAX=0.
      DO 5102 L=LLO,LHI
      XL=L
      PL=PJK*PARITY(L)
C     -----------------------TERM 1 -------------------
      PP=1.+EPSI*EPSF*PL
      PP=PP*PP
      IF (PP.LE.ZTOL)  GO TO 5200
      TJ=THRJ(XJF,XL,XJI,XKF,XKI-XKF,-XKI)
      TJ=TJ*TJ
      IF (TJ.LE.ZTOL)  GO TO 5200
      CALL IXQLF(L,MMIN,MMIN,0,INDEX,LMM,NQL)
      IF (INDEX.GT.0)  GO TO 5110
      IF (INDEX.EQ.-1) GO TO 5200
      IMSG=IMSG+1
      S=STAR
      GO TO 5200
 5110 TT=PP*TJ*RATES(INDEX)
      TMAX=DMAX1(ABS(TT),TMAX)
      SIG=SIG+TT
C     -----------------------TERM 2 -------------------
 5200 PP=(1.+EPSI*EPSF*PL)*(EPSF+EPSI*PL)
      IF (ABS(PP).LE.ZTOL)  GO TO 5300
      TJ=THRJ(XJF,XL,XJI,XKF,XKI-XKF,-XKI)*
     &             THRJ(XJF,XL,XJI,-XKF,XKF+XKI,-XKI)
      IF (ABS(TJ).LE.ZTOL)  GO TO 5300
      CALL IXQLF(L,MPLS,MMIN,1,INDEX,LMM,NQL)
      IF (INDEX.GT.0)  GO TO 5210
      IF (INDEX.EQ.-1)  GO TO 5300
      IMSG=IMSG+1
      S=STAR
C     ON HIGH PRNTLV WRITE MSG
      GO TO 5300
 5210 TT=2.*P2*PP*TJ*RATES(INDEX)
      TMAX=DMAX1(TMAX,ABS(TT))
      SIG=SIG+TT
C     -----------------------TERM 3 -------------------
 5300 PP=EPSF+EPSI*PL
      PP=PP*PP
      IF (PP.LE.ZTOL)  GO TO 5102
      TJ=THRJ(XJF,XL,XJI,-XKF,XKF+XKI,-XKI)
      TJ=TJ*TJ
       IF (TJ.LE.ZTOL)  GO TO 5102
      CALL IXQLF(L,MPLS,MPLS,0,INDEX,LMM,NQL)
      IF (INDEX.GT.0)  GO TO 5310
      IF (INDEX.EQ.-1)  GO TO 5102
      S=STAR
      IMSG=IMSG+1
C     ON 'PRNTLV' WRITEN  MSG
      GO TO 5102
 5310 TT=PP*TJ*RATES(INDEX)
      TMAX=DMAX1(ABS(TT),TMAX)
      SIG=SIG+TT
 5102 CONTINUE
      IF (SIG.GE.ZTOL*TMAX)  GO TO 5101
      IF (SIG.EQ.0.)  GO TO 5101
      SIG=0.
 5101 SIG=SIG*XNI*XNF*F(JF)
      RETURN
      END
      FUNCTION THRJ(F1,F2,F3,G1,G2,G3)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      SAVE MUNG,X,Y
      DIMENSION X(302),Y(302)
      DATA MUNG/0/,MXIX/302/
      IF (MUNG.EQ.21) GO TO 69
      MUNG = 21
      X(1) = 0.D0
      DO 100 I = 1, 201
      A = I
      X(I+1) = LOG(A) +X(I)
      Y(I+1) = LOG(A)
  100 CONTINUE
   69 IF(F1-ABS(G1)) 1,13,13
   13 IF(F2-ABS(G2))1,14,14
   14 IF(F3-ABS(G3))1,15,15
   15 SUM=F1+F2+F3
      NSUM=SUM+.001D0
      IF(SUM-NSUM)2,2,1
    1 THRJ=0.D0
      RETURN
    2 IF(ABS(G1+G2+G3)-1.D-08)3,3,1
    3 IF(F1+F2-F3)1,4,4
    4 IF(F1+F3-F2)1,5,5
    5 IF(F2+F3-F1)1,6,6
    6 J1=2.D0*F3+2.001D0
      J2=F1+F2-F3+1.001D0
      J3=F1-F2+F3+1.001D0
      J4=-F1+F2+F3+1.001D0
      J5=F1+F2+F3+2.001D0
      J6=F1+G1+1.001D0
      J7=F1-G1+1.001D0
      J8=F2+G2+1.001D0
      J9=F2-G2+1.001D0
      J10=F3+G3+1.001D0
      J11=F3-G3+1.001D0
      IF(J5.GT.MXIX) THEN
        WRITE(6,601) J5,MXIX
  601   FORMAT(' *** DIMENSION ERROR IN THRJ - INDEX.GT.MXIX',2I5)
        STOP
      ENDIF
      R=0.5D0*(Y(J1)+X(J2)+X(J3)+X(J4)-X(J5)
     1+X(J6)+X(J7)+X(J8)+X(J9)+X(J10)+X(J11))
      SUM=0.D0
      F=-1
      KZ=-1
    7 KZ=KZ+1
      F=-F
      J1=KZ+1
      J2=F1+F2-F3-KZ+1.001D0
      IF(J2)20,20,8
    8 J3=F1-G1-KZ+1.001D0
      IF(J3)20,20,9
    9 J4=F2+G2-KZ+1.001D0
      IF(J4)20,20,10
   10 J5=F3-F2+G1+KZ+1.001D0
      IF(J5)7,7,11
   11 J6=F3-F1-G2+KZ+1.001D0
      IF(J6)7,7,12
   12 JMAX=MAX(J1,J2,J3,J4,J5,J6)
      IF(JMAX.GT.MXIX) THEN
        WRITE(6,601) JMAX,MXIX
        STOP
      ENDIF
      S=-(X(J1)+X(J2)+X(J3)+X(J4)+X(J5)+X(J6))
      SUM=SUM+F*EXP(R+S)
      GO TO 7
   20 INT=ABS(F1-F2-G3)+0.0001D0
      VAL=((-1.D0)**INT)*SUM/SQRT(2.D0*F3+1.D0)
      IF(ABS(VAL).LE.1.D-6) VAL=0.D0
      THRJ=VAL
      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
      SUBROUTINE IXQLF(L,M1,M2,IC,IX,LMM,MX)
      DIMENSION LMM(3,MX)
      IF (IC.EQ.0 .OR. IC.EQ.1) GOTO 1000
      IX=-1
      RETURN
 1000 IX=0
      DO 2000 I=1,MX
      IF (LMM(1,I).NE.L) GO TO 2000
      IF (LMM(2,I).NE.M1) GO TO 2000
      IF (LMM(3,I).NE.M2) GO TO 2000
      IX=I
      RETURN
 2000 CONTINUE
      RETURN
      END
      SUBROUTINE EORDR(MXNRG,NNRG,E,NQL,QL)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C     ROUTINE TO ORDER QL VALUES ON ENERGY
      DIMENSION QL(MXNRG,NQL),E(MXNRG)
      IF (NNRG.LE.1) RETURN
      NM1=NNRG-1
      DO 1000 I=1,NM1
      IST=I+1
      DO 1000 J=IST,NNRG
      IF (E(I).LE.E(J)) GO TO 1000
      T=E(I)
      E(I)=E(J)
      E(J)=T
      DO 1001 L=1,NQL
      T=QL(I,L)
      QL(I,L)=QL(J,L)
 1001 QL(J,L)=T
 1000 CONTINUE
      RETURN
      END
      FUNCTION ENRG(J,K,A,B,C)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
C     NEAR SYMMETRIC TOP ENERGY LEVELS
C       K.GE.0 IS SYMMETRIC COMBINATION
C       K.LT.0 IS ANTI-SYMMETRIC COMBINATION
      HKK=(A+B)*0.5*(J*(J+1)-K*K)+C*(K*K)
      IF (IABS(K).NE.1) GO TO 1000
      SS=1.
      IF (K.LT.0)  SS=-1.
      HKK=HKK+SS*(A-B)*SQRT(FLOAT((J*(J+1)-K*(K-1))*
     &                           (J*(J+1)-(K-1)*(K-2))))
 1000 ENRG=HKK
      RETURN
      END
      FUNCTION INTEG(N,E,SIG,TH,T,ERR)
C
C     BOLTZMAN INTEGRATOR FOR SET OF ENERGIES IN 1/CM, E, AND CROSS
C       SECTIONS IN ANGSTROMS, SIG.  N IS NUMBER OF E, SIG VALUES.
C     T IS TEMPERATURE IN KELVINS, TH IS THRESHOLD ENERGY IN 1/CM
C     ERR RETURNS THE FRACTION FROM EXTRAPOLATION AT HIGH E.
C
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION INTEG,KBT
      DIMENSION  E(N),SIG(N)
C
C     CALCULATES NUMERICALLY INTEGRAL OF DE*E*EXP(-E/KT)*SIG(E)
C       IN EACH INTERVAL FITS LINE TO SIG(E)=B+SLOPE*E
C       AND ESTIMATES  CONTRIBUTION OF HIGH-ENERGY TAIL.
C
      ERR=0.
C  E IS IN 1/CM.  SIG IS IN A**2.
      KBT=208.4*T/300.
      XINT=0.
      EL=TH
      XL=-EL/KBT
      WL=DEXP(XL)
      SL=0.
      DO 1000 I=1,N
      ER=E(I)
      DE=ER-EL
C     SKIP UNTIL ABOVE THRESHOLD.
C     IF SEVERAL CARDS HAVE SAME ENERFY, THIS AVERAGES FIRST AND LAST.
      IF (DE.LE.0.0)  GO TO 1000
      SR=SIG(I)
      XR=-ER/KBT
      WR=DEXP(XR)
C     FIT STRAIGHT LINE TO (EL,SL), (ER,SR). . .
      SLOPE=(SR-SL)/(ER-EL)
      B=SL-EL*SLOPE
C     ADD INTEGRAL OVER THIS SEGMENT
      XINT=XINT+(KBT*KBT)*(B*(WR*(XR-1.)-WL*(XL-1.))
     1       -KBT*SLOPE*(WR*(XR*XR-2.*XR+2.)-WL*(XL*XL-2.*XL+2.)))
C  SET UP LEFT VALUES FOR NEXT POINT.
      EL=ER
      XL=XR
      WL=WR
      SL=SR
 1000 CONTINUE
C  COMPUT TAIL BY FITTING EXPONENTIAL TO LAST TWO POINTS.
      DEL=0.
      NM1=N-1
C     CHECK FOR DUPLICATE LAST ENERGIES.
 2300 IF (E(NM1).LT.E(N))  GO TO 2200
      NM1=NM1-1
      GO TO 2300
 2200 RAT=SIG(NM1)/SIG(N)
      IF (RAT.GE.1.) GO TO 2000
C   IF SIG(E) IS INCREASING, USE CONSTANT = AVERAGE OF N-1,N.
C     N.B. ONE COULD ALSO FIT A STRAIGHT LINE SIG(E)=A*E+B IN THIS CASE.
 2100 A=0.5*(SIG(NM1)+SIG(N))
      B=0.
      GO TO 2500
 2000 B=DLOG(RAT)/(E(N)-E(NM1))
      A = SIG(N)*DEXP( B*E(N) )
 2500 B=B+1./KBT
      DEL=A*(1.+E(N)*B)*DEXP(-E(N)*B)/(B*B)
      ERR=DEL/(XINT+DEL)
 3000 INTEG=XINT+DEL
      RETURN
      END