PROGRAM CALEGY
C
C    PROGRAM TO GET ENERGY FROM CATALOGUE FREQUENCIES
C                  AND COMPUTE PARTITION FUNCTION
C
      INTEGER NDEGY,NEGY,I,K,KK,NQN,IQFMM,IQFMU,IQFML,NQNP,IHF,
     &        NQNM,NQNX,IDGNL,IDGNU,IEGY,NQN,NTEMP,IQNX,ID,IQNX
      PARAMETER (NDEGY=20000)
      REAL*4 TEMP(57),QX(7),TF,STRMN,STRMX,STR,SONE,Q,QLOG,C,QOFF
      INTEGER IQNL(6),IQNU(6),IPTR(NDEGY),IDGN(NDEGY),IQN(0:6,NDEGY)
      CHARACTER FIL*60,NAM*15
      INTEGER*4 VAL,IVER,NLINE,ITAG
      REAL*8 EGY(NDEGY),EMIN,EMAX,FAC,TMC,CMC,ELOW,EHIGH,FRQ,FMAX
      COMMON /CMEGY/EGY,NEGY,IPTR,IDGN,IQN
      DATA TEMP/ 300.,225.,150.,75.,37.5,18.75,9.375,50*0./
      WRITE(*,*) 'ENTER CATALOG TAG'
      VAL=0
      READ(*,*) VAL
      Q=1.
      OPEN(9,FILE='calegy.lst',STATUS='UNKNOWN')
      CALL CATFIL(VAL,FIL)
      WRITE(9,'(A)') FIL
      CALL CATDIR(VAL,NAM,NLINE,QX,IVER)
      WRITE(*,*) 'file is :',NAM
      WRITE(9,*) 'catalog Q values'
      DO 3 I=1,7
  3       QX(I)=10.**QX(I)
      Q=QX(1)
      WRITE(9,'(4G18.7)') QX
      OPEN(50,FILE=FIL)
      FAC=0.434294481D0
      TMC=1.43878
      CMC=29979.2458D0
      TF=-TMC/300.
      IQNX=0
      NQNX=0
      NLINE=0
      NEGY=0
      STRMN=100.
      STRMX=-100.
      EMAX=0.
      FMAX=0.
      EMIN=1.E20
C    ASSEMBLE ENERGIES
C
  1   READ(50,900,END=10) FRQ,STR,ELOW,IDGNU,ITAG,IQFMU,
     1       IQFMM,IQFML,IQNU,IQNL
 900  FORMAT(F13.4,8X,F8.4,2X,F10.4,I3,I4,4X,3I1,12I2)
      NLINE=NLINE+1
      NQN=IQFML
      NQNP=MOD(IQFMU,5)
      STRMX=MAX(STRMX,STR)
      STRMN=MIN(STRMN,STR)
      IF(NQNP.EQ.NQN)  THEN
          I=1
      ELSE
          I=NQN
      ENDIF
      IHF=1-MOD(IQFMM,2)
      IDGNL=IDGNU*(IQNL(I)*2+IHF)/(IQNU(I)*2+IHF)
      IQNX=MAX(IQNX,IQNU(I),IQNL(I))
      EHIGH=ELOW+FRQ/CMC
      EMIN=MIN(EMIN,ELOW)
      FMAX=MAX(FMAX,FRQ)
      EMAX=MAX(EMAX,EHIGH)
      NQNX=MAX(NQNX,NQN)
      NQNM=MIN(NQNM,NQN)
      CALL EFIND(NQN,IQNU,EHIGH,IDGNU)
      CALL EFIND(NQN,IQNL,ELOW ,IDGNL)
      IF(IQNU(1).EQ.1) THEN
          SONE=24025.12*(10.**STR)/(FRQ*(EXP(TF*ELOW)-
     1             EXP(TF*EHIGH)))
          STR=SONE*Q
          WRITE(9,'(F13.4,2E15.7,12I3)') FRQ,SONE,STR,IQNU,IQNL
      ENDIF
      GO TO 1
 10   WRITE( 9,'(I6,A)') NLINE,' LINES FOUND',IQNX,' = MAX QN'
      WRITE( 9,'(A,2F10.4)') ' MIN MAX EGY',EMIN,EMAX,
     1                     ' MIN MAX STR',STRMN ,STRMX
      WRITE( 9,'(A,F13.4)') ' MAX FRQ',FMAX
C
C   FIND REDUNDANT QUANTUM NUMBERS
C
      IF(NQNM.LT.NQNX) THEN
          DO 20 IEGY=NEGY,2,-1
              KK=IPTR(IEGY)
              K =IPTR(IEGY-1)
              IF(IQN(0,K).GE.IQN(0,KK)) GO TO 20
              NQN=IQN(0,K)
              DO 21 I=1,NQN
                 IF(IQN(I,K).NE.IQN(I,KK)) GO TO 20
 21              CONTINUE
              IDGN(KK)=0
 20       CONTINUE
      ENDIF
C
C   COMPUTE PARTITION FUNCTION
C
      NTEMP=7
C     WRITE(*,'(7F10.3/A)') (TEMP(I),I=1,7),
C    1                    ' ENTER EXTRA TEMPERATURES '
C     DO 50 I=8,57
C         READ(*,*) TEMP(I)
C         IF(TEMP(I).LT.0.01) GO TO 51
C         NTEMP=NTEMP+1
C50       CONTINUE
C51   CONTINUE
      DO 55 I=1,NTEMP
          C=-TMC/TEMP(I)
          QOFF=FAC*C*EMIN
          Q=0.
          DO 52 K=1,NEGY
              KK=IPTR(K)
              ID=IDGN(KK)
              IF(ID.LE.0) GO TO 52
              Q=Q+ID*EXP(C*(EGY(KK)-EMIN))
  52      CONTINUE
          IF(Q.LT.1.E-38) THEN
              QLOG=-40.
          ELSE
              QLOG=LOG10(Q)+QOFF
          ENDIF
          IF(QLOG.GT.38.) THEN
              Q=1.E+38
          ELSE IF (QLOG.LT.-38.) THEN
              Q=1.E-38
          ELSE
              Q=10.**QLOG
          ENDIF
          IF(I.LE.7) THEN
              WRITE( 9,901) TEMP(I),Q,QLOG,QX(I)
          ELSE
              WRITE( 9,901) TEMP(I),Q,QLOG
          ENDIF
 55   CONTINUE
C
C   WRITE ENERGIES
C
C      NEGY=MIN(40,NEGY)
      DO 60 I=1,NEGY
          KK=IPTR(I)
          NQN=IQN(0,KK)
          WRITE( 9,'(F11.4,I5,1H:,6I3)') EGY(KK),IDGN(KK),
     1                               (IQN(K,KK),K=1,NQN)
 60   CONTINUE
      CLOSE(9)
      STOP
 901  FORMAT(F10.3,3G17.7)
      END
C
      SUBROUTINE EFIND(NQN,IQNN,ENOW,IDGNN)
      INTEGER NDEGY,I,NQN,IDGNN,INOW,KK,NQNM,LEFT,RIGHT,ICMP
      PARAMETER (NDEGY=20000)
      INTEGER IQNN(6),IPTR(NDEGY),IDGN(NDEGY),IQN(0:6,NDEGY),NEGY
      REAL*8 ENOW,EGY(NDEGY)
      COMMON /CMEGY/EGY,NEGY,IPTR,IDGN,IQN

      IF(NEGY.GE.NDEGY) RETURN
      LEFT=1
      RIGHT=NEGY
      DO WHILE (LEFT.LE.RIGHT)
          INOW=(LEFT+RIGHT)/2
          KK=IPTR(INOW)
          NQNM=MIN(NQN,IQN(0,KK))
          DO 11 I=1,NQNM
              ICMP=IQN(I,KK)-IQNN(I)
              IF(ICMP.NE.0) GO TO 12
  11          CONTINUE
          ICMP=IQN(0,KK)-NQN
C           RETURN IF ENERGY ALREADY IN LIST
          IF(ICMP.EQ.0) RETURN
  12      IF(ICMP.LT.0) THEN
              LEFT=INOW+1
          ELSE
              RIGHT=INOW-1
          ENDIF
      END DO
      INOW=RIGHT+1
      DO 20 I=NEGY,INOW,-1
20        IPTR(I+1)=IPTR(I)
      NEGY=NEGY+1
      IPTR(INOW)=NEGY
      EGY(NEGY)=ENOW
      IDGN(NEGY)=IDGNN
      DO 30 I=1,6
30        IQN(I,NEGY)=IQNN(I)
      IQN(0,NEGY)=NQN
      RETURN
      END