
      PROGRAM cap
C
C  CAP is a program written by Bonnie McBride at NASA Lewis 1992 
C
C  CAP stands for Coefficients and Properties.  It calculates
C     thermodynamic functions, including delta H and log K from
C     coeficients in the new form as produced by the PAC91 program.
C     In a way CAP is the reverse of PAC.
C
C  For delta H and log K calculations, a file of data in the same
C     format for the reference elements is required for I/O unit 13.
C     Either associations or OPEN statements need to be included.
C
C  For the standard input on I/O unit 5, a list of options and a
C     temperature schedule need to be included ahead of the coefficient
C     data.  This information is used for the entire run.  These two
C     records must each end with a slash (/).
C
C  Record one contains 5 possible character variables in single quotes
C     with space and/or commas between them.  They may be either all
C     upper or all lower case.  They contain the options as follows:
C
C     'logk'   Include delta H and log K values with rounded thermo
C              functions.
C     'mfig'   Do tables of thermo functions with many figures and no
C              delta H and log K values.
C     'joules'  Energy units for the output will be in joules or
C               kilojoules.  (JOULES is the default option)
C     'cal'  Energy units for the output will be in calories or Kcal.
C     'nodim'  Thermodynamic functions for a many-figured table will
C              be dimensionless.  (i.e. NODIM is only valid when used
C              with MFIG)
C
C  Record two gives the temperature schedule in the form of temperatures
C     and temperature intervals.  The even-numbered values are
C     intervals to be used between the odd-numbered temperature values.
C     The schedule will be used for the entire computer run.
C     For each set of coefficients, CAP automatically
C     limits the schedule according to the temperatures ranges given
C     with the coefficients.  Endpoints, 298.15, and transition points
C     are automatically added if they are within the requested range.
C     For transition points, values will be given for each phase at
C     the common point.  Up to 201 values may be given.  They are
C     read with a list-directed format.  That is, values are given
C     with either space or commas between them.
C     Even numbered interval values must be included even
C     if they are zero.  Two adjacent commas will indicate a zero value.
C     For example, to get the temperature set: 100,298.15,500,1000,1500,
C     2000,and 3000, the record might be as follows.
C       100,0,500,500. 2000,,3000./
C
C  The following is a sample set for the two header records:
C
C        'mfig' 'LOGK' 'joules' 'cal' /
C    200,100,500,500,3000,1000,6000. /
C
C      This should give both many-figured tables and rounded tables
C      with delta H and log K for both sets of enery units, joules and
C      calories.  The temperature schedule should be 200 K every 100 K
C      to 500 K, every 500 K to 3000 K, and every 1000 K to
C      6000 K with 298.15 inserted between 200 and 300 K.  (i.e. 200,
C      298.15,300,400,500,1000,1500,2000,2500,3000,4000,5000,6000.)
C
      LOGICAL TAB8,calclk,rec3
      DOUBLE PRECISION EXL,TEX,TT,tin(201),temps(202),wt,htf
      dimension sl(5)
C
      CHARACTER*4 ALOGK,rec1(4),ulogk,ucal,ujoule
      CHARACTER*4 CAL,gas,fmla(2,5)
      character*8 output(5),joules,mfig,nodim,unodim,umfig,code
      character*64 record
C
      CHARACTER*4 BLANK,DATE,EFAZ,FAZ,FORMLA
      CHARACTER*4 NAME,SYMBOL
      COMMON /CHAR4/ BLANK,DATE(2),EFAZ(100),FAZ,
     1 FORMLA(2,5),NAME(4,9),SYMBOL(2,100)
C
      DOUBLE PRECISION asindt,ATMWT,COEF,CPR,GHRT,htform,enth29,
     1 HHRT,R,RR,SR,T,tcoef,term,WEIGHT,WORD,TTRANS
      COMMON /REAL8/ asindt,ATMWT(100),COEF(10,20),CPR(202),
     1 GHRT(202),htform,enth29,HHRT(202),R,RR,
     2 SR(202),T(202),tcoef(2,20),term,WEIGHT,WORD(4),TTRANS(6)
C
      COMMON /REAL4/ EX(8,20),sel(5),sels(5)
C
      COMMON /INTGR/ IC80,itr,JF(5),
     1 numexp(20),NEL,NKIND,NKINDS,
     2 NMAT(100),NMLA(100),NOATMS,NT,NTMP,NFZ,INA,NNA
C
      LOGICAL hasind,refel,TOUT
      COMMON /LOGCL/ hasind,refel,TOUT(3)
C
      OPEN(5,file='cap.dat')
      OPEN(6,file='sampcap.out')
      OPEN(13,file='elements.dat')
C
      DATA gas/'G'/, nodim/'nodim'/, unodim/'NODIM'/
      DATA UMFIG/'MFIG'/, UCAL/'CAL'/, UJOULE/'JOULES'/, ULOGK/'LOGK'/
      data mfig/'mfig'/, cal/'cal'/, joules/'joules'/, alogk/'logk'/
C
C  INITIALIZE ONCE.
C
1     R=RR/4.184D0
      DO 10 I=1,202
      CPR(I)=0.
      temps(i)=0.d0
10    T(I)=0.
      NT=0
      do 15 i=1,5
      output(i)=blank
      formla(1,i)=blank
      sel(i) = 0.
15    continue
      DO 20 I=1,3
20    TOUT(I)=.FALSE.
      TAB8=.FALSE.
      calclk=.false.
C
C  READ AND WRITE CONTENTS OF INPUT RECORD with output information.
C
      READ (5,*,END=250) OUTPUT
      WRITE (6,*) 'OPTIONS:  ',OUTPUT
C
C     PROCESS OUTPUT RECORD
C
      DO 25  I=1,5
      IF (output(i).EQ.joules .or. output(i).eq.ujoule) THEN
      TOUT (2)=.TRUE.
      ELSEIF (output(i).EQ.CAL .or. output(i).eq.ucal) THEN
      TOUT (3)=.TRUE.
      ELSEIF (output(i).EQ.nodim .or. output(i).eq.unodim) THEN
      TOUT (1)=.TRUE.
      TAB8=.TRUE.
      ELSEIF (output(i).EQ.ALOGK .or. output(i).eq.ulogk) THEN
      calclk=.TRUE.
      ELSEIF (output(i).EQ.mfig .or. output(i).eq.umfig) THEN
      TAB8=.TRUE.
      ENDIF
25    CONTINUE
c
c  Read temperature schedule.
c
      nts=0
      read (5,*,end=250) tin
      do 35 i=1,201,2
      if (tin(i).eq.0.) go to 40
      i1=i+1
      nts=nts+1
      temps(nts)=tin(i)
      if (tin(i1).eq.0.)go to 35
30    tt=temps(nts) + tin(i1)
      if(tt.lt.tin(i+2)-.0001d0)then
         nts=nts+1
         temps(nts) = tt
         if(nts.le.202) go to 30
      endif
35    continue
40    do 45 i=1,nts
      if(temps(i).gt.298.15d0-.005) go to 60
      if (temps(i+1).lt.298.15d0+.005) go to 45
      it=i+1
      go to 50
45    continue
      go to 60
50    nts=nts+1
      do 55 i=nts,it+1,-1
      temps(i)=temps(i-1)
55    continue
      temps(it)=298.15d0
60    WRITE(6,65) (TEMPS(I),I=1,NTS)
65    format('0TEMPERATURE SCHEDULE:',//(6f12.3))
      WRITE (6,70)
70    FORMAT (//35x,'***** WARNING *****',//3x,'The temperature range of
     * the coefficients read in the standard input (I/O unit 5)',/,3x,
     * 'will be extended by 20 %.  E.g. a range of 300 to 5000 K will be
     * extended to'/,3x, '240 to 6000 K.  The extrapolated data may have
     * large errors.  The element data'/,3x, 'on I/O unit 13 will NOT be
     * extrapolated when calculating delta H and log K.'/'1')
      rec3=.false.
c
c   End Initialization.
c
80    nt=nts
      do 85 i=1,nt
      t(i)=temps(i)
85    continue
c
      DO 90 i=1,20
      DO 90 J=1,8
      COEF(J,I)=0.
      EX(J,I)=0.
90    CONTINUE
      IX=1
      ITR=0
      INA=1
      NNA=0
      if (rec3) go to 130
95    READ(5,100,end=260) (rec1(i),i=1,4),record
100   format(4a4,a64)
      READ(5,105,END=1000) nit,code,(fmla(1,I),fmla(2,I),SL(I),
     * I=1,5),ifz,wt,htf
105   FORMAT (i2,a8,5(2a1,f6.2),i2,f13.5,f15.3)
      if (formla(1,1).eq.blank) go to 130
      if (ifz.le.ifaz) go to 120
      do 110 i=1,5
      if (fmla(1,i).ne.formla(1,i)) go to 120
      if (fmla(2,i).ne.formla(2,i)) go to 120
      if (sel(i).ne.sl(i)) go to 120
110   continue
      go to 130
120   rec3=.true.
      go to 260
130   do 140 i=1,5
      formla(1,i)=fmla(1,i)
      formla(2,i)=fmla(2,i)
      sel(i)=sl(i)
140   continue
      rec3=.false.
      nna=nna+1
      do 145 i=1,4
      name(i,nna)=rec1(i)
145   continue
      numint=nit
      ifaz=ifz
      htform=htf
      weight=wt
c
      WRITE(6,150) (name(I,nna),I=1,4),record
150   format('0',4a4,a64)
      WRITE(6,105) NUMINT,code,(FORMLA(1,i),FORMLA(2,i),sel(i),
     * I=1,5),ifaz,weight,htform
      do 180 i=1,numint
      itr=itr+1
      read(5,160,end=1000) tcoef(1,itr),tcoef(2,itr),numexp(itr)
     * ,(ex(j,itr),j=1,8),enth29
160   format(1x,2f10.3,i2,8f5.1,f17.3)
      WRITE(6,160) TCOEF(1,ITR),TCOEF(2,ITR),NUMEXP(ITR)
     * ,(ex(j,itr),j=1,8),enth29
      read (5,170,end=1000) (coef(j,itr),j=1,10)
170   format(5d16.8)
      WRITE (6,170) (coef(j,itr),j=1,10)
180   continue
      ttrans(nna)=tcoef(2,itr)
      do 200 i=1,5
      if (sel(i).eq.0.) go to 95
      nkind=i
200   continue
      go to 95
C
250   WRITE (6,255)
255   FORMAT ('0ERROR IN INPUT, GO TO NEXT SPECIES, (MAIN PGM)')
      GO TO 1000
C
C  INSERT TRANSITION POINTS IN CTEM SCHEDULE
C
260   if(nna.le.1) go to 315
      it=1
      do 310 n=1,nna-1
      tt=ttrans(n)
      do 270 i=it,nt-1
      if(tt.gt.t(i)) go to 270
      ii=i
      ix=2
      if(tt.eq.t(i))ix=1
      go to 280
270   CONTINUE
      go to 310
280   it=ii
      NT=NT+IX
      do 300 i=nt,it,-1
300   t(i)=t(i-ix)
      t(it)=tt
      t(it+1)=tt
310   continue
315   IT=0
C
C   The next 2 statements extend the T range of the coefficients.
C
      tcoef(1,1)=.8d0*tcoef(1,1)
      tcoef(2,itr)=1.2d0*tcoef(2,itr)
      if(nt.le.0) go to 250
      L=1
320   IT=IT+1
322   TT=T(IT)
      if (tt.ge.tcoef(1,L)) go to 330
      nt=nt-1
      if(nt.le.it) go to 360
      do 324 i=it,nt
324   t(i)=t(i+1)
      go to 322
330   IF (TT.LT.TCoef(2,L)) GO TO 340
      if(dabs(tt-tcoef(2,l)).lt.0005 .and. tt.ne.t(it-1)) go to 340
      if (tt.gt.tcoef(2,l)+.001 .and. t(it-1).lt.tcoef(2,l)-.001)then
        nt=nt+1
        do 335 j=nt,it+1,-1
        t(j)=t(j-1)
335     continue
        tt=tcoef(2,l)
        t(it)=tt
        if (l.eq.itr) nt=it
        go to 340
      endif
      L=L+1
      IF (L.LE.itr) GO TO 330
      nt=it-1
      GO TO 360
340   CPR(IT)=0.
      HHRT(IT)=0.
      SR(IT)=0.
      GHRT(IT)=0.
      DO 350 J=1,numexp(L)
      EXL=EX(J,L)
      TEX=COEF(J,L)
      IF (EXL.NE.0.) TEX=TEX*TT**EXL
      IF (EXL.EQ.-1.) HHRT(IT)=HHRT(IT)+COEF(J,L)*DLOG(TT)/TT
      IF (EXL.NE.-1.) HHRT(IT)=HHRT(IT)+TEX/(EXL+1.D0)
      CPR(IT)=CPR(IT)+TEX
      IF (EXL.EQ.0.) SR(IT)=SR(IT)+COEF(J,L)*DLOG(TT)
      IF (EXL.NE.0.) SR(IT)=SR(IT)+TEX/EXL
350   CONTINUE
      HHRT(IT)=HHRT(IT)+(COEF(9,L))/TT
      SR(IT)=SR(IT)+COEF(10,L)
      GHRT(IT)=SR(IT)-HHRT(IT)
      IF (IT.LT.NT) GO TO 320
360   if (.not.tout(1).and..not.tout(3)) tout(2)=.true.
c
c  Store element numbers from BLOCK DATA in JF array.  If
c    species is a reference element, set relel = .true.
c
      refel=.false.
      do 400 k=1,nkind
      jf(k)=0
      do 380 i=1,100
      do 370 j=1,2
      if (SYMBOL(j,i).ne.formla(j,k)) go to 380
370   continue
      jf(k)=i
      if(nkind.gt.1 .or. htform.ne.0.) go to 400
      n=sel(1)+.001
      if (n.ne.nmla(i)) go to 400
      if (efaz(i).eq.gas .and. ifaz.ne.0) go to 400
      if (efaz(i).ne.gas .and. ifaz.eq.0) go to 400
      refel=.true.
      go to 400
380   continue
      WRITE (6,390) formla(1,k),formla(2,k)
390   format('0*** Element '2A1,' not found in BLOCK DATA')
      calclk = .false.
      go to 405
400   continue
c
405   if(refel.or.htform.ne.0.d0) go to 430
      if(t(1).gt.298.151 .or. t(nt).lt.298.149) go to 430
      do 410 i=1,nt
      if(dabs(t(i)-298.15d0).gt..0001)go to 410
      htform = hhrt(i)*rr*298.15d0
      go to 430
410   continue
C
430   term=(htform-enth29)/rr
      asindt=0.
      if (htform.ne.0. .or. refel) then
           hasind=.true.
           if(enth29.eq.0.) asindt=298.15d0
      else
           hasind=.false.
      endif
      IF (TAB8) CALL TABLES
      IF (calclk) CALL LOGK
      if(rec3) go to 80
C
1000  stop
      END
c
      SUBROUTINE LOGK
C
C  CALCULATE ROUNDED TABLES WITH DELTAH AND LOGK CALCULATIONS.
C
      DIMENSION PT(30),MARK(30),VFM2(15),idone(5)
C
      DOUBLE PRECISION HZERO,RT,RUSE(2)
      logical kdone
      DOUBLE PRECISION ALK(202),DH0,CP,S,H,HH,GH,DH(202),H0
      DOUBLE PRECISION HH0,DNMLA,EXL,TEX,TT,H298,dmla
C
      CHARACTER*4 CHH,CHDH,BZ,VFM2,FSTAR,F2B,ZERO,FI6,ename(2,30)
      CHARACTER*4 FA4,F3X,F8X,F9,F12,FP2,FP3,FP4,fm(2),fms(2),f312,f224
C
      CHARACTER*4 BLANK,DATE,EFAZ,FAZ,FORMLA
      CHARACTER*4 NAME,SYMBOL
      COMMON /CHAR4/ BLANK,DATE(2),EFAZ(100),FAZ,
     1 FORMLA(2,5),NAME(4,9),SYMBOL(2,100)
C
      DOUBLE PRECISION asindt,ATMWT,COEF,CPR,GHRT,htform,enth29,
     1 HHRT,R,RR,SR,T,tcoef,term,WEIGHT,WORD,TTRANS
      COMMON /REAL8/ asindt,ATMWT(100),COEF(10,20),CPR(202),
     1 GHRT(202),htform,enth29,HHRT(202),R,RR,
     2 SR(202),T(202),tcoef(2,20),term,WEIGHT,WORD(4),TTRANS(6)
C
      COMMON /REAL4/ EX(8,20),sel(5),sels(5)
C
      COMMON /INTGR/ IC80,itr,JF(5),
     1 numexp(20),NEL,NKIND,NKINDS,
     2 NMAT(100),NMLA(100),NOATMS,NT,NTMP,NFZ,INA,NNA
C
      LOGICAL hasind,refel,TOUT
      COMMON /LOGCL/ hasind,refel,TOUT(3)
C
      EQUIVALENCE (IT,TI)
C
      DATA  FSTAR/'2H *'/,F9/',F9'/,FP2/'.2'/, FI6/',I6'/,F3X/',3X'/
     1 ,FA4/',A4'/,FP4/'.4'/,FP3/'.3'/, ZERO/'0'/
     2 , F12/',F12'/, F8X/',8X'/, F2B/'2H  '/,f312/'3f12'/,f224/'2f24'/
C
      DATA VFM2/'(', '2H  ', ',F9', '.2', ',F9', '.3,', 'F12', '.3',
     1   ',F12', '.3', ',F12', '.3',  ',F12', '.4', ')'/
C
1     REWIND 13
      CHH=FSTAR
      CHDH=FSTAR
      RUSE(1)=RR
      RUSE(2)=R
      IK=1
      fms(1)=blank
      iend=0
      INIT=1
      NTX=nt
      if (hasind) then
         vfm2(7)=f312
      else
         vfm2(7)=f224
      endif
      VFM2( 9)=F12
      VFM2(10)=FP3
      DO 10 I=1,NT
      DH(I)=0.D0
10    ALK(I)=0.D0
      DH0=0.D0
C
C  CHDH = BLANK OR ZERO - USE CHARACTER FOR DELTA H AND LOGK COLUMNS
C  CHH = BLANK - USE BLANK FOR H COLUMN
C
C  SEARCH I/O 13 FOR ELEMENT DATA FOR THIS MOLECULE.
C
      DO 20 I=1,NKIND
      idone(i)=-1
20    continue
30    READ (13,40,END=210) nint,fm,fel
40    format(/,i2,8x,2a1,f6.2)
      do 50 i=1,nint
      read (13,45,end=1000) tcoef(1,i),tcoef(2,i),numexp(i),
     *  (ex(j,i),j=1,8),hzero,(coef(j,i),j=1,10)
45    format(1x,2f10.3,i2,8f5.1,f17.3,/(5d16.8))
50    continue
      DO 60 Ij=1,NKIND
      J=JF(Ij)
      IF (fm(1).NE.FORMLA(1,Ij)) GO TO 60
      IF (fm(2).NE.FORMLA(2,Ij)) GO TO 60
      ii=ij
      GO TO 70
60    CONTINUE
      GO TO 30
70    LL=1
      hzero=-hzero/rr
      IF (refel) GO TO 210
      DMLA=SEL(II)
      DNMLA=NMLA(J)
      if (idone(ii).lt.0) then
          idone(ii)=0
          DH0=DH0-HZERO*DMLA/DNMLA
          istart=0
      endif
      if (fm(1).ne.fms(1) .or. fm(2).ne.fms(2)) then
         if(fms(1).ne.blank) ntx=min(ntx,iend)
      endif
C
90    i1=idone(ii)+1
      L1=1
      DO 180 I=I1,NTX
      TT=t(i)
      HH=0.d0
      S=0.d0
      GH=0.d0
C
C   FIND INTERVAL FOR T
C
      do 100 LL=L1,nint
      if(TT.lt.tcoef(1,LL)-.0001 .or. TT.gt.tcoef(2,LL)+.0001)
     * go to 100
      if (istart.eq.0) then
         istart=i
         INIT=max(INIT,i)
      endif
      L=LL
      go to 110
100   continue
      go to 180
c
c           CALCULATE ELEMENT H AND G PROPERTIES
c
110   DO 120 J=1,numexp(L)
      EXL=EX(J,L)
      TEX=COEF(J,L)
      IF (EXL.NE.0.) TEX=TEX*TT**EXL
      IF (EXL.EQ.-1.) HH=HH+COEF(J,L)*DLOG(TT)/TT
      IF (EXL.NE.-1.) HH=HH+TEX/(EXL+1.D0)
      IF (EXL.EQ.0.) S=S+COEF(J,L)*DLOG(TT)
      IF (EXL.NE.0.) S=S+TEX/EXL
120   CONTINUE
      HH=HH+(COEF(9,L))/TT
      S=S+COEF(10,L)
      GH=S-HH
      iend=i
      idone(ii)=i
C
C  CALCULATE DELTA H AND DELTA F
C
170   DH(I)=DH(I)-HH*DMLA/DNMLA
      ALK(I)=ALK(I)-GH*DMLA/DNMLA
C
C  SET T INDEX IN MARK. IF THIS END NOT THE END OF THE DATA, IT
C     IS A TRANSITION POINT FOR THIS ELEMENT.
C  IK IS THE INDEX FOR THE TRANSITION POINTS.
C
      if(fm(1).eq.fms(1) .and. fm(2).eq.fms(2)) then
         MARK(IK)=I
         PT(IK)=Tcoef(1,L)
         ename(1,ik)=fm(1)
         ename(2,ik)=fm(2)
         IF (IK.LT.30) IK=IK+1
         fms(1)=blank
      endif
180   CONTINUE
190      fms(1)=fm(1)
         fms(2)=fm(2)
      go to 30
C
210   DH0=DH0/1000.D0
      IK=IK-1
      IF (hasind.and.asindt.eq.0.)DH0= DH0 + term/1000.D0
215   DO 220 I=1,NKIND
      IF (idone(I).lt.0) CHDH=BLANK
220   CONTINUE
      IF (refel) CHDH=ZERO
C
C  PRINT TABLES
C
      ntx = min(ntx,iend)
      kdone=.false.
      DO 530 NTABLE=1,2
      IF (.NOT.TOUT(NTABLE+1)) GO TO 530
      DH0=DH0*ruse(ntable)
c
      WRITE (6,230) ((NAME(I,J),I=1,3),J=1,nna)
230   FORMAT ('0','THERMODYNAMIC FUNCTIONS CALCULATED FROM COEFFICIENTS
     *FOR ',9A4,/60X,2(3X,6A4))
238   WRITE (6,240)
240   FORMAT ('0     T         CP        H-H298        S      -(G-H298)/
     1T      H        DELTA H      LOG K')
245   IF (NTABLE.EQ.1)WRITE(6,250)
      IF (NTABLE.EQ.2)WRITE(6,255)
250   FORMAT ('    DEG-K    J/MOL-K      KJ/MOL     J/MOL-K     J/MOL-K
     *     KJ/MOL      KJ/MOL'/)
255   FORMAT ('    DEG-K   CAL/MOL-K    KCAL/MOL   CAL/MOL-K   CAL/MOL-K
     *    KCAL/MOL    KCAL/MOL'/)
      IF (.NOT.HASIND.or.asindt.ne.0.0) GO TO 290
      H0=(term*RUSE(NTABLE))/1000.D0
      HH0=-enth29/1000.d0
      if(ntable.eq.2)hh0=hh0/4.184D0
      IF (CHDH.EQ.BLANK .OR. CHDH.EQ.ZERO) GO TO 280
      WRITE (6,275) HH0,H0,DH0
275   FORMAT (6X,'0',8X,'0.',F15.3,7X,'0.',8X,'INFINITE',F11.3,F12.3
     * ,4X,'INFINITE')
      GO TO 290
280   WRITE (6,285) HH0,H0,CHDH
285   FORMAT (6X,'0',8X,'0.',F15.3,7X,'0.',8X,'INFINITE',F11.3,
     * 8X,A4,4X,'INFINITE')
290   H=0.
      H298=htform
      if(ntable.eq.2)H298=H298/4.184d0
      IF (HASIND) GO TO 298
      CHH=BLANK
      CHDH=BLANK
C
298   DO 420 I=1,NT
      IT=T(I)
      VFM2(3)=FI6
      VFM2(4)=F3X
      IF (DMOD(T(I),1.D0).EQ.0.) GO TO 300
      TI=T(I)
      VFM2(3)=F9
      VFM2(4)=FP2
300   RT=RUSE(NTABLE)*T(I)
      H=HHRT(I)*RT/1000.D0
      CP=CPR(I)*RUSE(NTABLE)
      S=SR(I)*RUSE(NTABLE)
      HH=(HHRT(I)*RT-H298)/1000.D0
      GH=(GHRT(I)*RUSE(NTABLE)+H298/T(I))
C
302   VFM2(2)=F2B
       IF (IK.EQ.0) GO TO 320
       DO 310 IX=1,IK
       IF (MARK(IX).EQ.I) VFM2(2)=FSTAR
310   CONTINUE
320   BZ=ZERO
      IF (CHDH.EQ.ZERO)GO TO 350
      BZ=BLANK
      IF ((I.GT.NTX.AND.NTX.NE.0).OR.I.LT.INIT) GO TO 350
      if (kdone) then
          DH(I)=DH(I)/4.184D0
      else
          DH(I)=(DH(I)+HHRT(I))*RT/1000.D0
          ALK(I)=(GHRT(I)+ALK(I))/2.3025851D0
      endif
330   VFM2(11)=F12
      if ( hasind) then
         VFM2(12)=FP3
         VFM2(13)=F12
         VFM2(14)=FP4
         IF(VFM2(3).EQ.F9)WRITE (6,VFM2)TI,CP,HH,S,GH,H,DH(I),ALK(I)
         IF(VFM2(3).EQ.FI6)WRITE (6,VFM2)IT,CP,HH,S,GH,H,DH(I),ALK(I)
      else
         VFM2( 9)=F12
         VFM2(10)=FP3
         VFM2(12)=FP4
         IF(VFM2(3).EQ.F9)WRITE (6,VFM2)TI,CP,S,H,DH(I),ALK(I)
         IF(VFM2(3).EQ.FI6)WRITE (6,VFM2)IT,CP,S,H,DH(I),ALK(I)
      endif
      GO TO 400
350   VFM2(11)=F8X
      VFM2(12)=FA4
      if (hasind) then
         VFM2(13)=F8X
         VFM2(14)=FA4
         IF(VFM2(3).EQ.F9)WRITE (6,VFM2)TI,CP,HH,S,GH,H,BZ,BZ
         IF(VFM2(3).EQ.FI6)WRITE (6,VFM2)IT,CP,HH,S,GH,H,BZ,BZ
      else
         VFM2( 9)=F8X
         VFM2(10)=FA4
         IF(VFM2(3).EQ.F9)WRITE (6,VFM2)TI,CP,S,H,BZ,BZ
         IF(VFM2(3).EQ.FI6)WRITE (6,VFM2)IT,CP,S,H,BZ,BZ
      endif
C
400   IF (DMOD(T(I),500.D0).NE.0.0) GO TO 420
      IF (T(I).GT.10000.D0.AND.DMOD(T(I),2500.D0).NE.0.) GO TO 420
      WRITE (6,410)
410   FORMAT (1H )
420   CONTINUE
C
C     WRITE FOOTNOTE
C
      IF (IK.EQ.0) GO TO 510
      WRITE (6,440)
440   FORMAT ( '0*A Change in phase of an assigned reference element has
     1 occurred between this ',/
     *'   temperature and the preceding one,')
      WRITE (6,450) (ENAME(1,I),ename(2,i),PT(I),I=1,IK)
450   FORMAT (' ',37x,'for ',2A1,' at ',F8.3,' K')
      LINES=LINES+4
510   WRITE (6,520)
520   format ('1    ')
      kdone = .true.
530   continue
      REWIND 13
1000  RETURN
      END
      SUBROUTINE TABLES
C
C  WRITE OUTPUT FILE FOR THE FIRST 3 TABLES OF THERMODYNAMIC FUNCTIONS.
C
      REAL*8 RUSE(3),TI,HH29,GH29,HH,AR,ART,S,CP,GH,G,H,HR,h298
C
      CHARACTER*4 VFM(20),XX,F14,F16,P3,P6,P7
C
      CHARACTER*4 BLANK,DATE,EFAZ,FAZ,FORMLA
      CHARACTER*4 NAME,SYMBOL
      COMMON /CHAR4/ BLANK,DATE(2),EFAZ(100),FAZ,
     1 FORMLA(2,5),NAME(4,9),SYMBOL(2,100)
C
      DOUBLE PRECISION asindt,ATMWT,COEF,CPR,GHRT,htform,enth29,
     1 HHRT,R,RR,SR,T,tcoef,term,WEIGHT,WORD,TTRANS
      COMMON /REAL8/ asindt,ATMWT(100),COEF(10,20),CPR(202),
     1 GHRT(202),htform,enth29,HHRT(202),R,RR,
     2 SR(202),T(202),tcoef(2,20),term,WEIGHT,WORD(4),TTRANS(6)
C
      COMMON /REAL4/ EX(8,20),sel(5),sels(5)
C
      COMMON /INTGR/ IC80,itr,JF(5),
     1 numexp(20),NEL,NKIND,NKINDS,
     2 NMAT(100),NMLA(100),NOATMS,NT,NTMP,NFZ,INA,NNA
C
      LOGICAL hasind,refel,TOUT
      COMMON /LOGCL/ hasind,refel,TOUT(3)
C
      DATA XX/',14X'/, F16/',F16'/, P7/'.7'/, P3/'.3'/
      DATA P6/'.6'/,F14/',F14'/,VFM(1)/'(1H '/, VFM(2)/',F12'/
      DATA VFM(3)/'.2'/, VFM(4)/',F12'/, VFM(5)/'.5'/, VFM(20)/')'/
C
10    RUSE(1)=R
      RUSE(2)=RR
      RUSE(3)=R
      H298=htform/rr
C
C  INITIALIZE FORMAT FOR THE FIRST TABLE
C
      DO 30 I=6,14,2
      VFM(I)=F14
      VFM(I+1)=P7
30    CONTINUE
      VFM(16)=F16
      VFM(17)=P7
C
      if (asindt.ne.0. .or. .not.hasind) then
         vfm(6) = blank
         vfm(7) = xx
         vfm(12) = blank
         vfm(13) = xx
      endif
      if (.not.hasind) then
         VFM(8)=BLANK
         VFM(9)=XX
         VFM(14)=BLANK
         VFM(15)=XX
      endif
C
C  PRINT MFIG TABLES
C
60    DO 190 NTABLE=1,3
      IF (.NOT.TOUT(NTABLE)) GO TO 190
      WRITE (6, 65) ((NAME(I,J),I=1,3),J=1,nna)
65    FORMAT ('0',5X,'THERMODYNAMIC FUNCTIONS CALCULATED FROM COEFFICIEN
     *TS FOR ',6A4,/60X,2(3X,6A4))
C
      if (ntable.ne.1) then
        IF (VFM(7).EQ.P7) VFM(7)=P3
        IF (VFM(9).EQ.P7) VFM(9)=P3
        VFM(11)=P6
        IF (VFM(13).EQ.P7) VFM(13)=P3
        IF (VFM(15).EQ.P7) VFM(15)=P3
        VFM(16)=F16
        VFM(17)=P3
      endif
C
      IF (.not.HASIND) WRITE (6,70)
70      FORMAT ('0           NO HZERO VALUE IS AVAILABLE')
        HR=term*RUSE(NTABLE)
        IHT=asindt
C
        IF(NTABLE.eq.1) then
          if(hasind) WRITE(6,80)IHT,term
80        FORMAT('0          Assigned H(T)/R at',I4,' K =',F12.3,' K'/)
          WRITE (6,85)
85        FORMAT ('0',8X,'T',8X,'CP/R',8X,'(H-H0)/RT    (H-H298)/RT',6X,
     1    'S/R',8X,'-(G-H0)/RT   -(G-H298)/RT',8X,'H/RT',12X,'-G/RT'//)
C
        ELSEIF(NTABLE.EQ.2) then
          if(hasind) WRITE(6,90) IHT,HR
90        FORMAT('0Assigned H(T) at',I4,' K =',F12.3,' J/mol')
          WRITE (6,100)
100       FORMAT ('0',8X,'T',10X,'CP',11X,'H-H0',10X,'H-H298',8X,'S',11X
     1    ,'-(G-H0)',6X,'-(G-H298)',11X,'H',14X,'-G')
          WRITE (6,115)
115       FORMAT(7X,'DEG-K',5X,'J/MOL-K',8X,'J/MOL',10X,'J/MOL',6X,
     *    'J/MOL-K',9X,'J/MOL',9X,'J/MOL',11X,'J/MOL',11X,'J/MOL'/)
C
        ELSEIF(NTABLE.EQ.3) then
          if(hasind) WRITE(6,120) IHT,HR
120       FORMAT('0Assigned H(T) at',I4,' K =',F12.3,' cal/mol')
          WRITE (6,100)
          WRITE(6,130)
130    FORMAT(7X,'DEG-K',4X,'CAL/MOL-K',7X,'CAL/MOL',7X,'CAL/MOL',4X,
     * 'CAL/MOL-K',7X,'CAL/MOL',7X,'CAL/MOL',9X,'CAL/MOL',9X,'CAL/MOL'/)
        endif
C
      VFM(18)=VFM(16)
      VFM(19)=VFM(17)
C
      DO 180 I=1,NT
      TI=T(I)
      IF(DABS(T(I)-298.15D0).LT..02)TI=298.15D0
      CP=CPR(I)
      H=HHRT(I)
      S=SR(I)
      G=GHRT(I)
      IF (NTABLE.NE.1)GO TO 140
      AR=1./TI
      ART=1.
      GO TO 150
140   AR=RUSE(NTABLE)
      ART=AR*TI
      CP=CP*AR
      H=H*ART
      S=S*AR
      G=G*ART
C
150   HH=H-term*AR
      GH=G+term*AR
165   HH29= H-H298*AR
      GH29= G+H298*AR
      if (hasind) then
         IF (asindt.ne.0.) then
            WRITE (6,VFM) T(I),CP,hh29,s,gh29,h,g
         else
            WRITE (6,VFM) T(I),CP,hh,hh29,s,gh,gh29,h,g
         endif
      else
         WRITE (6,vfm) t(i),cp,s,h,g
      endif
180   CONTINUE
      WRITE (6,185)
185   format ('1')
190   CONTINUE
      RETURN
      END
c
      BLOCK DATA
C
      CHARACTER*4 BLANK,DATE,EFAZ,FAZ,FORMLA
      CHARACTER*4 NAME,SYMBOL
      COMMON /CHAR4/ BLANK,DATE(2),EFAZ(100),FAZ,
     1 FORMLA(2,5),NAME(4,9),SYMBOL(2,100)
C
      DOUBLE PRECISION asindt,ATMWT,COEF,CPR,GHRT,htform,enth29,
     1 HHRT,R,RR,SR,T,tcoef,term,WEIGHT,WORD,TTRANS
      COMMON /REAL8/ asindt,ATMWT(100),COEF(10,20),CPR(202),
     1 GHRT(202),htform,enth29,HHRT(202),R,RR,
     2 SR(202),T(202),tcoef(2,20),term,WEIGHT,WORD(4),TTRANS(6)
C
      COMMON /REAL4/ EX(8,20),sel(5),sels(5)
C
      COMMON /INTGR/ IC80,itr,JF(5),
     1 numexp(20),NEL,NKIND,NKINDS,
     2 NMAT(100),NMLA(100),NOATMS,NT,NTMP,NFZ,INA,NNA
C
      LOGICAL hasind,refel,TOUT
      COMMON /LOGCL/ hasind,refel,TOUT(3)
C
C    CONSTANTS FROM THE 1986 CODATA RECOMMENDED VALUES OF THE FUND.
C    PHYSICAL CONSTANTS; J OF RES NAT BUR STDS; COHEN AND TAYLOR;
C    VOL 92, P 85, 1987.
C
      DATA RR/8.314510D0/
      DATA SYMBOL/  'E',' ','H',' ','D',' ','H','E','L','I','B','E','B',
     1  ' ','C',' ','N',' ','O',' ','F',' ','N','E','N','A','M','G','A',
     2  'L','S','I','P',' ','S',' ','C','L','A','R','K',' ','C','A','S',
     3  'C','T','I','V',' ','C','R','M','N','F','E','C','O','N','I','C',
     4  'U','Z','N','G','A','G','E','A','S','S','E','B','R','K','R','R',
     5  'B','S','R','Y',' ','Z','R','N','B','M','O','T','C','R','U','R',
     6  'H','P','D','A','G','C','D','I','N','S','N','S','B','T','E','I',
     7  ' ','X','E','C','S','B','A','L','A','C','E','P','R','N','D','P',
     8  'M','S','M','E','U','G','D','T','B','D','Y','H','O','E','R','T',
     9  'M','Y','B','L','U','H','F','T','A','W',' ','R','E','O','S','I',
     $  'R','P','T','A','U','H','G','T','L','P','B','B','I','P','O','A',
     $  'T','R','N','F','R','R','A','A','C','T','H','P','A','U',' ','N',
     $  'P','P','U','A','M','C','M','B','K','C','F'/
C
C  ATOMIC WEIGHTS FROM: IUPAC COMMISSION ON ATOMIC WEIGHTS AND ISOTOPIC
C    ABUNDANCES 1985; PURE APPL. CHEM.; 1986, 58, 1677. (VALUES WERE
C    TAKEN FROM TRC WHICH REFERRED TO ABOVE.)
C    NOTE: THE ATOMIC WEIGHTS OF THE RADIOACTIVE ELEMENTS TC, PM, PO,
C          AT, RN, FR, RA, AC, NP, PU, AM, CM, BK AND CF WERE
C          UPDATED TO THE 1986 VALUES USING THE WEIGHT OF THE ISOTOPE
C          HAVING THE LONGEST HALF-LIFE AS REPORTED BY TRC.
C
      DATA ATMWT/5.48579903D-04,1.00794D0,2.014102D0,4.002602D0,6.941D0,
     1 9.012182D0,10.811D0,12.011D0,14.00674D0,15.9994D0,18.9984032D0,
     2 20.1797D0,22.989768D0,24.305D0,26.981539D0,28.0855D0,30.973762D0,
     3 32.066D0,35.4527D0,39.948D0,39.0983D0,40.078D0,44.95591D0,
     4 47.88D0, 50.9415D0,51.9961D0,54.93805D0,
     5 55.847D0,58.9332D0,58.69D0,63.546D0,65.39D0,69.723D0,72.61D0,
     6 74.92159D0,78.96D0,79.904D0,83.80D0,85.4678D0,87.62D0,88.90585D0,
     7 91.224D0,92.90638D0,95.94D0,97.9072D0,101.07D0,102.9055D0,
     $ 106.42D0,
     8 107.8682D0,112.411D0,114.82D0,118.710D0, 121.75D0,127.6D0,
     9 126.90447D0,131.29D0,132.90543D0,137.327D0,138.9055D0,140.115D0,
     $ 140.90765D0,144.9127D0,145.D0,150.36D0,151.965D0,157.25D0,
     $ 158.92534D0,
     $ 162.50D0,164.93032D0,167.26D0,168.93421D0,173.04D0,174.967D0,
     $ 178.49D0,180.9479D0,183.85D0,186.207D0,190.2D0,192.22D0,
     $ 195.08D0,196.96654D0,200.59D0,204.3833D0,207.2D0,208.98037D0,
     $ 208.9824D0, 209.9871D0,
     $ 222.0176D0,223.0197D0,226.0254D0,227.0278D0,232.0381D0,
     $ 231.03588D0,238.0289D0,237.0482D0,244.0642D0,243.0614D0,
     $ 247.0703D0,247.0703D0,251.0587D0/
C
      DATA EFAZ/4*'G',4*'S',4*'G',6*'S',2*'G',12*'S','L',3*'S','S','G',
     1 17*'S','G',25*'S','L',5*'S','G',12*'S'/, BLANK/'    '/
C
      DATA NMAT/0,1,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,
     1 21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,
     2 41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,
     3 61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,
     4 81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98/
C
      DATA NMLA/ 1,2,2,5*1,3*2,7*1,2,17*1,2,17*1,2,45*1/
      END
