CFTN77,Q,T,S
C$CDS ON
C$FILES(0,5)
C$EMA//
C$MSEG 6
C directives above used for HP100
C     INCLUDE 'CALPGM.INC'
      PROGRAM  CALCAT
C**********************************************************************
C
C   Copyright (C) 1989, California Institute of Technology
C   All rights reserved.  U. S. Government Sponsorship under
C   NASA Contract NAS7-918 is acknowledged.
C
C   Herbert M. Pickett, 20 Mar 1989
C
C   HMP, code added for binary read of variance
C   HMP, constained derivatives calculated within HAMX
C   HMP, IDPAR is now REAL*8
C   HMP, identical dipoles are summed on input
C   HMP, str output now will put out component dipoles
C   HMP, large uncertainties do not inhibit output to .CAT file
C
C THIS IS A GENERALIZED INTENSITY AND FREQUENCY CALCULATOR
C
C THE HAMILTONIAN IS ASSUMED TO BE ORGANIZED INTO BLOCKS
C     SUCH THAT ADJACENT SETS OF BLOCKS ARE CONNECTED BY TRANSITIONS
C
C***********************************************************************
      INTEGER   NDPAR,NTEMP,NDEGY,NBSAVX,NXDIP,LUINT,LUVAR,LUOUT,LUBIN,
     :          LUEGY,LUSTR,LUCAT,IFLG,NDIP,IQNFMT,NPAR,NPAR1,JNXT,
     :          INBLK,LBLK,I,ITD,IBLK,J,JBLK,NXT,NSIZD,IDGN,NBSAVM,
     :          NQN,NBKPJ,NSAV,NDM,ISIZ,IDMY,PCARD,SETOPT,SETBLK,KBGN,
     :          NPARTOT,K,JJ,JX,JSIZ,INTENS,IDF,JDGN,IJ,IGUP,NBSAV,
     :          IDAMAX,NPDIP
      REAL*8 ZERO,QROT,THRSH,THRSH1,THRSHF,FQMAX,FQMIN,FAC,CMC,DGNF,
     :          TMQ,STCOMP,STR,STRR,STRQ,TMC,TMQL,STARG,SCOMP,STRMN,ERR,
     :          CALERR,EGX,TE,DGN,FRQ,ELOW,STRLG,QLOG,EGYMIN,ONE
      INTEGER*4 HEAPLN,ISIZB,IOFF,IXDP,JXDP,ITAG,NLINE,VAR,S,NKNT,KNT,
     :          LDEDP,LEIG,IBUFOF,SIJ,SWK,NS
C standard dimensions:
      PARAMETER (NDPAR=999,NTEMP=7,NDEGY=450,NBSAVX=256,HEAPLN=560000)
C DOS-PC(16 bit compiler) dimensions:
C     PARAMETER (NDPAR=180,NTEMP=7,NDEGY=181,NBSAVX=99,HEAPLN=55000)
      PARAMETER (NXDIP=511)
      INTEGER   NFILE,LINT,LVAR,LOUT,LCAT,LSTR,LEGY,LBIN
      PARAMETER (NFILE=7,LINT=1,LVAR=2,LOUT=3,LCAT=4,LSTR=5,LEGY=6,
     :           LBIN=7)
      INTEGER*4 IDIP(0:NXDIP),IDIP1(0:1),NSBLK(NBSAVX)
      INTEGER   IQNI(12),NBLK(NBSAVX),IXBLK(NBSAVX)
      CHARACTER*2 IQN(12)
      CHARACTER*80 TITL
      CHARACTER*64 FNAME(NFILE)
      CHARACTER*4  EXT(NFILE)
      INTEGER IFDIAG,RQEXIT,GETPAR,I1,I0,ITMP,IFDUMP
      EXTERNAL      BRKQR
      INTEGER*2     CLEN5
      PARAMETER (I0=0,I1=1,CLEN5=5)
      LOGICAL   FIRST,PRDER,PREGY,PRFRQ,PREIG,PRSTR,PRIR,DIAG,DIAGD
      REAL*8  IDPAR(0:NDPAR),PAR(NDPAR),DERV(NDPAR),EGY(NDEGY),
     :          DIP(NXDIP),TEMP(NTEMP),QSUM(NTEMP),BIGERR
C      EMA  PAR,DERV,EGY,DIP,IDPAR,IDIP,NBLK,IXBLK
      REAL*8 HEAP(0:HEAPLN)
      COMMON HEAP
      SAVE
      DATA TEMP /300.D0,225.D0,150.D0,75.D0,37.5D0,18.75D0,9.375D0/
      DATA QSUM /7*0.D0/
      DATA IQN/12*'  '/
C
C >>>>>> BEGIN PROGRAM
C
C  following statement is for PRIME computer to enable ^P hot key   
C      CALL MKON$P('QUIT$',CLEN5,BRKQR)
      ONE=1
      ZERO=1.5E-38
      BIGERR=999.9999D0
      LUINT=LINT
      LUVAR=LVAR
      LUOUT=LOUT
      LUCAT=LCAT
      LUSTR=LSTR
      LUEGY=LEGY
      LUBIN=LBIN
      EXT(LINT)='int'
      EXT(LVAR)='var'
      EXT(LOUT)='out'
      EXT(LCAT)='cat'
      EXT(LSTR)='str'
      EXT(LEGY)='egy'
      EXT(LBIN)='unf'
      CALL FILGET(FNAME,EXT,NFILE)
      CALL OPENLU(LUOUT,FNAME(LOUT),'WP')
      CALL OPENLU(LUINT,FNAME(LINT),'R')
      NLINE=0
      EGYMIN=0.
C
C input read and echo
C
      READ (LUINT,430,END=425) TITL
      CALL CHTIME(TITL(51:80))
      WRITE (LUOUT,431) TITL
      WRITE (*,431) TITL
      PAR(1)=0
      PAR(2)=0
      PAR(3)=1000.
      PAR(4)=0
      PAR(5)=0
      PAR(6)=-100.0
      PAR(7)=-100.0
      PAR(8)=9999.99
      PAR(9)=300.
      I=9
      READ (LUINT,430,END=425) TITL
      IF(PCARD(TITL,PAR,I).EQ.0) GO TO 425
      IFLG=  PAR(1)
      ITAG=  PAR(2)
      QROT=  PAR(3)
      INBLK= PAR(4)
      LBLK=  PAR(5)
      THRSH= PAR(6)
      THRSH1=PAR(7)
      FQMAX= PAR(8)
      TMQ=   PAR(9)
      WRITE (LUOUT,450) ITAG,QROT,INBLK,LBLK,THRSH
     :              ,THRSH1,FQMAX,TMQ
      FAC=4.16231E-05/QROT
      CMC=29979.2458D0
      PRIR=IFLG.GE.1000
C     initialize intensities for infrared or microwave cases
      IF(PRIR) THEN
          FAC=FAC*CMC
          CMC=1
          IFLG=MOD(IFLG,1000)
          FQMIN=0.0000005
      ELSE
          FQMAX=FQMAX*1000.
          FQMIN=0.00005
      ENDIF
C   set up flags for different output options
      PRFRQ=IFLG.GE.100
      IFLG=MOD(IFLG,100)
      PRSTR=IFLG.GE.10
      STCOMP=1./ZERO
      NPDIP=1
      IF(PRSTR) THEN
          STCOMP=0.0000001
          IF(IFLG.GE.20) NPDIP=-1
          IFLG=MOD(IFLG,10)
      ENDIF
      IFDUMP=IFLG/5
      PRDER=IFLG.EQ.2.OR.IFLG.EQ.4
      PREIG=IFLG.GE.3
      PREGY=IFLG.NE.0
C  read dipole moments
      J=1
 10   IF(J.LT.NXDIP) THEN
          READ(LUINT,*,END=15,ERR=15) IDIP(J),DIP(J)
          DO 11 I=1,NDIP
              IF(IDIP(I).EQ.IDIP(J)) THEN
                  DIP(I)=DIP(I)+DIP(J)
                  WRITE(LUOUT,440) I,IDIP(I),DIP(I)
                  GO TO 10
              ENDIF
 11       CONTINUE
          WRITE(LUOUT,440) J,IDIP(J),DIP(J)
          NDIP=J
          J=J+1
          GO TO 10
      ENDIF
 15   IDIP(0)=NDIP
      IF(NPDIP.LE.0) NPDIP=NDIP
C     close .INT file
      CALL OPENLU(LUINT,FNAME(LINT),'C')
      CALL OPENLU(LUVAR,FNAME(LVAR),'R')
      CALL OPENLU(LUBIN,FNAME(LBIN),'RB')
      IF(LUBIN.EQ.0) LUBIN=-1
C read var file
      READ(LUVAR,430) TITL
      WRITE(LUOUT,431) '.VAR FILE TITLE :'//TITL
      NPARTOT=NDPAR
      READ(LUVAR,*) NPARTOT
C read option card
      I=SETOPT(LUVAR,IQNFMT,ITD,TITL)
      ITMP=10
      NQN=MOD(IQNFMT,ITMP)
C read parameters and variance
      IF(NPARTOT.GT.NDPAR) THEN
          WRITE(*,*) ' Warning: number of parameters less'
     :             , ' than requested',NPARTOT,NDPAR
          NPARTOT=NDPAR
      ENDIF
      VAR=HEAPLN
      NPAR=NPARTOT
      ITMP=GETPAR(LUVAR,VAR,HEAP,DERV,NPAR,IDPAR,PAR,LUBIN,LUOUT)
      NPARTOT=IDPAR(0)
C     close .VAR file
      CALL OPENLU(LUVAR,FNAME(LVAR),'C')
      CALL OPENLU(LUBIN,FNAME(LBIN),'C')
C     CALL PRVAR(LUOUT,NPAR,HEAP(VAR))
C set block structure
      NBKPJ=LBLK
      NDM=NDEGY
      KNT=SETBLK(LUOUT,IDPAR,PAR,NBKPJ,NDM)
      INBLK=NBKPJ*INBLK+1
      LBLK =NBKPJ*(LBLK+1)
C set up IFDIAG,NSAV  
      CALL SETINT(IFDIAG,NSAV,IDIP)
      NBSAV=(NSAV+1)*NBKPJ
      DIAG=IFDIAG.NE.0.AND.NDM.GT.1
      IF(IFDUMP.EQ.0) PREIG=PREIG.AND.DIAG
      DIAGD=DIAG.AND.(IFDUMP.EQ.0)
      IF (NBSAV.GE.NBSAVX) THEN
          WRITE(*,*) ' FATAL: NBSAV > NBSAVX',NBSAV,NBSAVX
          CALL OPENLU(LUOUT,FNAME(LOUT),'C')
          STOP
      ENDIF
      WRITE (LUOUT,460)
      KNT=0
C check that heap is within addressing constraints
      NPAR1=NPAR+1
      KNT=NDM
      NKNT=KNT*NDM
      LDEDP=0
      IF(DIAG) LDEDP=NKNT
      ISIZB=LDEDP+NPAR1*KNT
      IOFF=VAR
      SWK=IOFF-KNT
      S=SWK-NKNT*NPDIP
      NBSAVM=S/ISIZB
      IF(NBSAVM.LT.2) NBSAVM=2
      IF(NBSAVM.GT.NBSAV) NBSAVM=NBSAV
      IF(NBSAVM.LT.NBSAV) NBSAV=NBSAV+1
      KNT=NBSAVM*ISIZB-S
      IF(KNT.GT.0) THEN
          WRITE(*,*) ' Warning: HEALPLN should be larger by ',KNT
          STR=0.5*NPAR1*NBSAVM+0.5
          I=NPDIP
          IF(DIAG) I=I+NBSAVM
          NDM=(SQRT(STR*STR+IOFF*I) - STR)/I
          KNT=NDM
          NKNT=KNT*NDM
          LDEDP=0
          IF(DIAG) LDEDP=NKNT
          ISIZB=LDEDP+NPAR1*KNT
          SWK=IOFF-KNT
          S=SWK-NKNT*NPDIP
          WRITE(*,*) '        energy dimension truncated to ',NDM
      ENDIF
      K=-1
      I=RQEXIT(K)
C  if no diagonalization use S for working space
      LEIG=0
      IF(.NOT.DIAG) LEIG=S
C  set up intensity constants
      TMC=-1.43878
      TMQ=TMC/(TMQ*CMC)
C     TMQL= LOG10(EXP(1))*TMQ
      TMQL= 0.43429448D00*TMQ
      THRSH1=THRSH1-10.9542425
      STARG=THRSH-LOG10(FQMAX*FAC)
      SCOMP=-38.0
      STRMN= MAX (STARG,SCOMP)
      STRMN=10.**STRMN
      IF(PREGY)  CALL OPENLU(LUEGY,FNAME(LEGY),'WP')
      IF(PRSTR) CALL OPENLU(LUSTR,FNAME(LSTR),'WP')
      CALL OPENLU(LUCAT,FNAME(LCAT),'WP')
      DO 110 I=1,NBSAV
  110    IXBLK(I)=0
C********************************************************************
C
C START MAJOR LOOP OVER BLOCKS
C
      DO 400 IBLK=INBLK,LBLK
        FIRST=PRFRQ
        J=(IBLK-1)/NBKPJ
        IF(IBLK-J*NBKPJ.EQ.1) THEN
             IF(RQEXIT(0).NE.0) GO TO 499
C         ..brand old blocks
             IF(J.GT.NSAV) THEN
                 JBLK=(J-NSAV)*NBKPJ
                 DO 105 I=1,NBSAV
                     IF(IXBLK(I).LE.JBLK) IXBLK(I)=0
  105                CONTINUE
             ENDIF
             WRITE(*,495) J
        ENDIF
        IF(IXBLK(1).GT.0) THEN
C         move first block to unused position
            NXT=0
            IOFF=IBUFOF(NXT,NBSAVM,IXBLK,NBLK,NSBLK,ISIZB,HEAP)
            IXBLK(NXT)=IXBLK(1)
            NBLK(NXT)=NBLK(1)
            NSBLK(NXT)=NSBLK(1)
            IXBLK(1)=0
            NKNT=NSBLK(1)-1
            DO 130 KNT=0,NKNT
  130           HEAP(KNT+IOFF)=HEAP(KNT)
        ENDIF
C get energy and derivatives for new block
        CALL GETQN(IBLK,I0,IQNI,ISIZ,IDMY)
        IF (ISIZ.LE.0) GO TO 400
        IF (ISIZ.GT.NDM) GO TO 499
        NSIZD=ISIZ
        IOFF=ISIZ
        LDEDP=0
        IF(DIAG) LDEDP=IOFF*IOFF
        NSBLK(1)=LDEDP+NPAR1*IOFF
        NBLK(1)=ISIZ
        IXBLK(1)=IBLK
        CALL HAMX (IBLK,NSIZD,IDPAR,PAR,EGY,HEAP(LEIG),HEAP(LDEDP),
     :             IFDUMP)
        CALL DCOPY(ISIZ,EGY,I1,HEAP(NSBLK(1)-ISIZ),I1)
C print out energies and compute partition function
        IF(PRFRQ) WRITE (LUOUT,500) IBLK
        IOFF=-1
        LDEDP=LDEDP-1
        KBGN=1
        DO 190 I=1,ISIZ
          IF(IFDUMP.NE.0) KBGN=I
          CALL GETQN (IBLK,I,IQNI,IDGN,IDMY)
          IF(DIAGD.AND.IDGN.GT.0) THEN
C             check if eigenvector is zero
              IF(ABS(HEAP(IOFF+I)).LT.0.01) THEN
                  K=IDAMAX(ISIZ,HEAP(IOFF+1),I1)
                  IF(ABS(HEAP(IOFF+K)).LT.0.01) IDGN=0
              ENDIF
          ENDIF
          IF (IDGN.LE.0) GO TO 185
C         .. calculate uncertainty
          IXDP=LDEDP+I
          CALL DCOPY(NPAR,HEAP(IXDP),ISIZ,DERV,I1)
          ERR=MIN(BIGERR,CALERR(NPAR,HEAP(VAR),DERV)/CMC)
          EGX=EGY(I)/CMC
          IF (EGYMIN.GT.EGX) THEN
C             adjust energy zero
              TE=TMC*(EGYMIN-EGX)
              EGYMIN=EGX
              DO 160 K=1,NTEMP
  160             QSUM(K)=QSUM(K)*EXP(TE/TEMP(K))
          ENDIF
C         ..  accumulate partition function values
          TE=TMC*(EGX-EGYMIN)
          DGN=IDGN
          DO 180 K=1,NTEMP
  180         QSUM(K)=QSUM(K)+DGN*EXP(TE/TEMP(K))
          IF(PRFRQ)WRITE (LUOUT,510) I,IDGN,EGX,ERR,(IQNI(K),K=1,NQN)
          IF(PREGY)WRITE (LUEGY,510) IBLK,I,EGX,ERR,(IQNI(K),K=1,NQN)
          IF(PRDER)WRITE (LUEGY,596) (K,DERV(K),K=1,NPAR)
          IF(PREIG)WRITE (LUEGY,596) (K,HEAP(K+IOFF),K=KBGN,ISIZ)
C             subtract 1 from diagonal of eigenvector
  185     IF(DIAGD) HEAP(IOFF+I)=HEAP(IOFF+I)-ONE
  190     IOFF=IOFF+ISIZ
      JBLK=0
C...loop over previous blocks for intensity
  195     JJ=1
          JNXT=IBLK
C         find oldest block .GT. jblk
          DO 196 J=2,NBSAV
            JX=IXBLK(J)
            IF (JX.GT.JBLK.AND.JX.LT.JNXT) THEN
                JNXT=JX
                JJ=J
            ENDIF
  196       CONTINUE
          JBLK=JNXT
          JSIZ=NBLK(JJ)
          IOFF=IBUFOF(JJ,NBSAVM,IXBLK,NBLK,NSBLK,ISIZB,HEAP)
C
C  get intensity
C
          IDGN=0
          IDIP1(0)=1
          SIJ=S
          NS=ISIZ*JSIZ
          DO 200 I=1,NPDIP
            IF(NPDIP.EQ.1) THEN
              IDF=INTENS(IBLK,ISIZ,JBLK,JSIZ,IDIP,DIP,HEAP(S),NSIZD)
            ELSE
              IDIP1(1)=IDIP(I)
              IDF=INTENS(IBLK,ISIZ,JBLK,JSIZ,IDIP1,DIP(I),HEAP(SIJ),
     &                   NSIZD)
            ENDIF
            IF(IDF.NE.0) THEN
              IDGN=IDF
              IF (DIAGD) THEN
                CALL SIMT(ISIZ,JSIZ,HEAP(S),HEAP,HEAP(IOFF),HEAP(SWK))
              ENDIF
            ENDIF
            SIJ=SIJ+NS
  200       CONTINUE
          IF(IDGN.LE.0) GO TO 390
          DGNF=IDGN
C
          IF(DIAGD) IOFF=IOFF+JSIZ*JSIZ
          IOFF=IOFF-1
          IDF=-NDM
          IF (JBLK.EQ.IBLK) IDF=0
          DO 380 I=1,ISIZ
            CALL GETQN (IBLK,I,IQNI,IDGN,IDMY)
            IF (IDGN.LE.0) GO TO 380
            DGN=IDGN/DGNF
            DO 370 J=1,JSIZ
              IF ((I-J).LT.IDF) GO TO 380
              CALL GETQN (JBLK,J,IQNI(7),JDGN,IDMY)
              IF (JDGN.LE.0) GO TO 370
              IXDP=LDEDP+I
              JXDP=IOFF+J
              DO 340 K=1,NPAR
                DERV(K)=HEAP(IXDP)-HEAP(JXDP)
                IXDP=IXDP+ISIZ
  340           JXDP=JXDP+JSIZ
              FRQ=HEAP(IXDP)-HEAP(JXDP)
              IF (FRQ.LT.0.) THEN
                FRQ=-FRQ
                CALL QNFMT(IQNI(7),IQNI,NQN,IQN)
                IGUP=JDGN
                ELOW=EGY(I)
              ELSE
                CALL QNFMT(IQNI,IQNI(7),NQN,IQN)
                IGUP=IDGN
                ELOW=EGY(I)-FRQ
              ENDIF
              SIJ=S+I+ISIZ*(J-1)-1
              STRR=0.
              DO 350 IJ=1,NPDIP
                STR= HEAP( SIJ )
                STRR = STRR + STR
                STRQ = STR*STR
                IF (STRQ.GT.STCOMP) THEN
                  WRITE(LUSTR,530) FRQ,STR,IQNFMT,IQN,IJ
                ENDIF
                SIJ=SIJ+NS
  350           CONTINUE
              STRR = STRR * STRR
              IF (FRQ.GT.FQMAX) GO TO 370
              IF (FRQ.LT.FQMIN) GO TO 370
              IF (STRR.LT.STRMN) GO TO 370
              STR=DGN*STRR*FAC*FRQ*(1.D00-EXP(TMQ*FRQ))
              STRLG= LOG10(STR+ZERO)+TMQL*ELOW
              IF (STRLG.LT.THRSH) GO TO 370
              THRSHF=THRSH1+2.*LOG10(FRQ+ZERO)
              IF(ABS(THRSH-THRSHF).LT.4.)
     :            THRSHF=LOG10 (1.  +10.**(THRSH-THRSHF))+THRSHF
              IF (STRLG.LT.THRSHF) GO TO 370
C...calculate errors
              ERR=MIN(BIGERR,CALERR(NPAR,HEAP(VAR),DERV))
              ELOW=ELOW/CMC
              IF (FIRST) WRITE (LUOUT,520)
              FIRST=.FALSE.
              IF (PRFRQ) WRITE (LUOUT,550) FRQ,ERR,STRR,STRLG,ITD,ELOW,
     :                   IGUP,ITAG,IQNFMT,IQN
              NLINE=NLINE+1
              IF(PRIR) THEN
                  WRITE (LUCAT,545) FRQ,ERR,STRLG,ITD,ELOW,IGUP,ITAG,
     :                      IQNFMT,IQN
              ELSE
                  WRITE (LUCAT,540) FRQ,ERR,STRLG,ITD,ELOW,IGUP,ITAG,
     :                      IQNFMT,IQN
              ENDIF
  370         CONTINUE
  380       CONTINUE
  390     IF(JBLK.NE.IBLK) GO TO 195
  400   CONTINUE
C**************************************************************
C
C print partition function
C
  499 IF(IFDUMP.NE.0) THEN 
          WRITE(LUOUT,555)
          WRITE(*,555)
      ENDIF
      WRITE (LUOUT,560) NLINE,EGYMIN
      WRITE (*,560) NLINE,EGYMIN
      DO 405 I=1,NTEMP
          QLOG=-100.
          IF(QSUM(I).GT.ZERO) QLOG=LOG10(QSUM(I))
          WRITE (*,570)  I,TEMP(I),QSUM(I),QLOG
 405      WRITE (LUOUT,570)  I,TEMP(I),QSUM(I),QLOG
C
      IF(PREGY)  CALL OPENLU(LUEGY,FNAME(LEGY),'C')
      IF(PRSTR) CALL OPENLU(LUSTR,FNAME(LSTR),'C')
      CALL OPENLU(LUCAT,FNAME(LCAT),'C')
 600  CALL OPENLU(LUOUT,FNAME(LOUT),'C')
      STOP
C
 425  WRITE(*,*) ' ERROR IN READING  FIRST CARDS OF .INT'
      CALL OPENLU(LUINT,FNAME(LINT),'C')
      GO TO 600
C
C
  430 FORMAT (A)
  431 FORMAT (1X,A)
  440 FORMAT (I4,I10,E15.6)
  450 FORMAT('  ID=',I6,' QSPINROT=',F15.4/' MIN,'
     :,'MAX QN=',2I4/' MIN.LOG.STR='
     :,F10.3,' MIN.LOG.STR/(FRQ/300GHZ)**2=',F10.3,' MAX FREQ',F10.1,
     :  ' GHZ, TEMP=',F10.2)
  460 FORMAT ('  PARAMETER ID -- VALUE -- EST.ERROR')
  495 FORMAT (' STARTING QUANTUM ',I3)
  500 FORMAT (' ENERGIES FOR BLOCK NUMBER',I4/' INDEX-DEGEN-ENERGY'
     : ,'  -EST.ERROR   -   QUANTUM NUMBERS')
  510 FORMAT (2I5,2F18.6,6I3)
  520 FORMAT (' FREQUENCY-EST.ERROR.-LINE.STR. DIP**2-LGSTR.-ITD',
     : '-ELOW-GUP-I.D.-QNFORM-QUANTUM NUMBERS')
  530 FORMAT (F15.4,E15.6,I5,1X,12A2,I5)
  540 FORMAT (F13.4,2F8.4,I2,F10.4,I3,I7,I4,12A2)
  545 FORMAT (F13.6,F8.6,F8.4,I2,F10.4,I3,I7,I4,12A2)
  550 FORMAT (F13.4,F9.4,E13.6,F9.4,I2,F11.4,I4,I7,I4,12(1X,A2))
  555 FORMAT (' WARNING: THERE WAS NO DIAGONALIZATION')
  560 FORMAT (' NUMBER OF LINES = ',I6,/,
     :        ' TEMPERATURE - Q(SPIN-ROT.) RELATIVE TO MIN.EGY.=',F15.4)
  570 FORMAT (I3,F10.3,F15.4,F10.4)
  596 FORMAT(10X,I7,E13.6,I7,E13.6,I7,E13.6,I7,E13.6,
     :           I7,E13.6,I7,E13.6)
C
      END
      SUBROUTINE QNFMT(IQU,IQL,NQN,IQN)
C    subroutine to do quick conversion of quantum integers to characters
      INTEGER IQU(*),IQL(*),NQN
      CHARACTER*2 IQN(12),NNFMT(0:100),TQN
      INTEGER K,I,IX,IPCHR,INCHR
      PARAMETER (IPCHR=ICHAR('A')-10, INCHR=ICHAR('a')-1)
      SAVE NNFMT
      DATA NNFMT/' 0'
     1     ,' 1',' 2',' 3',' 4',' 5',' 6',' 7',' 8',' 9','10','11','12'
     2     ,'13','14','15','16','17','18','19','20','21','22','23','24'
     3     ,'25','26','27','28','29','30','31','32','33','34','35','36'
     4     ,'37','38','39','40','41','42','43','44','45','46','47','48'
     5     ,'49','50','51','52','53','54','55','56','57','58','59','60'
     6     ,'61','62','63','64','65','66','67','68','69','70','71','72'
     7     ,'73','74','75','76','77','78','79','80','81','82','83','84'
     8     ,'85','86','87','88','89','90','91','92','93','94','95','96'
     9     ,'97','98','99','**'/
      DO 310 K=1,NQN
          I=IQU(K)
          IF(I.LT.0) THEN
              I=-I 
              IF(I.LE.9) THEN
                  TQN=NNFMT(I)
                  TQN(1:1)='-'
              ELSE IF(I.LT.270) THEN
                  IX=I/10
                  I=I-10*IX
                  TQN=NNFMT(I)
                  TQN(1:1)=CHAR(IX+INCHR)
              ELSE
                  TQN=NNFMT(100)
              ENDIF
          ELSE IF(I.LT.100) THEN
              TQN=NNFMT(I)
          ELSE IF(I.LT.360) THEN
              IX=I/10
              I=I-IX*10
              TQN=NNFMT(I)
              TQN(1:1)=CHAR(IX+IPCHR)
          ELSE
              TQN=NNFMT(100)    
          ENDIF
          IQN(K)=TQN
          I=IQL(K)
          IF(I.LT.0) THEN
              I=-I
              IF(I.LE.9) THEN
                  TQN=NNFMT(I)
                  TQN(1:1)='-'
              ELSE IF(I.LT.270) THEN
                  IX=I/10
                  I=I-10*IX
                  TQN=NNFMT(I)
                  TQN(1:1)=CHAR(IX+INCHR)
              ELSE
                  TQN=NNFMT(100)
              ENDIF
          ELSE IF(I.LT.100) THEN
              TQN=NNFMT(I)
          ELSE IF(I.LT.360) THEN
              IX=I/10
              I=I-IX*10
              TQN=NNFMT(I)
              TQN(1:1)=CHAR(IX+IPCHR)
          ELSE
              TQN=NNFMT(100)    
          ENDIF
  310     IQN(K+6)=TQN
      RETURN
      END
      SUBROUTINE SIMT(ISIZ,JSIZ,S,T,U,WK)
C    subroutine to do similarity transform
C      S' = transpose(T+1) * S * (U+1)
      INTEGER ISIZ,JSIZ,I1,I,IJ,J,JT
      PARAMETER (I1=1)
      REAL*8 S(*),T(*),U(*),WK(*),DDOT
C      EMA S,T,U
C...do transform on the left
      IF(ISIZ.GT.I1) THEN
          IJ=1
          DO 270 J=1,JSIZ
            CALL DCOPY(ISIZ,S(IJ),I1,WK,I1)
            JT=1
            DO 260 I=1,ISIZ
              S(IJ)=S(IJ)+DDOT(ISIZ,WK,I1,T(JT),I1)
              JT=JT+ISIZ
  260         IJ=IJ+1
  270       CONTINUE
      ENDIF
C do transform on the right
      IF(JSIZ.GT.I1) THEN
          DO 230 I=1,ISIZ
            CALL DCOPY(JSIZ,S(I),ISIZ,WK,I1)
            JT=1
            IJ=I
            DO 220 J=1,JSIZ
              S(IJ)=S(IJ)+DDOT(JSIZ,WK,I1,U(JT),I1)
              JT=JT+JSIZ
  220         IJ=IJ+ISIZ
  230       CONTINUE
      ENDIF
      RETURN
      END
      FUNCTION IBUFOF(IBLK,NBSAVM,IXBLK,NBLK,NSBLK,ISIZB,HEAP)
      INTEGER NBSAVM,IBLK,IXBLK(*),NBLK(*),LREC,LU,JBLK,N,I,IOS
      PARAMETER (LREC=256)
      INTEGER*4 IBUFOF,NSBLK(*),ISIZB,LSIZB,MAXREC,NX,LSRC,LDEST
      REAL*8 HEAP(0:*)
      SAVE LU,LSIZB,MAXREC
      DATA LU /-1/
      IF(IBLK.LE.NBSAVM.AND.IBLK.NE.0) THEN
          IBUFOF=(IBLK-1)*ISIZB
          RETURN
      ENDIF
      JBLK=1
 1        JBLK=JBLK+1
          IF(IXBLK(JBLK).GT.0) GO TO 1
      IF(JBLK.GT.NBSAVM) THEN
C     ..... need to save block 2 in JBLK
          IF(LU.LT.0) THEN
              LU=LREC
              CALL OPENLU(LU,'Heap Scratch','S')
              LSIZB=(ISIZB-1)/LREC+1
              MAXREC=0
          ENDIF
          LDEST=(JBLK-NBSAVM-1)*LSIZB
          LSRC=ISIZB-1
  5       IF(LDEST.GT.MAXREC) THEN
C           .. fill in gaps with junk FIRST time
              MAXREC=MAXREC+1
              WRITE(LU,REC=MAXREC,ERR=90,IOSTAT=IOS) 
     :              (HEAP(LSRC+I),I=1,LREC)
              GO TO 5
          ENDIF      
          NX=NSBLK(2)
          NSBLK(JBLK)=NX
  2       IF(NX.GT.0) THEN
              LDEST=LDEST+1
              WRITE(LU,REC=LDEST,ERR=90,IOSTAT=IOS)
     :              (HEAP(LSRC+I),I=1,LREC)
              LSRC=LSRC+LREC
              NX=NX-LREC
              GO TO 2
          ENDIF
          IF(MAXREC.LT.LDEST) MAXREC=LDEST
          IXBLK(JBLK)=IXBLK(2)
          NBLK(JBLK)=NBLK(2)
          IXBLK(2)=0
          JBLK=2
      ENDIF
      IBUFOF=(JBLK-1)*ISIZB
      IF(IBLK.EQ.0) THEN
          IBLK=JBLK
      ELSE
C       read IBLK into JBLK
          LDEST=IBUFOF-1
          LSRC=(IBLK-NBSAVM-1)*LSIZB
          NX=NSBLK(IBLK)
          NSBLK(JBLK)=NX
  10      IF(NX.GT.0) THEN
              IF(NX.GT.LREC) THEN
                  N=LREC
              ELSE
                  N=NX
              ENDIF
              LSRC=LSRC+1
              READ(LU,REC=LSRC,ERR=91,IOSTAT=IOS)
     :            (HEAP(LDEST+I),I=1,N)
              LDEST=LDEST+LREC
              NX=NX-LREC
              GO TO 10
          ENDIF
          NBLK(JBLK)=NBLK(IBLK)
          IF(JBLK.NE.2) THEN
              IXBLK(JBLK)=IXBLK(IBLK)
              IXBLK(IBLK)=0
              IBLK=JBLK
          ENDIF
      ENDIF
      RETURN
 90   WRITE(*,*) ' Heap write ERROR, IOSTAT=',IOS
      STOP 'Aborted'
 91   WRITE(*,*) ' Heap read ERROR, IOSTAT=',IOS
      STOP 'Aborted'
      END