      SUBROUTINE TYPE64 (TIME,XIN,OUT,T,DTDT,PAR,INFO,ICNTRL,*)
C************************************************************************
C*    Copyright ASHRAE       A Toolkit for Primary HVAC System Energy
C*                           Calculation
C***********************************************************************
C*    SUBROUTINE:            TYPE64 (ENGINEID)
C*
C*    LANGUAGE:              FORTRAN 77
C*
C*    PURPOSE:               Parameter identification based on 
C*                           gas engine performances in full-load
C*                           regime.
C***********************************************************************
C*    INPUT VARIABLES:
C*    Fratio       Fuel/air ratio                                    (-) 
C*    xin(1)                                                         (-)
C*    Ta           Air temperature for specified working conditions  (K)
C*    xin(2)                                                        (øC) 
C*    MfrW         Water mass flow rate for specified working     (kg/s)
C*                 conditions
C*    xin(3)                                                     (kg/hr) 
C*    Twsu         Supply water temperature for specified working    (K)
C*                 conditions
C*    xin(4)                                                        (øC)
C*    Wsh          Shalft power for specified working conditions     (W)
C*    xin(5)                                                     (kJ/hr)
C*    Qw           Global water heat transfer for specified working  (W)
C*                 conditions
C*    xin(6)                                                     (kJ/hr)
C*    Qwenv        Water-environment heat transfer for specified     (W)
C*                 working conditions
C*    xin(7)                                                     (kJ/hr)
C*    MfrFuel      Gas mass flow rate for specified working       (kg/s)
C*                 conditions
C*    xin(8)                                                     (kg/hr)
C*    Wlo          Mechanical losses and auxiliary consumption for   (W)
C*                 specified working conditions
C*    xin(9)                                                     (kJ/hr)
C*    N            Rotation speed for specified working conditions (1/s)
C*    xin(10)                                                      (rpm)
C*    pa           Air pressure for specified working conditions    (Pa)
C*    xin(11)                                                      (atm)
C*
C*    OUTPUT VARIABLES
C*    Athroat    Nozzle throat area                               (m**2)
C*    out(1)                                                      (m**2)
C*    EffiInt    Internal efficiency                                 (-)
C*    out(2)                                                         (-)
C*    Tlo        Torque associated with the mechanical losses      (N*m)
C*                 and the auxiliary consumptions
C*    out(3)                                                       (N*m)
C*    AUgwNoEng  Gas-water heat transfer coefficient in nominal    (W/K)
C*                 conditions
C*    out(4)                                                  (kJ/hr/øC)
C*    MfrGasNom  Flue gas mass flow rate associated with the      (kg/s)
C*                 identified value of the gas-water heat transfer
C*                 coefficient
C*    out(5)                                                     (kg/hr)
C*    AUwenvEng  Water-environment heat transfer coefficient       (W/K)
C*    out(6)                                                  (kJ/hr/øC)
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*                 In  these cases, the routine stops running;
C*                 otherwise this variable is equal to 0.
C*    out(7)                                                         (-)
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*
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                            throttling 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 cylinders supply     (m**3/kg)
C     p3           Pressure at the cylinders supply                 (Pa)
C     pcritic      Critic pressure                                  (Pa)
C     p2           Pressure at nozzle throat                        (Pa)
C     v1           Air specific volume                         (m**3/kg)
C     Wpumping     Pumping losses                                    (W)
C     Win          Internal power                                    (W)
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     Qw           Net heat transfer brought by water                (W)
C     Qgw          Flue gas-water heat transfer                      (W)
C     Qwenv        Water-environment heat transfer                   (W)
C     Twex         Exhaust water temperature                         (K)
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     ErrRel       Relative error                                    (-)
C     hgex         Gas enthalpy at the exhaust of the    (J/kg flue gas)
C                  flue gas-water heat exchanger
C     CPgas        Mean specific heat of the combustion         (J/kg/K)
C                  products
C     Tgex         Flue gas temperature at the exhaust of the        (K)
C                  gas-water heat exchanger
C
C     Sum1,Sum2,Jm1,Dhgex,hgcal1,hgcal,Tgsup,hgex1 and Tgexp are 
C     variables used in the Newton-Raphson method.
C***********************************************************************
      
      INTEGER*4 INFO,INFO99
      DOUBLE PRECISION XIN,OUT,XIN99,OUT99
      
      REAL Kmolp(5)
      REAL Ifuel,MfrW,MfrFuel,MfrGasNom,N,Nmin,Nmax,i

      DIMENSION PAR(2),XIN(11),OUT(7),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) 

      INFO(6)=7
      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 11 (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)
      Wsh=SNGL(xin(5)/3.6)
      Qw=SNGL(xin(6)/3.6)
      Qwenv=SNGL(xin(7)/3.6)
      MfrFuel=SNGL(xin(8)/3600.)
      Wlo=SNGL(xin(9)/3.6)
      N=SNGL(xin(10)/60)
      pa=SNGL(xin(11)*101325)

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

      i=par(1)
      Vs=par(2)

C2*** The gaseous fuel used is methane

      Ifuel=4
      CALL FUEL (Ifuel,Cweight,FLHV,Tr,Cfuel,*1)
      CALL LINKCK('TYPE64','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 flue gas mass flow rate

      MfrGasNom=MfrFuel*(1+1/Fratio)

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

      VfrCyl=i*N*Vs

C2*** Calculate the specific volume at the cylinders supply

      vCyl=VfrCyl/MfrGasNom

C1*** Calculate the pressure at the cylinders supply

      p3=RAir*Ta/vCyl

C2*** Calculate the critic pressure and the pressure at nozzle
C2*** throad

      Gm1G=(GammaAir-1)/GammaAir
      pcritic=pa*(2/(GammaAir+1))**(1/Gm1G)
      p2=MAX(pcritic,p3)

C1*** Calculate the throttling area

      v1=RAir*Ta/pa
      Athroat=MfrGasNom*v1/(SQRT(2*CpAir*Ta)*SQRT((p2/pa)**
     &          (2/GammaAir)*(1-(p2/pa)**Gm1G)))

C1*** Calculate the pumping losses

      Wpumping=i*N*Vs*(pa-p3)

C1*** Calculate the internal power and the internal efficiency

      Win=Wsh+Wpumping+Wlo
      EffiInt=Win/(MfrFuel*FLHV)

C1*** Calculate the torque associated with the mechanical losses
C1*** and the auxiliary consumptions

      Tlo=Wlo/(2*Pi*N)

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 TYPE99 (TIME,XIN99,OUT99,T,DTDT,PAR99,INFO99,ICNTRL,*2)
      CALL LINKCK('TYPE64','TYPE99 ',1,99)
2     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/MfrGasNom

C1*** Calculate the flue gas temperature at the supply of the 
C1*** gas-water heat exchanger ( Newton-Raphson method )

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('TYPE64','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*** Calculate the flue gas-water heat transfer

      Qgw=Qw+Qwenv

C2*** Calculate the exhaust water temperature and the water 
C2*** temperature at the flue gas-water heat exchanger exhaust

      Twex=Twsu+Qw/(CpWat*MfrW)
      Twexs=Twsu+Qgw/(CpWat*MfrW)

C1*** Calculate the water-environment heat transfer coefficient

      AUwenvEng=CpWat*MfrW*LOG((Twexs-Ta)/(Twex-Ta))

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('TYPE64','ENTHALP',1,99)
51       CONTINUE
         hgex1=hgex1+Kmolp(J)*hpi
60    CONTINUE

C2*** The exhaust gas enthalpy is expressed in J/kg flue 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=MfrGasNom*CPgas
      Crw=MfrW*CpWat

C1*** Determine the value of ErrDetec

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

      Fct=Tgex-Tgsu+Qgw/Crgas

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)
      Dfct=1+Qgw/MfrW*(1/(hgex-hgsu)+(Tgsu-Tgex)*Dhgex/(hgsu-hgex)**2)
      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 coefficient

      AUgwNoEng=Qgw*LOG((Tgsu-Twexs)/(Tgex-Twsu))/
     &     ((Tgsu-Twexs)-(Tgex-Twsu))
 
90    CONTINUE


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

      out(1)=DBLE(Athroat)
      out(2)=DBLE(EffiInt)
      out(3)=DBLE(Tlo)
      out(4)=DBLE(AUgwNoEng*3.6)
      out(5)=DBLE(MfrGasNom*3600.)
      out(6)=DBLE(AUwenvEng*3.6)
      out(7)=DBLE(ErrDetec)    
            
      RETURN 1

      END


SUBROUTINE TYPE99 (TIME,XIN,OUT,T,DTDT,PAR,INFO,ICNTRL,*)
C************************************************************************
C*    Copyright ASHRAE       A Toolkit for Primary HVAC System Energy
C*                           Calculation
C***********************************************************************
C*    SUBROUTINE:            TYPE99 (COMBCH)
C*
C*    LANGUAGE:              FORTRAN 77
C*
C*    PURPOSE:               The purpose is to allow a very simplified
C*                           study of the combustion process that takes
C*                           place inside the combustion chamber.
C***********************************************************************
C*    INPUT VARIABLES
C*    Ifuel        Selection of the fuel                             (-)
C*                 If Ifuel
C*                    =1: Light fuel oil (liquid fuel)
C*                    =2: Heavy fuel oil (liquid fuel)
C*                    =3: Domestic gas oil (liquid fuel)
C*                    =4: Methane (gaseous fuel)
C*    xin(1)                                                         (-)
C*
C*    Choice       If Choice                                         (-)
C*                 =1: the fuel/air ratio is known
C*                 =2: the CO2 concentration in dry flue gas 
C*                     is known
C*                 =3: the O2 concentration in dry flue gas 
C*                     is known
C*    xin(2)                                                         (-)
C*   
C*    Val          This value is equal to the fuel/air ratio f ,     (-)
C*                 the CO2 concentration or the O2 concentration
C*                 according to the value of Choice
C*    xin(3)                                                         (-)
C*
C*    Tair         Air temperature                                   (K)
C*    xin(4)                                                        (øC)
C*
C*    Tfuel        Fuel temperature                                  (K)
C*    xin(5)                                                        (øC)
C*
C*
C*    OUTPUT VARIABLES
C*    Fratio       Actual fuel/air ratio                             (-)
C*    out(1)	                                                   (-)
C*
C*    Tadiab       Adiabatic temperature of the combustion products  (K)
C*    out(2)                                                        (øC)
C*
C*    Kmolp(2)     Number of kmol of O2                (kmol O2/kg fuel)
C*                 in the combustion products per kg of fuel
C*    out(3)                                           (kmol O2/kg fuel)
C*
C*    Kmolp(3)     Number of kmol of N2                (kmol N2/kg fuel)
C*                 in the combustion products per kg of fuel
C*    out(4)                                           (kmol N2/kg fuel)
C*
C*    Kmolp(4)     Number of kmol of CO2              (kmol CO2/kg fuel)
C*                 in the combustion products per kg of fuel
C*    out(5)                                          (kmol CO2/kg fuel)
C*
C*    Kmolp(5)     Number of kmol of HO2              (kmol HO2/kg fuel)
C*                 in the combustion products per kg of fuel
C*    out(6)                                          (kmol HO2/kg fuel)
C*
C*    hprod        Enthalpy of the combustion products       (J/kg fuel)
C*    out(7)                                                 (J/kg fuel)
C*
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:    The combustion process is assumed to be 
C                            adiabatic. The combustion reaction is
C                            assumed to be complete and dissociation is
C                            not taken into account.
C
C     DEVELOPER:             Philippe Ngendakumana
C                            Marc Grodent
C                            Jean-Pascal Bourdouxhe
C                            University of LiŠge, Belgium
C
C     DATE:                  March 1, 1995
C
C     SUBROUTINE CALLED:     ENTHALP
C                            FUEL
C                            LINKCK
C***********************************************************************
C     INTERNAL VARIABLES
C     Excess       Excess of air                                     (-)
C     Hweight      Weight of hydrogen in 1 kg of fuel               (kg)
C     ParC         Number of carbon atoms in equivalent hydrocarbon  (-)
C                  fuel
C     ParH         Number of hydrogen atoms in equivalent            (-)
C                  hydrocarbon fuel
C     xCO2         CO2 concentration in dry flue gas                 (-)
C     xO2          O2 concentration in dry flue gas                  (-)
C     O2st         Number of kmol of oxygen reacting with 1 kmol  (kmol)
C                  of fuel in a stoichiometric combustion process
C     Fratiost     Stoichiometric fuel/air ratio                     (-)
C     Qa           Heat transfer in the air preheater                (W)
C     Qf           Heat transfer in the fuel preheater               (W)
C     Qr           Heat transfer in the isothermal reactor           (W)
C     Qg           Heat transfer in the postheater                   (W)
C     MW           Array containing the molecular weights      (kg/kmol)
C                  of the species (H2,O2,N2,CO2,H2O)
C     ZC           Atomic number of carbon                           (-)
C     ParAir       Number of kmol of nitrogen per kmol of oxygen     (-)
C                  (air composition)
C     hO2          Enthalpy of oxygen per kmol of oxygen        (J/kmol)
C     hrO2         Enthalpy of oxygen per kmol of oxygen        (J/kmol)
C                  at the temperature Tr
C     hN2          Enthalpy of nitrogen per kmol of nitrogen    (J/kmol)
C     hrN2         Enthalpy of nitrogen per kmol of nitrogen    (J/kmol)
C                  at the temperature Tr
C     hKP          Enthalpy of species KP per kmol of           (J/kmol)
C                  species KP
C     hrKP         Enthalpy of species KP per kmol of           (J/kmol)
C                  species KP at the temperature Tr
C     Toler        Relative error tolerance                          (-)
C     ErrRel       Relative error                                    (-)
C     KP           Loop counter
C     IO           Integer replacing input Choice in the routine
C     Fct          Value of the function to be nullified     (J/kg fuel)
C
C     Tmin,Tmax,Qg1,Fct1,Qg2,Fct2,Fcti,Fctip1,Ti,Tip1,Tp and
C     T are variables used in the chord method
C***********************************************************************

      INTEGER*4 INFO
      DOUBLE PRECISION XIN,OUT
      
      REAL Kmolp(5),MW(5),Ifuel

      DIMENSION XIN(5),OUT(7),INFO(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
      
      INFO(6)=7

      DATA Toler/1E-5/
      DATA MW/2,32,28,44,18/
      DATA ZC/12/


C*** INPUTS 5 (converted in SI units)
C************
 
      Ifuel=SNGL(xin(1))
      Choice=SNGL(xin(2))
      Val=SNGL(xin(3))
      Tair=SNGL(xin(4)+273.15)
      Tfuel=SNGL(xin(5)+273.15)
  
      CALL FUEL (Ifuel,Cweight,FLHV,Tr,Cfuel,*1)
      CALL LINKCK('TYPE99','FUEL',1,99)
1     CONTINUE

      IO=INT(Choice)

C2*** Selection of the information available

      GO TO (10,20,30),IO

C2*** The fuel/air ratio is known

10    Fratio=Val
      GO TO 40

C2*** The CO2 concentration in dry flue gas is known

20    xCO2=Val
      GO TO 40

C2*** The O2 concentration in dry flue gas is known

30    xO2=Val

C2*** The only fuel constituents are carbon and hydrogen 

40    Hweight=1-Cweight

C2*** Calculate the number of carbon atoms in 
C2*** equivalent hydrocarbon fuel

      ParC=Cweight/ZC

C2*** Calculate the number of hydrogen atoms in 
C2*** equivalent hydrocarbon fuel

      ParH=Hweight

C2*** Calculate the number of kmol of nitrogen per kmol of oxygen
C2*** (air composition)

      ParAir=0.79/0.21

C2*** Calculate the number of kmol of oxygen reacting with 1 kmol
C2*** of fuel in a stoichiometric combustion process

      O2st=ParC+ParH/4

C2*** Calculate the stoichiometric fuel/air ratio

      Fratiost=1/(O2st*(MW(2)+ParAir*MW(3)))

C2*** Different cases according to the information available

      GO TO (50,60,70),IO

C2*** Calculate the excess of air

50    Excess=(Fratiost/Fratio)-1
      GO TO 80

60    Excess=(ParC-xCO2*(ParC+O2st*ParAir))/(xCO2*O2st*(1+ParAir))

C1*** Calculate the actual fuel/air ratio

      Fratio=Fratiost/(1+Excess)
      GO TO 80

70    Excess=-xO2*(ParC+O2st*ParAir)/(O2st*(xO2*(1+ParAir)-1))
      Fratio=Fratiost/(1+Excess)
    
C1*** Calculate the composition of the combustion products

80    Kmolp(2)=O2st*Excess
      Kmolp(3)=O2st*(1+Excess)*ParAir
      Kmolp(4)=ParC
      Kmolp(5)=ParH/2

C1*** Calculate the heat transfer in the fuel preheater

      Qf=Cfuel*(Tr-Tfuel)

C1*** Calculate the heat transfer in the air preheater

      CALL ENTHALP(Tair,2,hO2,*81)
      CALL LINKCK('TYPE99','ENTHALP',1,99)
81    CONTINUE
      CALL ENTHALP(Tair,3,hN2,*82)
82    CALL ENTHALP(Tr,2,hrO2,*83)
83    CALL ENTHALP(Tr,3,hrN2,*84)
84    Qa=O2st*(1+Excess)*((hrO2-hO2)+ParAir*(hrN2-hN2))

C1*** Calculate the heat transfer in the isothermal reactor

      Qr=FLHV

C1*** Calculate the adiabatic temperature by means of the chord method

      Tmin=700

C2*** For a given temperature calculate the heat transfer
C2*** in the postheater

      Qg1=0
      DO 90 KP=2,5
         CALL ENTHALP(Tmin,KP,hKP,*85)
85       CALL ENTHALP(Tr,KP,hrKP,*86)
86       Qg1=Qg1+(hKP-hrKP)*Kmolp(KP)
90    CONTINUE

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

      Fct1=(Qr-Qa-Qf)-Qg1

C2*** Same processes for another given temperature

      Tmax=3000
      Qg2=0
      DO 100 KP=2,5
         CALL ENTHALP(Tmax,KP,hKP,*91)
91       CALL ENTHALP(Tr,KP,hrKP,*92)
92       Qg2=Qg2+(hKP-hrKP)*Kmolp(KP)
100   CONTINUE
      Fct2=(Qr-Qa-Qf)-Qg2
      
      Fcti=Fct1
      Fctip1=Fct2
      Ti=Tmin
      Tip1=Tmax
      Tp=Tmax

C2*** New estimation of the adiabatic temperature

110   T=Ti+Fcti*(Ti-Tip1)/(Fctip1-Fcti)
      Qg=0
      DO 120 KP=2,5
         CALL ENTHALP(T,KP,hKP,*111)
111      CALL ENTHALP(Tr,KP,hrKP,*112)
112      Qg=Qg+(hKP-hrKP)*Kmolp(KP)
120   CONTINUE
      Fct=(Qr-Qa-Qf)-Qg

      ErrRel=ABS((T-Tp)/Tp)

C2*** If converged, leave loop

      IF (ErrRel.GT.Toler) THEN
         Fcti=Fctip1
         Fctip1=Fct
         Ti=Tip1
         Tip1=T
         Tp=T
         GO TO 110
      ENDIF

      Tadiab=T

C1*** Calculate the enthalpy of the combustion products
C1*** at the adiabatic temperature

      hprod=0
      DO 130 KP=2,5
        CALL ENTHALP(Tadiab,KP,hKP,*125)
125     hprod=hprod+Kmolp(KP)*hKP
130   CONTINUE

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

      out(1)=DBLE(Fratio)
      out(2)=DBLE(Tadiab-273.15)
      out(3)=DBLE(Kmolp(2))
      out(4)=DBLE(Kmolp(3))
      out(5)=DBLE(Kmolp(4))
      out(6)=DBLE(Kmolp(5))
      out(7)=DBLE(hprod)
 
      RETURN 1

      END

SUBROUTINE ENTHALP (Temp,I,Enthalpy,*)                                   
C************************************************************************
C*    Copyright ASHRAE       A Toolkit for Primary HVAC System Energy
C*                           Calculation
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:                  March 1, 1995
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***********************************************************************


      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 1                                                            
      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

SUBROUTINE FUEL (Ifuel,Cweight,FLHV,Tr,Cfuel,*)
C************************************************************************
C*    Copyright ASHRAE       A Toolkit for Primary HVAC System Energy
C*                           Calculation
C***********************************************************************
C*    SUBROUTINE:            FUEL
C*
C*    LANGUAGE:              FORTRAN 77
C*
C*    PURPOSE:               Selection of the properties of a given fuel
C***********************************************************************
C*    INPUT VARIABLE:
C*    Ifuel        Selection of the fuel                             (-)
C*                 If Ifuel
C*                    =1: Light fuel oil (liquid fuel)
C*                    =2: Heavy fuel oil (liquid fuel)
C*                    =3: Domestic gas oil (liquid fuel)
C*                    =4: Methane (gaseous fuel)
C***********************************************************************
C*    OUTPUT VARIABLES:
C*    Cweight      Weight of carbon in 1kg of fuel                  (kg)
C*    FLHV         The fuel lower heating value                   (J/kg)
C*    Tr           Reference temperature at which the FLHV is        (K)
C*                 evaluated
C*    Cfuel        Fuel specific heat                           (J/kg/K)
C***********************************************************************
C     DEVELOPER:             Jean-Pascal Bourdouxhe
C                            Marc Grodent
C                            University of LiŠge, Belgium
C
C     DATE:                  March 1, 1995
C***********************************************************************
    
      REAL Ifuel

      IF (Ifuel.EQ.1) THEN

         Cweight=0.88
         FLHV=40910E3
         Tr=298.15
         Cfuel=1840

      ENDIF

      IF (Ifuel.EQ.2) THEN

         Cweight=0.89
         FLHV=40430E3
         Tr=298.15
         Cfuel=1840

      ENDIF

      IF (Ifuel.EQ.3) THEN

         Cweight=0.86
         FLHV=42770E3
         Tr=298.15
         Cfuel=1880

      ENDIF

      IF (Ifuel.EQ.4) THEN

         Cweight=0.75
         FLHV=49997E3
         Tr=289.15
         Cfuel=2183

      ENDIF

      RETURN 1
      END

SUBROUTINE LINKCK(ENAME1,ENAME2,ILINK,LNKTYP)

C***************************************************************************
C  THIS SUBROUTINE WAS WRITTEN FOR TRNSYS 14.0 LINK CHECKING - THIS ROUTINE
C  IS CALLED BY OTHER SUBROUTINES WHEN AN UNLINKED SUBROUTINE HAS BEEN
C  FOUND.  LINKCK IS NEEDED IN ORDER TO AVOID PUTTING COMMON BLOCKS LUNITS
C  AND CONFIG IN THE TRNSYS TYPES - JWT -- 3/93
C***************************************************************************

      COMMON /LUNITS/ LUR,LUW,IFORM,LUK
      COMMON /CONFIG/ TRNEDT,PERCOM,HEADER,PRTLAB,LNKCHK,PRUNIT,IOCHEK,
     1                PRWARN
      COMMON /SIM/TIME0,TFINAL,DELT,IWARN


      CHARACTER*1 TRNEDT,PERCOM,HEADER,PRTLAB,LNKCHK,PRUNIT,IOCHEK,
     1            PRWARN
      CHARACTER*6 ENAME1,ENAME2
      INTEGER ILINK,LNKTYP

C     ILINK = 1 --> GENERATE AN ERROR MESSAGE AND STOP TRNSYS
C     ILINK = 2 --> GENERATE A WARNING BUT DON'T STOP TRNSYS
C     ILINK = 3 --> TRNSYS HAS FOUND AN UNLINKED TYPE - GENERATE AN ERROR AND
C                   STOP THE PROGRAM
C     ILINK = 4 --> WARN THE USER THAT A ROUTINE REQUIRES AN EXTERNAL FUNCTION
C     ENAME1    --> CALLING PROGRAM THAT NEEDED THE UNLINKED FILE
C     ENAME2    --> FILE THAT WAS NOT FOUND BY ENAME1 SUBROUTINE
C     LNKTYP    --> TYPE NUMBER THAT IS UNLINKED

      IF((LNKCHK.EQ.'Y').OR.(LNKCHK.EQ.'y')) THEN
        IF(ILINK.EQ.1) THEN
           WRITE(LUW,20) 104,ENAME1,ENAME2
           WRITE(LUW,15)
           CALL MYSTOP(104)
        ELSE IF(ILINK.EQ.2) THEN
           WRITE(LUW,20) 104,ENAME1,ENAME2
           IWARN=IWARN+1
        ELSE IF(ILINK.EQ.3) THEN
           WRITE(LUW,25) 105,LNKTYP,LNKTYP
           WRITE(LUW,15)
           CALL MYSTOP(105)
        ELSE IF(ILINK.EQ.4) THEN
           WRITE(LUW,35) LNKTYP,ENAME1,ENAME2
           IWARN=IWARN+1
        ELSE IF(ILINK.EQ.5) THEN
           WRITE(LUW,40) 105,LNKTYP,LNKTYP
           WRITE(LUW,15)
           CALL MYSTOP(105)
        ELSE
           WRITE(LUW,30) ENAME1
           IWARN=IWARN+1
        ENDIF
      ENDIF

15    FORMAT(//2X,47H*** SIMULATION TERMINATED WITH ERROR STATUS ***/)

20    FORMAT(//,1X,'***** ERROR *****',8X,'TRNSYS ERROR # ',I3,/1X,A6,
     1' REQUIRES THE FILE "',A6,'" WHICH WAS CALLED BUT NOT LINKED.',/1X
     1,'PLEASE LINK IN THE REQUIRED FILE AND RERUN THE SIMULATION.')

25    FORMAT(//,1X,'***** ERROR *****',8X,'TRNSYS ERROR # ',I3,/1X,
     1'TYPE ',I3,' WAS CALLED IN THE TRNSYS INPUT FILE BUT NOT LINKED.',
     1/1X,'LINK TYPE ',I3,' BEFORE RUNNING THIS SIMULATION.')

30    FORMAT(/1X,'*****WARNING*****',/1X,'THE LINKCK SUBROUTINE WAS CALL
     1ED WITH AN INVALID OPERAND.',/1X,'THE PROGRAM WHICH CALLED LINKCK
     1WITH THE IMPROPER OPERAND WAS ',A6,'.',/1X,'PLEASE MAKE SURE THAT
     1THE CALLING PROGRAM IS FIXED OR UNLINKED SUBROUTINES MAY ',/1X,'GO
     1 UNNOTICED.')

35    FORMAT(/1X,'*****WARNING*****',/1X,'UNIT ',I2,' ',A6,' REQUIRES TH
     1E SUBROUTINE ',A6,/1X,'MAKE SURE THAT THIS SUBROUTINE IS LINKED IN
     1 TO AVOID PROBLEMS. IT MAY ALREADY BE LINKED IN.',/)

40    FORMAT(//,1X,'***** ERROR *****',8X,'TRNSYS ERROR # ',I3,/1X,
     1'TYPE',I3,' WAS CALLED IN THE TRNSYS INPUT FILE BUT NOT LINKED.',
     1/1X,'A DUMMY TYPE SUBROUTINE WAS CALLED IN ITS PLACE. PLEASE LINK'
     1,/1X,'TYPE',I3,' BEFORE RUNNING THIS SIMULATION OR TURN OFF THE CH
     1ECK'/1X,'FOR UNLINKED SUBROUTINES OPTION IN THE CONFIGURATION FILE
     1.')

      RETURN
      END

