C     Five-step mechanism submitted to the 26th Symposium
C     version from March 1996
C     Citation:
C     J. Hewson and M. Bollig, 
C     Reduced Mechanisms for NOx Emissions from Hydrocarbon Diffusion Flames, 
C     Twenty-Sisxth Symposium (International) on Combustion,
C     1996
C
C     Global Reactions
C     I:        CH4 + 2 H + H2O -> 4H2 + CO             
C     II:       H2O + CO -> H2 + CO2                    
C     III:      2 H -> H2                               
C     IV:       O2 + 3 H2 -> 2 H2O + 2 H                
C     V:        H2 + 2 CO -> C2H2 + O2                  
C
C-----------------------------------------------------------------------
C ======= 5stepF.f =======
C-----------------------------------------------------------------------
      SUBROUTINE PRODRATES( CDOT, W, K, C, M, TEMP,
     &    PRESSURE )
C-----------------------------------------------------------------------
C     THIS SUBROUTINE COMPUTES RATES OF PRODUCTION CDOT
C     IN [KMOLE/(M^3S)]. THE PARAMETERS W ( REACTION RATE ),
C     K ( RATE COEFFICIENT ) AND M ( THIRD BODY CONCENTRATIONS ) ARE
C     JUST WORK SPACE FOR THIS FUNCTION.
C     C CONTAINS THE CONCENTRATIONS OF NON STEADY STATE SPECIES IN
C     [KMOLE/M^3] AND IS WORKSPACE FOR THE STEADY STATE CONCENTRATIONS,
C     WHICH ARE COMPUTED IN THIS FUNCTION.
C     TEMP IS THE TEMPERATURE IN [K] AND
C     PRESSURE IS THE PRESSURE IN [PA].
C     CALLED FUNCTIONS ARE 'GETLINDRATECOEFF', 'COMPSTEADYSTATES',
C     'CTCHZERO'
C-----------------------------------------------------------------------
      IMPLICIT NONE
      DOUBLE PRECISION CDOT(10), W(89), K(89),
     &    C(28), M(7), TEMP, PRESSURE
      INTEGER I
      DOUBLE PRECISION GETLINDRATECOEFF, LT
      DOUBLE PRECISION RGAS, RT
      PARAMETER ( RGAS = 8314.34 )
      DOUBLE PRECISION K34F0, K34FINF
      DOUBLE PRECISION FC34
      DOUBLE PRECISION K34B0, K34BINF
      DOUBLE PRECISION K36F0, K36FINF
      DOUBLE PRECISION FC36
      DOUBLE PRECISION K36B0, K36BINF
      DOUBLE PRECISION K51F0, K51FINF
      DOUBLE PRECISION FC51
      DOUBLE PRECISION K51B0, K51BINF
      DOUBLE PRECISION K58F0, K58FINF
      DOUBLE PRECISION FC58
      DOUBLE PRECISION K58B0, K58BINF
      DOUBLE PRECISION KN45F0, KN45FINF
      DOUBLE PRECISION FCN45
      DOUBLE PRECISION KN45B0, KN45BINF
      INCLUDE '5stepF.h'
C
      LT = DLOG( TEMP )
      RT = RGAS * TEMP 
C
      DO 100 I=1,89
        K(I)=A(I)*DEXP(N(I)*LT-E(I)/RT)
 100  CONTINUE
C
      K34F0 = 6.2570000000D+17 * DEXP(-1.8 * LT)
      K34FINF = 2.1080000000D+11
      FC34 =  0.577000 * DEXP ( -TEMP / 2370.00 )
      K(R34F) = GETLINDRATECOEFF( TEMP, PRESSURE, K34F0, K34FINF
     &    , FC34 )
      K34B0 = 6.5851871268D+22 * DEXP(-1.8 * LT - 439072242.3 / RT)
      K34BINF = 2.2185671189D+16 * DEXP(-439072242.3 / RT)
      K(R34B) = GETLINDRATECOEFF( TEMP, PRESSURE, K34B0, K34BINF
     &    , FC34 )
      K36F0 = 1.2720000000D+35 * DEXP(-7 * LT - 11560000 / RT)
      K36FINF = 1.8130000000D+10
      FC36 =  0.380000 * DEXP ( -TEMP / 73.0000 ) 
     &     + 0.620000 * DEXP ( -TEMP / 1180.00 )
      K(R36F) = GETLINDRATECOEFF( TEMP, PRESSURE, K36F0, K36FINF
     &    , FC36 )
      K36B0 = 5.0084110374D+41 * DEXP(-7 * LT - 384780198.3 / RT)
      K36BINF = 7.1385607004D+16 * DEXP(-373220198.3 / RT)
      K(R36B) = GETLINDRATECOEFF( TEMP, PRESSURE, K36B0, K36BINF
     &    , FC36 )
      K51F0 = 1.1870000000D+39 * DEXP(-7.5 * LT - 190400000 / RT)
      K51FINF = 2.0000000000D+14 * DEXP(-166290000 / RT)
      FC51 = 0.35
      K(R51F) = GETLINDRATECOEFF( TEMP, PRESSURE, K51F0, K51FINF
     &    , FC51 )
      K51B0 = 6.2450000000D+35 * DEXP(-7.5 * LT - 27500000 / RT)
      K51BINF = 1.0530000000D+11 * DEXP(-3390000 / RT)
      K(R51B) = GETLINDRATECOEFF( TEMP, PRESSURE, K51B0, K51BINF
     &    , FC51 )
      K58F0 = 1.0000000000D+13 * DEXP(-126000000 / RT)
      K58FINF = 1.3000000000D+13 * DEXP(-167000000 / RT)
      FC58 =  0.411 * DEXP( -73.4 / TEMP ) + DEXP( -TEMP / 422.8 )
      K(R58F) = GETLINDRATECOEFF( TEMP, PRESSURE, K58F0, K58FINF
     &    , FC58 )
      K58B0 = 1.5950000000D+10 * DEXP(27390000 / RT)
      K58BINF = 2.0730000000D+10 * DEXP(-13610000 / RT)
      K(R58B) = GETLINDRATECOEFF( TEMP, PRESSURE, K58B0, K58BINF
     &    , FC58 )
      KN45F0 = 3.9800000000D+11 * DEXP(-237000000 / RT)
      KN45FINF = 1.6000000000D+12 * DEXP(-262000000 / RT)
      FCN45 = 1.0
      K(RN45F) = GETLINDRATECOEFF( TEMP, PRESSURE, KN45F0, KN45FINF
     &    , FCN45 )
      KN45B0 = 4.7194802843D+06 * DEXP(-73764195.3 / RT)
      KN45BINF = 1.8972785063D+07 * DEXP(-98764195.3 / RT)
      K(RN45B) = GETLINDRATECOEFF( TEMP, PRESSURE, KN45B0, KN45BINF
     &    , FCN45 )
C

      M(MM1) = 3 * C(SCH4) + 6.5 * C(SH2O)
     &     + C(SH2) + 0.75 * C(SCO)
     &     + 0.4 * C(SO2) + 0.4 * C(SN2)
      M(MM6) = C(SCH4) + C(SH)
     &     + C(SH2O) + C(SH2)
     &     + C(SCO) + C(SCO2)
     &     + C(SO2) + C(SC2H2)
     &     + C(SN2) 
      M(MM3) = C(SCH4) + C(SH)
     &     + 10 * C(SH2O) + 2 * C(SH2)
     &     + C(SCO) + C(SCO2)
     &     + 2 * C(SO2) + C(SC2H2)
     &     + 2 * C(SN2) 
      M(MM8) = C(SCH4) + C(SH)
     &     + C(SH2O) + C(SH2)
     &     + C(SCO) + C(SCO2)
     &     + C(SO2) + C(SC2H2)
     &     + C(SN2) 
      M(MM7) = C(SCH4) + C(SH)
     &     + 18.6 * C(SH2O) + C(SH2)
     &     + C(SCO) + C(SCO2)
     &     + 1.5 * C(SO2) + C(SC2H2)
     &     + 1.5 * C(SN2) 
C
      CALL COMPSTEADYSTATES( K, C, M )
C
      W(R1F) = K(R1F) * C(SO2) * C(SH)
      W(R1B) = K(R1B) * C(SO) * C(SOH)
      W(R2F) = K(R2F) * C(SH2) * C(SO)
      W(R2B) = K(R2B) * C(SH) * C(SOH)
      W(R3F) = K(R3F) * C(SH2) * C(SOH)
      W(R3B) = K(R3B) * C(SH) * C(SH2O)
      W(R4F) = K(R4F) * C(SOH) * C(SOH)
      W(R4B) = K(R4B) * C(SO) * C(SH2O)
      W(R5F) = K(R5F) * C(SH) * C(SH) * M(MM1)
      W(R5B) = K(R5B) * C(SH2) * M(MM1)
      W(R6F) = K(R6F) * C(SO) * C(SO) * M(MM1)
      W(R6B) = K(R6B) * C(SO2) * M(MM1)
      W(R7F) = K(R7F) * C(SH) * C(SOH) * M(MM1)
      W(R7B) = K(R7B) * C(SH2O) * M(MM1)
      W(R8F) = K(R8F) * C(SH) * C(SO2) * M(MM1)
      W(R8B) = K(R8B) * C(SHO2) * M(MM1)
      W(R9) = K(R9) * C(SHO2) * C(SH)
      W(R10) = K(R10) * C(SHO2) * C(SH)
      W(R11) = K(R11) * C(SHO2) * C(SH)
      W(R12) = K(R12) * C(SHO2) * C(SO)
      W(R13F) = K(R13F) * C(SHO2) * C(SOH)
      W(R13B) = K(R13B) * C(SO2) * C(SH2O)
      W(R20F) = K(R20F) * C(SCO) * C(SOH)
      W(R20B) = K(R20B) * C(SH) * C(SCO2)
      W(R21) = K(R21) * C(SCO) * C(SHO2)
      W(R22) = K(R22) * C(SCO) * C(SO) * M(MM1)
      W(R25) = K(R25) * C(SCH) * C(SO2)
      W(R26) = K(R26) * C(SCH) * C(SCO2)
      W(R27F) = K(R27F) * C(SCH) * C(SH2O)
      W(R27B) = K(R27B) * C(SCH2OH)
      W(R28F) = K(R28F) * C(SCHO) * M(MM1)
      W(R28B) = K(R28B) * C(SH) * C(SCO) * M(MM1)
      W(R35F) = K(R35F) * C(S3XCH2) * C(SH)
      W(R35B) = K(R35B) * C(SH2) * C(SCH)
      W(R39) = K(R39) * C(S3XCH2) * C(SCH3)
      W(R40) = K(R40) * C(S3XCH2) * C(SO2)
      W(R41) = K(R41) * C(S3XCH2) * C(SO2)
      W(R42F) = K(R42F) * C(S1XCH2) * M(MM1)
      W(R42B) = K(R42B) * C(S3XCH2) * M(MM1)
      W(R44F) = K(R44F) * C(S1XCH2) * C(SH2)
      W(R44B) = K(R44B) * C(SH) * C(SCH3)
      W(R46) = K(R46) * C(SCH2O) * C(SH)
      W(R48) = K(R48) * C(SCH2O) * C(SOH)
      W(R53) = K(R53) * C(SCH3) * C(SO)
      W(R34F) = K(R34F) * C(SCH3) * C(SH)
      W(R34B) = K(R34B) * C(SCH4)
      W(R55) = K(R55) * C(SCH3) * C(SOH)
      W(R60) = K(R60) * C(SCH3) * C(SCH3)
      W(R36F) = K(R36F) * C(SCH3) * C(SCH3)
      W(R36B) = K(R36B) * C(SC2H6)
      W(R62) = K(R62) * C(SCH3O) * M(MM1)
      W(R63) = K(R63) * C(SCH3O) * C(SH)
      W(R69) = K(R69) * C(SCH2OH) * M(MM1)
      W(R70) = K(R70) * C(SCH2OH) * C(SH)
      W(R71) = K(R71) * C(SCH2OH) * C(SO2)
      W(R84F) = K(R84F) * C(SCH4) * C(SH)
      W(R84B) = K(R84B) * C(SCH3) * C(SH2)
      W(R86F) = K(R86F) * C(SCH4) * C(SOH)
      W(R86B) = K(R86B) * C(SCH3) * C(SH2O)
      W(R90F) = K(R90F) * C(SCH3OH)
      W(R90B) = K(R90B) * C(SOH) * C(SCH3)
      W(R91) = K(R91) * C(SCH3OH) * C(SH)
      W(R93) = K(R93) * C(SCH3OH) * C(SOH)
      W(R97) = K(R97) * C(SCH3O) * M(MM0)
      W(R105F) = K(R105F) * C(SHCCO) * C(SH)
      W(R105B) = K(R105B) * C(SCO) * C(S3XCH2)
      W(R109) = K(R109) * C(SC2H2) * C(SO2)
      W(R111) = K(R111) * C(SC2H2) * C(SO)
      W(R112) = K(R112) * C(SC2H2) * C(SO)
      W(R51F) = K(R51F) * C(SC2H3)
      W(R51B) = K(R51B) * C(SC2H2) * C(SH)
      W(R129F) = K(R129F) * C(SC2H4) * M(MM1)
      W(R129B) = K(R129B) * C(SH2) * C(SC2H2) * M(MM1)
      W(R131F) = K(R131F) * C(SC2H4) * C(SH)
      W(R131B) = K(R131B) * C(SH2) * C(SC2H3)
      W(R134F) = K(R134F) * C(SC2H4) * C(SOH)
      W(R134B) = K(R134B) * C(SH2O) * C(SC2H3)
      W(R58F) = K(R58F) * C(SC2H5)
      W(R58B) = K(R58B) * C(SC2H4) * C(SH)
      W(R164) = K(R164) * C(SC2H6) * C(SH)
      W(R166) = K(R166) * C(SC2H6) * C(SOH)
      W(R170F) = K(R170F) * C(SC2H6) * C(SCH3)
      W(R170B) = K(R170B) * C(SCH4) * C(SC2H5)

      REV36 = ( K(R164) * C(SH) + K(R166) * C(SOH)  ) / 
     &    CTCHZERO( K(R36B) + K(R164) * C(SH) + K(R166) * C(SOH) ) 

      W(RI) = K(RI) = W(R84F) + W(R86F) + W(R34B) 
     & - W(R84B) - W(R86B) - W(R34F) 
      W(RII) = K(RII) = W(R20F) - W(R20B) 
      W(RIII) = K(RIII) = W(R8F) - W(R8B) + W(R41) 
     & + W(R34F) - W(R34B) - W(R55) + W(R63) 
     & + W(R70) + W(R71) 
      W(RIV) = K(RIV) = W(R1F) - W(R1B) + W(R36F)*REV36 
     & + W(R60) + W(R25) - W(R35F) + W(R35B) 
     & - W(R44F) + W(R44B) 
      W(RV) = K(RV) = W(R36F)*REV36 + W(R60) - W(R109) 
     & - W(R111) - W(R112) 
 


      CDOT(SCH4) = - W(RI) 
      CDOT(SH) = - 2 * W(RI) - 2 * W(RIII) + 2 * W(RIV)
      CDOT(SH2O) = - W(RI) - W(RII) + 2 * W(RIV)
      CDOT(SH2) = 4 * W(RI) + W(RII) + W(RIII)
     &     - 3 * W(RIV) - W(RV) 
      CDOT(SCO) = W(RI) - W(RII) - 2 * W(RV)
      CDOT(SCO2) = W(RII) 
      CDOT(SO2) = - W(RIV) + W(RV)
      CDOT(SC2H2) = W(RV) 
      END
C
C
      DOUBLE PRECISION FUNCTION GETLINDRATECOEFF( TEMP, PRESSURE,
     &    K0, KINF, FC )
C
      IMPLICIT NONE
      DOUBLE PRECISION TEMP, PRESSURE, K0, KINF, FC
      DOUBLE PRECISION R
      PARAMETER (R = 8314.34)
C     [J / kmole K]
      DOUBLE PRECISION NTMP
      DOUBLE PRECISION KL
      DOUBLE PRECISION F
      DOUBLE PRECISION CONC
C
      CONC = PRESSURE / ( R * TEMP )
      NTMP = 0.75 - 1.27 * DLOG10( FC )
      K0 = K0 * CONC / KINF
      KL = K0 / ( 1.0 + K0 )
      F = DLOG10( K0 ) / NTMP
      F = FC ** ( 1.0 / ( F * F + 1.0 ) )
      GETLINDRATECOEFF = KINF * F * KL
      END
C
C
      DOUBLE PRECISION FUNCTION CTCHZERO( X )
C
      IMPLICIT NONE
      DOUBLE PRECISION X
      IF ( X.EQ.0.0D0 ) THEN
         X = 1.0D-20
      END IF
      CTCHZERO = X
      END
C
C
      DOUBLE PRECISION FUNCTION SOLVQDRT( A, B, C )
C-----------------------------------------------------------------------
C	SOLVES THE QUADRATIC EQUATION A x^2 + B x - C = 0
C-----------------------------------------------------------------------
C
      IMPLICIT NONE
C
      DOUBLE PRECISION RAD, A, B, C
C
      IF ( A.EQ.0.0D0 ) THEN
        IF ( B.EQ.0.0D0 ) THEN
          WRITE(*,*) '#WARNING: INVALID ARGUMENTS IN FUNCTION SOLVQDRT '
          SOLVQDRT = 0.0
        ELSE
          SOLVQDRT = C / B
        END IF
      ELSE
C
        B = B / A
        C = C / A
C
        RAD = 0.25D0 * B * B + C
        IF ( RAD.GE.0.0D0 ) THEN
          RAD = DSQRT( RAD ) + 0.5 * B
          IF ( RAD.NE.0.0D0 ) THEN
            SOLVQDRT = C / RAD
          ELSE
            SOLVQDRT = -0.5D0 * B + DSQRT( 0.25 * B * B + C )
          ENDIF
        ELSE
          SOLVQDRT = 0.0D0
        ENDIF
      END IF
      END
C
C
      SUBROUTINE GETMOLARMASS( MM )
C-----------------------------------------------------------------------
C	FILLS 'MM' WITH SPECIES MOLAR MASS IN KG/KMOLE
C-----------------------------------------------------------------------
C
      IMPLICIT NONE
C
      DOUBLE PRECISION MM(89)
      INCLUDE '6stepF.h'
C
      MM(SCH4) =  1.60430000e+01
      MM(SH) =  1.00800000e+00
      MM(SH2O) =  1.80160000e+01
      MM(SH2) =  2.01600000e+00
      MM(SCO) =  2.80110000e+01
      MM(SCO2) =  4.40110000e+01
      MM(SO2) =  3.20000000e+01
      MM(SC2H2) =  2.60380000e+01
      MM(SN2) =  2.80140000e+01
      MM(SOH) =  1.70080000e+01
      MM(SO) =  1.60000000e+01
      MM(SHO2) =  3.30080000e+01
      MM(SCH) =  1.30190000e+01
      MM(SCHO) =  2.90190000e+01
      MM(SCH2OH) =  3.10350000e+01
      MM(S3XCH2) =  1.40270000e+01
      MM(SCH3) =  1.50350000e+01
      MM(SC2H4) =  2.80540000e+01
      MM(S1XCH2) =  1.40270000e+01
      MM(SCH2O) =  3.00270000e+01
      MM(SCH3O) =  3.10300000e+01
      MM(SC2H6) =  3.00700000e+01
      MM(SCH3OH) =  3.20430000e+01
      MM(SHCCO) =  4.10300000e+01
      MM(SC2H3) =  2.70460000e+01
      MM(SC2H5) =  2.90620000e+01

      END
C
C
      SUBROUTINE GETNAMES( NAMES )
C-----------------------------------------------------------------------
C	FILLS 'NAMES' WITH SPECIES IDENTIFIER/KG
C-----------------------------------------------------------------------
C
      IMPLICIT NONE
C
      CHARACTER *20 NAMES(40)
      INCLUDE '6stepF.h'
C
      NAMES(SCH4)='CH4                 '
      NAMES(SH)='H                   '
      NAMES(SH2O)='H2O                 '
      NAMES(SH2)='H2                  '
      NAMES(SCO)='CO                  '
      NAMES(SCO2)='CO2                 '
      NAMES(SO2)='O2                  '
      NAMES(SC2H2)='C2H2                '
      NAMES(SN2)='N2                  '
      NAMES(SOH)='OH                  '
      NAMES(SO)='O                   '
      NAMES(SHO2)='HO2                 '
      NAMES(SCH)='CH                  '
      NAMES(SCHO)='CHO                 '
      NAMES(SCH2OH)='CH2OH               '
      NAMES(S3XCH2)='3-CH2               '
      NAMES(SCH3)='CH3                 '
      NAMES(SC2H4)='C2H4                '
      NAMES(S1XCH2)='1-CH2               '
      NAMES(SCH2O)='CH2O                '
      NAMES(SCH3O)='CH3O                '
      NAMES(SC2H6)='C2H6                '
      NAMES(SCH3OH)='CH3OH               '
      NAMES(SHCCO)='HCCO                '
      NAMES(SC2H3)='C2H3                '
      NAMES(SC2H5)='C2H5                '

C
      END
C
C
      SUBROUTINE COMPSTEADYSTATES( K, C, M )
C-----------------------------------------------------------------------
C     THIS SUBROUTINE COMPUTES THE STEADY STATE CONCENTRATIONS FROM
C     THE CONCENTRATIONS OF COMPUTED SPECIES AND RATE COEFFICIENTS.
C     CONCENTRATIONS OF COMPUTED SPECIES MAY NOT BE ALTERED.
C-----------------------------------------------------------------------
      IMPLICIT NONE
      DOUBLE PRECISION CTCHZERO
      DOUBLE PRECISION K(89), C(40), M(6)
      INCLUDE '6stepF.h'
      INCLUDE '6stepF.ss'
      END
C
C
      SUBROUTINE COMPTHERMODATA( H, CP, T )
C-----------------------------------------------------------------------
C     THIS FUNCTION COMPUTES ENTHALPY 'H' AND HEAT CAPACITY 'CP' AS
C     FUNCTION OF TEMPERATURE T FOR ALL NON STEADY STATE SPECIES
C     IN UNITS [J/KG] and [J/KG K], RESPECTIVELY.
C     THE PARAMETER H AND CP SHOULD PROVIDE WORKSPACE OF LENGTH 10
C-----------------------------------------------------------------------
      IMPLICIT NONE
      DOUBLE PRECISION H(10), CP(10), T
      INCLUDE '6stepF.h'
C
      IF (T.GT.1.0D3) THEN
      H(SCH4) =  5.18253444e+02 * (
     &      T * (  1.68347880e+00 + T * (  5.11861800e-03
     &      + T * ( -1.29170953e-06 + T * (  1.69639622e-10
     &      + T * ( -9.00684620e-15 ) ) ) ) ) -1.00807870e+04 )
      CP(SCH4) =  5.18253444e+02 * (
     &       1.68347880e+00 + T * (  1.02372360e-02 
     &      + T * ( -3.87512860e-06 + T * (  6.78558490e-10
     &      + T * ( -4.50342310e-14 ) ) ) ) )
      H(SH) =  8.24835317e+03 * (
     &      T * (  2.50000000e+00 + T * (  0.00000000e+00
     &      + T * (  0.00000000e+00 + T * (  0.00000000e+00
     &      + T * (  0.00000000e+00 ) ) ) ) ) +  2.54716270e+04 )
      CP(SH) =  8.24835317e+03 * (
     &       2.50000000e+00 + T * (  0.00000000e+00 
     &      + T * (  0.00000000e+00 + T * (  0.00000000e+00
     &      + T * (  0.00000000e+00 ) ) ) ) )
      H(SH2O) =  4.61497558e+02 * (
     &      T * (  2.61104720e+00 + T * (  1.57815650e-03
     &      + T * ( -3.09951460e-07 + T * (  3.33288450e-11
     &      + T * ( -1.49378702e-15 ) ) ) ) ) -2.98681670e+04 )
      CP(SH2O) =  4.61497558e+02 * (
     &       2.61104720e+00 + T * (  3.15631300e-03 
     &      + T * ( -9.29854380e-07 + T * (  1.33315380e-10
     &      + T * ( -7.46893510e-15 ) ) ) ) )
      H(SH2) =  4.12417659e+03 * (
     &      T * (  3.06670950e+00 + T * (  2.87368775e-04
     &      + T * (  4.64610633e-09 + T * ( -6.37087950e-12
     &      + T * (  5.81971480e-16 ) ) ) ) ) -8.65474120e+02 )
      CP(SH2) =  4.12417659e+03 * (
     &       3.06670950e+00 + T * (  5.74737550e-04 
     &      + T * (  1.39383190e-08 + T * ( -2.54835180e-11
     &      + T * (  2.90985740e-15 ) ) ) ) )
      H(SCO) =  2.96824105e+02 * (
     &      T * (  3.02507810e+00 + T * (  7.21344250e-04
     &      + T * ( -1.87694260e-07 + T * (  2.54645325e-11
     &      + T * ( -1.38219032e-15 ) ) ) ) ) -1.42683500e+04 )
      CP(SCO) =  2.96824105e+02 * (
     &       3.02507810e+00 + T * (  1.44268850e-03 
     &      + T * ( -5.63082780e-07 + T * (  1.01858130e-10
     &      + T * ( -6.91095160e-15 ) ) ) ) )
      H(SCO2) =  1.88915044e+02 * (
     &      T * (  4.45362280e+00 + T * (  1.57008435e-03
     &      + T * ( -4.26136833e-07 + T * (  5.98499175e-11
     &      + T * ( -3.33806640e-15 ) ) ) ) ) -4.89669610e+04 )
      CP(SCO2) =  1.88915044e+02 * (
     &       4.45362280e+00 + T * (  3.14016870e-03 
     &      + T * ( -1.27841050e-06 + T * (  2.39399670e-10
     &      + T * ( -1.66903320e-14 ) ) ) ) )
      H(SO2) =  2.59823125e+02 * (
     &      T * (  3.61221390e+00 + T * (  3.74265830e-04
     &      + T * ( -6.60688233e-08 + T * (  8.43725200e-12
     &      + T * ( -4.78147480e-16 ) ) ) ) ) -1.19781510e+03 )
      CP(SO2) =  2.59823125e+02 * (
     &       3.61221390e+00 + T * (  7.48531660e-04 
     &      + T * ( -1.98206470e-07 + T * (  3.37490080e-11
     &      + T * ( -2.39073740e-15 ) ) ) ) )
      H(SC2H2) =  3.19315616e+02 * (
     &      T * (  4.43677040e+00 + T * (  2.68801955e-03
     &      + T * ( -6.37605567e-07 + T * (  8.21594725e-11
     &      + T * ( -4.31341900e-15 ) ) ) ) ) +  2.56676640e+04 )
      CP(SC2H2) =  3.19315616e+02 * (
     &       4.43677040e+00 + T * (  5.37603910e-03 
     &      + T * ( -1.91281670e-06 + T * (  3.28637890e-10
     &      + T * ( -2.15670950e-14 ) ) ) ) )
      H(SN2) =  2.96792318e+02 * (
     &      T * (  2.85328990e+00 + T * (  8.01106400e-04
     &      + T * ( -2.09789643e-07 + T * (  2.86025550e-11
     &      + T * ( -1.56114930e-15 ) ) ) ) ) -8.90080930e+02 )
      CP(SN2) =  2.96792318e+02 * (
     &       2.85328990e+00 + T * (  1.60221280e-03 
     &      + T * ( -6.29368930e-07 + T * (  1.14410220e-10
     &      + T * ( -7.80574650e-15 ) ) ) ) )
      ELSE
      H(SCH4) =  5.18253444e+02 * (
     &      T * (  7.78741480e-01 + T * (  8.73834200e-03
     &      + T * ( -9.27803000e-06 + T * (  7.62427000e-09
     &      + T * ( -2.44786140e-12 ) ) ) ) ) -9.82522850e+03 )
      CP(SCH4) =  5.18253444e+02 * (
     &       7.78741480e-01 + T * (  1.74766840e-02 
     &      + T * ( -2.78340900e-05 + T * (  3.04970800e-08
     &      + T * ( -1.22393070e-11 ) ) ) ) )
      H(SH) =  8.24835317e+03 * (
     &      T * (  2.50000000e+00 + T * (  0.00000000e+00
     &      + T * (  0.00000000e+00 + T * (  0.00000000e+00
     &      + T * (  0.00000000e+00 ) ) ) ) ) +  2.54716270e+04 )
      CP(SH) =  8.24835317e+03 * (
     &       2.50000000e+00 + T * (  0.00000000e+00 
     &      + T * (  0.00000000e+00 + T * (  0.00000000e+00
     &      + T * (  0.00000000e+00 ) ) ) ) )
      H(SH2O) =  4.61497558e+02 * (
     &      T * (  4.16772340e+00 + T * ( -9.05748500e-04
     &      + T * (  1.98237627e-06 + T * ( -1.21730053e-09
     &      + T * (  3.05839820e-13 ) ) ) ) ) -3.02899690e+04 )
      CP(SH2O) =  4.61497558e+02 * (
     &       4.16772340e+00 + T * ( -1.81149700e-03 
     &      + T * (  5.94712880e-06 + T * ( -4.86920210e-09
     &      + T * (  1.52919910e-12 ) ) ) ) )
      H(SH2) =  4.12417659e+03 * (
     &      T * (  3.35535140e+00 + T * (  2.50680720e-04
     &      + T * ( -7.66896933e-08 + T * ( -1.19763310e-10
     &      + T * (  9.70451700e-14 ) ) ) ) ) -1.01916260e+03 )
      CP(SH2) =  4.12417659e+03 * (
     &       3.35535140e+00 + T * (  5.01361440e-04 
     &      + T * ( -2.30069080e-07 + T * ( -4.79053240e-10
     &      + T * (  4.85225850e-13 ) ) ) ) )
      H(SCO) =  2.96824105e+02 * (
     &      T * (  3.26245170e+00 + T * (  7.55970450e-04
     &      + T * ( -1.29391840e-06 + T * (  1.39548605e-09
     &      + T * ( -4.94990240e-13 ) ) ) ) ) -1.43105390e+04 )
      CP(SCO) =  2.96824105e+02 * (
     &       3.26245170e+00 + T * (  1.51194090e-03 
     &      + T * ( -3.88175520e-06 + T * (  5.58194420e-09
     &      + T * ( -2.47495120e-12 ) ) ) ) )
      H(SCO2) =  1.88915044e+02 * (
     &      T * (  2.27572470e+00 + T * (  4.96103615e-03
     &      + T * ( -3.46970433e-06 + T * (  1.71667170e-09
     &      + T * ( -4.23456020e-13 ) ) ) ) ) -4.83731410e+04 )
      CP(SCO2) =  1.88915044e+02 * (
     &       2.27572470e+00 + T * (  9.92207230e-03 
     &      + T * ( -1.04091130e-05 + T * (  6.86668680e-09
     &      + T * ( -2.11728010e-12 ) ) ) ) )
      H(SO2) =  2.59823125e+02 * (
     &      T * (  3.78371350e+00 + T * ( -1.51168170e-03
     &      + T * (  3.31642503e-06 + T * ( -2.45472752e-09
     &      + T * (  6.60636500e-13 ) ) ) ) ) -1.06381070e+03 )
      CP(SO2) =  2.59823125e+02 * (
     &       3.78371350e+00 + T * ( -3.02336340e-03 
     &      + T * (  9.94927510e-06 + T * ( -9.81891010e-09
     &      + T * (  3.30318250e-12 ) ) ) ) )
      H(SC2H2) =  3.19315616e+02 * (
     &      T * (  2.01356220e+00 + T * (  7.59522300e-03
     &      + T * ( -5.38772967e-06 + T * (  2.26974795e-09
     &      + T * ( -3.82549200e-13 ) ) ) ) ) +  2.61244430e+04 )
      CP(SC2H2) =  3.19315616e+02 * (
     &       2.01356220e+00 + T * (  1.51904460e-02 
     &      + T * ( -1.61631890e-05 + T * (  9.07899180e-09
     &      + T * ( -1.91274600e-12 ) ) ) ) )
      H(SN2) =  2.96792318e+02 * (
     &      T * (  3.70441770e+00 + T * ( -7.10937650e-04
     &      + T * (  9.55679733e-07 + T * ( -3.00722125e-10
     &      + T * ( -2.79093540e-15 ) ) ) ) ) -1.06407950e+03 )
      CP(SN2) =  2.96792318e+02 * (
     &       3.70441770e+00 + T * ( -1.42187530e-03 
     &      + T * (  2.86703920e-06 + T * ( -1.20288850e-09
     &      + T * ( -1.39546770e-14 ) ) ) ) )
      END IF
C
      END
