SUBROUTINE TYPE78 (TIME,XIN,OUT,T,DTDT,PAR,INFO,ICNTRL,*) C*********************************************************************** C* Copyright ASHRAE. Toolkit for HVAC System Energy Calculations C************************************************************************ C* C* ********************WARNING********************* C* THE SAMPLE INPUT AND OUTPUT VALUES LISTED BELOW WERE C* GENERATED USING THE TRNSYS DEBUG PROGRAM. THEY WERE NOT C* INCLUDED WITH THE ORIGINAL DOCUMENTATION AS WITH THE OTHER C* COMPONENTS!!! C* C************************************************************************ C* C* SUBROUTINE: DXDOE C* C* LANGUAGE: FORTRAN 77 C* C* PURPOSE: Model the performance of a direct expansion C* air conditioner using the performance at C* ARI rating conditions (95F outdoor, 80F C* indoor dry bulb, 67F indoor wet bulb, and C* 450 cfm/ton indoor airflow). The algorithm C* corrects for off-design operation using C* a bypass factor analysis and correction C* factors from the DOE 2.1 program. Gives C* total capacity, efficiency, and sensible C* heat ratio. C*********************************************************************** C* MUST USE SI UNITS! C*********************************************************************** C* INPUT VARIABLES DESCRIPTION(UNITS) SAMPLE VALUES C* XIN(1) MAir Dry air mass flow rate(kg/s) 5.0 C* XIN(2) TAirEnt Entering air dry bulb temperature(C) 25.0 C* XIN(3) WAirEnt Entering air humidity ratio(-) .01 C* XIN(4) TAirAmb Outdoor ambient dry bulb temperature(C) 35.0 C* C* OUTPUT VARIABLES C* OUT(1) QTot Total heat transfer rate(W) 67200.0 C* OUT(2) Cop Coefficient of Performance(-) 2.4940 C* OUT(3) Shr Sensible heat ratio(-) .6189 C* OUT(4) ErrStat Error status indicator,0=ok,1=error(-) 0.0 C* C* PARAMETERS C* PAR(1) CapRated Total capacity at ARI rating conditions(W) 70000.0 C* PAR(2) CopRated COP at ARI rating conditions(-) 2.5 C* PAR(3) ShrRated Sensible heat ratio at ARI rating conditions(-) .75 C* PAR(4) UACp Coefficient in bypass factor relationship(kg/s) 2.2 C* Bypass = EXP(-UACp/MAir) C* PAR(5) MRated Airflow at ARI rating point(kg/s) 5.0785 C*********************************************************************** C MAJOR RESTRICTIONS: Does not include the effect of indoor fan C power consumption or temperature rise. C C DEVELOPER: Michael J. Brandemuehl, PhD, PE C University of Colorado at Boulder C C Hugh Henderson & Kannan Rengarajan C Florida Solar Energy Center C C DATE: May 1, 1992 C C INCLUDE FILES: acdx.inc C prop.inc C SUBROUTINES CALLED: TDB_W C FUNCTIONS CALLED: ENTHSAT C HUMRATIO C SATPRESS C XITERATE C C REVISION HISTORY: None C C REFERENCE: DOE-2. 1982. Engineer's Manual, Version C 2.1A. Lawrence Berkeley Laboratory and C Los Alamos National Laboratory. NTIS, C Reference DE83004575. C************************************************************************ C INTERNAL VARIABLES: C mFract Ratio of airflow to rated airflow (-) C hDelta Enthalpy difference (J/kg) C twbAirEnt Entering air wet bulb temperature at rating (C) C hAirEnt Entering air enthalpy at rating (J/kg) C tAmbF Outdoor ambient dry bulb temperature at rating (F) C twbF Entering air wet bulb temperature at rating (F) C hAdp Apparatus dew point enthalpy (J/kg) C tAdp Apparatus dew point temperature (C) C wAdp Apparatus dew point humidity ratio (-) C hAdpEst Estimated apparatus dew point enthalpy (J/kg) C bypass Bypass factor (-) C error Deviation of dependent variable in iteration C iter Iteration counter C icvg Iteration convergence flag C F1,F2 Previous values of dependent variable in XITERATE C X1,X2 Previous values of independent variable in XITERATE C************************************************************************ DOUBLE PRECISION XIN,OUT INTEGER INFO(15),IOPT,NI,NP,ND INTEGER Patm,CpAir,CpWat,CpLiq,CpVap,DViscAir, & DViscLiq,KAir,KLiq,RhoLiq,Hfg,RAir, & TKelMult,TAbsAdd,PaMult,PAbsAdd,ERRSTAT REAL Prop(16),MFRACT,SHR,PAR, &CAP_T_COR,CAP_FLOW_COR, &EIR_T_COR,EIR_FLOW_COR DIMENSION XIN(4),OUT(4),PAR(5) CHARACTER*3 YCHECK(4),OCHECK(4) COMMON /LUNITS/LUR,LUW,IFORM,LUK PARAMETER (Patm = 1) PARAMETER (CpAir = 2) PARAMETER (CpWat = 3) PARAMETER (CpVap = 4) PARAMETER (CpLiq = 5) PARAMETER (DViscAir = 6) PARAMETER (DViscLiq = 7) PARAMETER (KAir = 8) PARAMETER (KLiq = 9) PARAMETER (RhoLiq = 10) PARAMETER (Hfg = 11) PARAMETER (RAir = 12) PARAMETER (TKelMult = 13) PARAMETER (TAbsAdd = 14) PARAMETER (PaMult = 15) PARAMETER (PAbsAdd = 16) PROP(PATM) = 101325.0 PROP(CPAIR) = 1006.0 PROP(CPWAT) = 4186.0 PROP(CPVAP) = 1805.0 PROP(CPLIQ) = 4186.0 PROP(DVISCAIR) = .0000182 PROP(DVISCLIQ) = .00144 PROP(KAIR) = .026 PROP(KLIQ) = .604 PROP(RHOLIQ) = 998.0 PROP(HFG) = 2501000.0 PROP(RAIR) = 287.055 PROP(TKELMULT) = 1.0 PROP(TABSADD) = 273.15 PROP(PAMULT) = 1.0 PROP(PABSADD) = 0.0 MAIR = XIN(1) TAIRENT = XIN(2) WAIRENT = XIN(3) TAIRAMB = XIN(4) CAPRATED = PAR(1) COPRATED = PAR(2) SHRRATED = PAR(3) UACP = PAR(4) MRATED = PAR(5) C* DOE 2.1 Correction factor statement functions CAP_T_COR = 0.6003404 + 0.00228726*twbF & - 0.0000128*twbF*twbf + 0.00138975*tAmbF & - 0.0000806*tambF*tAmbF + 0.00014125*twbF*tAmbF CAP_FLOW_COR = 0.8 + 0.2*mFract EIR_T_COR = -0.9617787 + 0.04817751*twbF & - 0.0002311*twbF*twbF + 0.00324392*tAmbF & + 0.00014876*tAmbF*tAmbF - 0.00029545*twbF*tAmbF EIR_FLOW_COR = 1.156 - 0.1816*mFract DATA itmax/20/ DATA YCHECK/'MF2','TE1','DM1','TE1'/ DATA OCHECK/'PW2','DM1','DM1','DM1'/ ErrStat = 0 IOPT = -1 NI = 4 !CORRECT NUMBER OF INPUTS NP = 5 !CORRECT NUMBER OF PARAMETERS ND = 0 !CORRECT NUMBER OF DERIVATIVES IF (INFO(7).EQ.-1) THEN CALL TYPECK(IOPT,INFO,NI,NP,ND) C CHECKS TO SEE IF USER'S INFO MATCHES CORRECT NUMBERS CALL RCHECK(INFO,YCHECK,OCHECK) C CHECKS TO SEE IF INPUT AND OUTPUT UNITS MATCH ENDIF C1*** Calculate psychrometric properties CALL TDB_W (Prop,TAirEnt,WAirEnt,rH,hAirEnt,twbAirEnt, & tDP,rhoD,rhoM,ErrStat) twbF = twbAirEnt*1.8+32. tAmbF = TAirAmb*1.8+32. C1*** Use DOE correction factor functions to determine capacity and C1*** efficiency at actual operating conditions mFract = MAir/MRated QTot = CapRated*CAP_T_COR*CAP_FLOW_COR Cop = CopRated/EIR_T_COR/EIR_FLOW_COR C1*** Calculate bypass factor bypass = EXP(-UACp/MAir) C1*** Iteratively solve for apparatus dewpoint that gives enthalpy C1*** at ADP, known from enthalpy difference and bypass factor hDelta = QTot/MAir hAdp = hAirEnt - hDelta/(1.-bypass) C1*** Initialize iteration loop parameters tAdp = TAirEnt - hDelta/Prop(CpAir) DO 100 iter = 1,itmax C1*** Calculate estimated bypass factor as enthalpy ratio involving C1*** entering, leaving, and apparatus dewpoint states hAdpEst = ENTHSAT(Prop,tAdp) C1*** Estimate new apparatus dewpoint temperature by comparing C1*** estimated bypass factor with known bypass factor error = hAdp - hAdpEst tAdp = XITERATE(tAdp,error,X1,F1,X2,F2,iter,icvg) C1*** If converged, leave iteration loop. IF (icvg .EQ. 1) GO TO 110 C1*** Apparatus dewpoint temperature not converged, repeat calculation C1*** with new estimate of temperature. 100 CONTINUE C1*** Apparatus dewpoint temperature has not converged after maximum C1*** iterations. Print error message, set return error flag, and RETURN WRITE(LUW,1001) itmax 1001 FORMAT(/1X,'*** ERROR IN SUBOUTINE DXDOE ***'/ & 1X,' Apparatus dewpoint temperature has not ' & 'converged after ',I2,' iterations'/) ErrStat = 1 110 CONTINUE C1*** Determine sensible heat ratio from ratio of latent-to-total C1*** enthalpy between entering and apparatus dewpoint conditions wAdp = HUMRATIO (Prop(Patm),SATPRESS(Prop,tAdp)) Shr = (ENTHALPY(Prop,TAirEnt,wAdp)-hAdp)/(hAirEnt-hAdp) 999 CONTINUE OUT(1) = QTOT OUT(2) = COP OUT(3) = SHR OUT(4) = ERRSTAT RETURN 1 END SUBROUTINE TDB_W (Prop,TDB,W,RH,H,TWB,TDP,RhoD,RhoM,ErrStat) C*********************************************************************** C* Copyright ASHRAE. Toolkit for HVAC System Energy Calculations C*********************************************************************** C* SUBROUTINE: TDB_W C* C* LANGUAGE: FORTRAN 77 C* C* PURPOSE: Calculate psychrometric properties of C* moist air given dry bulb temperature C* and humidity ratio. C*********************************************************************** C* INPUT VARIABLES C* TDB Dry bulb temperature (C) C* W Humidity ratio (-) C* C* OUTPUT VARIABLES C* RH Relative humidity (-) C* H Enthalpy of moist air (J/kg) C* TWB Wet bulb temperature (C) C* TDP Dewpoint temperature (C) C* RhoD Dry air density (kg/m3) C* RhoM Moist air density (kg/m3) C* ErrStat Error flag (0=ok, 1=error) (-) C* C* PROPERTIES C* Patm Atmospheric pressure (Pa) C*********************************************************************** C MAJOR RESTRICTIONS: None C C DEVELOPER: Shauna Gabel C Michael J. Brandemuehl, PhD, PE C University of Colorado at Boulder C C DATE: January 1, 1992 C C INCLUDE FILES: prop.inc C SUBROUTINES CALLED: None C FUNCTIONS CALLED: SATPRESS C RELHUM C WETBULB C ENTHALPY C DEWPOINT C RHODRY C RHOMOIST C C REVISION HISTORY: None C C REFERENCE: 1989 ASHRAE Handbook - Fundamentals C*********************************************************************** INTEGER Patm,CpAir,CpWat,CpLiq,CpVap,DViscAir, & DViscLiq,KAir,KLiq,RhoLiq,Hfg,RAir, & TKelMult,TAbsAdd,PaMult,PAbsAdd REAL Prop(16) PARAMETER (Patm = 1) PARAMETER (CpAir = 2) PARAMETER (CpWat = 3) PARAMETER (CpVap = 4) PARAMETER (CpLiq = 5) PARAMETER (DViscAir = 6) PARAMETER (DViscLiq = 7) PARAMETER (KAir = 8) PARAMETER (KLiq = 9) PARAMETER (RhoLiq = 10) PARAMETER (Hfg = 11) PARAMETER (RAir = 12) PARAMETER (TKelMult = 13) PARAMETER (TAbsAdd = 14) PARAMETER (PaMult = 15) PARAMETER (PAbsAdd = 16) INTEGER ErrStat ErrStat = 0 C1*** Calculate the saturation pressure at a given temperature. psat = SATPRESS (Prop,TDB) C1*** Calculate the relative humidity as a function of the partial C1*** pressure of water vapor and the saturation pressure of water vapor. RH = RELHUM (Prop(Patm),psat,W) C1*** Calculate enthalpy as a function of dry bulb temperature C1*** and humidity ratio. H = ENTHALPY (Prop,TDB,W) C1*** Calculate the wet bulb temperature as a function of enthalpy or C1*** dry bulb temperature and humidity ratio. TWB = WETBULB (Prop,TDB,W) C1*** Calculate the dewpoint temperature as a function of the humidity C1*** ratio. TDP = DEWPOINT (Prop,W) C1*** Calculate dry and moist air densities. RhoD = RHODRY (Prop,TDB,W) RhoM = RHOMOIST (RhoD,W) RETURN END REAL FUNCTION DEWPOINT (Prop,W) C*********************************************************************** C* Copyright ASHRAE. Toolkit for HVAC System Energy Calculations C*********************************************************************** C* FUNCTION: DEWPOINT C* C* LANGUAGE: FORTRAN 77 C* C* PURPOSE: Calculate the dewpoint temperature given C* humidity ratio C*********************************************************************** C* INPUT VARIABLES C* W Humidity ratio (-) C* C* OUTPUT VARIABLES C* DewPoint Dew point temperature of air (C) C* C* PROPERTIES C* Patm Atmospheric pressure (Pa) C*********************************************************************** C MAJOR RESTRICTIONS: None C C DEVELOPER: Michael J. Brandemuehl, PhD, PE C University of Colorado at Boulder C C DATE: January 1, 1992 C C INCLUDE FILES: None C SUBROUTINES CALLED: None C FUNCTIONS CALLED: None C C REVISION HISTORY: None C C REFERENCE: 1989 ASHRAE Handbook - Fundamentals C*********************************************************************** C INTERNAL VARIABLES: C pw Partial water vapor pressure (Pa) C small Small number C*********************************************************************** INTEGER Patm,CpAir,CpWat,CpLiq,CpVap,DViscAir, & DViscLiq,KAir,KLiq,RhoLiq,Hfg,RAir, & TKelMult,TAbsAdd,PaMult,PAbsAdd REAL Prop(16) PARAMETER (Patm = 1) PARAMETER (CpAir = 2) PARAMETER (CpWat = 3) PARAMETER (CpVap = 4) PARAMETER (CpLiq = 5) PARAMETER (DViscAir = 6) PARAMETER (DViscLiq = 7) PARAMETER (KAir = 8) PARAMETER (KLiq = 9) PARAMETER (RhoLiq = 10) PARAMETER (Hfg = 11) PARAMETER (RAir = 12) PARAMETER (TKelMult = 13) PARAMETER (TAbsAdd = 14) PARAMETER (PaMult = 15) PARAMETER (PAbsAdd = 16) DATA small/1.E-9/ C1*** Test for "dry" air IF (W .LT. small) THEN DewPoint = -999 ELSE C1*** Calculate the partial water vapor pressure as a function of C1*** humidity ratio. pw= Prop(Patm)*W/(.62198+W) C1*** Calculate dewpoint as saturation temperature at water vapor C1*** partial pressure DewPoint = SATTEMP(Prop,pw) ENDIF 999 RETURN END REAL FUNCTION DRYBULB (Prop,H,W) C*********************************************************************** C* Copyright ASHRAE. Toolkit for HVAC System Energy Calculations C*********************************************************************** C* FUNCTION: DRYBULB C* C* LANGUAGE: FORTRAN 77 C* C* PURPOSE: Calculate the dry bulb temperature of C* moist air from enthalpy and humidity. C*********************************************************************** C* INPUT VARIABLES: C* H Enthalpy (J/kg) C* W Humidity ratio (-) C* C* OUTPUT VARIABLES: C* Drybulb Dry bulb temperature (C) C* C* PROPERTIES: C* CpAir Specific heat of air (J/kg C) C* CpVap Specific heat of water vapor (J/kg C) C* Hfg Reference heat of vaporization of water (J/kg) C*********************************************************************** C MAJOR RESTRICTIONS: Uses perfect gas relationships C Fit for enthalpy of saturated water vapor C C DEVELOPER: Shauna Gabel C Michael J. Brandemuehl, PhD, PE C University of Colorado at Boulder C C DATE: January 1, 1992 C C INCLUDE FILES: PROP.INC C SUBROUTINES CALLED: None C FUNCTIONS CALLED: None C C REVISION HISTORY: None C C REFERENCE: 1989 ASHRAE Handbook - Fundamentals C*********************************************************************** INTEGER Patm,CpAir,CpWat,CpLiq,CpVap,DViscAir, & DViscLiq,KAir,KLiq,RhoLiq,Hfg,RAir, & TKelMult,TAbsAdd,PaMult,PAbsAdd REAL Prop(16) PARAMETER (Patm = 1) PARAMETER (CpAir = 2) PARAMETER (CpWat = 3) PARAMETER (CpVap = 4) PARAMETER (CpLiq = 5) PARAMETER (DViscAir = 6) PARAMETER (DViscLiq = 7) PARAMETER (KAir = 8) PARAMETER (KLiq = 9) PARAMETER (RhoLiq = 10) PARAMETER (Hfg = 11) PARAMETER (RAir = 12) PARAMETER (TKelMult = 13) PARAMETER (TAbsAdd = 14) PARAMETER (PaMult = 15) PARAMETER (PAbsAdd = 16) C1*** Calculate the dry bulb temperature as a function of enthalpy and C1*** humidity ratio. C2*** hDryAir = Prop(CpAir)*TDB C2*** hSatVap = Prop(Hfg) + Prop(CpVap)*TDB C2*** Enthalpy = hDryAir + W*hSatVap Drybulb = (H-Prop(Hfg)*W)/(Prop(CpAir)+Prop(CpVap)*W) RETURN END REAL FUNCTION ENTHALPY (Prop,TDB,W) C*********************************************************************** C* Copyright ASHRAE. Toolkit for HVAC System Energy Calculations C*********************************************************************** C* FUNCTION: ENTHALPY C* C* LANGUAGE: FORTRAN 77 C* C* PURPOSE: Calculate the enthalpy of moist air. C*********************************************************************** C* INPUT VARIABLES: C* TDB Dry bulb temperature (C) C* W Humidity ratio (-) C* C* OUTPUT VARIABLES: C* Enthalpy Enthalpy of moist air (J/kg) C* C* PROPERTIES: C* CpAir Specific heat of air (J/kg C) C* CpVap Specific heat of water vapor (J/kg C) C* Hfg Reference heat of vaporization of water (J/kg) C*********************************************************************** C MAJOR RESTRICTIONS Uses perfect gas relationships C Fit for enthalpy of saturated water vapor C C DEVELOPER: Shauna Gabel C Michael J. Brandemuehl, PhD, PE C University of Colorado at Boulder C C DATE: January 1, 1992 C C INCLUDE FILES: PROP.INC C SUBROUTINES CALLED: None C FUNCTIONS CALLED: None C C REVISION HISTORY: None C C REFERENCE: 1989 ASHRAE Handbook - Fundamentals C*********************************************************************** INTEGER Patm,CpAir,CpWat,CpLiq,CpVap,DViscAir, & DViscLiq,KAir,KLiq,RhoLiq,Hfg,RAir, & TKelMult,TAbsAdd,PaMult,PAbsAdd REAL Prop(16) PARAMETER (Patm = 1) PARAMETER (CpAir = 2) PARAMETER (CpWat = 3) PARAMETER (CpVap = 4) PARAMETER (CpLiq = 5) PARAMETER (DViscAir = 6) PARAMETER (DViscLiq = 7) PARAMETER (KAir = 8) PARAMETER (KLiq = 9) PARAMETER (RhoLiq = 10) PARAMETER (Hfg = 11) PARAMETER (RAir = 12) PARAMETER (TKelMult = 13) PARAMETER (TAbsAdd = 14) PARAMETER (PaMult = 15) PARAMETER (PAbsAdd = 16) C1*** Calculate the enthalpy as a function of dry bulb temperature and C1*** humidity ratio. hDryAir = Prop(CpAir)*TDB hSatVap = Prop(Hfg) + Prop(CpVap)*TDB Enthalpy = hDryAir + W*hSatVap RETURN END REAL FUNCTION ENTHSAT (Prop,TDB) C*********************************************************************** C* Copyright ASHRAE. Toolkit for HVAC System Energy Calculations C*********************************************************************** C* FUNCTION: ENTHSAT C* C* LANGUAGE: FORTRAN 77 C* C* PURPOSE: Calculate the enthalpy at saturation C* for given dry bulb temperature C*********************************************************************** C* INPUT VARIABLES C* TDB Dry bulb temperature (C) C* C* OUTPUT VARIABLES C* EnthSat Enthalpy at saturation (J/kg) C* C* PROPERTIES C* Patm Atmospheric pressure (Pa) C*********************************************************************** C MAJOR RESTRICTIONS: None C C DEVELOPER: Shauna Gabel C Michael J. Brandemuehl, PhD, PE C University of Colorado at Boulder C C DATE: January 1, 1992 C C INCLUDE FILES: PROP.INC C SUBROUTINES CALLED: None C FUNCTIONS CALLED: SATPRESS C HUMRATIO C ENTHALPY C C REVISION HISTORY: None C C REFERENCE: 1989 ASHRAE Handbook - Fundamentals C*********************************************************************** C INTERNAL VARIABLES: C psat Saturated water vapor pressure (Pa) C w Humidity ratio (-) C*********************************************************************** INTEGER Patm,CpAir,CpWat,CpLiq,CpVap,DViscAir, & DViscLiq,KAir,KLiq,RhoLiq,Hfg,RAir, & TKelMult,TAbsAdd,PaMult,PAbsAdd REAL Prop(16) PARAMETER (Patm = 1) PARAMETER (CpAir = 2) PARAMETER (CpWat = 3) PARAMETER (CpVap = 4) PARAMETER (CpLiq = 5) PARAMETER (DViscAir = 6) PARAMETER (DViscLiq = 7) PARAMETER (KAir = 8) PARAMETER (KLiq = 9) PARAMETER (RhoLiq = 10) PARAMETER (Hfg = 11) PARAMETER (RAir = 12) PARAMETER (TKelMult = 13) PARAMETER (TAbsAdd = 14) PARAMETER (PaMult = 15) PARAMETER (PAbsAdd = 16) C1*** Calculate the saturation pressure at the given temperature. psat = SATPRESS (Prop,TDB) C1*** Calculate the humidity ratio from the saturation pressure w = HUMRATIO (Prop(Patm),psat) C1*** Calculate the enthalpy as a function of dry bulb temperature C1*** and humidity ratio. ENTHSAT = ENTHALPY (Prop,TDB,w) RETURN END REAL FUNCTION HUMRATIO (Patm,Pw) C*********************************************************************** C* Copyright ASHRAE. Toolkit for HVAC System Energy Calculations C*********************************************************************** C* FUNCTION: HUMRATIO C* C* LANGUAGE: FORTRAN 77 C* C* PURPOSE: Calculate the humidity ratio from water C* vapor pressure and atmospheric pressure C*********************************************************************** C* INPUT VARIABLES C* Patm Atmospheric pressure (Pa) C* Pw Partial water vapor pressure (Pa) C* C* OUTPUT VARIABLES C* HumRatio Humidity ratio (-) C*********************************************************************** C MAJOR RESRICTIONS: None C C DEVELOPER: Shauna Gabel C Michael J. Brandemuehl, PhD, PE C University of Colorado at Boulder C C DATE: January 1, 1992 C C INCLUDE FILES: None C SUBROUTINES CALLED: None C FUNCTIONS CALLED: None C C REVISION HISTORY: None C C REFERENCE: 1989 ASHRAE Handbook - Fundamentals C************************************************************************ C1*** Calculate the humidity ratio. HumRatio = 0.62198*Pw/(Patm-Pw) RETURN END REAL FUNCTION HUMTH (Prop,TDB,H) C*********************************************************************** C* Copyright ASHRAE. Toolkit for HVAC System Energy Calculations C*********************************************************************** C* FUNCTION: HUMTH C* C* LANGUAGE: FORTRAN 77 C* C* PURPOSE: Calculate the humidity ratio of moist air C* from dry bulb temperature and enthalpy. C*********************************************************************** C* INPUT VARIABLES: C* H Enthalpy (J/kg) C* TDB Dry bulb temperature (C) C* C* OUTPUT VARIABLES: C* HumTH Humidity ratio (-) C* C* PROPERTIES: C* CpAir Specific heat of air (J/kg C) C* CpVap Specific heat of water vapor (J/kg C) C* Hfg Reference heat of vaporization of water (J/kg) C*********************************************************************** C MAJOR RESTRICTIONS: Uses perfect gas relationships C Fit for enthalpy of saturated water vapor C C DEVELOPER: Shauna Gabel C Michael J. Brandemuehl, PhD, PE C University of Colorado at Boulder C C DATE: January 1, 1992 C C INCLUDE FILES: prop.inc C SUBROUTINES CALLED: None C FUNCTIONS CALLED: None C C REVISION HISTORY: None C C REFERENCE: 1989 ASHRAE Handbook - Fundamentals C*********************************************************************** INTEGER Patm,CpAir,CpWat,CpLiq,CpVap,DViscAir, & DViscLiq,KAir,KLiq,RhoLiq,Hfg,RAir, & TKelMult,TAbsAdd,PaMult,PAbsAdd REAL Prop(16) PARAMETER (Patm = 1) PARAMETER (CpAir = 2) PARAMETER (CpWat = 3) PARAMETER (CpVap = 4) PARAMETER (CpLiq = 5) PARAMETER (DViscAir = 6) PARAMETER (DViscLiq = 7) PARAMETER (KAir = 8) PARAMETER (KLiq = 9) PARAMETER (RhoLiq = 10) PARAMETER (Hfg = 11) PARAMETER (RAir = 12) PARAMETER (TKelMult = 13) PARAMETER (TAbsAdd = 14) PARAMETER (PaMult = 15) PARAMETER (PAbsAdd = 16) C1*** Calculate humidity ratio from dry bulb temperature and enthalpy C2*** hDryAir = Prop(CpAir)*TDB C2*** hSatVap = Prop(Hfg) + Prop(CpVap)*TDB C2*** Enthalpy = hDryAir + W*hSatVap HumTH = (H-Prop(CpAir)*TDB)/(Prop(Hfg)+Prop(CpVap)*TDB) RETURN END REAL FUNCTION RELHUM (Patm,Psat,HumRatio) C*********************************************************************** C* Copyright ASHRAE. Toolkit for HVAC System Energy Calculations C*********************************************************************** C* FUNCTION: RELHUM C* C* LANGUAGE: FORTRAN 77 C* C* PURPOSE: Calculate the relative humidity from C* saturation and atmospheric pressures C*********************************************************************** C* INPUT VARIABLES C* Patm Atmospheric pressure (Pa) C* Psat Saturation pressure (Pa) C* HumRatio Humidity ratio (-) C* C* OUTPUT VARIABLES C* RelHum Relative humidity (-) C****** ****************************************************************** C MAJOR RESTRICTIONS: None C C DEVELOPER: Shauna Gabel C Michael J. Brandemuehl, PhD, PE C University of Colorado at Boulder C C DATE: January 1, 1992 C C INCLUDE FILES: None C SUBROUTINES CALLED: None C FUNCTIONS CALLED: None C C REVISION HISTORY: None C C REFERENCE: 1989 ASHRAE Handbook - Fundamentals C*********************************************************************** C INTERNAL VARIABLES: C pw Partial water vapor pressure (Pa) C*********************************************************************** C1*** Calculate the partial water vapor pressure as a function of C1*** humidity ratio. pw = Patm*HumRatio/(.62198+HumRatio) C1*** Calculate the relative humidity as a function of partial water C1*** vapor pressure and water vapor pressure at saturation. RelHum = pw/Psat RETURN END REAL FUNCTION RHODRY (Prop,TDB,W) C*********************************************************************** C* Copyright ASHRAE. Toolkit for HVAC System Energy Calculations C*********************************************************************** C* FUNCTION: RHODRY C* C* LANGUAGE: FORTRAN 77 C* C* PURPOSE: Calculate dry air density. C*********************************************************************** C* INPUT VARIABLES C* TDB Dry bulb temperature (C) C* W Humidity ratio (-) C* C* OUTPUT VARIABLES C* RhoDry Density of dry air (kg/m3) C* C* PROPERTIES C* Patm Atmospheric pressure (Pa) C* RAir Gas constant for air (J/kg C) C* TAbsAdd Additive constant to convert user T to absolute T C*********************************************************************** C MAJOR RESTRICTIONS: Perfect gas relationships C C DEVELOPER: Shauna Gabel C Michael J. Brandemuehl, PhD, PE C University of Colorado at Boulder C C DATE: January 1, 1992 C C INCLUDE FILES: prop.inc C SUBROUTINES CALLED: None C FUNCTIONS CALLED: None C C REVISION HISTORY: None C C REFERENCE: 1989 ASHRAE Handbook - Fundamentals C*********************************************************************** C INTERNAL VARIABLES: C pAir Partial pressure of dry air (Pa) C*********************************************************************** INTEGER Patm,CpAir,CpWat,CpLiq,CpVap,DViscAir, & DViscLiq,KAir,KLiq,RhoLiq,Hfg,RAir, & TKelMult,TAbsAdd,PaMult,PAbsAdd REAL Prop(16) PARAMETER (Patm = 1) PARAMETER (CpAir = 2) PARAMETER (CpWat = 3) PARAMETER (CpVap = 4) PARAMETER (CpLiq = 5) PARAMETER (DViscAir = 6) PARAMETER (DViscLiq = 7) PARAMETER (KAir = 8) PARAMETER (KLiq = 9) PARAMETER (RhoLiq = 10) PARAMETER (Hfg = 11) PARAMETER (RAir = 12) PARAMETER (TKelMult = 13) PARAMETER (TAbsAdd = 14) PARAMETER (PaMult = 15) PARAMETER (PAbsAdd = 16) C1*** Calculate the dry air density from perfect gas laws. pAir = 0.62198*Prop(Patm)/(0.62198+W) RhoDry = pAir/Prop(RAir)/(TDB+Prop(TAbsAdd)) RETURN END REAL FUNCTION RHOMOIST (RhoDry,W) C*********************************************************************** C* Copyright ASHRAE. Toolkit for HVAC System Energy Calculations C*********************************************************************** C* FUNCTION: RHOMOIST C* C* LANGUAGE: FORTRAN 77 C* C* PURPOSE: Calculate moist air density from dry air C* density and humidity ratio C*********************************************************************** C* INPUT VARIABLES: C* RhoDry Dry air density (kg/m3) C* W Humidity ratio (-) C* C* OUTPUT VARIABLES: C* RhoMoist Density of dry air (kg/m3) C*********************************************************************** C MAJOR RESTRICTIONS: None C C DEVELOPER: Shauna Gabel C Michael J. Brandemuehl, PhD, PE C University of Colorado at Boulder C C DATE: January 1, 1992 C C INCLUDE FILES: None C SUBROUTINES CALLED: None C FUNCTIONS CALLED: None C C REVISION HISTORY: None C C REFERENCE: 1989 ASHRAE Handbook - Fundamentals C************************************************************************ C1*** Calculate the moist air density RhoMoist = RhoDry*(1.+W) RETURN END REAL FUNCTION SATPRESS (Prop,T) C*********************************************************************** C* Copyright ASHRAE. Toolkit for HVAC System Energy Calculations C*********************************************************************** C* SUBROUTINE: SATPRESS C* C* LANGUAGE: FORTRAN 77 C* C* PURPOSE: Calculate saturation pressure of water C* vapor as a function of temperature C*********************************************************************** C* INPUT VARIABLES C* T Temperature (C) C* C* OUTPUT VARIABLES C* SatPress Saturation pressure (Pa) C* C* PROPERTIES C* TKelMult Multiplying factor to convert user T to Kelvin C* TAbsAdd Additive factor to convert user T to absolute T C* tKel = Prop(TKelMult) * (T + Prop(TAbsAdd)) C* PaMult Multiplying factor to convert user P to Pascals C* PAbsAdd Additive factor to convert user P to absolute P C* Pa = Prop(PaMult) * (P + Prop(PAbsAdd)) C*********************************************************************** C MAJOR RESTRICTIONS: 173.16 K <= Temp <= 473.15 K C C DEVELOPER: Shauna Gabel C Michael J. Brandemuehl, PhD, PE C University of Colorado at Boulder C C DATE: January 1, 1992 C C INCLUDE FILES: prop.inc C SUBROUTINES CALLED: None C FUNCTIONS CALLED: None C C REVISION HISTORY: None C C REFERENCE: 1989 ASHRAE Handbook - Fundamentals C C Hyland, R.W., and A. Wexler. 1983. C Formulations for the thermodynamic C properties of the saturated phases of H2O C from 173.15 K to 473.15 K. ASHRAE C Transactions, Vol. 89, No. 2A, pp. 500-519 C*********************************************************************** C INTERNAL VARIABLES: C tKel Temperature in Kelvin (K) C pascals Saturation pressure (Pa) C*********************************************************************** INTEGER Patm,CpAir,CpWat,CpLiq,CpVap,DViscAir, & DViscLiq,KAir,KLiq,RhoLiq,Hfg,RAir, & TKelMult,TAbsAdd,PaMult,PAbsAdd REAL Prop(16) PARAMETER (Patm = 1) PARAMETER (CpAir = 2) PARAMETER (CpWat = 3) PARAMETER (CpVap = 4) PARAMETER (CpLiq = 5) PARAMETER (DViscAir = 6) PARAMETER (DViscLiq = 7) PARAMETER (KAir = 8) PARAMETER (KLiq = 9) PARAMETER (RhoLiq = 10) PARAMETER (Hfg = 11) PARAMETER (RAir = 12) PARAMETER (TKelMult = 13) PARAMETER (TAbsAdd = 14) PARAMETER (PaMult = 15) PARAMETER (PAbsAdd = 16) DATA C1/-5674.5359/,C2/6.3925247/,C3/-0.9677843E-2/ DATA C4/0.62215701E-6/,C5/0.20747825E-8/,C6/-0.9484024E-12/ DATA C7/4.1635019/,C8/-5800.2206/,C9/1.3914993/,C10/-0.048640239/ DATA C11/0.41764768E-4/,C12/-0.14452093E-7/,C13/6.5459673/ C1*** Convert temperature from user units to Kelvin. tKel = Prop(TKelMult)*(T+Prop(TAbsAdd)) C1*** If below freezing, calculate saturation pressure over ice. IF (tKel .LT. 273.15) THEN pascals = EXP(C1/tKel+C2+C3*tKel+C4*tKel**2+C5*tKel**3+C6* & tKel**4+C7*ALOG(tKel)) C1*** If above freezing, calculate saturation pressure over liquid water. ELSE IF (tKel .GE. 273.15) THEN pascals = EXP(C8/tKel+C9+C10*tKel+C11*tKel**2+C12*tKel**3+C13 & *ALOG(tKel)) ENDIF C1*** Convert pressure from Pascals to user units SatPress = pascals/Prop(PaMult) - Prop(PAbsAdd) RETURN END REAL FUNCTION SATTEMP (Prop,P) C*********************************************************************** C* Copyright ASHRAE. Toolkit for HVAC System Energy Calculations C*********************************************************************** C* FUNCTION: SATTEMP C* C* LANGUAGE: FORTRAN 77 C* C* PURPOSE: Calculate the saturation (boiling) C* temperature of water given pressure C*********************************************************************** C* INPUT VARIABLES C* P Pressure (Pa) C* C* OUTPUT VARIABLES C* SatTemp Saturation temperature of water vapor (C) C*********************************************************************** C MAJOR RESTRICTIONS: None C C DEVELOPER: Shauna Gabel C Michael J. Brandemuehl, PhD, PE C University of Colorado at Boulder C C DATE: January 1, 1992 C C INCLUDE FILES: prop.inc C SUBROUTINES CALLED: None C FUNCTIONS CALLED: SATPRESS C XITERATE C C REVISION HISTORY: None C C REFERENCE: 1989 ASHRAE Handbook - Fundamentals C*********************************************************************** C INTERNAL VARIABLES: C tSat Water temperature guess (C) C pSat Pressure corresponding to temp. guess (Pa) C error Deviation of dependent variable in iteration C iter Iteration counter C icvg Iteration convergence flag C F1,F2 Previous values of dependent variable in XITERATE C X1,X2 Previous values of independent variable in XITERATE C*********************************************************************** DATA itmax/50/ C1*** Use an iterative process to determine the saturation temperature C1*** at a given pressure using a correlation of saturated water vapor C1*** pressure as a function of temperature C1*** Initial guess of boiling temperature tSat = 100. C1*** Iterate to find the saturation temperature C1*** of water given the total pressure C2*** Set iteration loop parameters DO 100 iter = 1,itmax C1*** Calculate saturation pressure for estimated boiling temperature pSat = SATPRESS(Prop,tSat) C1*** Compare with specified pressure and update estimate of temperature error = P - pSat tSat = XITERATE (tSat,error,X1,F1,X2,F2,iter,icvg) C2*** If converged leave loop iteration IF (icvg .EQ. 1) GO TO 110 C2*** Water temperature not converged, repeat calculations with new C2*** estimate of water temperature 100 CONTINUE C1*** Saturation temperature has not converged after maximum specified C1*** iterations. Print error message, set return error flag, and RETURN WRITE(LUW,1001) itmax 1001 FORMAT(/1X,'*** ERROR IN FUNCTION SatTemp ***'/ & 1X,' Saturation temperature has not ' & 'converged after ',I2,' iterations'/) 110 SatTemp = tSat RETURN END REAL FUNCTION TAIRSAT (Prop,HSat) C*********************************************************************** C* Copyright ASHRAE. Toolkit for HVAC System Energy Calculations C*********************************************************************** C* FUNCTION: ENTHSAT C* C* LANGUAGE: FORTRAN 77 C* C* PURPOSE: Calculate the dry bulb temperature given C* enthalpy at saturation. C*********************************************************************** C* INPUT VARIABLES: C* HSat Enthalpy at saturation (J/kg) C* C* OUTPUT VARIABLES: C* TAirSat Dry bulb temperature (C) C*********************************************************************** C MAJOR RESTRICTIONS: None C C DEVELOPER: Michael J. Brandemuehl, PhD, PE C University of Colorado at Boulder C C DATE: January 1, 1992 C C INCLUDE FILES: prop.inc C SUBROUTINES CALLED: None C FUNCTIONS CALLED: ENTHSAT C C REVISION HISTORY: None C C REFERENCE: 1989 ASHRAE Handbook - Fundamentals C*********************************************************************** C INTERNAL VARIABLES: C error Deviation of dependent variable in iteration C iter Iteration counter C icvg Iteration convergence flag C F1,F2 Previous values of dependent variable in XITERATE C X1,X2 Previous values of independent variable in XITERATE C*********************************************************************** DATA itmax/20/,tSat/50./ C1*** Estimate saturation temperature if reasonable value not available IF(tSat .LT. -200. .OR. tSat .GT. 1000.) tSat = 50. C1*** Calculate saturation temperature by iteration using function to C1*** calculate saturation enthalpy from temperature DO 100 iter=1,itmax error = HSat - ENTHSAT(Prop,tSat) tSat = XITERATE(tSat,error,X1,F1,X2,F2,iter,icvg) C1*** If converged, leave iteration loop. IF (icvg .EQ. 1) GO TO 110 C1*** Temperature not converged, repeat calculation with new C1*** estimate of temperature. 100 CONTINUE C1*** Temperature has not converged after maximum specified C1*** iterations. Print error message and RETURN WRITE(LUW,1001) itmax 1001 FORMAT(/1X,'*** ERROR IN FUNCTION TAIRSAT ***'/ & 1X,' Temperature has not ' & 'converged after ',I2,' iterations'/) 110 CONTINUE TAirSat = tSat RETURN END REAL FUNCTION WETBULB (Prop,TDB,W) C*********************************************************************** C* Copyright ASHRAE. Toolkit for HVAC System Energy Calculations C*********************************************************************** C* FUNCTION: WETBULB C* C* LANGUAGE: FORTRAN 77 C* C* PURPOSE: Calculate wet bulb temperature from dry C* bulb temperature and humidity ratio C*********************************************************************** C* INPUT VARIABLES C* TDB Dry bulb temperature (C) C* W Humidity ratio of air (-) C* C* OUTPUT VARIABLES C* WetBulb Wet bulb temperature (C) C* C* PROPERTIES: C* Patm Atmospheric pressure (Pa) C* Hfg Latent heat of vaporization of water (J/kg) C* CpAir Specific heat of air (J/kg C) C* CpVap Specific heat of water vapor (J/kg C) C* CpWat Specific heat of water (J/kg C) C*********************************************************************** C MAJOR RESTRICTIONS: None C C DEVELOPER: Shauna Gabel C Michael J. Brandemuehl, PhD, PE C University of Colorado at Boulder C C DATE: January 1, 1992 C C INCLUDE FILES: prop.inc C SUBROUTINES CALLED: None C FUNCTIONS CALLED: SATPRESS C HUMRATIO C SATTEMP C XITERATE C C REVISION HISTORY: None C C REFERENCE: 1989 ASHRAE Handbook - Fundamentals C*********************************************************************** C INTERNAL VARIABLES: C tBoil Boiling temperature of water at given pressure (C) C psatStar Saturation pressure at wet bulb temperature (C) C wStar Humidity ratio as a function of PsatStar (-) C newW Humidity ratio calculated with wet bulb guess (-) C error Deviation of dependent variable in iteration C iter Iteration counter C icvg Iteration convergence flag C F1,F2 Previous values of dependent variable in XITERATE C X1,X2 Previous values of independent variable in XITERATE C*********************************************************************** INTEGER Patm,CpAir,CpWat,CpLiq,CpVap,DViscAir, & DViscLiq,KAir,KLiq,RhoLiq,Hfg,RAir, & TKelMult,TAbsAdd,PaMult,PAbsAdd REAL Prop(16) PARAMETER (Patm = 1) PARAMETER (CpAir = 2) PARAMETER (CpWat = 3) PARAMETER (CpVap = 4) PARAMETER (CpLiq = 5) PARAMETER (DViscAir = 6) PARAMETER (DViscLiq = 7) PARAMETER (KAir = 8) PARAMETER (KLiq = 9) PARAMETER (RhoLiq = 10) PARAMETER (Hfg = 11) PARAMETER (RAir = 12) PARAMETER (TKelMult = 13) PARAMETER (TAbsAdd = 14) PARAMETER (PaMult = 15) PARAMETER (PAbsAdd = 16) REAL newW DATA itmax/20/ C1*** Initial temperature guess tBoil = SATTEMP (Prop,Prop(Patm)) WetBulb = MAX( MIN(WetBulb,TDB,(tBoil-0.1)), 0.) C1*** Begin iteration loop DO 100 iter = 1,itmax IF (WetBulb .GE. (tBoil-0.09) ) WETBULB = tBoil-0.1 C1*** Determine the saturation pressure for wet bulb temperature psatStar = SATPRESS (Prop,WetBulb) C1*** Determine humidity ratio for given saturation pressure wStar = HUMRATIO (Prop(Patm),psatStar) C1*** Calculate new humidity ratio and determine difference from known C1*** humidity ratio newW = ((Prop(Hfg)-(Prop(CpWat)-Prop(CpVap))*WetBulb)*wStar- & Prop(CpAir)*(TDB-WetBulb))/(Prop(Hfg)+Prop(CpVap)*TDB & -Prop(CpWat)*WetBulb) C1*** Check error, if not satisfied, calculate new guess and iterate error = W-newW WetBulb = XITERATE(WetBulb,error,X1,F1,X2,F2,iter,icvg) C1*** If converged, leave iteration loop. IF (icvg .EQ. 1) GO TO 900 C1*** Wet bulb temperature not converged, repeat calculation with new C1*** estimate of wet bulb temperature. 100 CONTINUE C1*** Wet bulb temperature has not converged after maximum specified C1*** iterations. Print error message, set return error flag, and RETURN WRITE(LUW,1009) itmax 1009 FORMAT(/1X,'*** ERROR IN FUNCTION WetBulb ***'/ & 1X,' Wet bulb temperature has not ' & 'converged after ',I2,' iterations'/) 900 IF (WetBulb .GT. TDB) WetBulb = TDB 999 RETURN END REAL FUNCTION XITERATE (X0,F0,X1,F1,X2,F2,ICount,ICvg) C*********************************************************************** C* Copyright ASHRAE. Toolkit for HVAC System Energy Calculations C*********************************************************************** C* C* SUBROUTINE: XITERATE C* C* LANGUAGE: FORTRAN 77 C* C* PURPOSE: Iterately solves for the value of X which C* satisfies F(X)=0. Given Xi,F(Xi) pairs, C* the subroutine tests for convergence and C* provides a new guess for the value of the C* independent variable X. C*********************************************************************** C* INPUT VARIABLES C* F0 Current value of the function F(X) C* X0 Current value of X C* F1,F2 Two previous values of F(Xi) C* X1,X2 Two previous values of X C* C* NOTE: F1,X1,F2,X2 MUST BE STORED AND SAVED IN CALLING C* ROUTINE. THEY NEED NO INITIALIZATION C* C* ICount Number of iterations C* C* OUTPUT VARIABLES C* XIterate New estimate of X for F(X)=0 C* ICvg Convergence flag ICvg = 0: Not converged C* ICvg = 1: Converged C*********************************************************************** C DEVELOPER: Michael J. Brandemuehl, PhD, PE C University of Colorado at Boulder C C DATE: January 1, 1992 C C INCLUDE FILES: None C SUBROUTINES CALLED: None C FUNCTIONS CALLED: None C C REVISION HISTORY: None C C REFERENCE: None C*********************************************************************** C INTERNAL VARIABLES C small Small number used in place of zero C mode Number of points used in fit C mode = 1: Use XPerburb to get new X C mode = 2: Linear equation to get new X C mode > 2: Quadratic equation to get new X C coef(i) Coefficients for quadratic fit C F(X) = coef(1) + coef(2)*X + coef(3)*X*X C check Term under radical in quadratic solution C FiQ,XiQ Double precision values of Fi,Xi C slope Slope for linear fit C tolRel Relative error tolerance C xPerturb Perturbation applied to X to initialize iteration C*********************************************************************** DOUBLE PRECISION coef(3),check,F0Q,F1Q,F2Q,X0Q,X1Q,X2Q DATA tolRel/1.E-5/,xPerturb/0.1/,small/1.E-9/ C1*** Check for convergence by comparing change in X IF ((ABS(X0-X1) .LT. tolRel*MAX(ABS(X0),small) .AND. & ICount .NE. 1) .OR. F0 .EQ. 0.) THEN XIterate = X0 ICvg=1 RETURN ENDIF C1*** Not converged. C2*** If after the second iteration there are enough previous points to C2 fit a quadratic for the new X. If the quadratic fit is not C2 applicable, mode will be set to 1 or 2 and a new X will be C2 determined by incrementing X from xPerturb or from a linear fit. ICvg=0 mode=ICount 10 IF (mode .EQ. 1) THEN C1*** New guess is specified by xPerturb IF (ABS(X0) .GT. small) THEN XIterate = X0*(1.+xPerturb) ELSE XIterate = xPerturb ENDIF ELSEIF (mode .EQ. 2) THEN C1*** New guess calculated from LINEAR FIT of most recent two points SLOPE=(F1-F0)/(X1-X0) IF(slope.EQ.0) THEN mode=1 GO TO 10 ENDIF XIterate=X0-F0/SLOPE ELSE C1*** New guess calculated from QUADRATIC FIT C1*** If two Xi are equal, set mode for linear fit and return to top IF (X0 .EQ. X1) THEN X1=X2 F1=F2 mode=2 GO TO 10 ELSEIF (X0 .EQ. X2) THEN mode=2 GO TO 10 ENDIF C1*** Determine quadratic coefficients from the three data points C1*** using double precision. F2Q=F2 F1Q=F1 F0Q=F0 X2Q=X2 X1Q=X1 X0Q=X0 coef(3)=((F2Q-F0Q)/(X2Q-X0Q)-(F1Q-F0Q)/(X1Q-X0Q))/(X2Q-X1Q) coef(2)=(F1Q-F0Q)/(X1Q-X0Q)-(X1Q+X0Q)*coef(3) coef(1)=F0-(coef(2)+coef(3)*X0Q)*X0Q C1*** If points are colinear, set mode for linear fit and return to top IF (ABS(coef(3)) .LT. 1.D-10) THEN mode=2 GO TO 10 ENDIF C1*** Check for precision. If the coefficients do not accurately C1*** predict the given data points due to round-off errors, set C1*** mode for a linear fit and return to top. IF (ABS((coef(1)+(coef(2)+coef(3)*X1Q)*X1Q-F1Q)/F1Q) .GT. & 1.D-4) THEN mode=2 GO TO 10 ENDIF C1*** Check for imaginary roots. If no real roots, set mode to C1*** estimate new X by simply incrementing by xPerturb check=coef(2)**2-4*coef(1)*coef(3) IF (check .LT. 0) THEN C1*** Imaginary roots -- go back to linear fit mode=2 GO TO 10 ELSEIF (check .GT. 0) THEN C1*** Real unequal roots -- determine root nearest to most recent guess XIterate=(-coef(2)+SQRT(check))/coef(3)/2 xOther=-XIterate-coef(2)/coef(3) IF (ABS(XIterate-X0) .GT. ABS(xOther-X0)) XIterate=xOther ELSE C1*** Real Equal Roots -- one solution XIterate=-coef(2)/coef(3)/2 ENDIF ENDIF C1*** Set previous variable values for the next iteration IF (mode .LT. 3) THEN C1*** No valid previous points to eliminate. X2=X1 F2=F1 X1=X0 F1=F0 ELSE C1*** Eliminate one previous point based on sign and magnitude of F(X) C2*** Keep the current point and eliminate one of the previous ones. IF (F1*F0 .GT. 0 .AND. F2*F0 .GT. 0) THEN C2*** All previous points of same sign. Eliminate one with biggest F(X) IF (ABS(F2) .GT. ABS(F1)) THEN X2=X1 F2=F1 ENDIF ELSE C1*** Points of different sign. C1*** Eliminate the previous one with the same sign as current F(X). IF (F2*F0 .GT. 0) THEN X2=X1 F2=F1 ENDIF ENDIF X1=X0 F1=F0 ENDIF RETURN END C*********************************************************************** C* FILE: PROP.INC C* C* This file assigns a numbers to air and water property names to be C* used in the "Prop" array. C*********************************************************************** C DEVELOPER: Inger Andresen C Michael J. Brandemuehl, PhD, PE C C DATE: July 1, 1991 C C FILES REQUIRED: None C*********************************************************************** C INTERNAL VARIABLES: C Patm Atmospheric pressure (Pa) C CpAir Specific heat of dry air (J/kg C) C CpLiq Specific heat of liquid water (J/kg C) C CpVap Specific heat of saturated water vapor (J/kg C) C DViscAir Air dynamic viscosity (kg/m s) C DViscLiq Liquid dynamic viscosity (kg/m s) C KAir Air thermal conductivity (W/m C) C KLiq Liquid thermal conductivity (W/m C) C RhoLiq Liquid density (kg/m3) C Hfg Latent heat of vaporization of water (J/kg) C RAir Gas constant for air (J/kg C) C TKelMult Multiplying factor to convert user T to Kelvin C TAbsAdd Additive factor to convert user P to Kelvin C tKel = Prop(TKelMult)*T + Prop(TKelAdd) C PaMult Multiplying factor to convert user P to Pascals C PAbsAdd Additive factor to convert user P to Pascals C Pa = Prop(PaMult)*P + Prop(PaAdd) C*********************************************************************** C C INTEGER Patm,CpAir,CpWat,CpLiq,CpVap,DViscAir, C & DViscLiq,KAir,KLiq,RhoLiq,Hfg,RAir, C & TKelMult,TAbsAdd,PaMult,PAbsAdd C REAL Prop(16) C C PARAMETER (Patm = 1) C PARAMETER (CpAir = 2) C PARAMETER (CpWat = 3) C PARAMETER (CpVap = 4) C PARAMETER (CpLiq = 5) C PARAMETER (DViscAir = 6) C PARAMETER (DViscLiq = 7) C PARAMETER (KAir = 8) C PARAMETER (KLiq = 9) C PARAMETER (RhoLiq = 10) C PARAMETER (Hfg = 11) C PARAMETER (RAir = 12) C PARAMETER (TKelMult = 13) C PARAMETER (TAbsAdd = 14) C PARAMETER (PaMult = 15) C PARAMETER (PAbsAdd = 16)