SUBROUTINE TYPE69 (TIME,XIN,OUT,T,DTDT,PAR,INFO,ICNTRL,*) C*********************************************************************** C* Copyright ASHRAE. Toolkit for HVAC System Energy Calculations C* Revised for TRNSYS by N.BLAIR and R.Schwarz C* C*********************************************************************** C* C* WARNING!!! THE ASHRAE TOOLKIT DIDN'T CONTAIN SAMPLE INPUT AND C* OUTPUT! THE SAMPLE VALUES WERE GENERATED, USING THE C* TRNSYS PROGRAM DEBUG!!!!! C* C*********************************************************************** C* SUBROUTINE: WETCOIL C* C* LANGUAGE: FORTRAN 77 C* C* PURPOSE: Calculate the performance of a cooling C* coil when the external fin surface is C* complete wet. Results include C* outlet air temperature and humidity, C* outlet water temperature, sensible and C* total cooling capacities, and the wet C* fraction of the air-side surface area. C*********************************************************************** C* INPUT VARIABLES DESCRIPTION(UNITS) SAMPLE VALUES C* XIN(1) MLiq Liquid mass flow rate(kg/s) 3.600 C* XIN(2) TLiqEnt Entering water temperature(C) 12.000 C* XIN(3) MAir Dry air mass flow rate(kg/s) 9.500 C* XIN(4) TAirEnt Entering air dry bulb temperature(C) 25.000 C* XIN(5) WAirEnt Entering air humidity ratio(-) 0.008 C* C* XIN(6) UAIntTot Internal overall heat transfer coefficient(W/m2 C) 1000 C* XIN(7) UAExtTot External overall heat transfer coefficient(W/m2 C) 700 C* XIN(8) ConfigHX Heat exchanger configuration(-) 1.0 C* 1 - Counterflow C* 2 - Parallel flow C* 3 - Cross flow, both streams unmixed C* 4 - Cross flow, both streams mixed C* 5 - Cross flow, stream 1 unmixed C* 6 - Cross flow, stream 2 unmixed C* C* OUTPUT VARIABLES C* OUT(1) TLiqLvg Leaving water temperature(C) 12.190 C* OUT(2) TAirLvg Leaving air dry bulb temperature(C) 24.700 C* OUT(3) WAirLvg Leaving air humidity ratio(-) 0.008 C* OUT(4) QTot Total heat transfer rate(W) 2827.000 C* OUT(5) QSen Sensible heat transfer rate(W) 2827.000 C* OUT(6) FWet Fraction of surface area wet(-) 1.0 C* OUT(7) TSurfEnt Surface temperature at air entrance(C) 14.780 C* OUT(8) ErrStat Error status indicator,0=ok,1=error(-) 0.0 C* C* PROPERTIES C* CpLiq Specific heat of liquid (J/kg C) C* CpAir Specific heat of dry air (J/kg C) C*********************************************************************** C MAJOR RESTRICTIONS: Models coil as counterflow heat exchanger C Approximates saturated air enthalpy as C a linear function of temperature 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: HEATEX C WCOILOUT C FUNCTIONS CALLED: ENTHALPY C ENTHSAT C TAIRSAT C C REVISION HISTORY: None C C REFERENCE: Elmahdy, A.H. and Mitalas, G.P. 1977. C "A Simple Model for Cooling and C Dehumidifying Coils for Use In Calculating C Energy Requirements for Buildings," C ASHRAE Transactions,Vol.83 Part 2, C pp. 103-117. C C TRNSYS. 1990. A Transient System C Simulation Program: Reference Manual. C Solar Energy Laboratory, Univ. Wisconsin- C Madison, pp. 4.6.8-1 - 4.6.8-12. C C Threlkeld, J.L. 1970. Thermal C Environmental Engineering, 2nd Edition, C Englewood Cliffs: Prentice-Hall,Inc. C pp. 254-270. C*********************************************************************** C INTERNAL VARIABLES: C extResist Air-side resistance to heat transfer (m2 C/W) C intResist Liquid-side resistance to heat transfer (m2 C/W) C tDewEnt Entering air dew point (C) C uaH Overall enthalpy heat transfer coefficient (kg/s) C capAirWet Air-side capacity rate (kg/s) C capLiqWet Liquid-side capacity rate (kg/s) C resistRatio Ratio of resistances (-) C hAirLvg Outlet air enthalpy C hLiqEntSat Saturated enthalpy of air at (J/kg) C entering water temperature C hLiqLvgSat Saturated enthalpy of air at exit (J/kg) C water temperature C hSurfEntSat Saturated enthalpy of air at (J/kg) C entering surface temperature C hSurfLvgSat Saturated enthalpy of air at exit (J/kg) C surface temperature C cpSat Coefficient for equation below (J/kg C) C EnthSat1-EnthSat2 = cpSat*(TSat1-TSat2) C (all water and surface temperatures are C related to saturated air enthalpies for C wet surface heat transfer calculations) C************************************************************************ DOUBLE PRECISION XIN, OUT REAL intResist,MLIQ,SMALL INTEGER ErrStat, INFO, IOPT, NI, NP, ND DIMENSION XIN(8),OUT(8),INFO(15) CHARACTER*3 YCHECK(8), OCHECK(8) COMMON /LUNITS/LUR,LUW,IFORM,LUK DATA small/1.E-9/ DATA PATM/101325.0/,CPAIR/1006.0/,CPWAT/4186.0/,HFG/2501000.0/ &,RAIR/287.055/,TKELMULT/1.0/,TABSADD/273.15/,PAMULT/1.0/, &PABSADD/0.0/,CPLIQ/4186.0/ DATA YCHECK/'MF2','TE1','MF2','TE1','DM1','NAV','NAV','DM1'/ DATA OCHECK/'TE1','TE1','DM1','PW2','PW2','DM1','TE1','DM1'/ MLIQ = XIN(1) TLIQENT = XIN(2) MAIR = XIN(3) TAIRENT = XIN(4) WAIRENT = XIN(5) UAINTTOT = XIN(6) UAEXTTOT = XIN(7) CONFIGHX = XIN(8) IOPT = -1 NI = 8 !CORRECT NUMBER OF INPUTS NP = 0 !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 THE USER'S INPUT MATCHES THE CORRECT NUMBERS CALL RCHECK(INFO,YCHECK,OCHECK) C CHECKS TO SEE IN INPUT AND OUTPUT UNITS MATCH ENDIF FWet = 1. extResist = 1./UAExtTot intResist = 1./UAIntTot C1*** Calculate enthalpies of entering air and water hAirEnt = ENTHALPY(CPAIR,CPVAP,HFG,TAirEnt,WAirEnt) hLiqEntSat = ENTHSAT(PATM,CPAIR,CPVAP,HFG,TKELMULT, & TABSADD,PAMULT,PABSADD,TLiqEnt) C1*** Estimate cpSat using entering air dewpoint and water temperature tDewEnt = DEWPOINT(PATM,TKELMULT,TABSADD,PAMULT, & PABSADD,WAirEnt) cpSat = (ENTHSAT(PATM,CPAIR,CPVAP,HFG,TKELMULT, & TABSADD,PAMULT,PABSADD,tDewEnt)-hLiqEntSat) & /(tDewEnt-TLiqEnt) C1*** Enthalpy-based heat transfer calculations C2*** Heat transfer in a wet coil is calculated based on enthalpy C2*** rather than temperature to include latent effects. Air enthalpies C2*** are evaluated using conventional psychrometric equations. The C2*** corresponding enthalpies of the coil and water are related to C2*** that of the air through "fictitious enthalpies," defined as the C2*** enthalpy of saturated air at the temperature of the coil or water. C2 C2*** While heat transfer rates are commonly expressed as the product C2*** of an overall heat transfer coefficient, UA, and a temperature C2*** difference, the use of enthalpy-based heat transfer calculations C2*** requires an enthalpy-based heat transfer coefficient, UAH. C2 C2*** q = UAH * (H1-H2) C2 C2*** where UAH = UA / cp C2*** UA = conventional heat transfer coefficient C2*** cp = specific heat across enthalpy difference C2 C2*** When using fictitious enthalpies, a corresponding fictitious C2*** specific heat must be defined. C2 C2*** EnthSat1-EnthSat2 = cpSat * (Temp1-Temp2) C2 C2*** UAH can be calculated from a combination of series or parallel C2*** enthalpy resistances, similar to thermal resistances modified for C2*** enthalpy as above. Enthalpy capacity rates relate heat transfer C2*** to the enthalpy change of a fluid between inlet and outlet. C2 C2*** q = CapH * (HAirLvg - HAirEnt) C2 C2*** On the air side, enthalpy capacity rate is the air mass flow rate. C2*** On the water side, the enthalpy capacity rate is based on the C2*** enthalpy of saturated air at the water temperature. C1*** Determine air and water enthalpy outlet conditions by modeling C1*** coil as counterflow enthalpy heat exchanger uaH = 1./(cpSat*intResist+CpAir*extResist) capAirWet = MAir capLiqWet = MLiq * (CpLiq/cpSat) CALL HEATEX(capAirWet,hAirEnt,capLiqWet,hLiqEntSat,uaH, & ConfigHX,hAirLvg,hLiqLvgSat) C1*** Calculate entering and leaving external surface conditions from C1*** air and water conditions and the ratio of resistances resistRatio = (intResist)/(intResist + & CpAir/cpSat*extResist) hSurfEntSat = hLiqLvgSat + resistRatio*(hAirEnt-hLiqLvgSat) hSurfLvgSat = hLiqEntSat + resistRatio*(hAirLvg-hLiqEntSat) TSurfEnt = TAIRSAT(PATM,CPAIR,CPVAP,HFG,TKELMULT,TABSADD, & PAMULT,PABSADD,hSurfEntSat) C1*** Calculate outlet air temperature and humidity from enthalpies and C1*** surface conditions. QTot = MAir*(hAirEnt-hAirLvg) TLiqLvg = TLiqEnt+QTot/MAX(MLiq,small)/CpLiq CALL WCOILOUT (PATM,CPAIR,CPVAP,HFG,TKELMULT,TABSADD, & PAMULT,PABSADD,MAir,TAirEnt,WAirEnt, & hAirEnt,hAirLvg,UAExtTot,TAirLvg, & WAirLvg,QSen,ErrStat) 999 CONTINUE OUT(1) = TLIQLVG OUT(2) = TAIRLVG OUT(3) = WAIRLVG OUT(4) = QTOT OUT(5) = QSEN OUT(6) = FWET OUT(7) = TSURFENT OUT(8) = ERRSTAT RETURN 1 END SUBROUTINE HEATEX (Cap1,In1,Cap2,In2,UA,ConfigHX,Out1,Out2) C*********************************************************************** C* Copyright ASHRAE. Toolkit for HVAC System Energy Calculations C*********************************************************************** C* SUBROUTINE: HEATEX C* C* LANGUAGE: FORTRAN 77 C* C* PURPOSE: Calculate the outlet states of a simple C* heat exchanger using the effectiveness-Ntu C* method of analysis. C*********************************************************************** C* INPUT VARIABLES C* Cap1 Capacity rate of stream 1 (W/C) C* In1 Inlet state of stream 1 (C) C* Cap2 Capacity rate of stream 2 (W/C) C* In2 Inlet state of stream 2 (C) C* UA Overall heat transfer coefficient (W/C) C* ConfigHX Heat exchanger configuration (-) C* 1 - Counterflow C* 2 - Parallel flow C* 3 - Cross flow, both streams unmixed C* 4 - Cross flow, both streams mixed C* 5 - Cross flow, stream 1 unmixed C* 6 - Cross flow, stream 2 unmixed C* C* OUTPUT VARIABLES C* Out1 Outlet state of stream 1 (C) C* Out2 Outlet state of stream 2 (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: None C SUBROUTINES CALLED: None C FUNCTIONS CALLED: None C C REVISION HISTORY: None C C REFERENCE: Kays, W.M. and A.L. London. 1964. C Compact Heat Exchangers, 2nd Ed., McGraw- C Hill: New York. C*********************************************************************** C* INTERNAL VARIABLES: C* cMin Minimum capacity rate of the streams (W/C) C* cMax Maximum capacity rate of the streams (W/C) C* cRatio Ratio of minimum to maximum capacity rate C* ntu Number of transfer units (-) C* effectiveness Heat exchanger effectiveness (-) C* qMax Maximum heat transfer possible (W) C*********************************************************************** REAL ntu,qMax,In1,In2,large DATA small/1.E-15/, large/1.E15/ C1*** Ntu and Cmin/Cmax (cRatio) calculations cMin = MIN(Cap1,Cap2) cMax = MAX(Cap1,Cap2) IF( cMax .EQ. 0.) THEN cRatio = 1. ELSE cRatio = cMin/cMax ENDIF IF( cMin .EQ. 0.) THEN ntu = large ELSE ntu = ua/cMin ENDIF C1*** Calculate effectiveness for special limiting cases mode = NINT(ConfigHX) IF(ntu .LE. 0) THEN effectiveness = 0. ELSE IF(cRatio .LT. small) THEN C2*** Cmin/Cmax = 0 and effectiveness is independent of configuration effectiveness = 1 - EXP(-ntu) C1*** Calculate effectiveness depending on heat exchanger configuration ELSE IF (mode .EQ. 1) THEN C2*** Counterflow IF (ABS(cRatio-1.) .LT. small) THEN effectiveness = ntu/(ntu+1.) ELSE e=EXP(-ntu*(1-cRatio)) effectiveness = (1-e)/(1-cRatio*e) ENDIF ELSE IF (mode .EQ. 2) THEN C2*** Parallel flow effectiveness = (1-EXP(-ntu*(1+cRatio)))/(1+cRatio) ELSE IF (mode .EQ. 3) THEN C2*** Cross flow, both streams unmixed eta = ntu**(-0.22) effectiveness = 1 - EXP((EXP(-ntu*cRatio*eta)-1)/(cRatio*eta)) ELSE IF (mode .EQ. 4) THEN C2*** Cross flow, both streams mixed effectiveness = ((1/(1-EXP(-ntu)))+ & (cRatio/(1-EXP(-ntu*cRatio)))-(1/(-ntu)))**(-1) ELSE C2*** One stream is mixed and one is unmixed. Determine whether the C2*** minimum or maximum capacity rate stream is mixed. IF ( (ABS(Cap1-cMin).LT.small .AND. mode.EQ.5) .OR. & (ABS(Cap2-cMin).LT.small .AND. mode.EQ.6) ) THEN C2*** Cross flow, stream with minimum capacity rate unmixed effectiveness = (1-EXP(-cRatio*(1-EXP(-ntu))))/cRatio ELSE C2*** Cross flow, stream with maximum capacity rate unmixed effectiveness = 1-EXP(-(1-EXP(-ntu*cRatio))/cRatio) ENDIF ENDIF C1*** Determine leaving conditions for the two streams qMax = MAX(cMin,small)*(In1-In2) Out1 = In1 - effectiveness*qMax/MAX(Cap1,small) Out2 = In2 + effectiveness*qMax/MAX(Cap2,small) RETURN END SUBROUTINE WCOILOUT (PATM,CPAIR,CPVAP,HFG,TKELMULT, & TABSADD,PAMULT,PABSADD,MAir, & TAirEnt,WAirEnt,HAirEnt,HAirLvg, & UAExt,TAirLvg,WAirLvg,QSen,ErrStat) C*********************************************************************** C* Copyright ASHRAE. Toolkit for HVAC System Energy Calculations C*********************************************************************** C* SUBROUTINE: WCOILOUT C* C* LANGUAGE: FORTRAN 77 C* C* PURPOSE: Calculate the leaving air temperature, C* the leaving air humidity ratio and the C* sensible cooling capacity of wet cooling C* coil. C*********************************************************************** C* INPUT VARIABLES C* MAir Dry air mass flow rate (kg/s) C* TAirEnt Entering air dry bulb temperature (C) C* WAirEnt Entering air humidity ratio (-) C* HAirEnt Entering air enthalpy (J/kg) C* HAirLvg Leaving air enthalpy (J/kg) C* UAExt Heat transfer coefficient for external surface (W/C) C* C* OUTPUT VARIABLES C* TAirLvg Leaving air dry bulb temperature (C) C* WAirLvg Leaving air humidity ratio (-) C* Qsen Sensible heat transfer rate (W) C* ErrStat Error status indicator, 0 = ok, 1 = error (-) C* C* PROPERTIES C* Patm Atmospheric pressure (-) C* CpAir Specific heat of air (J/kg C) C* CpVap Specific heat of water vapor (J/kg C) C*********************************************************************** C MAJOR RESTRICTIONS: Assumes condensate at uniform temperature. 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: TAIRSAT C HUMTH C C C REVISION HISTORY: None C C REFERENCE: Elmahdy, A.H. and Mitalas, G.P. 1977. C "A Simple Model for Cooling and C Dehumidifying Coils for Use In Calculating C Energy Requirements for Buildings," C ASHRAE Transactions,Vol.83 Part 2, C pp. 103-117. C*********************************************************************** C INTERNAL VARIABLES: C capAir Air capacity rate (W/C) C ntu Number of heat transfer units (-) C effectiveness Heat exchanger effectiveness (-) C hCondSat Saturated air enthalpy at temperature of C condensate (J/kg) C tempCond Temperature of condensate (C) C*********************************************************************** INTEGER ErrStat REAL ntu,MAir,HFG,CPAIR,PATM DATA small/1.E-9/ ErrStat = 0 C1*** Determine the temperature effectiveness, assuming the temperature C1*** of the condensate is constant (Cmin/Cmax = 0) and the specific heat C1*** of moist air is constant capAir = MAir*(CpAir+WAirEnt*CpVap) ntu = UAExt/MAX(capAir,small) effectiveness = 1 - EXP(-ntu) C1*** Calculate coil surface enthalpy and temperature at the exit C1*** of the wet part of the coil using the effectiveness relation effectiveness = MAX(effectiveness,small) hCondSat = HAirEnt-(HAirEnt-HAirLvg)/effectiveness C1*** Calculate condensate temperature as the saturation temperature C1*** at given saturation enthalpy tempCond = TAIRSAT(PATM,CPAIR,CPVAP,HFG,TKELMULT, & TABSADD,PAMULT,PABSADD,hCondSat) C1*** Calculate exit air conditions and sensible heat transfer IF (tempCond .LT. DEWPOINT(PATM,TKELMULT,TABSADD, & PAMULT,PABSADD,WAirEnt)) THEN TAirLvg = TAirEnt-(TAirEnt-tempCond)*effectiveness WAirLvg = HUMTH(CPAIR,CPVAP,HFG,TAirLvg,HAirLvg) ELSE WAirLvg = WAirEnt TAirLvg = DRYBULB(CPAIR,CPVAP,HFG,HAirLvg,WAirLvg) ENDIF Qsen = capAir*(TAirEnt-TAirLvg) RETURN END REAL FUNCTION DEWPOINT (PATM,TKELMULT,TABSADD, & PAMULT,PABSADD,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: SATTEMP 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*********************************************************************** 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= Patm*W/(.62198+W) C1*** Calculate dewpoint as saturation temperature at water vapor C1*** partial pressure DewPoint = SATTEMP(TKELMULT,TABSADD,PAMULT,PABSADD,pw) ENDIF 999 RETURN END REAL FUNCTION DRYBULB (CPAIR,CPVAP,HFG,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*********************************************************************** REAL CPAIR,CPVAP,HFG 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-Hfg*W)/(CpAir+CpVap*W) RETURN END REAL FUNCTION ENTHALPY (CPAIR,CPVAP,HFG,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*********************************************************************** REAL CPAIR,CPVAP,HFG C1*** Calculate the enthalpy as a function of dry bulb temperature and C1*** humidity ratio. hDryAir = CpAir*TDB hSatVap = Hfg + CpVap*TDB Enthalpy = hDryAir + W*hSatVap RETURN END REAL FUNCTION ENTHSAT (PATM,CPAIR,CPVAP,HFG,TKELMULT, & TABSADD,PAMULT,PABSADD,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*********************************************************************** REAL TKELMULT,TABSADD,PAMULT,PABSADD,PATM,CPAIR,CPVAP,HFG C1*** Calculate the saturation pressure at the given temperature. psat = SATPRESS (TKELMULT,TABSADD,PAMULT,PABSADD,TDB) C1*** Calculate the humidity ratio from the saturation pressure w = HUMRATIO (Patm,psat) C1*** Calculate the enthalpy as a function of dry bulb temperature C1*** and humidity ratio. ENTHSAT = ENTHALPY (CPAIR,CPVAP,HFG,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 (CPAIR,CPVAP,HFG,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*********************************************************************** 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-CpAir*TDB)/(Hfg+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 (PATM,TABSADD,RAIR,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*********************************************************************** C1*** Calculate the dry air density from perfect gas laws. pAir = 0.62198*Patm/(0.62198+W) RhoDry = pAir/RAir/(TDB+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 (TKELMULT,TABSADD,PAMULT,PABSADD,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*********************************************************************** REAL TKELMULT,TABSADD,PAMULT,PBSADD 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 = TKelMult*(T+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/PaMult - PAbsAdd RETURN END REAL FUNCTION SATTEMP (TKELMULT,TABSADD,PAMULT,PABSADD,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(TKELMULT,TABSADD,PAMULT,PABSADD,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 (PATM,CPAIR,CPVAP,HFG, & TKELMULT,TABSADD,PAMULT,PABSADD,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(PATM,CPAIR,CPVAP,HFG, & TKELMULT,TABSADD,PAMULT,PABSADD,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 (PATM,HFG,CPAIR,CPVAP,CPWAT, & TKELMULT,TABSADD,PAMULT,PABSADD,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*********************************************************************** REAL newW,HFG,PATM,PAMULT,CPWAT,CPVAP,CPAIR DATA itmax/20/ C1*** Initial temperature guess tBoil = SATTEMP (TKELMULT,TABSADD,PAMULT,PABSADD,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 (TKELMULT,TABSADD, & PAMULT,PABSADD,WetBulb) C1*** Determine humidity ratio for given saturation pressure wStar = HUMRATIO (Patm,psatStar) C1*** Calculate new humidity ratio and determine difference from known C1*** humidity ratio newW = ((Hfg-(CpWat-CpVap)*WetBulb)*wStar- & CpAir*(TDB-WetBulb))/(Hfg+CpVap*TDB & -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