      SUBROUTINE THERMO(A,B,C,LINEAR,SYM,WT,VIBS,NVIBS,ESCF)
      implicit double precision (A-H,O-Z)
      dimension VIBS(*)
      logical LINEAR
      character KEYWRD*241, KOMENT*81, TITLE*81, TMPKEY*241
      include 'MOP6FILE'
      common /KEYWRD/ KEYWRD
      common /TITLES/ KOMENT,TITLE
C
C
C   THERMO CALCULATES THE VARIOUS THERMODYNAMIC QUANTITIES FOR A
C   SPECIFIED TEMPERATURE GIVEN THE VIBRATIONAL FREQUENCIES, MOMENTS OF
C   INERTIA, MOLECULAR WEIGHT AND SYMMETRY NUMBER.
C
C   REFERENCE: G.HERZBERG MOLECULAR SPECTRA AND MOLECULAR STRUCTURE
C              VOL 2, CHAP. 5
C
C   ----    TABLE OF SYMMETRY NUMBERS    ----
C
C        C1 CI CS     1      D2 D2D D2H  4       C(INF)V   1
C        C2 C2V C2H   2      D3 D3D D3H  6       D(INF)H   2
C        C3 C3V C3H   3      D4 D4D D4H  8       T TD     12
C        C4 C4V C4H   4      D6 D6D D6H  12      OH       24
C        C6 C6V C6H   6      S6          3
C
C
C   PROGRAM LIMITATIONS:  THE EQUATIONS USED ARE APPROPRIATE TO THE
C   HIGH TEMPERATURE LIMIT AND WILL BEGIN TO BE INADEQUATE AT TEMPERA-
C   TURES BELOW ABOUT 100 K.  SECONDLY THIS PROGRAM IS ONLY APPROPRIATE
C   IN THE CASE OF MOLECULES IN WHICH THERE IS NO FREE ROTATION
C
C
C
C
*******************************************************************
*
*  THE FOLLOWING CONSTANTS ARE NOW DEFINED:
*          PI  = CIRCUMFERENCE TO DIAMETER OF A CIRCLE
*          R   = GAS CONSTANT IN CALORIES/MOLE
*          H   = PLANCK'S CONSTANT IN ERG-SECONDS
*          AK  = BOLTZMANN CONSTANT IN ERG/DEGREE
*          AC  = SPEED OF LIGHT IN CM/SEC
*******************************************************************
      save PI, R, H, AK, AC
      dimension TRANGE(300)
      data PI /3.14159D0 /
      data R/1.98726D0/
      data H/6.626D-27/
      data AK/1.3807D-16/
      data AC/2.99776D+10/
*******************************************************************
      IT1=200
      IT2=400
      ISTEP=10
      TMPKEY=KEYWRD
      I=INDEX(TMPKEY,'THERMO(')
      if (I.NE.0) then
C
C   ERASE ALL TEXT FROM TMPKEY EXCEPT THERMO DATA
C
         TMPKEY(:I)=' '
         TMPKEY(INDEX(TMPKEY,')'):)=' '
         IT1 = READA(TMPKEY,I)
         if (IT1.LT.100) then
            write(MIOB,'(//10X,''TEMPERATURE RANGE STARTS TOO LOW,'',
     %        '' LOWER BOUND IS RESET TO 30K'')')
            IT1=100
         end if
         I=INDEX(TMPKEY,',')
         if (I.NE.0) then
            TMPKEY(I:I)=' '
            IT2 = READA(TMPKEY,I)
            if (IT2.LT.IT1) then
               IT2=IT1+200
               ISTEP=10
               go to 10
            end if
            I=INDEX(TMPKEY,',')
            if (I.NE.0) then
               TMPKEY(I:I)=' '
               ISTEP = READA(TMPKEY,I)
               if (ISTEP.LT.1) ISTEP=1
            else
               ISTEP=(IT2-IT1)/20
               if (ISTEP.EQ.0)                    ISTEP=1
               if (ISTEP.GE.2.AND. ISTEP.LT.5)    ISTEP=2
               if (ISTEP.GE.5.AND. ISTEP.LT.10)   ISTEP=5
               if (ISTEP.GE.10.AND. ISTEP.LT.20)  ISTEP=10
               if (ISTEP.GT.20.AND. ISTEP.LT.50)  ISTEP=20
               if (ISTEP.GT.50.AND. ISTEP.LT.100) ISTEP=50
               if (ISTEP.GT.100)                  ISTEP=100
            end if
         else
            IT2=IT1+200
         end if
      end if
   10 continue
      write(MIOB,'(//,A)') TITLE
      write(MIOB,'(A)')    KOMENT
      if (LINEAR) then
         write(MIOB,'(//10X,''MOLECULE IS LINEAR'')')
      else
         write(MIOB,'(//10X,''MOLECULE IS NOT LINEAR'')')
      end if
      write(MIOB,
     %  '(/10X,''THERE ARE'',I3,'' GENUINE VIBRATIONS IN THIS '',
     %         ''SYSTEM'')')  NVIBS
      write(MIOB,20)
   20 format(10X,'THIS THERMODYNAMICS CALCULATION IS LIMITED TO',/
     110X,'MOLECULES WHICH HAVE NO INTERNAL ROTATIONS'//)
      write(MIOB,'(//20X,''CALCULATED THERMODYNAMIC PROPERTIES'')')
      write(MIOB,'(42X,''*'')')
      write(MIOB,'(''   TEMP. (K)   PARTITION FUNCTION   H.O.F.'',
     %                  ''    ENTHALPY   HEAT CAPACITY  ENTROPY'')')
      write(MIOB,'(  ''                                    KCAL/MOL'',
     %                  ''   CAL/MOLE    CAL/K/MOL   CAL/K/MOL''/)')
      do 30 I=1,NVIBS
   30 VIBS(I)=ABS(VIBS(I))
      ILIM=1
      do 40 ITEMP=IT1,IT2,ISTEP
         ILIM=ILIM+1
   40 TRANGE(ILIM)=ITEMP
      TRANGE(1)=298.D0
      do 80 IR=1,ILIM
         ITEMP=TRANGE(IR)
         T=ITEMP
C   ***   INITIALISE SOME VARIABLES   ***
         C1=H*AC/AK/T
         QV=1.0D0
         HV=0.0D0
         E0=0.0D0
         CPV=0.0D0
         SV1=0.0D0
         SV2=0.0D0
C   ***   CONSTRUCT THE FREQUENCY DEPENDENT PARTS OF PARTITION FUNCTION
         do 50 I=1,NVIBS
            WI=VIBS(I)
            EWJ=EXP(-WI*C1)
            QV=QV/(1-EWJ)
            HV=HV+WI*EWJ/(1-EWJ)
            E0=E0+WI
            CPV=CPV+WI*WI*EWJ/(1-EWJ)/(1-EWJ)
            SV1=SV1+LOG(1.0D0-EWJ)
   50    SV2=SV2+WI*EWJ/(1-EWJ)
C   ***   FINISH CALCULATION OF VIBRATIONAL PARTS   ***
         HV=HV*R*H*AC/AK
         E0=E0*1.4295D0
         CPV=CPV*R*C1*C1
         SV=SV2*R*C1-R*SV1
C   ***   NOW CALCULATE THE ROTATIONAL PARTS  (FIRST LINEAR MOLECULES
         if (.NOT.LINEAR) go to 60
         QR=1/(C1*A*SYM)
         HR=R*T
         CPR=R
         SR=R*(LOG(T*AK/(H*AC*A*SYM)))+R
         go to 70
   60    QR=SQRT(PI/(A*B*C*C1*C1*C1))/SYM
         HR=3.0D0*R*T/2.0D0
         CPR=3.0D0*R/2.0D0
         SR=0.5D0*R*(3.D0*LOG(T*AK/(H*AC))
     1-2.D0*LOG(SYM)+LOG(PI/(A*B*C))+3.D0)
   70    continue
C   ***   CALCULATE INTERNAL CONTRIBUTIONS   ***
         QINT=QV*QR
         HINT=HV+HR
         CPINT=CPV+CPR
         SINT=SV+SR
C   ***   CONSTRUCT TRANSLATION CONTRIBUTIONS   ***
         QTR=(SQRT(2.D0*PI*WT*T*AK*1.6606D-24)/H)**3
         HTR=5.0D0*R*T/2.0D0
         CPTR=5.0D0*R/2.0D0
         STR=2.2868D0*(5.0D0*LOG10(T)+3.0D0*LOG10(WT))-2.3135D0
C   ***   CONSTRUCT TOTALS   ***
         CPTOT=CPTR+CPINT
         STOT=STR+SINT
         HTOT=HTR+HINT
C   ***   OUTPUT SECTION   ***
         if (IR.EQ.1) then
            H298=HTOT
         else
            write(MIOB,'(/,I7,''  VIB.'',G18.4
     %           ,13X,3F11.5        )') ITEMP,QV,  HV,  CPV,  SV
            write(MIOB,'(7X,''  ROT.'',G13.3
     %           ,16X,3F11.3        )')       QR,  HR,  CPR,  SR
            write(MIOB,'(7X,''  INT.'',G13.3
     %           ,16X,3F11.3        )')     QINT,HINT,CPINT,SINT
            write(MIOB,'(7X,''  TRA.'',G13.3
     %           ,16X,3F11.3)')
     %                                       QTR, HTR, CPTR, STR
            write(MIOB,'(7X,''  TOT.'',13X,F17.3,F11.4,2F11.4)')
     %                  ESCF+(HTOT-H298)/1000.D0,HTOT,CPTOT,STOT
         end if
   80 continue
      write(MIOB,
     %  '(/3X,'' * NOTE: HEATS OF FORMATION ARE RELATIVE TO THE''/
     %    12X,'' ELEMENTS IN THEIR STANDARD STATE AT 298K'')')
      end
