[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[TRNSYS-users] linking subroutines in trnsys16



Dear all,

I tried to use a type that calls different subroutines. I get the following error, indicating that it did not link the subroutines well. I do not see the solution to this error.
I attached the fortran files as well.

*** Fatal Error at time   :         0.025000
   Generated by Unit     : Not applicable or not available
   Generated by Type     : Not applicable or not available
TRNSYS Message 104 : The TRNSYS processor has reported that a subroutine was called that has not been found in the available TRNSYS libraries.
   Reported information  : Reported by LINKCK



***** ERROR *****        TRNSYS ERROR # 104
TYPE905 FUEL REQUIRES THE FILE "ENTHALP " WHICH WAS CALLED BUT NOT LINKED.
PLEASE LINK IN THE REQUIRED FILE AND RERUN THE SIMULATION.




Kind regards,
Leen



--
Ir.Jeroen Van der Veken
Afdeling Bouwfysica
Katholieke Universiteit Leuven
Kasteelpark Arenberg 40
3001 Heverlee
T: +32 16 32 13 47
F: +32 16 32 19 80
@: jeroen.vanderveken@bwk.kuleuven.be

NEW MAILADRES WITHOUT .AC !!
      SUBROUTINE TYPE905 (TIME,XIN,OUT,T,DTDT,PAR,INFO,ICNTRL,*)
	!DEC$ATTRIBUTES DLLEXPORT :: TYPE 905
C************************************************************************
C*    Copyright ASHRAE       A Toolkit for Primary HVAC System Energy
C*                           Calculation
C***********************************************************************
C*    SUBROUTINE:            TYPE905 (ENGIFLSI)
C*
C*    LANGUAGE:              FORTRAN 77
C*
C*    PURPOSE:               Calculates the shaft power when the gas
C*                           engine is running in steady-state regime.
C***********************************************************************
C*    INPUT VARIABLES:
C*    Fratio       Fuel/air ratio                                    (-) 
C*    xin(1)                                                         (-)
C*    Ta           Air temperature                                   (K)
C*    xin(2)                                                        (øC)
C*    MfrW         Water mass flow rate                           (kg/s)
C*    xin(3)                                                     (kg/hr)
C*    Twsu         Supply water temperature                          (K)
C*    xin(4)                                                        (øC)
C*    pa           Air pressure                                     (Pa)
C*    xin(5)                                                       (atm) 
C*    N            Rotation speed                                  (1/s)
C*    xin(6)                                                       (rpm)
C*
C*    OUTPUT VARIABLES
C*    CPgas        Mean specific heat of the combustion         (J/kg/K)
C*                 products
C*    out(1)                                                   (J/kg/øC)
C*    Twex         Exhaust water temperature                         (K)
C*    out(2)                                                        (øC)   
C*    Tgex         Flue gas temperature at the exhaust of the        (K)
C*                 gas-water heat exchanger
C*    out(3)                                                        (øC)
C*    Wsh          Shaft power                                       (W)
C*    out(4)                                                     (kJ/hr)
C*    MfrFuel      Gas mass flow rate                             (kg/s)
C*    out(5)                                                     (kg/hr)
C*    MfrGas       Flue gas mass flow rate                        (kg/s)
C*    out(6)                                                     (kg/hr)
C*    Effic        Gas engine efficiency                             (-)
C*    out(7)                                                         (-)
C*    ErrDetec     = 1: the ratio of water capacity flow rate to     (-)
C*                      flue gas capacity flow rate is too small ( <1 )
C*                 = 2: the rotation speed for specified working 
C*                      conditions is lower than the minimum rotation
C*                      speed (Nmin)
C*                 = 3: the rotation speed for specified working 
C*                      conditions is greater than the maximum
C*                      rotation speed (Nmax)
C*                 = 4: the routine does not converge
C*                 In  these cases, the routine stops running;
C*                 otherwise this variable is equal to 0.
C*    out(8)                                                         (-)
C*
C*    PARAMETERS
C*    i            Intermittency factor                              (-)
C*    par(1)                                                         (-)
C*    Vs           Swept volume corresponding to all the          (m**3)
C*                 cylinders
C*    par(2)                                                      (m**3)
C*    Athroat      Nozzle throat area                             (m**2)
C*    par(3)                                                      (m**2)
C*    EffiInt      Internal efficiency                               (-)
C*    par(4)                                                         (-)
C*    Tlo          Torque associated with the mechanical losses    (N*m)
C*                 and the auxiliary consumptions
C*    par(5)                                                       (N*m)
C*    AUgwNoEng    Gas-water heat transfer coefficient in nominal  (W/K)
C*                 conditions
C*    par(6)                                                  (kJ/hr/øC)
C*    MfrGasNom    Flue gas mass flow rate in nominal conditions  (kg/s)
C*    par(7)                                                     (kg/hr)
C*    AUwenvEng    Water-environment heat transfer coefficient     (W/K)
C*    par(8)                                                  (kJ/hr/øC)
C*
C*    AIR PROPERTIES
C*    CpAir        Air specific heat                            (J/kg/K)
C*    RAir         Air constant                                 (J/kg/K)
C*    GammaAir     Air isentropic coefficient                        (-)
C*
C*    WATER PROPERTIES
C*    CpWat        Specific heat of liquid water                (J/kg/K)
C*
C*    FUEL PROPERTIES
C*    Cweight      Weight of carbon in 1kg of fuel                  (kg)
C*    FLHV         Fuel lower heating value                       (J/kg)
C*    Tr           Reference temperature at which the FLHV is 
C*                 evaluated                                         (K)
C*    Cfuel        Fuel specific heat                           (J/kg/K)
C***********************************************************************
C     MAJOR RESTRICTIONS:    It is assumed that the water-environment
C                            heat transfer coefficient as well as the
C                            nozzle throat area, the internal efficiency
C                            ,the fuel/air ratio and the torque
C                            associated with the mechanical losses and
C                            the auxiliary consumptions are constant.
C                            Air-fuel mixing properties are the same as 
C                            for pure air.
C                            The gas-water heat transfer coefficient is 
C                            function of the flue gas mass flow rate.
C
C     DEVELOPER:             Jean Lebrun
C                            Marc Grodent
C                            Jean-Pascal Bourdouxhe
C                            Mark Nott
C                            University of Li?ge, Belgium
C
C     DATE:                  March 1, 1995
C
C     SUBROUTINES CALLED:    TYPE99 (COMBCH)
C                            ENTHALP
C                            FUEL
C                            LINKCK 
C***********************************************************************
C     INTERNAL VARIABLES
C     Nmin         Minimum rotation speed                          (1/s)
C     Nmax         Maximum rotation speed                          (1/s)
C     VfrCyl       Volume flow rate corresponding to all the    (m**3/s)
C                  cylinders
C     vCyl         Specific volume at the cylinder supply      (m**3/kg)
C     p3           Pressure at the cylinder supply                  (Pa)
C     pcritic      Critical pressure                                (Pa)
C     v1           Air specific volume                         (m**3/kg)
C     Wpumping     Pumping losses                                    (W)
C     Win          Internal power                                    (W)
C     Wlo          Power associated with the mechanical losses       (W)
C                  and the auxiliary consumptions
C     hg0f         Flue gas enthalpy at the exhaust of the    (J/kg gas)
C                  adiabatic combustion chamber
C     hg0          Flue gas enthalpy at the exhaust of   (J/kg flue gas)
C                  the adiabatic combustion chamber
C     hgsu         Flue gas enthalpy at the supply of    (J/kg flue gas)
C                  the gas-water heat exchanger
C     Tg0          Flue gas temperature at the exhaust of the        (K)
C                  adiabatic combustion chamber
C     Tgsu         Flue gas temperature at the supply of the         (K)
C                  gas-water heat exchanger
C     Twexs        Water temperature at the flue gas-water heat      (K)
C                  exchanger exhaust
C     TolRel       Relative error tolerance                          (-)
C     Crgas        Capacity flow rate of the combustion products   (W/K)
C     Crw          Water capacity flow rate                        (W/K)
C     Fct          Value of the function to be nullified             (K)
C     Dfct         Value of the first derivative                     (-)
C     Effgw        Effectiveness of the gas-water heat exchanger     (-)
C     ErrRel       Relative error                                    (-)
C     hgex         Gas enthalpy at the exhaust of the    (J/kg flue gas)
C                  flue gas-water heat exchanger
C     Qgw          Flue gas-water heat transfer                      (W)
C     Qwenv        Water-environment heat transfer                   (W)
C     AUgwEng      Gas-water heat transfer coefficient             (W/K)
C
C     Sum1,Sum2,Jm1,Dhgex,DCPgas,Dcrgas,Deffgw,hgcal1,hgcal,Tgsup,hgex1
C     and Tgexp are variables used in the Newton-Raphson method.
C***********************************************************************
      INCLUDE 'c:/Trnsys15/Include/param.inc'

      INTEGER*4 INFO,INFO99
      DOUBLE PRECISION XIN,OUT,XIN99,OUT99
      
      REAL Kmolp(5)
      REAL Ifuel,MfrW,MfrFuel,MfrGas,N,Nmin,Nmax,i,MfrGasNom

      DIMENSION PAR(8),XIN(6),OUT(8),INFO(15),
     &          XIN99(5),OUT99(7),INFO99(15)
      COMMON /LUNITS/ LUR,LUW,IFORM,LUK
      COMMON /SIM/ TIME0,TFINAL,DELT,IWARN
      COMMON /STORE/ NSTORE,IAV,S(5000)
      COMMON /CONFIG/ TRNEDT,PERCOM,HEADER,PRTLAB,LNKCHK,PRUNIT,IOCHEK,
     &                PRWARN
      
      COMMON/COMCP/PFCP(5,10) 

      ! Set the version information for TRNSYS
		if (INFO(7) == -2) then
			INFO(12) = 15
			return 1
		endif
	
	INFO(6)=8
      INFO99(6)=7

      DATA TolRel,Nmin,Nmax,CpWat,Pi/1E-05,8,85,4187,3.14159265359/
      DATA CpAir,RAir,GammaAir/1005,287.06,1.4/   

C*** INPUTS 6 (converted in SI units)
C************

      Fratio=SNGL(xin(1))
      Ta=SNGL(xin(2)+273.15)
      Mfrw=SNGL(xin(3)/3600.)
      Twsu=SNGL(xin(4)+273.15)
      pa=SNGL(xin(5)*101325)
      N=SNGL(xin(6)/60)

C*** PARAMETERS 8 (converted in SI units)
C****************

      i=par(1)
      Vs=par(2)
      Athroat=par(3)
      EffiInt=par(4)
      Tlo=par(5)
      AUgwNoEng=par(6)/3.6
      MfrGasNom=par(7)/3600.
      AUwenvEng=par(8)/3.6

C2*** The gaseous fuel used is methane

      Ifuel=4
      CALL FUEL (Ifuel,Cweight,FLHV,Tr,Cfuel,*1)
      CALL LINKCK('TYPE905','FUEL',1,99)
1     CONTINUE

C1*** Test on the value of the rotation speed

      IF (N.LT.Nmin) THEN
      ErrDetec=2
      GOTO 90
      ELSE
      IF (N.GT.Nmax) THEN
      ErrDetec=3
      GOTO 90
      ENDIF
      ENDIF

C1*** Calculate the critic pressure at the nozzle throat

      Gm1G=(GammaAir-1)/GammaAir
      pcritic=pa*(2/(GammaAir+1))**(1/Gm1g)

C2*** Calculate the volume flow rate corresponding to all
C2*** the cylinders

      VfrCyl=i*N*Vs

C2*** Calculate the pressure at the cylinder supply if we
C2*** assumed to be in sonic regime at the nozzle throat

      v1=Rair*Ta/pa
      MfrGas=Athroat/v1*SQRT(2*CpAir*Ta)*SQRT((pcritic/pa)**
     &       (2/GammaAir)*(1-(pcritic/pa)**Gm1G))
      vCyl=VfrCyl/MfrGas
      p3=RAir*Ta/vCyl

C2*** Compare the pressure at the cylinder supply with the 
C2*** critic pressure

      IF (p3.GT.pcritic) THEN

C2*** No sonic regime at the nozzle throat; calculate the pressure
C2*** at the cylinder supply by means of the Newton-Raphson method

C2*** First guess of the value of the pressure at the cylinder supply

      p3=0.9*pa

5     CONTINUE

C2*** Calculate the function to be nullified

      Fct=Athroat/v1*SQRT(2*CpAir*Ta)*SQRT((p3/pa)**
     &    (2/GammaAir)*(1-(p3/pa)**Gm1G))-VfrCyl*p3/(RAir*
     &     Ta)

C2*** Calculate the value of the first derivative

      prate=p3/pa
      Den=SQRT(prate**(2/GammaAir)*(1-prate**Gm1G))
      DFct=Athroat/v1*SQRT(2*CpAir*Ta)*(2/GammaAir*prate
     &     **((2-GammaAir)/GammaAir)-(GammaAir+1)/
     &     GammaAir*prate**(1/GammaAir))/(2*pa*Den)-
     &     VfrCyl/(RAir*Ta)

C2*** A new estimated value is calculated

      p3p=p3
      p3=p3-Fct/DFct
      ErrRel=ABS((p3-p3p)/p3p)

C2*** If converged, leave the loop

      IF (ErrRel.GT.TolRel) GOTO 5

      IF (p3.LT.pcritic) THEN
      ErrDetec=4
      GOTO 90
      ENDIF

C2*** Calculate the flue gas mass flow rate

      vCyl=RAir*Ta/p3
      MfrGas=VfrCyl/vCyl
      ENDIF

C1*** Calculate the gas mass flow rate

      MfrFuel=MfrGas*(Fratio/(1+Fratio))

C2*** Calculate the internal power

      Win=EffiInt*MfrFuel*FLHV

C2*** Calculate the pumping loss

      Wpumping=i*N*Vs*(pa-p3)
     
C2*** Calculate the mechanical losses and the auxiliary consumptions

      Wlo=Tlo*2*Pi*N

C1*** Calculate the shaft power

      Wsh=Win-Wpumping-Wlo
     
C1*** Calculate the gas-water heat transfer coefficient

      AUgwEng=AUgwNoEng*(MfrGas/MfrGasNom)**0.65

C1*** Calculate the adiabatic temperature, the fuel/air ratio as well as
C1*** the enthalpy (expressed in J/kg fuel) and composition of the 
C1*** combustion products

      xin99(1)=DBLE(Ifuel)
      xin99(2)=1
      xin99(3)=DBLE(Fratio)
      xin99(4)=DBLE(Ta-273.15)
      xin99(5)=DBLE(Ta-273.15)
      CALL TYPE908 (TIME,XIN99,OUT99,T,DTDT,PAR99,INFO99,ICNTRL,*7)
      CALL LINKCK('TYPE905','TYPE908 ',1,99)
7     CONTINUE
      Fratio=SNGL(out99(1))
      Tg0=SNGL(out99(2)+273.15)
      Kmolp(2)=SNGL(out99(3))
      Kmolp(3)=SNGL(out99(4))
      Kmolp(4)=SNGL(out99(5))
      Kmolp(5)=SNGL(out99(6))
      hg0f=SNGL(out99(7))

C2*** The flue gas enthalpy at the exhaust of the adiabatic 
C2*** combustion chamber is expressed in J/kg (flue gas)

      hg0=hg0f/(1+1/Fratio)

C1*** Calculate the flue gas enthalpy at the supply of the gas-water
C1*** heat exchanger 

      hgsu=hg0-Wsh/MfrGas

C1*** Calculate the flue gas temperature at the supply of the 
C1*** gas-water heat exchanger

C2*** First guess of the flue gas temperature at the heat exchanger
C2*** supply

      Tgsu=Tg0/2

10    hgcal1=0
      DO 20 J=2,5
        CALL ENTHALP (Tgsu,J,hpi,*11)
        CALL LINKCK('TYPE905','ENTHALP',1,99)
11      CONTINUE
        hgcal1=hgcal1+Kmolp(J)*hpi
20    CONTINUE
      hgcal=hgcal1/(1+1/Fratio)

C2*** Calculate the function to nullify

      Fct=hgcal-hgsu

C2*** Calculate the value of the first derivative

      Sum1=0
      DO 30 K=1,5 
      Sum2=0
      DO 40 J=1,10
      Sum2=Sum2+PFCP(K,J)*Tgsu**(J-1)
40    CONTINUE
      Sum1=Sum1+Kmolp(K)*Sum2
30    CONTINUE
      DFct=Sum1/(1+1/Fratio)

C2*** A new estimated value is calculated

      Tgsup=Tgsu
      Tgsu=Tgsu-Fct/DFct
      ErrRel=ABS((Tgsu-Tgsup)/Tgsup)

C2*** If converged, then leave the loop

      IF (ErrRel.GT.TolRel) GOTO 10
      
C2*** First guess of the exhaust flue gas temperature

      Tgex=Tgsu/2

C1*** Calculate the exhaust flue gas enthalpy (expressed in J/kg fuel)

50    hgex1=0
      DO 60 J=2,5
         CALL ENTHALP (Tgex,J,hpi,*51)
         CALL LINKCK('TYPE905','ENTHALP',1,99)
51       CONTINUE
         hgex1=hgex1+Kmolp(J)*hpi
60    CONTINUE

C2*** The exhaust flue gas enthalpy is expressed in J/kg gas

      hgex=hgex1/(1+1/Fratio)

C1*** Calculate the flue gas mean specific heat

      CPgas=(hgsu-hgex)/(Tgsu-Tgex)

C1*** Calculate a new estimated value of the exhaust flue gas
C1*** temperature by using the Newton-Raphson method

C2*** Calculate the value of the function to be nullified

      Crgas=MfrGas*CPgas
      Crw=MfrW*CpWat

C1*** Determine the value of ErrDetec

      IF (Crgas.GT.Crw) THEN
        ErrDetec=1
        GOTO 90
      ELSE
        ErrDetec=0
      ENDIF

      par1=EXP(-AUgwEng*(1/Crgas-1/Crw))
      Effgw=(1-par1)/(1-Crgas*par1/Crw)
      Fct=Effgw*(Tgsu-Twsu)-Tgsu+Tgex

C2*** Calculate the value of the first derivative

      Sum1=0
      DO 70 K=2,5
         Sum2=0
         DO 80 J=1,10
            Jm1=J-1
            Sum2=Sum2+PFCP(K,J)*Tgex**Jm1
80       CONTINUE
         Sum1=Sum1+Sum2*Kmolp(K)
70    CONTINUE
      Dhgex=Sum1/(1+1/Fratio)
      DCPgas=(hgsu-hgex-Dhgex*(Tgsu-Tgex))/(Tgsu-Tgex)**2
      DCrgas=MfrGas*DCPgas
      DEffgw=(AUgwEng*DCrgas*par1*(1/Crw-1/Crgas)/Crgas+DCrgas*par1*
     &        (1-par1)/Crw)/(1-(Crgas/Crw)*par1)**2
      Dfct=(Tgsu-Twsu)*DEffgw+1
      Tgexp=Tgex

C2*** The new estimated value is calculated

      Tgex=Tgex-Fct/Dfct
      ErrRel=ABS((Tgex-Tgexp)/Tgexp)

C2*** If converged, leave loop 

      IF (ErrRel.GT.TolRel) GO TO 50

C1*** Calculate the gas-water heat transfer
      
      Qgw=MfrGas*(hgsu-hgex)

C1*** Calculate the exhaust water temperature

      Twexs=Twsu+Qgw/(MfrW*CpWat)
      Twex=Ta+(Twexs-Ta)/EXP(AUwenvEng/(MfrW*CpWat))

C1*** Calculate the water-environment heat transfer

      Qwenv=MfrW*CpWat*(Twexs-Twex)

C1*** Calculate the gas engine efficiency

      Effic=Wsh/(MfrFuel*FLHV)

90    CONTINUE


C*** OUTPUTS 8 (converted in TRNSYS units)
C*************

      out(1)=DBLE(CPgas)
      out(2)=DBLE(Twex-273.15)
      out(3)=DBLE(Tgex-273.15)
      out(4)=DBLE(Wsh*3.6)
      out(5)=DBLE(MfrFuel*3600.)
      out(6)=DBLE(MfrGas*3600.)
      out(7)=DBLE(Effic)
      out(8)=DBLE(ErrDetec)    
            
      RETURN 1

      END

      SUBROUTINE ENTHALP (Temp,I,Enthalpy,*)                                   
C***********************************************************************
C*    SUBROUTINE:            ENTHALP
C*    
C*    LANGUAGE:              FORTRAN 77
C*   
C*    PURPOSE:               Calculate the enthalpy (J/kmol) of each
C*                           species (H2,O2,N2,CO2,H2O) at a given
C*                           temperature
C***********************************************************************
C*    INPUT VARIABLES
C*    Temp         Temperature at which enthalpy must be calculated  (K)
C*    I            Selection of the species to be considered         (-)
C*                 I=1: H2
C*                 I=2: O2
C*                 I=3: N2
C*                 I=4: CO2
C*                 I=5: H2O
C*
C*    OUTPUT VARIABLES
C*    Enthalpy     Enthalpy of the species                      (J/kmol)
C***********************************************************************
C     DEVELOPER:             Philippe Ngendakumana
C                            Marc Grodent
C                            University of Li?ge, Belgium
C
C     DATE:                  November 8, 1993
C
C     REFERENCE:             A. Brohmer and P. Kreuter
C                            FEV Motorentechnik GmbH & Co KG
C                            Aachen, Germany
C***********************************************************************
C     INTERNAL VARIABLES
C     PFCP         Array containing the coefficients used     (J/kmol/K)
C                  in the polynomial expressions     
C     Tref         Array containing the temperatures at which        (K)
C                  the reference enthalpies are calculated
C     href         Array containing the reference enthalpies    (J/kmol)
C     h            Enthalpy of species I                        (J/kmol)
C     J            Loop counter 
C***********************************************************************
!export this subroutine for its use in external DLLs.
!DEC$ATTRIBUTES DLLEXPORT :: ENTHALP

      COMMON/COMCP/PFCP(5,10)                                           
      COMMON/THREF/Tref(5),href(5)

      h=href(I)                                                         
      Enthalpy=0                                                            
      DO 10 J=1,10                                                      
      h=h+((PFCP(I,J)*Temp**J)-(PFCP(I,J)*Tref(I)**J))/J                
  10  CONTINUE                                                          
      Enthalpy=h                                                             

      RETURN                                                            
      END                                                               

      BLOCK DATA                                                        

      COMMON/COMCP/PFCP(5,10)                                           
      COMMON/THREF/Tref(5),href(5)
                                                                        
C1*** Coefficients are given for H2

      DATA PFCP(1,1),PFCP(1,2),PFCP(1,3),                              
     $PFCP(1,4),PFCP(1,5),PFCP(1,6),PFCP(1,7),                          
     $PFCP(1,8),PFCP(1,9),PFCP(1,10)/                                   
     $ 2.12183E+04, 4.90483E+01,-1.18908E-01, 1.50167E-04,              
     $-1.07285E-07, 4.66644E-11,-1.26418E-14, 2.08562E-18,              
     $-1.91864E-22, 7.54661E-27/                                        

C1*** Coefficients are given for O2

      DATA PFCP(2,1),PFCP(2,2),PFCP(2,3),                              
     $PFCP(2,4),PFCP(2,5),PFCP(2,6),PFCP(2,7),                          
     $PFCP(2,8),PFCP(2,9),PFCP(2,10)/                                   
     $ 3.12398E+04,-2.51025E+01, 9.50643E-02,-1.29283E-04,              
     $ 9.56020E-08,-4.25012E-11, 1.16866E-14,-1.94778E-18,              
     $ 1.80410E-22,-7.12717E-27/                                        

C1*** Coefficients are given for N2

      DATA PFCP(3,1),PFCP(3,2),PFCP(3,3),                              
     $PFCP(3,4),PFCP(3,5),PFCP(3,6),PFCP(3,7),                          
     $PFCP(3,8),PFCP(3,9),PFCP(3,10)/                                   
     $ 3.10052E+04,-1.65866E+01, 4.37297E-02,-4.10720E-05,              
     $ 2.08732E-08,-6.27548E-12, 1.11654E-15,-1.08777E-19,              
     $ 4.47487E-24, 0.E0         /                                      

C1*** Coefficients are given for CO2

      DATA PFCP(4,1),PFCP(4,2),PFCP(4,3),                              
     $PFCP(4,4),PFCP(4,5),PFCP(4,6),PFCP(4,7),                          
     $PFCP(4,8),PFCP(4,9),PFCP(4,10)/                                   
     $ 1.89318E+04, 8.20742E+01,-8.47204E-02, 5.92177E-05,              
     $-2.92546E-08, 1.01523E-11,-2.39525E-15, 3.62658E-19,              
     $-3.15882E-23, 1.19863E-27/                                        

C1*** Coefficients are given for H2O

       DATA PFCP(5,1),PFCP(5,2),PFCP(5,3),                              
     $PFCP(5,4),PFCP(5,5),PFCP(5,6),PFCP(5,7),                          
     $PFCP(5,8),PFCP(5,9),PFCP(5,10)/                                   
     $ 3.42084E+04,-1.04650E+01, 3.61342E-02,-2.73709E-05,              
     $ 1.12406E-08,-2.93883E-12, 5.25323E-16,-6.54907E-20,              
     $ 5.27765E-24,-2.04468E-28/                                        

C1*** Reference values are given for H2

      DATA Tref(1),href(1)/2.E3,6.144129E7/                                      

C1*** Reference values are given for O2

      DATA Tref(2),Href(2)/2.E3,6.7926643E7/                                     

C1*** Reference values are given for N2

      DATA Tref(3),Href(3)/2.E3,6.485353E7/                                      

C1*** Reference values are given for CO2

      DATA Tref(4),Href(4)/2.E3,-2.9253172E8/                                    

C1*** Reference values are given for H2O

      DATA Tref(5),Href(5)/2.E3,-1.5643141E8/
                                    
      END